Optical distortion correction in Optical Coherence Tomography for quantitative ocular anterior segment by three-dimensional imaging

Sergio Ortiz; Damian Siedlecki; Ireneusz Grulkowski; Laura Remon; Daniel Pascual; Maciej Wojtkowski; Susana Marcos

doi:10.1364/OE.18.002782

1. Introduction

Anterior segment imaging has received a strong attention in the last few years. The availability of different imaging modalities and the increased resolution and visibility of the ocular structures has allowed several of the new developments to have an impact on the clinical practice. However, the increased quality of the images obtained from these methods has not been in general paralleled by the quantitative analysis of the raw data. Quantitative biometry is essential in many areas in physiological optics (i.e. customized eye modeling [1] understanding of the structural changes in the crystalline lens during accommodation [2] and aging [3], crystalline lens gradient index of refraction reconstruction, etc…) and ophthalmology (i.e. changes in the anterior and posterior corneal topography after refractive surgery [4],[5] and the accuracy [6] and evaluation of cataract surgery [7],[8], among others). Optical Coherence Tomography (OCT) is one of the most versatile and promising imaging techniques of the anterior segment of the eye, due to its high resolution, and the possibility of imaging the anterior and posterior cornea, the iris and the crystalline lens over a large cross-sectional area. Recent studies [9] have shown the possibility of increasing the range of OCT imaging from the anterior cornea to the posterior surface of the crystalline lens. OCT presents advantages over other techniques such as Purkinje (limited to the estimation of lens alignment, keratometry and phakometry) and Scheimpflug imaging (which has lower resolution, and the posterior lens visibility is limited by the pupil aperture). However, unlike exceptionally in OCT, Purkinje and Scheimpflug imaging have been thoroughly validated and corrected from the optical distortions produced by the preceding ocular surfaces and the accuracy of the results demonstrated using model eyes with known geometry [5],[7],[10–12].

In general, OCT is subject to two types of distortion: fan and optical distortions. Fan (field) distortion is related to the rastering of the surfaces with optical scanners. If the fan of principal rays is not perfectly flat with respect to the optical axis of the system, then the first surface it is not well reproduced and even a flat surface becomes curved on the image. This effect can be characterized and minimized numerically by application of correction algorithms based on ray propagation, Snell law and Fermat principle. We have addressed this effect in a previous publication [13].

Optical distortion is an effect produced while imaging a sample embedded in a medium different from vacuum (air), and therefore the rays are bended by the effect of refraction by preceding surfaces. This produces a very significant loss of the geometrical features of the observed surface, and therefore its geometrical parameters, such as curvature and thickness cannot be retrieved directly from the image. Imaging the inner layers in the eye is therefore subject to this type of distortion, as each layer differs in shape and structure [14] and will show a different index of refraction. Apart from the strongest refractive effect at the air-anterior cornea (or tear film) interface, the major effects in a normal eye will happen in the posterior cornea-aqueous humor, aqueous humor-anterior lens, and posterior lens-vitreous humor interfaces.

Since OCT was introduced in 1991 by Huang and his group [15] and [16], great advances have been produced in the area, and several instruments have become commercially available. Most of the research efforts have been directed to increase speed, allowing real-time data acquisition in living eyes (i.e. Fourier domain [17–19]) and to increase resolution by using new light sources (i.e. femto-second lasers or supercontinuum) achieving axial resolutions of down to 2-5 μm [20–22] or even to the submicron scale [23]. However, while OCT is capable of producing images of very high quality, its use to provide quantitative estimates of the ocular surfaces geometry has been scarce, and mostly limited to axial biometry. Developments in quantitative anterior segment OCT corneal topography are recent [24–26, 13] and primarily limited to anterior corneal surface.

The optical surfaces such as those in the eye are generally more complicated because they do not exhibit rotational symmetry and are not centered around an optical axis. As a consequence, the technique of distortion correction based on 2-D ray tracing seems to be insufficient. Therefore, a complete optical distortion correction requires a 3-D approach. Providing accurate topographies of the internal surfaces of the eye requires correction of the optical distortion, and this problem, to our knowledge, has not been addressed before in 3-D. Only a few studies have addressed the problem of optical distortion, and have proposed methods for 2-D corrections for OCT images [27, 28], based on the Fermat’s principle. While this is an excellent first approach to the problem and undoubtedly better than the standard correction consisting of dividing the OCT optical distances by the refractive index, we will show that the errors of minimizing parameters that define the surfaces are better estimated correcting for 3-D instead of doing it only on cross-sectional images. Podoleanu et al. [27] showed corrections of cross-sectional images of an intraocular lens (IOL), an intralipid drop and the anterior chamber angle in a living human eye. A limitation of the work is the use of analytical fits of circles (according to Gullstrand model) to the surfaces, and therefore ignoring local irregularities of the surfaces. Westphal et al. [28] recognized the potential improvement of applying the optical correction in 3-D, instead of 2-D, since the normal has three components, as well as the refracted cosine director.

In this paper we present a numerical method to correct for the optical distortion of 3-D data from OCT systems. Experimental data on lenses (artificial cornea and IOL) with known geometrical parameters and in real eyes (porcine eye in vitro and human eye in vivo) were collected with two calibrated laboratory OCT set-ups, one time-domain OCT (TD-OCT) system and the other a spectral OCT (sOCT) system. Previously to the optical correction, the raw images were corrected from fan distortion [13]. The images are also segmented using custom-developed algorithms and fitted by a Delaunay surface decomposition [29], which is used to perform 3-D Ray Tracing. The obtained geometrical parameters on the artificial cornea and IOL were validated against a non-contact profilometer, while a porcine and a human cornea were validated against a Scheimpflug imaging system.

