This paper proposes an automatic point spread function (PSF) estimation method to de-blur out-of-focus optical coherence tomography (OCT) images. The method utilizes Richardson-Lucy deconvolution algorithm to deconvolve noisy defocused images with a family of Gaussian PSFs with different beam spot sizes. Then, the best beam spot size is automatically estimated based on the discontinuity of information entropy of recovered images. Therefore, it is not required a prior knowledge of the parameters or PSF of OCT system for de-convoluting image. The model does not account for the diffraction and the coherent scattering of light by the sample. A series of experiments are performed on digital phantoms, a custom-built phantom doped with microspheres, fresh onion as well as the human fingertip in vivo to show the performance of the proposed method. The method may also be useful in combining with other deconvolution algorithms for PSF estimation and image recovery.
Optical coherence tomography (OCT) is a non-invasive, micrometer resolution, three-dimensional optical imaging technique capable of capturing backscattering photons from within optical scattering media such as biological tissues [1, 2]. One of the advantages of OCT is that the axial and the lateral resolutions are decoupled from each other. The axial resolution of the OCT system is determined by the spectral bandwidth of light source and is defined as the half of the coherence length of the light source which can be expressed by
The lateral resolution of the OCT system is determined by the objective lens that is used to deliver/focus the light onto the sample and is expressed byFig. 1 for the cases of low (left figure) and high NAs (right figure). From the illustrations, we can appreciate that only the portion of the 3D image within DOF exhibits the desired lateral resolution, while that in the de-focal plane is blurred. There exist several dynamic focusing methods in the literature to mitigate the compromise between lateral resolution and DOF. In general, these methods can be divided into hardware-based and software-based (digital) approaches.
Yasuno et al.  applied a numerical deconvolution method to cancel lateral defocus in Fourier-domain OCT using a depth-dependent lateral PSF and applied the method to OCT images of in vivo human anterior eye segment. Kulkarni et al.  proposed a linear shift-invariant system model to describe coherent light specimen interactions in OCT images and applied an iterative deconvolution algorithm to enhance the sharpness of the images of biological structures. A constrained iterative restoration algorithm was utilized to estimate the impulse response of the system model using a prior knowledge of the properties of that response . The algorithm was applied to a defocused OCT image of fresh onion sample and the deconvolved OCT image revealed polygonal cellular structure in the specimen. Ralston et al.  proposed a method for reducing transverse blurring and improving the transverse resolution in OCT using Gaussian beam deconvolution. First, the direct inverse problem was investigated using two types of regularization, truncated singular value decomposition and Tikhonov. Then, an iterative expectation-maximization algorithm, the Richardson-Lucy algorithm, with a beam-width-dependent iteration scheme was developed. Using a dynamically iterative Richardson-Lucy algorithm reduced the transverse blurring by providing an improvement in the transverse PSF for sparse scattering samples in regions up to two times larger than the confocal region of the lens. Also, Wiener and Richardson-Lucy OCT image restoration algorithms were implemented and compared with each other where the images deconvolved using Richardson-Lucy had a better contrast and quality . Besides, JM Schmitt  derived CLEAN deconvolution kernel to provide an effective way of improving resolution and contrast in optical coherence tomography. Rolland et al.  developed a Gabor-based fusion technique based on the concept of the inverse local Fourier transform and the Gabor’s signal expansion to produce a high lateral resolution image throughout the depth of imaging.
Since these previous digital methods require a prior knowledge of the optical parameters or PSF of the optical system, an initial calibration experiment is required. Tomlins et al.  and  made PSF phantoms doped with a small quantity of calibrated spherical micro-particles to estimate the Gaussian PSF parameters and researched the potential PSF degrading effects of optical dispersion and spherical aberration. Tomlins et al.  also adopted femtosecond laser subsurface micro-inscription techniques to fabricate an OCT test artifact for validating the resolution performance of a commercial OCT system. The method  and  were generalized by Woolliams et al.  to estimate 2D PSF at multiple locations across the B-scans to account for the variation in PSF as a function of both depth and lateral position. Agrawal et al.  present a PSF measurement approach using a nanoparticle-embedded phantom, towards the development of standardized test methods. However, the optical parameters or the PSF of OCT system needed to be calculated by performing a separate experiment on a specialized phantom to calibrate the parameters and the parameters may not be valid for complex heterogeneous structures.
In this paper, by using the forward OCT model previously developed by Ralston et al , we propose an automatic PSF estimation method based on discontinuity of information entropy and Richardson-Lucy deconvolution algorithm to improve the image contrast and quality. The advantage of this method is that there is no need to implement a new optical hardware or to measure the optical parameters of the OCT system, and importantly it is relatively in-sensitive to the phase-instability of the measurements. Since the Richardson-Lucy algorithm tends to concentrate energies near boundaries, it provides a good approximation to cellular boundaries and subcellular features and is more robust to noise. Therefore, combining PSF estimation method and Richardson-Lucy algorithm can estimate PSF of OCT at different defocal planes automatically and objectively while decreasing the transverse blurring.
2. Image recovering principle and algorithms
In Gaussian optics model, the lateral PSF of OCT is approximated by :
In OCT system, the final image at each depth location is the convolution of the lateral PSF and the sample profile which can be formulated by23]. By estimating the best beam spot size at each depth location, the lateral resolution at out-of-focus depth locations can be improved by solving an inverse problem.
2.1 Image recovering flow diagram
The flow diagram of the proposed method for recovering defocused FD-OCT images is given in Fig. 2 .
Step1: The non-linear spectral interferogram captured in an A-scan is rescaled into linear k-space followed by Fast Fourier Transform (FFT) to become A-line complex information (phase and amplitude) along the z axis, whose amplitude represents the structural information. The 3D structure of a sample is acquired by two-dimensional scanning of the probe beam by a paired galvanometer mirror.
Step2 and Step3: The 3D structure image is resampled along the depth (in the z direction) to obtain a sequence of en-face images I(x,y) at different depth locations followed by a Gaussian low-pass filter. Here, the Gaussian LPF acts to reduce the noise levels in the image, including the system noise and the spurious effects from speckle, so that the recovered image becomes smooth. It should be noted that applying a Gaussian LPF will reduce the lateral PSF of OCT system as follows:
Step4 and Step5: Richardson-Lucy deconvolution method requires a prior knowledge about the lateral PSF, , which is a function of the actual beam spot size . In order to estimate the best choice of at each depth location, the filtered en-face images are deconvolved using Richardson-Lucy algorithm by a Gaussian family PSFs with stepped beam spot size and an evaluation function defined in subsection 2.3 is analyzed. The evaluation functions tends to have a discontinuity around , so, PSF(The best choice of ) is acquired by searching a peak of the Laplacian of evaluation function .
Richardson-Lucy is an iterative method based on the maximum-likelihood estimation for Poisson statistical data, which is updated by
Step6: Each original en-face image can be recovered using Richardson-Lucy algorithm and estimated lateral PSF ().
2.2 Discontinuity of recovered image around
The output of OCT imaging system is described mathematically by convolving its PSF with an input image. Let’s assume that a sufficiently small object can be modeled as a Dirac delta function. Then, the output of the imaging system is its lateral PSF as follows:Equation (12) is obtained from Eq. (10) and Eq. (11) and can be solved by using 2-dimentional convolution theorem as Eqs. (12) and (13).
It can be observed that the recovered image will have a discontinuity at Wz = Wza.
According to sampling property of delta function, any input image structure can be modeled as weighted superposition of delta functions.
2.3 Evaluation function for PSF estimation
Information entropy of the recovered image and the blurry out-of-focus raw image can be expressed by
The evaluation function is defined byFig. 3 .
2.4 Image recovery
If the object en-face image is located in the focal plane (Wz = W0), the received image will be given by
Assuming that the object en-face image O(x,y) is located in the out-of-focus plane (Wz = Wza), the received image will be given by
The focused image I0(x,y) can be recovered from out-of-focus image I(x,y) by filtering the out-of-focus image with the compensation function g(x,y) as follows
3. Simulation studies
In order to show the performance of the proposed method, we performed a series of simulations on digital phantoms. First, a Dirac delta function digital phantom was generated and manually defocused at different Wza. For each Wza scenario, Richardson-Lucy algorithm was performed to recover the image on a family of Wz values and the discontinuity (peak of the Laplacian) of the evaluation function was measured to estimate the beam spot size. Figure 4(a) shows a 3D plot of the evaluation function for different Wzas where the x axis is Wz, the y axis is Wza and z axis (intensity) is the evaluation function. It can be observed that the discontinuity of the evaluation function is at Wz = Wza. Also, the Laplacian of the evaluation function is shown in Fig. 4(b) where x axis is Wz, y axis is Wza and z axis is the Laplacian of the evaluation function. The peak values are observed around Wz = Wza, which was expected.
Figure 5(a) shows another digital phantom to simulate the cellular structures of onion. The image was manually defocused at = 3(Fig. 5(b)) and the proposed method was performed to recover the image where the recovered image is shown in Fig. 5(f). In order to study the effectiveness of the proposed method in suppressing noise, the image was degraded by Poisson and speckle noise. Figure 5(c-e) show defocused and degraded images by speckle noise, Poisson noise, and speckle and Poisson noise, respectively. Figure 5(g-i) show the recovered images from different noise conditions. The multiplicative speckle noise was a zero-mean uniformly distributed random noise with variance 0.05.
The experiment was repeated for different defocused depth locations () and the difference between the actual and estimated a number of beam spot size (estimation error) was measured at different noise conditions. Figure 6 shows the error in estimating actual under different noise conditions. It can be observed that the algorithm can effectively suppress the noise effect and the results are comparable with the experiment when no noise is present. The proposed method was specifically more effective when the image was slightly defocused ( [4 9]). In general, the estimation error was acceptable and the proposed method could effectively estimate the beam spot size at the present of different noise conditions.
4. Experimental results
The OCT system used to achieve 3D volume data is shown in Fig. 7 , which is similar to one of previous studies . Briefly, the system used a super luminescent diode light source, with a central wavelength of 1300 nm and a bandwidth of 80 nm that provided a 12 µm axial resolution in air. The light is split into two paths in a fiber based Michelson interferometer. One beam is coupled onto a stationary reference mirror and the second beam is focused with a collimating lens and a focal lens (NA = 0.05) at the sample with a theoretical lateral resolution of ~16 μm and a depth of focus of ~350 μm. The focal spot on the sample was scanned using a pair of galvanometer mirror mounted in the sample arm. The spectral interferogram between the reference light and the light backscattered from the sample was sent to a homebuilt spectrometer via an optical circulator. The spectrometer consisted of a transmission grating (1200 lines/mm), a camera lens with a focal length of 100 mm, and 1024 element line scan InGaAs detector (capable of a line rate of 47 kHz), with a wavelength coverage of 200 nm. A pilot laser is also coupled into the interferometer for beam guidance.
In order to show the performance of the proposed method in real OCT images, three sets of experiments were performed on: 1) a phantom doped with microspheres that simulate the PSF of the system, 2) a fresh onion to show the efficiency of the proposed method on the still sample, and 3) in vivo tissue sample to demonstrate whether the proposed method has the potential for in vivo imaging applications. For the results presented below, the processing time of one frame image (512*512) was 1.23s for a HP dv6 laptop computer (processor: Intel i3 2.4GHz; RAM: 6G) when number of iterations was 10. Also note that during the image reconstruction process, there was no threshold applied to the images. Only when the OCT images were displayed, a fixed threshold (equal to the noise level in the raw OCT image) was applied to better show the final images.
First, a phantom was made from gelatin, doped with microspheres, having a predominant diameter of 300 nm. Since the microspheres are small relative to the system resolutions, they can be approximated by delta functions and the phantom can be modeled as a superposition of delta functions. 3D OCT volume imaging was obtained from the phantom and the focus of the system was set to the center of the phantom.
Figure 8(a) shows one B-scan of the phantom where the center of the image is in focus but the upper and lower portion of the image has a degraded lateral resolution and the microspheres look blurry. Then, the proposed automatic PSF estimation and image recovery (see Eq. (22) and (23)) method was applied on to the 3D volume data at different depth locations and the recovered image of the same B-scan is shown in Fig. 8(b). Note that the blurred image in the out-of-focus area is now focused after performing the proposed algorithm on the expérimental data. In addition, we can also recover the raw B-scan image (Fig. 8(a)) to object profile (Fig. 8(c)), which are delta functions when microspheres are sufficiently small, by deconvolution (see Eq. (7)) using estimated beam spot size at different depth locations. Note that some connected microspheres in raw B-scan image are seperated and the image quality is further improved.
Then, an experiment was performed on a fresh onion and 3D OCT images with 256 A-scans and 256 B scans were acquired from the focus area and 1.5 mm away from the focal point. Figure 9(a) shows one en face image of the onion in the focal depth and the proposed method is applied to improve the contrast and resolution of the image in the presence of noise. Figure 9(b) shows the recovered image from Fig. 9(a). Then, while the probe beam was kept unchanged, the sample was shifted downwards by 1.5 mm. Figure 9(c) is the resulted image that corresponded to the same depth location of sample shown in Fig. 9(a). The proposed method was applied to recover that image (Fig. 9(c)), and the result is shown in Fig. 9(d). Note the similarity between Fig. 9(d) and Fig. 9(b) shows onion structures at the same depth location.
The recovered images of onion layers at different depths are shown in Fig. 10 where Figs. 10(a-d) show the defocused onion layers at 1500 µm, 1604 µm, 1725 µm and 1823 µm, respectively. Figures 10(e-h) show the recovered images at their corresponding depth locations in the onion.
To show the discontinuity during the recovering process, Fig. 11 gives the 2D plots of PSF estimation curve for recovering from Fig. 10(a) to Fig. 10(e), where the changes of evaluation function and Laplacian of the evaluation function with the beam spot size Wz are displayed. It can be observed that the discontinuity of the evaluation function and peak value of the Laplacian of the evaluation function occur when Wz = 4.0, meaning that the estimated beam size was 4.0 for this en face image.
The theoretical beam spot size is compared with the estimated beam size in the onion at different depths and is shown in Fig. 12 where the x axis is the depth in um and y axis is the beam spot size in µm. It can be observed that the estimated beam spot sizes agree with the expected theoretical values. The maximum estimation error is 4.2 µm and the maximum relative error is 12.8% for a depth location between 1.5 mm to 2.3 mm.
Lastly, the preliminary experiments were performed on the OCT images that were acquired from the human fingertip using the system described in Fig. 7 to demonstrate whether the proposed method is applicable in vivo. For this set of experiments, the system was running at 120 frames per second. 3D OCT images consisted of 256 A-scans and 200 B scans, covering an area of 1.5x1.2 mm2 on the skin surface. The data was first obtained from the focus region. And then with the other system setup fixed, the finger was moved 500 µm away from the focal point, and another 3D image was acquired. Because of the system spatial resolution was about 16 μm, the image quality enhancement for the regions within dermis was not easily appreciated. For this reason, we decided to show here the ability of proposed method to provide enhanced imaging quality for the de-focused sweat glands. Figure 13(a) shows one typical en-face image located at the focus region, which was dissected from the 3D volume data set. This en-face image represented a plane approximately cut through the upward oriented sweat glands with 90 degrees, in which the sweat glands (shown as the bright spots) and the skin ridges are clearly seen. Figure 13(b) shows one defocused en-face image at the approximately the same location as in Fig. 13(a). The proposed method was then used to de-blur the de-focused image. The result is given in Fig. 13(c), where it can be seen that the recovered image has qualitatively the same quality as the focused image. Quantitatively, the estimated entropies for the three images were 6.9699 (Fig. 13(a)), 7.0128 (Fig. 13(b)), and 6.9356 (Fig. 13(c)), respectively. Thus, these preliminary results demonstrate that the proposed approach is efficient in de-blurring the in vivo OCT images. Further refinement of the current model for in vivo imaging applications are still in progress.
We proposed an automatic PSF estimation method for recovering defocused OCT images based on the discontinuity of information entropy. A family of beam spot sizes were used to recover the original image using Richardson-Lucy algorithm and the maxima in the Laplacian of the introduced evaluation function was utilized to find the best PSF. The performance of the proposed method was first tested on digital phantoms and then on a custom-built phantom, onion cells as well as on human finger in vivo. The quantitative performance of the method was quantized by knowing the parameters of the system and comparing the estimated parameters with their real values. It was observed that the proposed method is effective to recover defocused images in the presence of noise and extend depth of imaging range to 2.8mm. This method can be used in recovering defocused OCT images where the PSF is not known, especially when the index of refraction at different depths of the sample is unknown and not fixed. The proposed method used an approximated OCT forward model when NA of OCT system is relatively small. In future, the combination of the proposed PSF estimation method with accurate OCT forward model, such as that described in interferometric synthetic aperture microscopy , will be explored to account for the diffraction and the coherent scattering of light by the sample, so that a better representation of the sample’s microstructure outside the focus region is presented. The proposed image recovery approach assumes that the PSF exhibits linear shift invariance throughout an en face plane, which was achieved by confining the beam scanning to the central region of the objective lens.
The proposed PSF estimation method can be combined with other deconvolution algorithms such as wiener filters for PSF estimation and image recovery. Furthermore, the proposed PSF estimation and image deconvolution method is not limited to OCT and may be useful for image defocusing in other confocal optical system.
6. T. Xie, S. Guo, Z. Chen, D. Mukai, and M. Brenner, “GRIN lens rod based probe for endoscopic spectral domain optical coherence tomography with fast dynamic focus tracking,” Opt. Express 14(8), 3238–3246 (2006). [CrossRef] [PubMed]
7. Z. H. Ding, H. W. Ren, Y. H. Zhao, J. S. Nelson, and Z. P. Chen, “High-resolution optical coherence tomography over a large depth range with an axicon lens,” Opt. Lett. 27(4), 243–245 (2002). [CrossRef] [PubMed]
8. J. Holmes and S. Hattersley, “Image blending and speckle noise reduction in multi-beam OCT,” Optical Coherence Tomography and Coherence Domain Optical Methods in Biomedicine XIII, Proc. of SPIE Vol. 7168, 71681N (2009).
9. T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, “Interferometric synthetic aperture microscopy,” Nat. Phys. 3(2), 129–134 (2007). [CrossRef]
10. L. Yu, B. Rao, J. Zhang, J. Su, Q. Wang, S. Guo, and Z. Chen, “Improved lateral resolution in optical coherence tomography by digital focusing using two-dimensional numerical diffraction method,” Opt. Express 15(12), 7634–7641 (2007). [CrossRef] [PubMed]
11. Y. Yasuno, J. I. Sugisaka, Y. Sando, Y. Nakamura, S. Makita, M. Itoh, and T. Yatagai, “Non-iterative numerical method for laterally superresolving Fourier domain optical coherence tomography,” Opt. Express 14(3), 1006–1020 (2006). [CrossRef] [PubMed]
12. M. D. Kulkarni, C. W. Thomas, and J. A. Izatt, “Image enhancement in optical coherence tomography using deconvolution,” Electron. Lett. 33(16), 1365–1367 (1997). [CrossRef]
13. R. K. Wang, “Resolution improved optical coherence-gating tomography for imaging biological tissue,” J. Mod. Opt. 46, 1905–1913 (1999).
14. T. S. Ralston, D. L. Marks, F. Kamalabadi, and S. A. Boppart, “Deconvolution Methods for Mitigation of Transverse Blurring in Optical Coherence Tomography,” IEEE Trans. Image Process. 14(9), 1254–1264 (2005). [CrossRef] [PubMed]
16. J. M. Schmitt, “Restoration of optical coherence images of living tissue using the CLEAN algorithm,” J. Biomed. Opt. 3(1), 66–75 (1998). [CrossRef]
18. P. H. Tomlins, P. Woolliams, M. Tedaldi, A. Beaumont, and C. Hart, “Measurement of the three-dimensional point-spread function in an optical coherence tomography imaging system,” Proc. SPIE 6847, 68472Q, 68472Q-8 (2008). [CrossRef]
19. P. H. Tomlins, R. A. Ferguson, C. Hart, and P. D. Woolliams, “Point-spread function phantoms for optical coherence tomography,” NPL Report OP 2 (National Physical Laboratory, pp: 1754–2944, (2009).
21. P. H. Tomlins, G. N. Smith, P. D. Woolliams, J. Rasakanthan, and K. Sugden, “Femtosecond laser micro-inscription of optical coherence tomography resolution test artifacts,” Biomed. Opt. Express 2(5), 1319–1327 (2011). [CrossRef] [PubMed]
22. A. Agrawal, T. J. Pfefer, N. Gilani, and R. Drezek, “Three-dimensional characterization of optical coherence tomography point spread functions with a nanoparticle-embedded phantom,” Opt. Lett. 35(13), 2269–2271 (2010). [CrossRef] [PubMed]