3 D refraction correction and extraction of clinical parameters from spectral domain optical coherence tomography of the cornea

Capable of three-dimensional imaging of the cornea with micrometer-scale resolution, spectral domain-optical coherence tomography (SDOCT) offers potential advantages over Placido ring and Scheimpflug photography based systems for accurate extraction of quantitative keratometric parameters. In this work, an SDOCT scanning protocol and motion correction algorithm were implemented to minimize the effects of patient motion during data acquisition. Procedures are described for correction of image data artifacts resulting from 3D refraction of SDOCT light in the cornea and from non-idealities of the scanning system geometry performed as a pre-requisite for accurate parameter extraction. Zernike polynomial 3D reconstruction and a recursive half searching algorithm (RHSA) were implemented to extract clinical keratometric parameters including anterior and posterior radii of curvature, central cornea optical power, central corneal thickness, and thickness maps of the cornea. Accuracy and repeatability of the extracted parameters obtained using a commercial 859nm SDOCT retinal imaging system with a corneal adapter were assessed using a rigid gas permeable (RGP) contact lens as a phantom target. Extraction of these parameters was performed in vivo in 3 patients and compared to commercial Placido topography and Scheimpflug photography systems. The repeatability of SDOCT central corneal power measured in vivo was 0.18 Diopters, and the difference observed between the systems averaged 0.1 Diopters between SDOCT and Scheimpflug photography, and 0.6 Diopters between SDOCT and Placido topography. © 2010 Optical Society of America OCIS code: (170.4500) Optical coherence tomography; (170.4460) Ophthalmic optics and devices; (170.3880) Medical and biological imaging; (190.5330) Photorefractive optics References and links 1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254(5035), 4 (1991). 2. M. Wojtkowski, V. Srinivasan, T. Ko, J. Fujimoto, A. Kowalczyk, and J. Duker, “Ultrahigh-resolution, highspeed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12(11), 2404–2422 (2004). 3. A. F. Fercher, C. K. Hitzenberger, G. Kamp, and S. Y. El-Zaiat, “Measurement of intraocular distances by backscattering spectral interferometry,” Opt. Commun. 17, 6 (1995). 4. M. Pircher, E. Götzinger, R. Leitgeb, H. Sattmann, O. Findl, and C. Hitzenberger, “Imaging of polarization properties of human retina in vivo with phase resolved transversal PS-OCT,” Opt. Express 12(24), 5940 5951 (2004). 5. B. Cense, T. C. Chen, B. H. Park, M. C. Pierce, and J. F. de Boer, “In vivo depth-resolved birefringence measurements of the human retinal nerve fiber layer by polarization-sensitive optical coherence tomography,” Opt. Lett. 27(18), 1610 (2002). 6. B. Hofer, B. Považay, B. Hermann, A. Unterhuber, G. Matz, and W. Drexler, “Dispersion encoded full range frequency domain optical coherence tomography,” Opt. Express 17(1), 7 – 24 (2009). (C) 2010 OSA 26 April 2010 / Vol. 18, No. 9 / OPTICS EXPRESS 8923 #122547 $15.00 USD Received 11 Jan 2010; revised 6 Apr 2010; accepted 6 Apr 2010; published 14 Apr 2010 7. M. Yamanari, Y. Lim, S. Makita, and Y. Yasuno, “Visualization of phase retardation of deep posterior eye by polarization-sensitive swept-source optical coherence tomography with 1-μm probe,” Opt. Express 17(15), 12385 (2009). 8. R. Zawadzki, S. Jones, S. Olivier, M. Zhao, B. Bower, J. Izatt, S. Choi, S. Laut, and J. Werner, “Adaptive-optics optical coherence tomography for high-resolution and high-speed 3D retinal in vivo imaging,” Opt. Express 13(21), 8532 – 8546 (2005). 9. M. Zhao, and J. A. Izatt, “Single-camera sequential-scan-based polarization-sensitive SDOCT for retinal imaging,” Opt. Lett. 34(2), 205 (2009). 10. L. An, and R. K. Wang, “In vivo volumetric imaging of vascular perfusion within human retina and choroids with optical micro-angiography,” Opt. Express 16(15), 11438 – 11452 (2008). 11. R. J. Zawadzki, C. Leisser, R. Leitgeb, M. Pircher, and A. F. Fercher, “Three-dimensional ophthalmic optical coherence tomography with a refraction correction algorithm,” Proc. SPIE 5140, 8 (2003). 12. V. Westphal, A. Rollins, S. Radhakrishnan, and J. Izatt, “Correction of geometric and refractive image distortions in optical coherence tomography applying Fermat’s principle,” Opt. Express 10, 8 (2002). 13. M. V. Sarunic, B. E. Applegate, and J. A. Izatt, “Real-time quadrature projection complex conjugate resolved Fourier domain optical coherence tomography,” Opt. Lett. 31(16), 2426 (2006). 14. M. V. Sarunic, S. Asrani, and J. A. Izatt, “Imaging the Ocular Anterior Segment With Real-Time, Full-Range Fourier-Domain Optical Coherence Tomography,” Arch. Ophthalmol. 126(4), 6 (2008). 15. M. Tang, Y. Li, M. Avila, and D. Huang, “Measuring total corneal power before and after laser in situ keratomileusis with high-speed optical coherence tomography,” J. Cataract Refract. Surg. 32(11), 18438 (2006). 16. S. Radhakrishnan, A. M. Rollins, J. E. Roth, S. Yazdanfar, V. Westphal, D. S. Bardenstein, and J. A. Izatt, “Real-time optical coherence tomography of the anterior segment at 1310 nm,” Arch. Ophthalmol. 119, 7 (2001). 17. Y. Wang, J. S. Nelson, Z. Chen, B. J. Reiser, R. S. Chuck, and R. S. Windeler, “Optimal wavelength for ultrahigh-resolution optical coherence tomography,” Opt. Express 11(12), 1411 (2003). 18. Y. Yasuno, M. Yamanari, K. Kawana, T. Oshika, and M. Miura, “Investigation of post-glaucoma-surgery structures by three-dimensional and polarization sensitive anterior eye segment optical coherence tomography,” Opt. Express 17(5), 3980 (2009). 19. M. Yamanari, S. Makita, and Y. Yasuno, “Polarization-sensitive swept-source optical coherence tomography with continuous source polarization modulation,” Opt. Express 16(8), 5892 (2008). 20. M. Pircher, E. Götzinger, R. Leitgeb, A. F. Fercher, and C. Hitzenberger, and C. K. Hitzenberger, “Measurement and imaging of water concentration in human cornea with differential absorption optical coherence tomography,” Opt. Express 11(18), 2190 (2003). 21. C. Kerbage, H. Lim, W. Sun, M. Mujat, and J. F. de Boer, “Large depth-high resolution full 3D imaging of the anterior segments of the eye using high speed optical frequency domain imaging,” Opt. Express 15(12), 7117 (2007). 22. M. Akiba, N. Maeda, K. Yumikake, T. Soma, K. Nishida, Y. Tano, and K. P. Chan, “Ultrahigh-resolution imaging of human donor cornea using full-field optical coherence tomography,” J. Biomed. Opt. 12(4), 4 (2007). 23. K. Grieve, A. Dubois, M. Simonutti, M. Paques, J. Sahel, J.-F. Le Gargasson, and C. Boccara, “In vivo anterior segment imaging in the rat eye with high speed white light full-field optical coherence tomography,” Opt. Express 13(16), 6286 (2005). 24. A. N.Kuo, M. Zhao, and J. A. Izatt, “Corneal aberration measurement with three-dimensional refraction correction for high-speed spectral domain optical coherence tomography,” ARVO E-Abstract, 2 (2009). 25. M. Tang, Y. Li, and D. Huang, “Corneal topography and power measurement with optical coherence tomography,” ARVO E-Abstract 5791(A93), 2 (2009). 26. M. Gora, K. Karnowski, M. Szkulmowski, B. J. Kaluzny, R. Huber, A. Kowalczyk, and M. Wojtkowski, “Ultra high-speed swept source OCT imaging of the anterior segment of human eye at 200 kHz with adjustable imaging range,” Opt. Express 17(17), 14880 (2009). 27. A. N. S. Institute, “Corneal topography systems standard terminology, requirements,” ANSI 80, 23–2008 (2008). 28. V. A. D. P. Sicam, M. Dubbelman, and R. G. L. van der Heijde, “Spherical aberration of the anterior and posterior surfaces of the human cornea,” J. Opt. Soc. Am. A 23(3), 544 (2006). 29. R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. 66(3), 207 (1976). 30. V. A. D. P. Sicam, J. Coppens, T. J. T. P. van den Berg, and R. G. L. van der Heijde, “Corneal surface reconstruction algorithm that uses Zernike polynomial representation,” J. Opt. Soc. Am. A 21, 7 (2004). 31. C. K. Hitzenberger, “Optical measurement of the axial eye length by laser doppler interferometry,” Invest. Ophthalmol. Vis. Sci. 32, 9 (1991). 32. W. Drexler, C. K. Hitzenberger, A. Baumgartner, O. Findl, H. Sattmann, and A. F. Fercher, “Investigation of disperison effects in ocular media by multiple wavelength partial coherence interferometry,” Exp. Eye Res. 66(1), 25 (1998). 33. R. C. Lin, M. A. Shure, A. M. Rollins, J. A. Izatt, and D. Huang, “Group index of the human cornea at 1.3-mm wavelength obtained in vitro by optical coherence domain reflectometry,” Opt. Lett. 29(1), 83 (2004). 34. I. Grulkowski, M. Gora, M. Szkulmowski, I. Gorczynska, D. Szlag, S. Marcos, A. Kowalczyk, and M. Wojtkowski, “Anterior segment imaging with Spectral OCT system using a high-speed CMOS camera,” Opt. Express 17(6), 4842 (2009). (C) 2010 OSA 26 April 2010 / Vol. 18, No. 9 / OPTICS EXPRESS 8924 #122547 $15.00 USD Received 11 Jan 2010; revised 6 Apr 2010; accepted 6 Apr 2010; published 14 Apr 2010 35. T. O. Salmon, and L. N. Thibos, “Videokeratoscope-line-of-sight misalignment and its effect on measurements of corneal and internal ocular aberrations,” J. Opt. Soc. Am. A 19(4), 657 (2002).