2. Methods

2.1. Theoretical considerations and computational implementation of optical distortion correction

Any low-coherence interferometry system records optical paths. The easiest and most frequent method for reconstructing the position of the surfaces and interface distances is just by dividing the recorded optical paths by the value of refractive index without taking into account the refraction effect on the optical surfaces [30, 31]. Although this approach can be justified for the paraxial region and coaxial surfaces, where the axial ray goes through the optical system without being refracted, it cannot be applied for outer (peripheral) regions.

The scanning rays in an OCT system are generally characterized by angular coordinates $(θ, φ)$ of the scanner. Even though the OCT may be calibrated and the angular coordinates can be converted directly to spatial $(x, y)$ coordinates, such calibration is valid only for the first surface, as it is the only one that is reached by parallel beams (assuming the fan distortion can be negligible or has been corrected [13]). For this reason, a simple division of the optical path by the value of refractive index of the media, which is the simplest and the most common way to correct the optical distortion, is insufficient for retrieval of the correct physical distances and dimensions of the sample.

Figure 1 presents an illustration of the sample arm of the OCT measurement configuration, along with a curved sample (of refractive index $n^{1}$ , in blue) immersed in an external medium (in white) of refractive index $n^{0}$ . The elevation and azimuth angles $θ_{i}$ _and $φ_{i}$ _{introduced to the galvanometric scanners define the} i-th ray. Once this ray has been propagated through the lens, the outcoming ray defines a reference point ${\bar{R}}_{i}^{0}$ and a cosine director ${\hat{T}}_{i}^{0}$ , where $L_{i}^{0}$ is the distance to meet the surface at the point ${\bar{R}}_{i}^{1}$ . The normal at the meeting point is denoted by ${\hat{N}}_{i}^{0}$ . Then, the ray is subject to refraction and it changes its directional cosine ${\hat{T}}_{i}^{1}$ . If the ray continues propagating through the second medium (in blue), it meets the second interface surface at point ${\bar{R}}_{i}^{2}$ travelling a geometrical distance denoted by $L_{i}^{1}$ . Therefore every pair of angles of the galvanometric scanner $(θ_{i}, φ_{i})$ can be related to a sampling ray and its equation in a generalized vectorial form:

{\bar{R}}_{i}^{j + 1} = {\bar{R}}_{i}^{j} + L_{i}^{j} {\hat{T}}_{i}^{j} j = 0, ..., m

where the subscript i is the number of the ray, and the superscript jdenotes the number of the interface.

Fig. 1 Illustration of a scanning system plus a collimation-condensing lens in an anterior segment OCT system. See text for details.

Download Full Size | PDF

