Improved lateral resolution in optical coherence tomography by digital focusing using two-dimensional numerical diffraction method

: This paper proposes a non-iterative, two-dimensional numerical method to alleviate the compromise between the lateral resolution and wide depth measurement range in optical coherence tomography (OCT). A two-dimensional scalar diffraction model was developed to simulate the wave propagation process from out-of-focus scatterers within the short coherence gate of the OCT system. High-resolution details can be recovered from outside the depth-of-field region with minimum loss of lateral resolution. Experiments were performed to demonstrate the effectiveness of the proposed method.


Introduction
Optical coherence tomography (OCT) [1] has recently been developed for diverse areas of medical imaging including in vivo imaging of human retina, skin, and internal body tissues. In standard time-domain OCT (TDOCT), a Michelson-type interferometer is illuminated by a femtosecond laser or superluminescent LED, and an optical-delay-line is used for depth scanning. However, this embodiment of OCT needs specific optical and mechanical designs for scanning and is mainly limited by its relatively slow imaging speed. An alternative approach records the interferometric signal between light from a reference and the backscattered light from the sample point in frequency domain rather than time domain. This is referred to as Fourier domain OCT (FDOCT) [2][3] or spectral domain OCT. The spectral information discrimination in FDOCT is accomplished either by using a dispersive spectrometer [2][3][4][5][6] in the detection arm or rapidly scanning a swept laser source [7][8][9][10][11]. FDOCT has attracted more attention recently because of its high sensitivity and imaging speed compared to TDOCT [12][13][14].
As a coherent cross-sectional imaging technique, OCT is capable of penetrating ~3 millimeters into highly scattering biological tissues, and axial resolution of a few μ m is provided by the low coherence nature of the light source. However, in OCT, a high lateral resolution and wide depth of field (DOF) are always exclusive. Although a higher effective numerical aperture enhances lateral resolution, it narrows the DOF which is inversely proportional to the square of the effective numerical aperture (NA 2 ) of the optical system. Only a very small range around the DOF will exhibit the desired lateral resolution of the system, and the OCT image in the out-of-focus range is blurred laterally. A typical DOF for a small NA system will be several hundred microns, which is about 10 orders smaller than the scanning depth range of an OCT system. Adaptive optics [15], axicon lens [16], dynamic focusing or focus tracking [17][18] is used for maintaining high lateral resolution over a large imaging depth. However, each technique requires special hardware in the system design and may limit the scanning speed and its application in real-time. Inverse scattering [19] and deconvolution algorithms [20][21][22][23][24] have also been designed to improve lateral resolution. Of these methods, most are nonlinear and do not use the phase information of the OCT image. Others [19,24] have taken account of the phase information; however, simulations or experimental studies were mainly limited to one lateral dimension. Inverse scattering has been recently extended for two-dimensional studies [25,26].
In this paper, a novel non-iterative, two-dimensional numerical method for lateral resolution improvement is proposed based on a two-dimensional scalar diffraction model, and the compromise between the lateral resolution and depth measurement range in OCT is alleviated. Numerical diffraction algorithms can be used to calculate back the original detail of the sample without moving the focal plane, thus numerically canceling the lateral defocus. High-resolution details are recovered from outside of the depth-of-field region with minimal loss of lateral resolution. The diffraction model is relatively simple to implement because it only involves propagation of wave between different planes.

