Open Access
Add to BibSonomyAbstract
The resolution in optical coherence tomography imaging is an important parameter which determines the size of the smallest features that can be visualized. Sparse sampling approaches have shown considerable promise in producing high resolution OCT images with fewer camera pixels, reducing both the cost and the complexity of an imaging system. In this paper, we propose a non-local approach to the reconstruction of high resolution OCT images from sparsely sampled measurements. An iterative strategy is introduced for minimizing a homotopic, non-local regularized functional in the spatial domain, subject to data fidelity constraints in the k-space domain. The novel algorithm was tested on human retinal, corneal, and limbus images, acquired in-vivo, demonstrating the effectiveness of the proposed approach in generating high resolution reconstructions from a limited number of camera pixels.
© 2012 Optical Society of America
1. Introduction
Optical coherence tomography (OCT) is an emerging biomedical optical imaging technique that allows for high-resolution, cross-sectional tomographic imaging of the microstructure of biological systems [1]. OCT was first introduced by Huang et al. [2], and was initially applied to retinal imaging in ophthalmology. Since its introduction in the early 1990s, OCT has become a powerful method for imaging the internal structure of biological systems and is particularly widely used in obtaining high-resolution images of the retina and the anterior segment of the eye [3–5]. To date, among other biomedical and clinical applications, OCT has its largest clinical impact in ophthalmology.
Spectral domain OCT (SD-OCT) is an alternative to the original time domain OCT (TD-OCT) technique [6]. SD-OCT uses a broad bandwidth light source at the input and a combination of a high resolution spectrometer and a linear array camera at the detection end of the system to achieve high axial resolution imaging over a large scanning range. The axial imaging resolution in SD-OCT, as in the case of TD-OCT, is determined by the central wavelength and spectral bandwidth of the light source, while the scanning range is determined by the number of illuminated pixels of the camera [7]. SD-OCT has significant advantages over TD-OCT in terms of improved SNR and orders of magnitude higher data acquisition rate [8]. For clinical applications such as surgical guidance and intervention, real-time image tracking, large imaging depth, and high axial resolution are all necessary. A large scanning range requires cameras with larger number of pixels and high resolution spectrometers in order to capture spectra at a high sampling rate, because conventional image reconstruction algorithms require spectral domain sampling beyond the Nyquist rate to achieve a certain imaging depth [9]. In other words, the camera has to have enough pixels to guarantee that at least two data points are sampled within one period of the spectral interferogram, implying CCD or CMOS cameras which are complicated and expensive, with the cost increasing monotonically with the number of pixels, while the smaller size of the pixels poses design and production difficulties and results in higher cost of the spectrometers.
To address the issues associated with large camera arrays, high sampling rates, and long acquisition times, recent studies in sparse sampling and reconstruction [10, 11] have shown that certain types of images can be reconstructed with a high level of accuracy and resolution from data sampled below the Nyquist rate. Much of the early literature on the topic revolved around magnetic resonance imaging (MRI), with the first being a L1 based minimization method [12]. Subsequently, a homotopic L0 minimization approach [13] was shown to provide superior reconstruction results compared to L1 for spine and wrist MRI images, however the effectiveness of this approach may be limited, given its use of the finite differential transform, which is poorly suited for handling the fine details and variations found in other MRI images [14]. Liang et al. [15] explored the use of the non-local total variation regularization framework for MRI reconstruction, where regional characteristics between non-local sites is used as a penalty constraint, and found that such an approach allowed for improved handling of fine details and variations. More recently, the combination of a regional differential sparsifying transform and a homotopic L0 framework was proposed by Wong et al. [14], which was able to lower computational complexity while still maintaining structural fidelity.
More recently, the concept of sparse sampling and reconstruction has begun to be explored for high resolution OCT imagery. Both Mohan et al. [16] and Liu et al. [17] explored the use of L1 minimization methods to reconstruct OCT imagery with highly undersampled k-space data. However, these early works have been evaluated using either simulated signals or objects such as onions, which do not necessarily reflect the characteristics of living tissue. One of the key applications of ophthalmic OCT is the imaging of retinal, corneal, and limbus tissues, which are characterized by complex structural characteristics that are important for clinical observation and diagnosis. Given the complexities of the human retina and limbus, the classical L1 based minimization approach may not provide an accurate reconstruction given sparser sampling, and may lead to blurring and loss of fine structure. As such, the integration of a non-local strategy for sparse OCT reconstruction is worth investigating for improved image quality.
In this paper, a non-local strategy for the sparse reconstruction of OCT imagery is proposed. A homotopic non-local regularization framework is introduced to better handle fine details in the underlying data, hence providing improved high-resolution reconstructions at sparser sampling rates. The framework is implemented via an iterative reconstruction strategy. First, a homotopic non-local regularization step is performed in the spatial domain. Second, data consistency constraints are reinforced in the k-space domain. This novel non-local sparse reconstruction is specifically applied to human retinal, corneal, and limbus OCT data, with results demonstrating that high-quality, high resolution OCT images can be reconstructed using the proposed method.
The organization of this paper is as follows. Section 2 describes the theory and implementation of the proposed method. Experimental results using OCT retinal, corneal, and limbus data are presented and discussed in Section 3, with conclusions drawn in Section 4.
2. Methods
The main principle behind the proposed approach is the integration of non-local regularization [18] within a homotopic L0 minimization framework for reconstructing OCT images from sparsely-sampled measurements acquired using a spectral domain OCT (SD-OCT) imaging system.
In SD-OCT, the light from a broadband light source is split between the reference and sample arms of an interferometer. Light reflected from various depths within the imaged object in the sample arm of the interferometer interferes with light reflected from a stationary mirror in the reference arm. The interference pattern is decomposed into different frequencies by the use of a high resolution spectrometer and each frequency is detected by an individual pixel of a linear array camera attached to the spectrometer and mapped into k-space. Denoting the sample as f(x), and the measurements in the k-space domain as F(k), the relationship between f(x) and F(k) is formulated as
where 𝔽−1 is the inverse Fourier operator.According to the Nyquist criterion [17], the number of samples K in the k-space domain should obey
where xmax is the maximum imaging depth, and Δk is the spectral bandwidth, which is inversely proportional to spatial resolution. A large value for K is necessary to achieve a high imaging resolution and large imaging depth, which means that a large total number of CCD camera pixels is required. As such, strategies are desired that obtain high resolution OCT imagery without requiring a large number of camera pixels.2.1. Lp minimization for sparse reconstruction
In SD-OCT, the number of sampling points is related to the number of illuminated camera pixels, which is dependent on the design of the spectrometer. Such CCD or CMOS cameras and associated electronics are usually costly and cumbersome, making an increase in the number of sampled points technologically very challenging [17]. Instead, our goal is to measure only a fraction of the spectral coefficients, hence fewer pixels required in the camera array, and still be able to get an accurate, high resolution image. Consequently, an undersampled reconstruction fu can be expressed as
where Φ is a measurement operator defining which of the spectral sites are measured (non-measured k-space indices will be assigned to 0 by Φ). Our goal is to reconstruct f(x) from a sparse sampling of F(k). With only few measurements made, not all sites are measured, making this problem essentially an ill-posed inverse problem [19], and may have multiple solutions. A common method to solve this inverse problem is through L2 minimization, as where f̂(x) and F̂(k) are the estimated signals in the spatial and k-space domains, respectively. However, solving the problem in this manner results in significant artifacts such as aliasing and blurring in practical situations [13].Recently, scientists realized that many signals, such as OCT imagery, possessing an inherent sparsity in some sparsifying transformed domain Ψ. According to the emerging theory of compressive sensing [10, 11], a better estimate of f(x) can be obtained by maximizing the sparsity of the signal in the transformed domain and enforcing data fidelity in the k-space domain. This can be formulated as a constrained L0 minimization problem,
Two commonly used transformation operators are the finite differential transformation and the wavelet transformation [12].Unfortunately, solving the L0 problem is essentially NP hard, intractable in practice [20]. To address this problem, pioneering work was done by Candes [10] and Donoho [11], showing that under certain conditions, L0 and L1 minimization lead to the same solutions. Theoretically, solving the L1 problem can get exactly the same solution as solving the L0 problem, at the cost of a substantial increase in the number of measurements required [13].
Typically, a total variation (TV) penalty [18] is employed in the L1 minimization framework, which is known to have an edge-preservation effect: Both Mohan et al. [16] and Liu et al. [17] explored this L1 minimization framework, which represents the state-of-the-art in sparse OCT reconstruction, and their results show that OCT imagery can be reconstructed in a meaningful manner using sparsely sampled measurements in k-space.2.2. Homotopic, non-local regularized reconstruction
There are two fundamental limitations of the reconstruction framework described in Eq. (7) for the purpose of reconstructing retinal, corneal, and limbus OCT imagery [15]. First, the use of a TV penalty can lead to the loss of fine structures at low sampling rates [21]. Since retinal, corneal, and limbus OCT images are characterized by complex morphological details and structural variations, the TV approach may not be well-suited for reconstructing these data. Second, as already mentioned, the number of measurements required for solving the L1 problem can be noticeably higher than that for the L0 problem, therefore requiring a greater number of camera pixels in a SD-OCT system. Hence, an alternative reconstruction strategy that addresses these two problems is desired.
To achieve the theoretical capability of the constrained L0 minimization approach without a drastic increment in the number of measurements, Trzasko et al. [13] introduced a homotopic L0 minimization framework, which uses an increasing approximation framework to get close to the L0 minimization problem:
where ρ is the homotopic approximation of the L0 norm, which approaches the L0 norm as σ approaches 0. To account for the unavoidable noise in the measurements, the data fidelity constraint is typically replaced by a L2 norm constraint: Experimental results showed that this strategy is able to approach the capabilities of the L0 minimization approach [13], thus addressing the problem pertaining to the number of measurements required.To address the problem of detail loss from using total variation, we instead integrate the concept of non-local regularization into the homotopic L0 minimization framework, as it is more is well-suited for handling fine structural details. The proposed homotopic, non-local regularization framework is formulated as:
where η denotes the homotopic non-local regularization functional, and σ controls the approximation degree: where Ω denotes the set of all pixels in the image, 𝒩(x) denote the neighborhood search space for pixel x, N(x) and N(y) denotes the neighborhoods similarity space around x and y, σ controls the decay of the exponential function, and w(x,y,σ) is defined as To reduce the computation complexity, the neighborhood search space is limited to a search window around the pixel to be estimated, so that the regularization by η(|Ψf(x)|, σ) efficiently suppresses erroneous artifacts without destroying the structure of the image. At the same time, we use a data fidelity term to ensure that the reconstructed image complies with the original signal while accounting for some level of noise.2.3. Implementation
An iterative strategy is introduced for solving Eq. (10), with an overview of the proposed sparse OCT reconstruction framework shown in Algorithm 1. Specifically, considering that the data are sampled line-by-line, we assert the non-local regularization constraint over the set of measurements in the spatial domain. The proposed iterative strategy consists of two steps: first, non-local regularization is enforced in the spatial domain, and as such we assert sparsity (i.e., Ψ) in the regional differential sparsifying transform domain in [14], which has shown strong detail preservation in complex tissues. Many other sparsifying transforms may be used in medical imaging reconstruction, nevertheless, a comprehensive analysis of different sparsifying transforms is beyond the scope of this paper. Then data fidelity is enforced in the k-space domain to correct estimation errors. This two-step process is iteratively repeated, with the approximation degree σ used to control the homotopic non-local regularization functional η decreasing as the number of iterations increases, until convergence. Empirically, all experimental applications of this proposed method converged successfully. The pseudocode for the proposed homotopic, non-local sparse reconstruction algorithm is shown in Algorithm 1.
Algorithm 1. Homotopic, non-local sparse reconstruction for OCT imagery
3. Results and discussions
To evaluate the effectiveness of the proposed homotopic non-local regularized (NLR) method, we applied it to a series of OCT images acquired in-vivo from the human i) retinal fovea (Figs. 1, 2), ii) cornea (Figs. 3, 4), and iii) limbus (Figs. 5, 6).
Fig. 1 Reconstruction results from 50% of the acquired human retinal fovea data. The colored boxes mark sections that are enlarged in Fig. 2. The reconstruction results using 100% of the samples are provided as a reference.
Fig. 2 Zoomed-in regions from Fig. 1. The fine details and the structural detail of the individual layers are better maintained in the reconstructed results produced using NLR compared to L2 and L1.
Fig. 3 Reconstruction results from 50% of the acquired human corneal data. The colored boxes mark sections that are enlarged in Fig. 4.
Fig. 4 Zoomed-in regions from Fig. 3. The structural detail of the individual layers as well as the keratocyte cells are better maintained in the reconstructed results produced using NLR compared to L2 and L1.
Fig. 5 Reconstruction results from 50% of the acquired human limbus data. The colored boxes mark sections that are enlarged in Fig. 6.
Fig. 6 Zoomed-in regions from Fig. 5. The structural detail of the individual layers as well as the fine structural details are better maintained in the reconstructed results produced using NLR compared to L2 and L1.
The images were acquired with a research grade, high-speed, SD-OCT system [22], operating at 1060nm wavelength, that utilizes a super-luminescent diode (λ-c = 1020nm, δλ = 110nm, Pout = 10mW) and a 47kHz InGaAs linear array, 1024 pixel camera (SUI, Goodrich) interfaced with a high performance spectrometer. The SD-OCT system provides 3μm axial and 15μm lateral resolution in the human cornea and limbus, and 6μm axial resolution in the human retina, and 100dB SNR for 1.3mW of optical power incident on the sample. The OCT images were acquired from healthy subjects using an imaging procedure carried out in accordance with University of Waterloo research ethics regulations, at 1000 A-scans. To reconstruct OCT imagery from only a percentage of camera pixels, SD-OCT spectral data is randomly sampled from the camera array using a uniformly distributed pseudo-random mask. The obtained spectral data was then used to populate the k-space grid according to the known functional dependency of wavenumber on pixel index [17].
The tested sparse sampling approaches were implemented in MATLAB and tested on an AMD Athlon II X3 3.10 GHz machine with 3.25GB of RAM. For comparison purposes, the proposed method was compared with two other sparse reconstruction techniques: the baseline L2 minimization approach described in Section 2.1, and the L1 minimization approach with total variation, as proposed by Mohan et al. [16] and Liu et al. [17] for sparse OCT reconstruction. To perform a comprehensive and systematic assessment of the reconstruction performance of the different methods, the signal to noise ratio (SNR) was computed for reconstructed data using between 30% and 70% of the camera pixels. The SNR metric was computed as follows:
and (MSE) was defined as mean squared error between original image and reconstructed image: where f(x) is original image, f̂(x) is reconstructed image, and N is the number of pixels in each image.Furthermore, a qualitative visual assessment is performed on the reconstructed data to investigate the reconstruction performance and the preservation of details achieved using the tested methods.
For implementation purposes, we restrict the size of the search window (𝒩) to 7 × 7, and the size of a neighborhood (N) to 3×3. If N is the number of pixels of the image, then the final computational complexity of the algorithm is 9 × 49 × N.
Figures 1, 3, and 5 show examples for each of the three types of in-vivo human OCT imaging data, each reconstructed using the three reconstruction methods from spectral data acquired using 50% of the camera pixels. To visualize the improvements obtained from using the proposed method, regions of interest (ROI) extracted from the reconstructed images (shown as colored boxes in Figs. 1, 3, and 5) are shown for better visual comparison in Figs. 2, 4, and 6. The human retina (Fig. 1) contains a number of morphological details such as individual cellular and plexiform layers, cross-sections of retinal capillaries (circular black features in the retina), as well as large choroidal blood vessels (pale circular or elongated features below the retinal pigmented epithelium — the last thick black line). We can see very clearly that the L2 method leads to considerable blur and artifacts, making it difficult to see any of the underlying structure and detail in the reconstructed retinal image except for the very highly reflective (black) lines corresponding to the inner and outer photoreceptor junctions and the retinal pigmented epithelium. The rest of the retinal layers, along with the inner retina vasculature, cannot be vizualized because of the algorithm induced blur. The L1 method results in noticeably better image quality as compared to that produced using the L2 method, although the contrast of the individual retinal layers is not as good as in the original image. The proposed novel NLR algorithm results in overall higher contrast of the individual retinal layers as compared to L1, with features such as the thin, highly reflective retinal layers and the cross-sections of the retinal capillaries appearing sharp and are distinctly visible (Fig. 2), appearing similar to the reconstructed retinal image from 100% of the acquired samples.
The human cornea (Fig. 3) contains a number of morphological features of different size and optical properties. It has four distinct layers: the epithelium, Bowman’s membrane, stroma, and Descemett’s-endothelial complex, each of which are clearly visible in the image reconstructed from 100% of the acquired samples. The small black dots in the corneal stroma correspond to reflections from keratocyte cells. As observed in Fig. 3 and Fig. 4, the L1 method results in an image where most of the layers and some of the keratocytes are still visible, however, the overall contrast of the image is drastically lower as compared to the image reconstructed from 100% of the acquired samples. The NLR approach result in significantly better reconstruction of the corneal morphological details, as well as higher image contrast as compared to the L2 and L1 methods. Once again, the image reconstructed using the NLR approach is closer to the image reconstructed from 100% of the acquired samples.
The human limbus is a transitional region of tissue, located between the cornea and the sclera, that contains a number of blood and lymph vessels, as well as the transition of the corneal epithelium into the scleral conjunctiva. When applied to OCT images of the human limbus (Fig. 5 and Fig. 6), the L2 method generates artifacts and results in overall poor contrast and visibility of the tissue morphology. The L1 and the NLR approaches generate images of significantly higher quality than L2 with fairly good contrast and preservation of the morphological details. As expected, the overall image contrast of the image reconstructed using the NLR method is better than that produced using the L1 method, and closer to the quality of the image reconstructed from 100% of the acquired samples.
Results from a quantitative comparison of the image SNR for the sparse sampling approaches are presented in Fig. 7, as a function of the percentage of camera pixels used for reconstruction. It can be observed that the proposed NLR method produces reconstructed OCT data with higher SNR values for all levels of camera pixel undersampling when compared to the other two methods.
Fig. 7 SNR of the reconstructed retinal, corneal and limbus data as a function of the percentage of camera pixels used for sparse reconstruction. Observe the substantial improvement in all cases of the proposed NLR method over the L2 and L1 reconstruction methods.
Therefore, based on both quantitative SNR analysis as well as qualitative visual assessment, we can conclude that the proposed NLR method demonstrates improved reconstruction performance and visual quality compared to the other tested methods when using a reduced number of camera pixels.
4. Conclusions
A novel approach for the sparse reconstruction of OCT data using a homotopic, non-local regularized approach was presented. The reconstruction algorithm is compared with conventional L2 and L1 total variation reconstruction methods. Results show that the proposed approach is able to achieve a significantly higher signal-to-noise ratio and visual quality using a lower number of camera pixels, thus illustrating the potential for obtaining high resolution images with lower equipment costs and reduced imaging times. In future work, a comprehensive analysis of different sparsifying transforms (wavelets, curvelets, etc.) will be investigated on different types of OCT imagery to identify the optimal transform for reconstructing OCT imagery from sparse measurements.
Acknowledgments
This project was funded in part by the Natural Sciences and Engineering Research Council of Canada, the Canadian Institutes for Health Research, the Chinese Council Scholarship, and the University of Waterloo.
References and links
1. S. A. Boppart, “Optical coherence tomography: technology and applications for neuroimaging,” Psychophysiology 40, 529–541 (2003). [CrossRef] [PubMed]
2. 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, 1178–1181 (1991). [CrossRef] [PubMed]
3. B. E. Bouma and G. J. Tearney, Handbook of Optical Coherence Tomography (Informa Healthcare, 2001).
4. M. Wojtkowski, V. Srinivasan, J. G. Fujimoto, T. Ko, J. S. Schuman, A. Kowalczyk, and J. S. Duker, “Three-dimensional retinal imaging with high-speed ultrahigh-resolution optical coherence tomography,” Ophthalmology 112, 1734–1746 (2005). [CrossRef] [PubMed]
5. W. Drexler and J. G. Fujimoto, Optical Coherence Tomography (Springer, 2008). [CrossRef]
6. M. Wojtkowski, A. Kowalczyk, R. Leitgeb, and A. F. Fercher, “Full range complex spectral optical coherence tomography technique in eye imaging,” Opt. Lett. 27, 1415–1417 (2002). [CrossRef]
7. R. Leitgeb, W. Drexler, A. Unterhuber, B. Hermann, T. Bajraszewski, T. Le, A. Stingl, and A. Fercher, “Ultrahigh resolution Fourier domain optical coherence tomography,” Opt. Express 12, 2156–2165 (2004). [CrossRef] [PubMed]
8. R. Leitgeb, C. Hitzenberger, and A. Fercher, “Performance of Fourier domain vs. time domain optical coherence tomography,” Opt. Express 11, 889–894 (2003). [CrossRef] [PubMed]
9. A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography—principles and applications,” Rep. Prog. Phys. 66, 239–303 (2003). [CrossRef]
10. E. Candës, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52, 489–509 (2006). [CrossRef]
11. D. Donoho, “Compressive sensing,” IEEE Trans. Inf. Theory 52, 1289–1306 (2006). [CrossRef]
12. M. Lustig, D. Donoho, and J. Pauly, “Sparse MRI: the application of compressed sesing for rapid MR imaging,” Magn. Reson. Med. 58, 1182–1195 (2007). [CrossRef] [PubMed]
13. J. Trzasko and A. Manduca, “Highly undersampled magnetic resonance image reconstruction via homotopic l0-minimization,” IEEE Trans. Med. Imag. 28, 106–121 (2009). [CrossRef]
14. A. Wong, A. Mishra, D. Clausi, and P. Fieguth, “Sparse reconstruction of breast MRI using homotopic L0 minimization in a regional sparsified domain,” Biomed. Eng. IEEE Trans, 1–10 (2010).
15. D. Liang, H. Wang, and L. Ying, “SENSE reconstruction with nonlocal TV regularization,” Proc. IEEE Eng. Med. Biol. Soc., 1032–1035 (2009).
16. N. Mohan, I. Stojanovic, W. C. Karl, B. E. A. Saleh, and M. C. Teich, “Compressed sensing in optical coherence tomography,” Proc. SPIE 7570, 75700L (2010). [CrossRef]
17. X. Liu and J. U. Kang, “Compressive SD-OCT: the application of compressed sensing in spectral domain optical coherence tomography,” Opt. Express 18, 22010–22019 (2010). [CrossRef] [PubMed]
18. G. Gilboa and S. Osher, “Nonlocal operators with applications to image processing,” Tech. Rep. CAM Report 07-23, Univ. California, Los Angeles, 2007.
19. P. Fieguth, Statistical Image Processing and Multidimensional Modeling (Springer, 2010).
20. B. K. Natarajan, “Sparse approximate solutions to linear systems,” SIAM J. Comput. 24, 227–234 (1995). [CrossRef]
21. W. Guo and F. Huang, “Adaptive total variation based filtering for MRI images with spatially inhomogeneous noise and artifacts,” Int. Sym. Biomed Imag , 101–104 (2009).
22. P. Puvanathasan, P. Forbes, Z. Ren, D. Malchow, S. Boyd, and K. Bizheva, “High-speed, high-resolution Fourier-domain optical coherence tomography system for retinal imaging in the 1060 nm wavelength region,” Opt. Lett. 33, 2479–2481 (2008). [PubMed]