However, in the OCT image representation the geometrical distances cannot be seen because each single A-Scan represents the optical path of one particular ray. In general the total path $D_{i}$ of a particular ray (i) can be expressed as the summation of the partial optical paths of every media passed by the ray, assuming that the media are homogeneous the optical distance is simply the multiplication of the geometrical path and the value of the refractive index (see Eq. (2).

D_{i}^{j} = \sum_{j = 0}^{m} n^{j} L_{i}^{j}

where

n^{j}

is the refractive index of the medium between j-th and (

j + 1

) interface.

This way any OCT A-Scan can be described by two angular coordinates $(θ_{i}, φ_{i})$ and one spatial coordinate $D_{i}$ .

In order to correct for the effect of optical distortion in general (fan distortion and refraction), it is required to propagate the optical distances obtained by the OCT along the directional cosines of the ray, taking into account refraction effects that happen at each interface between two media of different refractive indices. Starting from the calibrated zero position ${\bar{R}}_{i}^{0}$ and the initial directional cosine ${\hat{T}}_{i}^{0}$ the consecutive coordinates (spatial ${\bar{R}}_{i}^{1}$ and directional ${\hat{T}}_{i}^{1}$ ) can be estimated considering that the position of this point ( ${\bar{R}}_{i}^{1}$ ) is actually a geometrical distance $L_{i}^{0}$ measured along the ray (see Eq. (1).). Once the ray goes through the first surface, it undergoes refraction if the consecutive medium has a refractive index different from air. In order to calculate the direction of ray propagation after refraction, the vectorial refraction equation needs to be applied [36]:

{\hat{T}}_{i}^{j + 1} = (\frac{n^{j}}{n^{j + 1}}) {\hat{T}}_{i}^{j} - {(\frac{n^{j}}{n^{j + 1}}) {\hat{T}}_{i}^{j} {\hat{N}}_{i}^{j} - \sqrt{1 - {(\frac{n^{j}}{n^{j + 1}})}^{2} [1 - {({\hat{T}}_{i}^{j} {\hat{N}}_{i}^{j})}^{2}]}} {\hat{N}}_{i}^{j}

In order to do this, it is required to know the cosine director ${\hat{T}}_{i}^{j}$ of the incoming ray; the refractive index of every media in the ray path ${n_{i + 1}, n_{i}}$ , as well as the normal vector ${\hat{N}}_{i}^{j}$ to the surface at point of refraction.

In order to apply this procedure computationally we have performed the calculation of the 3-D normal vectors by performing the cross products of the 3-D gradients for vertical and horizontal components. The gradients are obtained (Eq. (4) by convolving a gradient mask vector in the horizontal direction ( $M_{h}$ ) and in the vertical direction ( $M_{v}$ ) with the component of the detected points of every surface ${X^{j}, Y^{j}, Z^{j}}$ arranged in matrices. Notice that the arrangement of points corresponding to horizontal, vertical and axial coordinates is denoted in capital letters. The results of the operations are six matrices corresponding to the horizontal, vertical and axial coordinates for the horizontal ${G x_{h}^{j}, G y_{h}^{j}, G z_{h}^{j}}$ and vertical ${G x_{v}^{j}, G y_{v}^{j}, G z_{v}^{j}}$ components of the gradients.

{\begin{matrix} {G x_{h}^{j}, G y_{h}^{j}, G z_{h}^{j}} = M_{h}^{} \otimes {X^{j}, Y^{j}, Z^{j}} = (\begin{matrix} \frac{1}{2} & 0 & - \frac{1}{2} \end{matrix}) \otimes {X^{j}, Y^{j}, Z^{j}} \\ {G x_{v}^{j}, G y_{v}^{j}, G z_{v}^{j}} = M_{v}^{} \otimes {X^{j}, Y^{j}, Z^{j}} = (\begin{matrix} - 1 / 2 \\ 0 \\ 1 / 2 \end{matrix}) \otimes {X^{j}, Y^{j}, Z^{j}} \end{matrix}

In order to avoid the edge effects of filtering, the boundaries of the matrices are expanded by bicubic local extrapolation. This way the gradients for the three components are obtained in matricial form. Once the vertical and horizontal gradients are obtained, the resultant matrices are cropped by their boundaries to prevent them to be affected by filtering operations and to retrieve its original size.

The local gradient of a point (an element in matrix representation) is obtained by arranging again the matrices by the concatenation operation for the particular element of the gradient matrices. Let us to assume that the i-th ray corresponds to the element (m,n) in matricial notation. Then the horizontal and gradients vectors ( ${\vec{G}}_{h}^{j}$ and ${\vec{G}}_{v}^{j}$ ) can be expressed by:

\begin{matrix} {\vec{G}}_{h}^{j} (i) = [G x_{h}^{j} (m, n), G y_{h}^{j} (m, n), G z_{h}^{j} (m, n)] \\ {\vec{G}}_{v}^{j} (i) = [G x_{v}^{j} (m, n), G y_{v}^{j} (m, n), G z_{v}^{j} (m, n)] \end{matrix}

The cross-product of these gradients (

{\vec{G}}_{h}^{j}

and

{\vec{G}}_{v}^{j}

), for each point of these matrices, provides non-normalized normal vectors

{\vec{N}}_{i}^{j}

perpendicular to the surface. The vectors can be normalized as follows:

{\hat{N}}_{i}^{j} = \frac{{\bar{G}}_{h}^{j} (i) \times {\bar{G}}_{v}^{j} (i)}{‖ {\bar{G}}_{h}^{j} (i) \times {\bar{G}}_{v}^{j} (i) ‖}

As explained above, the point which met the surface acts as the reference for the refracted ray $R_{i}^{j + 1}$ , while the new cosine director ${\hat{T}}_{i}^{j + 1}$ is calculated by the refraction formula. In a last step, as the measured distance between layers is an optical path, it needs to be divided by the refractive index to retrieve the real geometrical aspects of the surface. In such a way, the points under the layer can be shifted following the new ray equation. This processing should be repeated for every layer until all the points are corrected. Figure 2 illustrates the steps of the 3-D correction from the original surfaces to the reconstructed surface after correction.

Fig. 2 a) Original segmented surfaces of a human cornea. b) Anterior surface plus incoming rays in blue c) Normals to the surface the surface in red. d) Refracted rays in red and reconstruction of the second surface.

Download Full Size | PDF

In this work, two more methods to correct the optical distortion have been implemented for comparison purposes with the 3-D correction previously described. The most common method consists of dividing the optical path by the refractive index of every layer. The so called 2-D Ray tracing method consists of applying the Snell law to every segmented curve of each collected B-Scan. This way the images are corrected from refraction in the very same plane as they are collected. Finally, the resultant layers or surfaces obtained by the three methods are fitted to 3-D conics in order to get relevant parameters (radius or conic constant) of the surface.

In order to predict the magnitude of the effect of the optical distortion in the eye, we simulated OCT imaging of the anterior chamber of a simple 4-surface computer model (with conic anterior and posterior cornea and lens surface) with realistic amounts of lens tilt and decentration [7,10] and foveal misalignment [12]. The light was propagated across surfaces by exact 3-D Ray Tracing, the points of intersection between the surface and the beams recorded, and the optical distances evaluated. The resultant four surfaces were corrected from optical distortion by three approaches: (1) the most common method of dividing the optical path by the refractive index, (2) by performing a 2-D Ray tracing of every cross section (B-Scan) and (3) the method of 3-D correction of optical distortions.

2.2. Experimental setups and protocols

Experimental optical coherence tomography images were obtained using two laboratory anterior segment OCT systems set up at the Visual Optics and Biophotonics Lab (Instituto de Optica, CSIC, Madrid, Spain). One is a classical time domain OCT (TD-OCT) system with a confocal channel, described in detail before [13]. The other system is a high-speed spectral OCT (sOCT) using custom designed software, optics and spectrometer with CMOS camera (Basler Sprint), constructed in collaboration with Nicolaus Copernicus University, Torun, Poland based on commercially available Copernius OCT device (Poland) [9]. The parameters of light sources and resolutions of both systems can be found in our earlier publications.

These two devices were used to collect the 3-D data from different samples. In vitro samples (imaged with the time-domain OCT system) included (1) IOL made of PMMA, with spherical surfaces of known radii of curvature (as provided by the manufacturer, AJL Ophthalmics, Spain, and measured using non-contact profilometry); (2) a custom-made contact lens made of PMMA with spherical anterior and posterior surfaces, used in previous studies as an artificial cornea in physical model eyes [7]; (3) a freshly enucleated porcine cornea from a local slaughterhouse. In vivo measurements were performed on one normal human eye with the spectral OCT system.