Principle
We start by briefly reviewing the illumination and detection of an OCT system. After passing through an objective of the OCT system, a Gaussian probe beam is focused onto a sample, as shown in Fig. 1(a). In an ideal situation, the sample is placed within the DOF of the probe beam. However, in some cases, if the sample is located outside the DOF, the probe beam expands to illuminate a larger area of the sample. Because of the confocal detection scheme of the OCT system, we can consider that the detector is placed in the focal plane of the probe beam and OCT detects the backscattered light, as shown in Fig. 1(b). Note that scatterers at different depths can be differentiated by the short coherence gate of the OCT system. Since the expanded probe beam is illuminating a larger area of the sample, multiple outof-focus scatterers from the same depth will contribute to a single detection (along one Aline). The same scatterer will contribute to different detections as the probe beam is scanning different A-lines across the sample. Interestingly, this is exactly the physical model of diffraction. Optically, it can be considered that the detector (at different lateral positions of the focal plane) records a two-dimensional diffraction pattern from the sample plane, which is located within the short coherence gate. Hence, it is reasonable to use the scalar diffraction model to simulate the wave propagation process for direct scattering [27] from these out-offocus scatterers. This is based on the assumption that only ballistic and quasi-ballistic backscattered photons contribute to OCT imaging because of coherence gating, and contribution from multiple scattering will tend to be delayed and thus may fall outside the coherence gate.
From Fourier optics [28,29], if ( , ;0) E x y is the en-face wave field distribution at plane 0 z = , the corresponding angular spectrum of the field at this plane can be obtained by taking the Fourier transform: where x k and y k are corresponding spatial frequencies of x and y. The field ( , ;0) E x y can be rewritten as the inverse Fourier transform of its angular spectrum, The complex-exponential function exp[ ( )] x y i k x k y + may be regarded as a projection, onto the plane 0 z = , of a wave propagating with a wave vector ( , , ) x y z k k k , where . Thus, the field ( , ;0) E x y can be viewed as a superposition of many wave components propagating in different directions in space and with complex amplitude of each component equal to ( , ;0) x y S k k . After propagating along the z axis to a new plane, the new angular spectrum, ( , ; ) x y S k k z , at plane z can be calculated from Thus, the complex field distribution of any plane perpendicular to the propagating z axis can be calculated from Fourier theory as The lateral resolution of the system is determined by the effective numerical aperture of the objective as in Fig. 1. It is required that the actual lateral sampling step (or scanning step) should be smaller than the theoretical lateral resolution of the system. The scanning step functions as the pixel size of the focal plane, and the pixel size of the reconstructed field from the angular spectrum method is always the same as that of the focal plane. Thus, the proposed method could provide spatially invariant lateral resolution. The schematic of the FDOCT system is shown in Fig. 2. Low-coherence light having a 1310 nm center wavelength with a full width at half maximum of 95 nm was coupled into the source arm of a fiber-based Michelson interferometer. Back-reflected lights from the reference and sample arms were guided into a spectrometer. The dispersed spectrum by a diffraction grating (500g/mm) was sampled by the spectrometer with a 1×1024 InGaAs detector array (SU1024-1.7T, Sensors Unlimited) at 7.7 kHz. The wavelength range on the array was 130 nm, corresponding to a spectral resolution of 0.13 nm and an imaging depth of 3.6 mm in air. Figure 2 also shows the flow diagram of the whole process for lateral resolution improvement based on an FDOCT system. After the captured spectral interferogram in a linescan camera is rescaled as evenly k-spaced, a Fast Fourier transform is followed to achieve the complex A-line information along the z axis. The whole three-dimensional (3-D) complex volume of the sample is obtained by two-dimensional scanning of a galvo-system. Since the sample is located outside the DOF range of the probe beam, the OCT system suffers greatly from lateral resolution degradation. Note here that the 3-D out-of-focus complex volume needs to be re-sampled in the z direction to achieve a sequence of en-face (x-y) images ( , ; ) i I x y z . Then the angular spectrum ( , ; ) x y i S k k z of each en-face image was calculated by Eq. (1), and the new angular spectrum after propagating a distance z was calculated by multiplying ( , ; )  To better understand the process, Figure 3(a) shows one evenly k-spaced spectral interferogram after interpolation when the galvo-system is scanning a point ( , ) m n x y on the sample surface. It is then inverse Fourier transformed and filtered to extract a specific layer ( , ; ) m n I x y k whose amplitude is proportional to the spectral density of the light source, as shown in Fig. 3(d). Here, we have assumed that the spectral interferogram is interpolated with several times of 1024 points so that each layer may cross several pixels in Fig. 3(c). By two-dimensionally scanning the galvo-system and repeating the process, both the full-field layer information ( , ; ) μ m underneath the onion surface. Since the onion was placed outside the DOF, the probe Gaussian beam was expanded to illuminate a larger area and multiple scatterers in the lateral directions contributed to a single A-line detection. This resulted in a great degradation of the lateral resolution, as is clearly shown in Figs. 4(c) and 4(d). Figure 4(e) shows the digitally focused B-scan images at the same cross-section positions as in Fig. 4(c), and the corrected en-face image of Fig. 4(d) is shown in Fig. 4(f) where the scatterers are now greatly sharpened to small bright points. Note that only those scatterers whose centers are located exactly on the observed B-scan image will shrink to bright points, but those scatterers centered on neighboring B-scan images tend to be eliminated or removed. This clearly demonstrates the advantage of a two-dimensional algorithm over one-dimensional algorithms.
In order to laterally differentiate more details of onion cells, an objective of 10 mm focal length (20×, Nachet) was used to improve the lateral resolution of the system to ~3.5 μ m, but the DOF was narrowed to 14.4 μ m. Figure 5(a) shows onion cells under a microscope. Figure  5   In our experiment, both the amplitude and phase information were used to calculate the wave propagation between different planes. A microscope cover glass was placed on top of the sample, and the top surface of the glass was used as a reference to eliminate the effect of phase fluctuations. For in-vivo imaging applications, it is challenging to maintain phase stability over a whole two-transverse-dimensional scan. The problem of phase stability might be minimized or eventually solved with the development of ultrahigh speed swept laser sources [11]. It is also noteworthy from the above examples that outside the DOF, not only the lateral resolution decreases, but the illumination amplitude attenuates and affects the signal to noise ratio. Although the proposed method improves the lateral resolution and concentrates the energy within a specific out-of-focus layer, it does not correct the attenuation of illumination amplitude outside the DOF along the probe beam. Some fringe artifacts exist in the reconstructed images, and that might be because of the residual phase instability, existence of multiple scattering or because of the limited axial resolution of the system so that structures of adjacent layers overlap together. Multiple scattering is a problem that exists for all OCT systems and will also affect the performance of the proposed method. However, as long as the full-field phase information of the sample can be correctly extracted, the proposed method should work for lateral resolution improvement. Speckle and the diffraction of the Gaussian beam illuminating the out-of-focus sample plane are currently not taken into account in the model. If the Gaussian beam propagation were considered, the digital focusing effects outside the DOF might be further improved.

Conclusion
In conclusion, we have shown that the wave scattering process from out-of-focus scatterers in OCT can be considered as a two-dimensional scalar diffraction model. The lateral resolution in either x or y direction can be improved by use of a non-iterative numerical diffraction algorithm, and high-resolution details can be reconstructed from outside the depth-of-field region without any special hardware in system design. Although a spectrometer-based system is considered in the paper, the proposed method is also applicable to swept-source or full-field OCT systems.
This work was supported by research grants from the National Institutes of Health (EB-00293, NCI-91717, and RR-01192), Air Force Office of Scientific Research (FA9550-04-1-0101), and the Beckman Laser Institute Endowment.