Introduction
Spectral domain optical coherence tomography (SDOCT) has recently emerged as a powerful new tool for noninvasive ophthalmic imaging [1][2][3][4][5][6][7][8][9][10].Accurate determinations of corneal curvature, refractive power, and wave aberrations are important in clinical diagnosis and management.Unlike Placido ring based systems which only image the anterior corneal surface, optical coherence tomography (OCT) offers the ability to image both the anterior and posterior refractive surfaces of the cornea for the determination of clinical parameters.Initial reports of real time anterior segment OCT utilized illumination light in the 1310 nm region [11][12][13][14][15][16][17][18][19][20][21].Some groups employed sources centered at 760nm [22] and 650nm [23] respectively.However more recently there have appeared several spectral domain OCT (SD-OCT) systems utilizing the commonly used retinal 840nm and 859nm light sources specifically for high resolution corneal imaging [24][25][26].These OCT systems offer higher resolution than that found in photographic tomographic imaging techniques such as Scheimpflug photography, which should result in more accurate modeling of the imaged structures.
All OCT imaging systems are subject to the effects of refraction at surfaces corresponding to interfaces between regions of differing refractive index within a sample.Additionally, departures from ideal optical performance of OCT scanners (such as non-telecentricity) will generate additional geometric image distortions in anterior segment OCT imaging [11][12]15].To produce accurate volumetric representations of the cornea, the effects of the nontelecentric SDOCT scanning source and refraction in three-dimensions (3D) of this source light through the cornea must be determined and applied to correct the raw images for accurate clinical parameter extraction.Prior methods to correct for refraction in anterior segment OCT images have been described [12,[14][15].However these methods were limited to 2D processing, which assumes that refraction occurs only within the plane of individually acquired B-scans (defined as sets of A-scans acquired along a straight line comprising a cross-sectional image of the sample).For a curved 3D structure such as the cornea, 2D refraction correction is correct only if the sample is rotationally conically symmetric about an axis passing through the apex of the sample, and if the acquired B-scan data all pass exactly through the apex point.The first condition is rarely true for realistic samples such as the human cornea, especially if they have been modified surgically.The second condition is true only for idealized radial scan patterns, which may not be optimal because they oversample the central region and undersample the outer region and because unavoidable patient motion or operator misalignment of the OCT system make it difficult to acquire perfectly centered radial scans.The only published report of 3D refraction correction [11] does not provide a description of the method used and only shows in vitro correction of a model eye.
In this paper, we describe algorithms for three-dimensional correction of scanner nontelecentricity and for refraction at the corneal surface in in vivo anterior segment OCT imaging.Motivated by sag error analysis, a radial scanning pattern is employed to mitigate errors due to axial patient motion during volume scan acquisition.Algorithms for extraction of clinical corneal parameters from corrected image data including anterior and posterior radii of curvature, central corneal optical power, and thickness maps of the cornea are described.

Requirements for quantitative corneal biometry: sag error analysis
OCT imaging has had a significant impact on ophthalmic diagnostics through its ability to elucidate morphological correlates of disease processes by allowing for qualitative imaging of gross anatomical changes as well as quantitative measurement of corneal and retinal layer thicknesses.Thickness measurements of conformal layers are relatively straightforward to perform and are immune to patient motion on time scales long compared to individual B-scans, since they correspond to monitoring the relative motion of multiple layer boundaries within individual images.For quantitative corneal biometry, however, absolute rather than relative surface measurements are required to calculate clinically relevant refractive parameters of these high quality optical surfaces.Despite the resolution and depth-resolved imaging advantages of OCT compared to Placido ring and Scheimpflug photography-based approaches, even high-speed Fourier domain OCT approaches have not yet overcome the effective resolution loss due to patient motion during the seconds-long acquisition times required for volumetric OCT data set acquisition.
To illustrate the demanding requirements for performing quantitative corneal biometry with OCT, we present the following simple sag error analysis.The keratometric equation commonly used to calculate corneal power is given by [15]: Here, P k is the power of the cornea expressed in Diopters, n k is the keratomeric refractive index usually given as 1.3375 [27], n 0 is the refractive index of air (1.0), and R a is the radius of curvature of the anterior corneal surface.Figure 1.illustrates the principle of the sag error analysis.The z axis represents the axial direction of the sag error or motion, and r is the radial (or lateral) distance from the z axis.Assuming that we wish to use OCT to measure corneal power with an accuracy at least matching the resolution of clinical refraction measurements (0.25D), the change in corneal power as a function of a change in radius of curvature (error) is derived as as the desired corneal power resolution accuracy and using a nominally average R a = 7.79mm [28], then 46µm d = z and 234µm dr = at a lateral position of r = 1.5mm , corresponding to the edge of a 3mm diameter fitting area.Thus, any errors in OCT measurement of the corneal surface location resulting from mis-calibrations of the scanning apparatus or patient motion must be less than 46µm axially and 234µm laterally on the anterior corneal surface (where the requirement is more stringent since the index of refraction contrast is higher).In this study, these requirements were met by a combination of careful OCT scanner alignment and calibration, in combination with the use of a radial scan pattern-based patient motion correction algorithm.Fig. 1.Patient bulk motion and sag error analysis.Here the z axis represents the axial direction of motion, and r is the radial (lateral) distance from the z axis.At a lateral position of r= 1.5mm corresponding to a 3mm diameter fitting area, axial motion (dz) of more than 46 µm and lateral motion (dr) of more than 234 µm will change the calculated power of an average cornea with an anterior radius of curvature Ra = 7.79mm by more than 0.25D.

Non-telecentric scan correction
Anterior segment OCT scanners described to date [11,[14][15] have utilized a nominally telecentric scanning geometry, here defined as comprised of mutually parallel A-scans which are normal to the corneal surface at its apex.However, due to small optical misalignments and aberrations in real scanners, even small departures from this ideal may detract from the stringent accuracy required for corneal biometry.To correct for non-telecentric scanning geometric image distortion effects, we extend the 2D non-telecentric correction method previously described by Westphal et al. [12,16] to the 3D case.The extension is described in Eq. ( 2 Here x, y and ρ are the original physical target coordinates, and x', y' and ρ' are the uncorrected raw positions caused by nontelecentric scanning.D is the pivot scanning distance, z is the axial coordinate which is illustrated in Fig. 2, and φ is the scanning angle.To correct for the effects of corneal surface refraction in 3D, a generalized vector implementation of Snell's law was utilized as illustrated in Fig. 3.The unit normal vector on the corneal epithelium at point C can be expressed as:

Three dimensional refraction correction using vector Snell's law
The non-telecentric incident light MC → is a unit vector denoted as (a,b,c) and the unit vector of the refracted ray In any realistic biological or industrial sample, the magnetic current densities M are zero.
According to Maxwell's equations, the electric magnetic field boundary condition can be represented: .
Here, in E and ref E are the incident and refracted electric fields individually, and × is the vector product.After straightforward vector manipulation, our derivation of Snell's law for use in three-dimensional OCT applications can be represented in the following equations: OPL is the optical path length of the refracted ray.Other definitions are illustrated in Fig. 3. To obtain 3D refraction correction, we first obtain the corrected z axis physical position using vector projection of the optical path length.The lateral physical position is then obtained using Eq. ( 7) based on the retrieved z axial value.This procedure may optionally be iterated for all layers to correct the entire volume data.Any correction for non-telecentric illumination is included in the incident vector IN V ( , , ) a b c = .Thus by using this vector form of Snell's law, we can correct refraction issues in one step of computation.For the telecentric case, IN V ( , , ) a b c = becomes (0, 0,1) . Hence Eq. ( 7) can be further simplified in this case.

Radial scanning pattern for patient motion correction with Zernike interpolation
The two most common scan patterns employed in ophthalmic OCT imaging are raster scanning (stacking laterally displaced B-scans to form a box of images) and radial scanning (with each B-scan is centered on a rotational axis).The advantage of the raster scanning is that it provides complete and evenly sampled x, y, and z data throughout the volume (each Bscan creates the x-axis, the stacking of B-scans creates the y-axis, and the individual points in an A-scan define the z axis).This is ideal for a Cartesian implementation of 3D refraction correction.However, since the SNR of the OCT image falls off with distance from the corneal apex, B-scans obtained near the periphery of raster scans suffer from low SNR throughout their length which can complicate image processing.While the radial scanning pattern has the advantage of high SNR for at least a portion of all B-scan frames (surrounding the apex area), radial scan pattern data is unevenly sampled with denser sampling centrally and lower sampling peripherally.For a Cartesian implementation as in Eq. ( 7), this requires that the data first be interpolated and then resampled evenly.A very significant advantage of the radial scan pattern is for patient motion correction.Assuming that the radial B-scan pattern is centered on a high SNR point near the corneal apex, all B-scans in the data set then share a common physical point which is re-acquired at the B-scan rate (typically ~20Hz).The offset of these sequentially acquired common A-scans from each other is then be used to axially re-align all A-scans in each B-scan.
To re-sample motion-corrected radial OCT volumes to a uniform Cartesian distribution suitable for 3D refraction correction using Eq. ( 7), we employ Zernike fitting as illustrated in Fig. 4. Zernike polynomials are a set of complete orthogonal polynomials defined on a unit circle, which mathematically form a complete compact basis to represent a 3D surface in polar coordinates.This is useful for representing structures which have an axis of rotation, such as the cornea.Additionally Zernike polynomials can be used to estimate the spherical aberration of the corneal surface.These polynomials are represented in polar coordinates by [29][30]: The indices n and m are the radial degree and the azimuthal frequency respectively; i is a mode-ordering number; R is used to normalize the measurement data within a dimensionless unit circle; r is the normalized radius; and i c is the Zernike elevation interpolation coefficient.( , ) A Rr θ is the Zernike 3D representation of the fitting surfaces.Another issue with the newly generated, 3D refraction corrected volumetric data is that the refraction corrected pixels are no longer evenly sampled which is illustrated conceptually in Fig. 5(a).The complete endothelial surface must be reconstructed from this uneven distribution to obtain a corneal thickness map (as each epithelial normal may no longer intersect a physical endothelial pixel) or to generate wave front aberration analyses of the endothelium.To reconstruct an evenly-sampled endothelial surface, Zernike polynomial 3D interpolation was again employed to obtain the even-grid refraction corrected endothelial surfaces.

Recursive half searching algorithm for pachymetry map
After reconstructing the endothelium, the corneal pachymetry (thickness) was measured along the epithelial surface normal vector, which is represented by: It is difficult to derive an analytical solution for a pachymetry map between Eq. ( 8)b) and Eq.(9).A numerical solution called the recursive half searching algorithm (RHSA) was used to solve this problem.The idea is to search for the intersection between the epithelial surface normal equation and the endothelium.The starting guess point was the refraction uncorrected endothelium ( , , ) x y z .The new coordinates of ' ' ' ( , , ) x y z of the first step were estimated by manipulating Eq. ( 9) into the form: ε is the stop condition which we have defined as 0.5µm in our case.Step 3 will be followed if it does not reach the convergent accuracy. ( The above search method recursively halves the difference of the search distance.Hence it converges rapidly, typically only requiring 7 steps to define the intersection point.The searching scheme is illustrated in Fig. 6.Fig. 6.Principle of searching intersection points between the surface normal equation of epithelial and refraction-corrected endothelial surface.The top surface is epithelium and the second is endothelium.Z ҩ : Obtained using Eq. ( 9) after half distance recursion; Z Ҫ : retrieved using Zernike coefficients in Eq. ( 10).

Central corneal power calculation
Best spherical fitting by least squares fitting of the epithelial and endothelial surfaces using the following equation was employed to extract the radii of curvature [28].
Here c is the curvature of this curve, 0 0 ( , ) r z is the coordinate of the vertex and k is the asphericity of this curve.
The central corneal power for OCT instruments was then obtained using the following thin lens equation: In the above equation, the unit of radius is in meters and the optical power has units of Diopters (inverse meters).R a is the radius of the anterior surface, and R p is the radius of the posterior surface as derived above.The refractive index of the cornea was taken to be 1.385 ± 0.002 as measured empirically in vitro at 859 nm using three human donor corneas in our lab.This value is close to previous measurement results [31][32][33].The refractive index of aqueous was taken from the literature as 1.336 [15].P t is then the overall power.k P was obtained using the simulated keratometric refractive index 1.3375 [27].
This keratometric equation is used by topography instruments to determine corneal power.

Experimental methods
To implement and test the algorithms described in the previous section, we used an anterior segment OCT system consisting of a commercially available spectral domain OCT retinal imager with a corneal adapter (Bioptigen Inc., Research Triangle Park, NC).The bandwidth of the OCT spectrometer is 53.9nm centered at 859nm.The axial resolution is 6.0µm in air and 4.3µm in cornea.The sensitivity is 104db.To perform scanner calibration, SDOCT raster volume scans comprising 50 laterally displaced B-scans of 1000 A-scans each were acquired of a dual axis linear scale square calibration target with 25µm markings (Edmund Optics).Three raster scanning field of views (FOVs) of 6mm × 6mm, 7mm × 6mm and 8mm × 8mm were scanned to obtain the individual pivot scanning distance (D, Eq. ( 2) using least square curve fitting.The three pivot scanning distances were averaged to give the value reported for D in the results section.
To test the artifact removal and 3D refraction correction algorithms, a rigid gas permeable (RGP) contact lens (Conforma; Norfolk, VA) of known refractive index, curvature, central thickness, and power served as the phantom target.The refractive index of the contact lens was 1.466.This contact lens was mounted horizontally (simulating patient corneal orientation in clinical use) onto a stage for imaging and its epithelial and endothelial surfaces were imaged three times using a radial scanning pattern.Each radial scan included 50 B-frames consisting of 1000 A-scans per frame.The radial scanning range was 6mm and the scan integration time of each A-scan was 50µs.
For in vivo imaging, volumetric corneal images were obtained from the eyes of three healthy volunteers under an Investigational Review Board approved protocol.The scanning protocol reflected that for phantom imaging: 50 radial B-scans centered on the visual axis, each B-scan consisting of 1000 A-scans spanning 6 mm.For comparison, each volunteer's eye was also imaged using a clinical topography unit (Atlas; Carl Zeiss Meditec; Dublin, CA) and a clinical Scheimpflug camera (Oculus; Pentacam; Lynnwood, WA) using standard clinical protocols.For repeatability information, one subject was imaged three times on all three instruments.
All SDOCT volume data was processed offline using custom software to manually segment the anterior and posterior surfaces and to apply the correction algorithms as described in the previous section.Axial motion artifact correction was performed by aligning each B-scan image with respect to the axis of scan rotation.The radii of curvature for the anterior and posterior surfaces of the target structures were then calculated using central best spherical fitting (6mm diameter for the phantom and 3mm diameter for each in vivo cornea).A total central power (P t ) was then calculated from these radii using the thin lens equation.

Nontelecentric calibration
Figure 7 (a) shows an axially summed volume projection of a 50 B-scan raster volume of the flat calibration test target.Figure 7 (b) shows the nontelecentric field curvature image which is formed by overlaying 50 B-frames in the (r,z) plane.The top surface of the overlaid images represents the field curve distortion caused by non-telecentric illumination of the flat test object (Fig. 7(b)).A least square circular curve fitting of the field curvature was used to obtain the pivot scanning distance D. The fitting radius is equal to the pivot scanning distance D. The fitted curvature is overlaid on the plot in red color.To retrieve the pivot scanning distance D, three raster scanning field of views (FOVs) of 6mm×6mm, 7mm×7mm and 8mm×8mm were employed to obtain the pivot scanning distance D. The average result was 205.6mm in our system setup which is quite close to the ideal telecentric illumination condition.This non-telecentric quantification was used to correct all subsequent data by employing Eq. (1).In this nearly ideal telecentric illumination condition, the influence of nontelecentric effects on the overall corneal power was very limited; for example, it was less than 0.1D for the corneal power measurements reported in the following sections.

Contact lens phantom results
Figure 8(a) is a photograph of the rigid gas-permeable contact lens phantom used for a first test of the described algorithms applied on SDOCT acquisitions of this phantom.Figure 8(b) is the pachymetry map (thickness map) of the contact lens obtained using the RHSA algorithm.The central thickness was 136.2µm.This corresponds closely to the manufacturer's specified value of 130±20µm.The manufacturer's value for the base curvature of the contact lens was 7.6mm.Our SDOCT measured result was 7.63±0.02mm.Knowing that this contact lens had a refractive power of −3.00 D in air, the expected curvature of the anterior surface was calculated using the thin lens equation and was determined to be nominally 7.99mm.The experimental value obtained using SDOCT was 8.04±0.03mm.Additionally we measured the same surface via a commercial interferometric optical surface profiler (Zygo TM ), which reported the value to be 8.10mm.Note that the field of view of the Zygo TM was only 200µm×200µm.These data indicate that the SDOCT derived values agree well with the expected values.The measurement results for this phantom contact lens are listed in Table 1.The error for the total power measurement was 0.18D, less than the target 0.25D resolution for a clinical prescription.

In vivo corneal clinical parameter extraction
Clinical biometric parameters were extracted from SDOCT corneal images of three healthy volunteers using the methods described above.A typical single corneal B-scan image from a radially acquired volumetric data set is illustrated in Fig. 9 (a).For comparison of the extracted clinical parameters, each volunteer's eye was also imaged using a clinical topography unit (Atlas; Carl Zeiss Meditec; Dublin, CA) and a clinical Scheimpflug camera (Oculus; Pentacam; Lynnwood, WA). Figure 9 (b) shows a plot illustrating the patient bulk motion prior to patient motion correction via B-scan alignment to the center A scans.To estimate the patient axial bulk motion, we imaged a patient who was very experienced in retinal imaging.The obtained maximum patient motion was 156µm in the axial direction.The thickness mapping (pachymetry map) of an in vivo corneal volume obtained using our RHSA algorithm is illustrated in Fig. 9 (c).Given the segmented corneal images (50 frames x 1000 points x 2 surfaces), a 2.13 GHz dual core CPU with 4Gb RAM takes 59.83 seconds on average to artifact correct, interpolate, and extract the described clinical parameters in Matlab (MathWorks; Natick, MA).The central corneal power for the commercial instruments was obtained from the on-board analysis provided by each instrument.The results are listed in the Table 2 for comparison.We also tested the in vivo repeatability of our SDOCT system by measuring one patient's eye in triplicate on each instrument.The standard deviations of the measurement results by each device are listed in the following Table 3.

Discussion
The contact lens phantom measurement indicated that the error for the total power measured by SDOCT was 0.18D, less than the 0.25D resolution of a clinical prescription.This implies that our 3D refraction correction algorithm works well for a stationary object.For the in vivo measurements, if we consider corneal power as determined by both the anterior and posterior surface as in Eq. ( 14), the central corneal power calculated by SDOCT, Scheimpflug photography, and single surface power from topography were also very similar.The difference observed between the systems averaged 0.1D between SDOCT and Scheimpflug photography, and 0.6D between SDOCT and topography.Our SDOCT system implementation has a comparable performance with the Scheimpflug camera in clinical use.Although we are not aware of the exact algorithm used by the Pentacam for calculating corneal power, presumably information from both surfaces is used in its calculation which may explain why our SDOCT corneal powers are closer to Pentacam values.In contrast, topography uses only the curvature of the anterior surface and a simulated keratometric refractive index (1.3375) to compensate for its inability to directly measure the posterior surface contribution.Hence it would be reasonable to expect that differences between SDOCT and topography would be larger than that of SDOCT versus Scheimpflug photography.Additionally, we used all of the data within the central 3mm for 3D best spherical fitting.Topographic methods interpolate central information as the rings typically do not extend completely to the middle.However the central areas are most important part for clinical applications.
The repeatability results illustrate that each system has good repeatability.The standard deviations are all less than the 0.25D clinical resolution.We have shown with sag error analysis in Fig. 1 the estimated influence of bulk motion on corneal power calculations.The sag error analysis shows that the axial bulk motion must be less than 46µm and the lateral bulk motion must be less than 234µm if we are to obtain an accuracy of at least 0.25D.In the actual imaging procedure, it was common to have at least 150µm axial motion for imaged subjects.Figure 9 (b) showed the axial bulk motion of a patient experienced in imaging.While lateral curvature may change due to irregularities in corneal shape, presumably in a radial scan pattern, the center A-scan should be capturing the same region of the cornea in the same location in space.Despite this, the center A-scan varied axially by a maximum of 156µm in this particular patient.With an axial change of 150µm, the corresponding error in corneal power would be 0.82D; with axial motion of 200 µm, 1.1D of error would be expected.Potential solutions to remove bulk motion would be full field imaging, ultra high speed imaging [26,34] or better bulk motion correction algorithms.Other potential sources of error include acquisition artifacts (decentration of the radial scan from the corneal apex) and segmentation errors.
The 3D refraction correction algorithm is a more generalized method than that of the specific 2D case.The 2D case is restricted to rotationally symmetric surfaces, such as roughly spherical normal corneas.Pathologic corneas with localized disturbances and global asymmetry are unlikely to meet this restriction, and refraction of light is unlikely to occur within the plane of the B-scan as the 2D case is usually implemented.This means that 2D refraction correction would be unable to recover truly accurate spatial information from pathologic corneas.Additionally, the generalized 3D algorithm can degrade into the 2D version for a rotationally symmetric surface.Then the 3D equation affects only the pixels in the local plane of the particular scan.For example in the x-z plane / 0 , and this is exactly the 2D refraction correction method for the non-telecentric illumination case: While our artifact removal algorithms (for non-telecentricity and refraction) and the RHSA for pachymetry are general to any cornea imaged by SDOCT, it should be noted that the described clinical parameter extraction for refractive power and central radii of curvature apply only to normal, approximately spherical corneas.Other specific clinical applications, such as screening for corneas with local curvature disturbances, will require specific analyses for the desired pathology.These custom analyses should be applied after first using the acquisition artifact removal processes described in this paper to create spatially accurate representations of the cornea from SDOCT data.The visual axis of the patient may not necessarily coincide with the axis defined by the anatomical corneal apex [35].The induced power error is estimated around 0.067D using a low order spherical fitting algorithm across a 3mm diameter zone of the central cornea.This is not a clinically significant effect for this low order aberration.Examination of higher order aberrations or combining measurements from instruments with differing reference axes would require formal alignment of the different axes.

Conclusions
In conclusion, we present here the first in vivo 3D algorithms correcting for non-telecentric scan deformation and refraction correction in SDOCT.Zernike polynomial 3D reconstruction and a novel numerical recursive half searching algorithm (RHSA) were developed to compute the thickness map of corneal volumes.Both a phantom contact lens and healthy volunteers were imaged and analyzed using the described protocols.The error for the total power in a static imaging test of the phantom was 0.18 D, which is less than the 0.25D clinical resolution requirement.This indicates that our methods accurately extracted the physical parameters and power of a stationary contact lens phantom from SDOCT images.In vivo clinical parameters were extracted from SDOCT images of three subjects and the results compared with clinical Placido topography and Scheimpflug camera systems.The overall optical power difference was 0.1 Diopters between SDOCT and Scheimpflug photography and 0.6 Diopters between SDOCT and Placido topography.

Fig. 3 .
Fig. 3. Schematic of the principle of 3D refraction correction.θ1: Angle between surface normal of epithelium and z-axis; θin: incidence angle; θRef: refraction angle; OPL: optical path length.∆z: Projection of OPL on z-axis.VIN = (a,b,c) is the incidence unit vector;(x0,y0,z0) is the coordinate of the epithelium; (x,y,z) is the coordinate of the refraction corrected endothelium.

Fig. 4 .
Fig. 4. Demonstration of the re-sampling issue in the radial scanning pattern.(a) Radial scanning pattern centered on the apex.(b) Radial scanning pattern non-even sampling issue.The sampling points are much sparser in the outer circles than that in the inner circles.(c) Zernike polynomial even-grid interpolation for both the epithelium and endothelium surfaces based on the noneven-grid radial sampling points.

Figure 4 (
Figure 4(a) illustrates the radial scanning pattern around the apex.Figure 4(b) demonstrates the non-even sampling issue with radial scan patterns.The sampling points are much denser in the inner circles than in the outer circles.To resample this uneven distribution into an evenly distributed data set, Zernike interpolation based on the radially acquired data is employed.Another issue with the newly generated, 3D refraction corrected volumetric data is that the refraction corrected pixels are no longer evenly sampled which is illustrated conceptually in Fig.5(a).

Fig. 5 .
Fig. 5. Illustration of Zernike 3D interpolation.(a) The refraction corrected voxels (e.g., voxels on the endothelial surface) are no longer on an even sampling grid after 3D Snell's law correction; (b) Regenerated evenly sampled voxels after Zernike 3D interpolation.Elevation value of arbitrary points is obtained using Zernike interpolation coefficients.

Fig. 7 .
Fig. 7. nontelecentric calibration (a) Surface volume projection of 50-frames raster scanning Bscan image with FOV of 6mm ×6mm.(b) Deformed field curvature in blue color of the top surface of the calibration target caused by nontelecentric illumination.Least square curve fitting in red color is used to find imaging scanning pivot distance D. The original segmented deformed FOV is 6mm ×6mm in this figure.All the units are in millimeters.

Fig. 8 .
Fig. 8. Demonstration of power measurement accuracy in a rigid gas-permeable contact lens phantom.(a) Phantom contact lens.(b) Pachymetry map (thickness map) of the contact lens using RHSA algorithm.(c) 3D surface rendering after 3D refraction correction.Axes are in millimeters.

Fig. 9 .
Fig. 9.In vivo anterior segment imaging results (a) In vivo corneal B-scan image which is saturated at the apex.(b) Patient bulk motion before aligning centered on the apex, which indicates maximum patient motion of 156µm in the axial direction during the duration of volume acquisition.50 anterior surfaces are plotted in this figure.The units are in millimeters.(c) Pachymetry map (thickness map) in micrometers of an in vivo corneal volume using our recursive half searching algorithm.The central thickness is 578.5µm.The image range is 6mm × 6mm.The color scale reflects those used clinically. ):