The correction was tested experimentally on a PMMA spherical biconvex IOL and a PMMA spherical convex-concave cornea. The radii of curvature of all the surfaces were measured with a microscopy-based non-contact profilometer (Sensofar Plμ 2300) which allows surface topographies with an axial resolution of 0.15 μm using a 20X microscope objective. The OCT measurement was done using the Time-domain OCT, optimized for fan distortion minimization [13]. The central wavelength of this OCT is 819 nm. Measurements were taken on a grid of 30x30 scanning angles with an axial depth of 4 mm (and therefore the 3-D set consisted of 30 B-scans, each of 30 A-Scans). The samples were measured in air. The confocal channel was used to align the sample with respect to the optical axis of the sample arm. The 3-D data cloud after denoising, detection and segmentation algorithms was fitted by a spherical using a least–mean-squares algorithm. The index of refraction used in the calculations was 1.4856 for PMMA.

The correction was also tested experimentally on a porcine eye. As the TD-OCT is a slow system the measurements must be done on in vitro samples. Besides, the introduction of in vitro measurements allows comparing OCT and Scheimpflug at the same conditions without the movements inherent to living tissue. The sample was placed in an eye-holder without any additional pressure control. The experimental procedure consisted of a first measurement with a commercial Scheimpflug imaging rotating camera (Pentacam, Oculus), using 25 cross-sections and its built-in software for 3-D reconstruction, and immediately a 3-D measurement using the TD-OCT, with similar procedures to those explained before. For this experiment, the grid size was distributed regularly in 55 A-Scans x 55 B-scans, which resulted in a tested area of 7.35 x 7.35 mm. Fan distortion correction, and denoising, segmentation, Delaunay decomposition used for ray tracing and fitting to a sphere were performed as described above, prior to optical distortion correction. The corneal refractive index used in the reconstruction was 1.376.

Additionally, 3-D measurements on the left eye of one normal subject were obtained using the laboratory-based spectral OCT. The participating subject was 28 years of age. The study was approved by Institutional Review Boards and followed the tenets of the Declaration of Helsinki. The subject signed a consent form and was aware of the nature of the study. We used 100 horizontal B-scans equally spaced along 10 mm. Each B-Scan was (450x2048 pixels) with an averaging of the 3 rows, resulting in 150x2048 pixels across a 15 x 7 mm region, with a total measurement time of less than 2 s and 800 µW corneal power exposure at a central wavelength of 840 nm. Surface elevations of the anterior and posterior cornea were fitted by conical surfaces on an optical zone of 9-mm diameter. The image processing procedures were performed as in previous examples with lenses and in vitro porcine eye. As in previous experiments, the data were analyzed using the conventional processing (division by the refractive index of the cornea), and after application of optical distortion correction. The index of refraction used for the cornea was 1.376 and 1.336 for the aqueous humor. The same subject was measured with the Pentacam Scheimpflug imaging system, using conventional clinical protocols (corneal reflex centration, 50 meridians, chin and forehead rest stabilization).

2.3 Signal processing

The 3-D images consist of a stack of images (B-Scans) being horizontal cross-sections of the object at different vertical positions. Custom algorithms were developed to reduce noise and improve the detectability of the edge detection algorithms. The routines were written in Matlab® and developed to meet the particular detection architecture of the OCT system (a linear CMOS camera in the sOCT set up and a photodetector in a balanced configuration in the time-domain OCT set-up). Some studies in the literature analyze the sources of noise in the images [32–34]. If we assume that the noise is random, it can be characterized by its bandwidth, which is directly related to the amount of energy that reaches the detector. In OCT, the energy coming from the sample arm changes for every A-Scan, so that the data processing should be performed for every A-Scan. In the first denoising step a rotational kernel transform (RKT) was used [35]. Since the denoising processing in OCT images is one dimensional only, the actual vertical kernel was used. Alternatively, a wavelet low-pass filtering processing (from Matlab Wavelet Toolbox catalogue) was used, obtaining similar detection performance. The images were further improved using a contrast enhancement algorithm consisting of the average of a number of points along the A-Scan, which increased the signal in the peaks and averaged noise, which was then subtracted as an offset. Figure 3(a) shows an original image from an artificial eye as captured by the TD-OCT. Figure 3(b) shows the result after applying the denoising and contrast enhancement algorithms. The structural information was distinguished from noise using a simple statistical algorithm. The mean and standard deviation of every A-Scan was estimated within the noise (non-image) region, and then all points above the mean and the standard deviation (multiplied by a value of the threshold) were saved. The number of significant points is further reduced by selecting those belonging to the actual surfaces of the optical components, by finding the separation between consecutive peaks, as we assume that the maximum of the autocorrelation distribution represents the surface where a scattering event occurs. The accuracy of the detection depends on the stability of superluminescent diode (SLED), on the mechanical stability of the components, and the resolution of the acquisition system. The procedures were applied to every single A-Scan, to create a mesh of points representing the signal from the imaged structures. The classification of peaks and the detection of the maximum value in each A-Scan allow identifying the structures that maximally reflect the light in each A-Scan. Apart from reducing the number of points in the mesh (by a factor of 5-20 per A-Scan) it also allows defining a set of consecutive of layers in the sample.

Fig. 3 a) Cross section of an artificial eye composed by a plastic cornea of PMMA and an IOL obtained by the TD-OCT. b) Resultant image after denoising and contrast enhancement (the additional, lower contrast lines are ghost images and easily removed in the segmentation process).

Download Full Size | PDF

Additionally, spurious points were eliminated using a neighborhood algorithm, which allows identifying a layer if there is continuity in at least one of the directions of the tensor mask (which is 3-D convolved with the tensor of data). A tensor data of the size $(m + 2) \times (n + 2) \times (p + 2)$ was used, where m, n and p are dimensions of the data tensor corresponding to the number of images (m), number of pixels in the horizontal direction (n) and in the vertical direction of the image (p) and a tensor mask $M_{T}$ of size 3x3x3 Eq. (7).

M_{T} = [\begin{matrix} 1 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{matrix}] \times [\begin{matrix} 1 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{matrix}] \times [\begin{matrix} 1 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{matrix}]

The convolution provides a set of numbers between 0 and 24 indicating the degree of points detected in the mask neighborhood. The optimal threshold was determined and the algorithm applied to segment the different layers 3-dimensionally (Fig. 4 ).

Fig. 4 a) Movie (Media 3) with the data without filtering with the 3-D neighborhood algorithm. b) Movie (Media 4) showing the data after spurious points removal. c) Movie (Media 5) showing the points and the segmented layers.

Download Full Size | PDF

Once the 3-D cloud of points is obtained, the data need to be corrected for fan distortion. The characterization of the fan distortion is important not only to correct the position of the points, but also to characterize the directional cosine of the propagation vector for every position of the scanner mirrors. Otherwise, the correction of the optical distortion must be done as if the ray that reaches the surface was collimated, resulting in improper propagation vector and therefore in a wrong location of the optically distorted point in space. An algorithm for fan distortion correction proposed in a previous study [13] to retrieve the anterior surface quantitative topographical data was applied prior to the optical distortion correction algorithm developed in the current study.

3. Results

3.1. Computer simulations

We simulated the effects of optical distortion on OCT anterior segment images and studied the errors introduced when applying the most commonly used methods for processing the images: (1) simple division of the optical path by the index of refraction for each A-Scan; (2) 2-D correction of the optical distortion only; (3) 3-D correction of optical distortions. In all cases the fan distortion correction, arising from the scanning architecture of the OCT system was applied [13]. The parameters of the elements of the OCT systems used in the experiment (galvanometer scanner and achromatic lens) were introduced in the simulation. The optical system of the model eye consisted of four surfaces simulating cornea and lens. The nominal geometrical features for the ocular surfaces (radii, conic constant and refractive index) of the computer eye model used in the simulations are summarized in Table 1 . The data for refractive indices of the cornea, lens and ocular media are taken from the literature [37–40], where no gradient index but an equivalent index is assumed for the crystalline lens. The simulation was performed for a square grid of 100 x 100 pixels with a 6-mm side.

Table 1. Nominal data for computer eye model used in computer simulations

View Table | View all tables in this article

Table 2 shows the results of radii of curvature of the lens surfaces and the thickness from simulated OCT surfaces by exact 3-D Ray Tracing with the values of the Table 1. The data were obtained from the simulated OCT images and reconstructed using the standard processing (division by refractive index), correction of 2-D data, and correction of 3-D data.

Table 2. Retrieved data from the simulated OCT images, with three different type of data processing

View Table | View all tables in this article

For the posterior corneal radius of curvature, the difference between the nominal value and the retrieved value correcting by the refractive only is 3.7% (equivalent to 0.23 diopters, D). Applying the optical correction in 2-D the difference is smaller 1.4% (equivalent to 0.09 D). The difference between the nominal value and that obtained after 3-D correction is negligible. For the anterior lens radius of curvature, the difference between the nominal value and the refractive index correction is 8.9% and 3.7% applying the 2-D correction (equivalent to 0.38 and 0.17 D respectively). For the posterior lens surface radius of curvature the discrepancies are 59.7% and 12.2% (equivalent to 3.4 and 1.0 D respectively). In both cases, the difference between the nominal radius of curvature and the value obtained after application of the 3-D correction is lower than 0.3%. For the conic constant the discrepancies with respect to the nominal value when applying only the refractive index correction range from 16.9% for the posterior cornea to 207% for the posterior lens surface. Applying the 3-D correction produces discrepancies of less than 1% in all surfaces.

3.2 Experimental data

In the next stage of our study we acquired images of an IOL and artificial cornea of known parameters. Table 3 shows the measurements of the surfaces radii of curvature on the IOL and on the PMMA cornea (reference values from non-contact profilometry and estimates from the TD- OCT). Figure 5 illustrates the effect of the correction. The anterior surface radius of the IOL from OCT shows a difference of 1.0% from the reference value. The posterior radius of curvature of the IOL from TD-OCT (with no correction) differs by 66.7% from the reference value, while after 3-D correction the difference is only 1%. For the PMMA the anterior cornea radius of curvature differs 1.0% from the value of reference. Non-contact profilometry of the posterior PMMA corneal surface was not due to the limitations of the profilometer in measuring convex surfaces. Nominal values were available from the manufacturer, and those differed less than 1.0% from the estimates from TD-OCT after 3-D correction.

Table 3. IOL and artificial cornea dimensions: values from manufacturer, and estimates from non-contact profilometry and TD-OCT measurements

View Table | View all tables in this article

Fig. 5 a) 2-D representation of all cross sections obtained from the IOL using the TD-OCT. b) 3-D Representation of the surfaces extracted from the IOL, the anterior and the corrected posterior surfaces are wire represented while the raw data without correction is represented by the red discrete points.

Download Full Size | PDF

In the second stage, experiments were performed on a porcine eye. Table 4 shows results of its anterior and posterior corneal curvature. Data from the Pentacam Scheimpflug imaging (using the commercial software), as well as the values obtained after standard correction (division by the index of refraction), and the optical distortion correction of TD-OCT images are included. Figure 6 shows a cross-section of the porcine eye as obtained from the OCT, after segmentation, and the corrected 3-D result. The agreement between measurements on the anterior cornea obtained with Scheimpflug and OCT confirms the appropriate fan distortion correction of the OCT data. Correction of the optical distortion on the posterior corneal changes the estimate by 6% (the difference is equivalent to 0.39 D) with respect to estimates obtained only dividing by the refractive index. Differences in the measured posterior corneal surface radius of curvature between Scheimpflug and corrected OCT are about 1.0%.

Table 4. Anterior and posterior corneal dimensions of an enucleated porcine cornea: estimates from Scheimpflug imaging and TD-OCT measurements

View Table | View all tables in this article

Fig. 6 a) 2-D representation of a cross-section obtained from the porcine cornea using the TD-OCT. b) The same cross section corrected using a 2-D ray tracing algorithm. c) 3-D representation of the anterior and posterior surfaces extracted from the porcine cornea. The raw data are represented by a mesh and the corrected posterior surface by solid color.

Download Full Size | PDF

Finally, the correction algorithms were tested on the in vivo images of the human eye. In the case of in vivo OCT imaging there is a good correspondence between Scheimpflug and OCT anterior corneal radii of curvature (after fan distortion correction): (7.46 and 7.66 mm for the radius of curvature, and 0.42 and 0.54 for the asphericity). Table 5 shows results for the posterior cornea of a human eye in vivo (data from Scheimpflug imaging, as well as from sOCT after 2-D optical distortion correction). The difference in the radius of curvature from Scheimpflug and 3-D corrected sOCT is 2.8%, while the correction of the optical distortion changes the estimate of the posterior radius by 6.2% (with respect to data only divided by the refractive index), equivalent to 0.37 D. The largest differences are observed in the conic constant, where the application of the 3-D correction distortion algorithm modifies the estimate by 73% with respect to the simplest corrections (also increasing the agreement with respect to the Scheimpflug estimate). Figure 7 illustrates the differences of correcting by dividing by the refractive index and the 3-D distortion correction in the living human eye. Note for example the differences in the iris diameter between the refractive-index corrected image (a), and that of the image corrected from optical distortion in 3-D (7.5% smaller).

Table 5. Posterior cornea dimensions of a human eye in vivo: estimates from Scheimpflug imaging and sOCT.

View Table | View all tables in this article

Fig. 7 a) Movie (Media 1) of the 3-D cloud of points by correcting by the refractive index, b) (Media 2) Correcting in 3-D.

Download Full Size | PDF

4. Conclusions

We have proposed a method for 3-dimensional correction of optical distortion in ocular anterior segment optical coherence tomography images, and demonstrated that after application of fan and optical distortion correction, reliable quantitative reconstruction of the ocular surface internal topography can be achieved.

The method has been validated using simulations of OCT imaging in a realistic computer eye model, in vitro on intraocular lenses and a PMMA artificial cornea of known parameters (from profilometry), and in real eyes (a porcine eye in vitro and a human eye in vivo) – where comparisons have been established with Scheimpflug imaging. Application of 3-D optical distortion correction increases the accuracy of the radii of curvature estimates by significant amounts, with respect to the simple division of the optical path by the index of refraction and over 2-D correction. Differences between the retrieved 3-D corrected and nominal values are less than 1% in most cases. The correction is particularly important for the surface asphericities, where the simulations predict discrepancies as large as 200% for the posterior lens, and large differences with respect to applying the conventional correction have been also found experimentally in the posterior cornea. Experiments have been performed using two different laboratory-developed OCT systems (time-domain and spectral).

Correcting for optical distortion (along with for fan distortion) is therefore critical to reliably use OCT for quantitative estimates of ocular surfaces geometry, particularly that of internal surfaces. This correction is important in numerous applications that rely on quantitative analysis of the ocular surfaces (corneal changes after refractive surgery, cataract surgery and IOL implantation, corneal biomechanics, etc…) as well as those that involve quantitative measurements in the crystalline lens. While the current study has primarily addressed the cornea (in the examples in real eyes) the technique would be applied similarly (and subsequently) for the reconstruction of the crystalline lens surfaces. An additional complication (with respect to the examples of artificial eyes in vitro presented here) is the presence of a gradient index distribution in the crystalline lens.

Two previous studies in the literature have discussed the presence of optical distortion and addressed it in 2-dimensions. We have demonstrated higher accuracies when the correction is applied in 3-D. The importance of the correction has been demonstrated using laboratory-developed instruments, but it is inherent to the technique, and therefore has potential of being applicable to standard OCT commercial devices.

Acknowledgements

The authors acknowledge funding from MICINNFIS2008-02065, PETRIPET-2006-0478 (Spanish Ministry of Science and Innovation) to SM; EURYI-05-102-ES (EURHORCs-ESF), to SM; and EURYI-01/2008-PL to MW; CSIC JAE Program to DS and LR; Alcon, Spain Scholarship to SO. We would like to acknowledge support of N. Chía, C. Teijeiro and P. Pérez form the Visual Optics and Biophotonics group; and D. Szlag, M. Szkulmowski and M. Gora from the NCU Medical Physics Group.

References and links

1. P. Rosales and S. Marcos, “Customized computer models of eyes with intraocular lenses,” Opt. Express 15(5), 2204–2218 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-5-2204. [CrossRef] [PubMed]

2. P. Rosales, M. Dubbelman, S. Marcos, and G.L. van der Heijde, “Crystalline lens radii of curvature from Purkinje and Scheimpflug imaging,” J. Vis. 6(10), 1057–1067 (2006). [CrossRef] [PubMed]

3. M. Dubbelman, G. L. Van der Heijde, and H. A. Weeber, “Change in shape of the aging human crystalline lens with accommodation,” Vision Res. 45(1), 117–132 (2005). [CrossRef]

4. E. Moreno-Barriuso, J. M. Lloves, S. Marcos, R. Navarro, L. Llorente, and S. Barbero, “Ocular aberrations before and after myopic corneal refractive surgery: LASIK-induced changes measured with Laser Ray Tracing,” Invest. Ophthalmol. Vis. Sci. 42(6), 1396–1403 (2001). [PubMed]

5. A. Pérez-Escudero, C. Dorronsoro, L. Sawides, L. Remón, J. Merayo-Lloves, and S. Marcos, “Minor Influence of Myopic Laser in Situ Keratomileusis on the Posterior Corneal Surface,” Invest. Ophthalmol. Vis. Sci. 50(9), 4146–4154 (2009). [CrossRef] [PubMed]

6. S. Norrby, “Sources of error in intraocular lens power calculation,” J. Cataract Refract. Surg. 34(3), 368–376 (2008). [CrossRef] [PubMed]

7. A. de Castro, P. Rosales, and S. Marcos, “Tilt and decentration of intraocular lenses in vivo from Purkinje and Scheimpflug imaging. Validation study,” J. Cataract Refract. Surg. 33(3), 418–429 (2007). [CrossRef] [PubMed]

8. S. Marcos, P. Rosales, L. Llorente, and I. Jiménez-Alfaro, “Change in corneal aberrations after cataract surgery with two types of aspherical intraocular lenses,” J. Cataract Refract. Surg. 33(2), 217–226 (2007). [CrossRef] [PubMed]

9. 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–4858 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-6-4842. [CrossRef] [PubMed]

10. P. Rosales and S. Marcos, “Phakometry and lens tilt and decentration using a custom-developed Purkinje imaging apparatus: validation and measurements,” J. Opt. Soc. Am. A 23(3), 509–520 (2006). [CrossRef]

11. M. Dubbelman, H. A. Weeber, R. G. van der Heijde, and H. J. Völker-Dieben, “Radius and asphericity of the posterior corneal surface determined by corrected Scheimpflug photography,” Acta Ophthalmol. Scand. 80(4), 379–383 (2002). [CrossRef] [PubMed]

12. P. Rosales and S. Marcos, “Pentacam Scheimpflug Quantative Imaging of the crystalline lens and intraocular lens,” J. Refract. Surg. 25, 421–428 (2009). [CrossRef]

13. S. Ortiz, D. Siedlecki, L. Remon, and S. Marcos, “Optical Coherence Tomography for quantitative surface topography,” Appl. Opt. 48(35), 6708–6715 (2009). [CrossRef] [PubMed]

14. F. A. Jakobiec, Ocular anatomy, embryology, and teratology (Harper & Row, Philadelphia, 1982).

15. 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), 1178–1181 (1991). [CrossRef] [PubMed]

16. C. K. Hitzenberger, W. Drexler, and A. F. Fercher, “Measurement of corneal thickness by laser Doppler interferometry,” Invest. Ophthalmol. Vis. Sci. 33(1), 98–103 (1992). [PubMed]

17. M. Choma, M. Sarunic, C. Yang, and J. Izatt, “Sensitivity advantage of swept source and Fourier domain optical coherence tomography,” Opt. Express 11(18), 2183–2189 (2003), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-11-18-2183. [CrossRef] [PubMed]

18. R. Leitgeb, C. K. Hitzenberger, and A. F. Fercher, “Performance of Fourier domain vs. time domain optical coherence tomography,” Opt. Express 11(8), 889–894 (2003). [CrossRef] [PubMed]

19. M. Wojtkowski, T. Bajraszewski, P. Targowski, and A. Kowalczyk, “Real-time in vivo imaging by high-speed spectral optical coherence tomography,” Opt. Lett. 28(19), 1745–1747 (2003). [CrossRef] [PubMed]

20. K. Bizheva, B. Povazay, B. Hermann, H. Sattmann, W. Drexler, M. Mei, R. Holzwarth, T. Hoelzenbein, V. Wacheck, and H. Pehamberger, “Compact, broad-bandwidth fiber laser for sub-2-microm axial resolution optical coherence tomography in the 1300-nm wavelength region,” Opt. Lett. 28(9), 707–709 (2003). [CrossRef] [PubMed]

21. A. Unterhuber, B. Považay, B. Hermann, H. Sattmann, W. Drexler, V. Yakovlev, G. Tempea, C. Schubert, E. M. Anger, P. K. Ahnelt, M. Stur, J. E. Morgan, A. Cowey, G. Jung, T. Le, and A. Stingl, “Compact, low-cost Ti:Al2O3 laser for in vivo ultrahigh-resolution optical coherence tomography,” Opt. Lett. 28(11), 905–907 (2003). [CrossRef] [PubMed]

22. H. Lim, Y. Jiang, Y. Wang, Y.-C. Huang, Z. Chen, and F. W. Wise, “Ultrahigh-resolution optical coherence tomography with a fiber laser source at 1 µm,” Opt. Lett. 30(10), 1171–1173 (2005). [CrossRef] [PubMed]

23. B. Povazay, K. Bizheva, A. Unterhuber, B. Hermann, H. Sattmann, A. F. Fercher, W. Drexler, A. Apolonski, W. J. Wadsworth, J. C. Knight, P. St. J. Russell, M. Vetterlein, and E. Scherzer, “Submicrometer axial resolution optical coherence tomography,” Opt. Lett. 27(20), 1800–1802 (2002). [CrossRef]

24. T. Simpson and D. Fonn, “Optical Coherence Tomography of the Anterior Segment,” Ocul. Surf. 6(3), 117–127 (2008). [PubMed]

25. 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(8), 1179–1185 (2001). [PubMed]

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–14894 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-17-14880. [CrossRef] [PubMed]

27. A. Podoleanu, I. Charalambous, L. Plesea, A. Dogariu, and R. Rosen, “Correction of distortions in optical coherence tomography imaging of the eye,” Phys. Med. Biol. 49(7), 1277–1294 (2004). [CrossRef] [PubMed]

28. V. Westphal, A. M. Rollins, S. Radhakrishnan, and J. A. Izatt, “Correction of geometric and refractive image distortions in optical coherence tomography applying Fermat’s principle,” Opt. Express 10(9), 397–404 (2002), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-10-9-397. [PubMed]

29. S. Ortiz, D. Siedlecki, L. Remon, and S. Marcos, “Three-dimensional ray tracing on Delaunay-based reconstructed surfaces,” Appl. Opt. 48(20), 3886–3893 (2009). [CrossRef] [PubMed]

30. X. Wang, C. Zhang, L. Zhang, L. Xue, and J. Tian, “Simultaneous refractive index and thickness measurements of bio tissue by optical coherence tomography,” J. Biomed. Opt. 7(4), 628–632 (2002). [CrossRef] [PubMed]

31. E. Kim, K. Ehrmann, S. Uhlhorn, D. Borja, and J.-M. Parel, “Automated Analysis of OCT Images of the Crystalline Lens,” in: Ophthalmic Technologies XIX, edited by F. Manns, P. G. Söderberg, A. Ho Proc. of SPIE Vol. 7163, 716313 (2009).

32. I. Takada, “Noise in Optical Low-Coherence Reflectrometry,” IEEE J. Quantum Electron. 34(7), 1098–1108 (1998). [CrossRef]

33. A. G. Podoleanu, “Unbalanced versus Balanced Operation in an Optical Coherence Tomography System,” Appl. Opt. 39(1), 173–182 (2000). [CrossRef]

34. R. C. Haskell, D. Liao, A. E. Pivonka, T. L. Bell, B. R. Haberle, B. M. Hoeling, and D. C. Petersen, “Role of beat noise in limiting the sensitivity of optical coherence tomography,” J. Opt. Soc. Am. A 23(11), 2747–2755 (2006). [CrossRef]

35. J. Rogowska and M. E. Brezinski, “Image processing techniques for noise removal, enhancement and segmentation of cartilage OCT images,” Phys. Med. Biol. 47(4), 641–655 (2002). [CrossRef] [PubMed]

36. J. B. Keller and H. B. Keller, “Determination of reflected and transmitted fields by geometrical optics,” J. Opt. Soc. Am. 40(1), 48–52 (1950). [CrossRef]

37. J. E. Weddell, J. A. Alvarado, and M. J. Hogan, Histology of the human eye (W.B. Saunders and Co, 1971).

38. M. J. Stafford, The histology and biology of the lens (Bausch & Lomb, 2001).

39. G. Smith, “The optical properties of the crystalline lens and their significance,” Clin. Exp. Optom. 86(1), 3–18 (2003). [CrossRef] [PubMed]

40. B. A. Moffat, D. A. Atchison, and J. M. Pope, “Age-related changes in refractive index distribution and power of the human lens as measured by magnetic resonance micro-imaging in vitro,” Vision Res. 42(13), 1683–1693 (2002). [CrossRef] [PubMed]

	Radius [mm]	Conic Constant		Refractive Index	Thickness [mm]
Anterior Cornea	7.87	0.82	Air	1	-
Posterior Cornea	6.4	0.9	Cornea	1.376	0.56
Anterior Lens	11.76	−4	Aqueous humor	1.336	1.56
Posterior Lens	−5.96	−3	Lens	1.399	5.6

	Refractive Index Correction	2-D optical distortion correction	3-D optical distortion correction
Cornea
Anterior Radius [mm]	7.88	7.88	7.88
Posterior Radius [mm]	6.64	6.49	6.39
Anterior conic constant	0.83	0.83	0.83
Posterior conic constant	0.77	0.93	0.91
Lens
Anterior Radius [mm]	12.80	12.20	11.77
Posterior Radius [mm]	9.55	6.69	5.96
Anterior conic constant	3.22	2.98	3.96
Posterior conic constant	1.45	1.09	2.99

	Nominal value	Non Contact profilometry[mm]	OCT without optical distortion correction [mm]	OCT with optical distortion correction [mm]
Intraocular lens
Anterior[mm]	17.29	17.66	17.82	17.82
Posterior[mm]	−11.04	−10.70	−6.42	−10.82
Thickness [mm]	0.93	0.89	1.37	0.92
PMMA Cornea
Anterior[mm]	7.8	7.72	7.77	7.77
Posterior[mm]	6.48	-	10.29	6.35
Thickness [mm]	0.55	-	0.87	0.6

	Scheimpflug imaging	OCT refractive index corrected	OCT optical distortion corrected
Anterior[mm]	8.62	8.63	8.63
Posterior[mm]	7.37	7.89	7.33
Thickness [mm]	0.87	0.90	0.90

	Scheimpflug	OCT Refractive Index Corrected	OCT 2-D distortion corrected	OCT 3-D distortion corrected
Radius[mm]	6.14	6.70	6.28	6.31
Conic Constant	−0.40	−0.28	−0.26	−0.45
Thickness [mm]	0.59	0.59	0.59	0.58

Optical distortion correction in Optical Coherence Tomography for quantitative ocular anterior segment by three-dimensional imaging

Abstract

1. Introduction

2. Methods

2.1. Theoretical considerations and computational implementation of optical distortion correction

2.2. Experimental setups and protocols

2.3 Signal processing

3. Results

3.1. Computer simulations

3.2 Experimental data

4. Conclusions

Acknowledgements

References and links

Supplementary Material (5)

Cited By

Figures (7)

Tables (5)

Equations (7)

Optics Express