1. Introduction

Optical coherence tomography (OCT) is an imaging modality for obtaining tomography images of subsurface structures particularly for biological tissues [2, 3]. Since it was invented in 1991 [1], the technique has evolved from its early version in time-domain OCT (TD-OCT) to the prevalent frequency domain OCT (FD-OCT) implementation [4]. FD-OCT makes use of the property that the depth profile of the subsurface interfaces are encoded by Fourier Transform in the spectral interferogram [5]. Then, the depth profile can be recovered by leveraging on the fast Fourier transform. However, this relation is only valid when the spectral interferogram is expressed in wavenumber or in frequency domain. In practice, the spectral interferogram is measured by an optical spectrum meter in the wavelength or the lambda domain. Thus, interpolation and resampling are usually needed to transform the frequency measurements from the wavelength domain to the frequency domain before fast Fourier transform is performed. Such a process results in resolution degradation and limited imaging depth of FD-OCT method [6].

In order to address this issue, one attempt is to use advanced interpolation methods [6], and the other attempt is to develop reconstruction methods which are alternative to Fourier transform and directly based on the frequency measurements in the wavelength domain [7]. The use of non-uniform Fourier transform via the Vandermonde matrix [7,8] is a result of the second attempt. The Vandermonde matrix is formed by the Fourier kernel on a non-uniform wavenumber grid which is converted from the uniform wavelength grid captured from an optical spectrum meter. The depth profile is reconstructed by applying the inverse of the Vandermonde matrix on the frequency measurements in the same non-uniform wavenumber grid. Considering the numerical instabilities of calculating the inverse of the Vandermonde matrix, a forward transform is adopted in later literatures [9].

An alternative approach of the second attempt makes use of algebraic reconstruction technique (ART) [10] to reconstruct the depth profile from the frequency measurements in wavelength domain. ART relates the frequency measurements directly to the depth profile in an algebraic format like Ax=b, where b denotes the frequency measurements in wavelength domain, vector x denotes the discretized depth profile of the subsurface interfaces in an axial scan, and A is the measurement matrix which also contains the information of the white light used. The basic principle of ART involves solving the algebraic inverse problem that is to reconstruct x from b. Compared to other OCT methods, the ART based OCT method gives directly the depth profile without using any resampling or interpolation. High resolution can be achieved by using small discretization interval in spatial coordinates, which results in a high dimension of vector x. Another merit point of ART is that the imaging depth would be free from Shannon Nyquist sampling theorem by using advanced signal processing techniques such as the compressed sensing technique [?]. However, ART method is not widely adopted because of the high computational load particularly when ones attempt to have high resolution and large imaging depth.

This paper proposes an ART based OCT method by exploring the sparsity of the depth profile to reduce the computational load. The idea is drawn from the observation that the depth profile detectable by OCT are either of a layer structure [11] or of discrete scatterers [12] or of mixture of both. This implies that vector x in the above mentioned algebraic representation is essentially sparse. Consequently, this presents an opportunity to use a greatly reduced dimension of b in ART, and in turn, the computational load of ART can be greatly reduced. Recent advances in compressed sensing [?] has demonstrated the feasibility of solving sparse unknowns from measurements of reduced dimension. The ART method proposed in this paper is based on the ℓ₁-optimization technique used in compressed sensing. It inherits the merits of ART in delivering high resolution of OCT, while overcoming the limitations in computational load.

In this paper, we relate the wavelength domain OCT measurements to the depth profile in an algebraic form by discretizing the depth profile in spatial domain. Then ℓ₁ norm of x is minimized to reconstruct the sparse depth profile including the positions of the subsurface interfaces and the reflectivity at the surfaces. By this way, the problem of axial scan reconstruction of the wavelength domain OCT is formulated into a convex optimization problem that can be solved by interior point based algorithm. The simulation and experiment results have shown that the proposed method delivers higher resolution and longer imaging depth compared to the FD-OCT method. The performance of the proposed method is also demonstrated by imaging an onion skin.

2. Principle of Algebraic OCT Based on ℓ₁-optimization

The proposed method is based on the frequency measurement of the FD-OCT setup which is shown in Fig. 1. A white light beam is tapped from a fiber-based low coherence light source and it is incident through a fiber collimator onto a beam splitter which splits the light into two arms of beam. One arm of beam penetrates the sample along the axial direction and it is reflected at subsurface interfaces in the sample. Another beam would reflect off the reference mirror that is stationary at a predetermined position. Eventually, light reflected from both arms mutually interferes at the beam splitter and the interfered light is transmitted through the fiber collimator and a fiber circulator to an optical spectrum analyzer (OSA). At OSA, the spectral interferogram is obtained, based on which an axial scan of the depth profile is reconstructed. The OCT image of the sample is obtained by assembling all the axial scan reconstruction results with the assistance of a motion stage used to scan the sample laterally over the axial beam.

 

Fig. 1 Fourier domain OCT setup.

Download Full Size | PDF

For convenience of analysis, we set up a z-coordinate for axial scan as shown in Fig. 1 where the origin is set on the zero optical path difference with respect to the reference mirror, M1. Assuming that absorption is weak or negligible, the frequency measurements captured from OSA are in the wavelength domain given as follows [4]:

The first and second terms in Eq. (1) represent the spectrum of the reflection from the reference mirror and the mutual interference among subsurface features in the sample, respectively. The third term contains the useful subsurface information a(z) of the sample. The second term and the third term in Eq. (1) can be made non-overlapped in the reconstructed axial scan by positioning appropriately the reference mirror. The basic problem of OCT is to reconstruct a(z) from the known frequency measurements of S(λ) and G(λ).

In the conventional FD-OCT method [4], the frequency measurements are converted to the frequency domain by using κ = 2π/λ, then they are interpolated and resampled on a uniformly spaced frequency grid with the same number of measurements. Thereafter, the depth profile is obtained from the fast Fourier transform (FFT) of the resampled measurements. This principle clearly shows that the resolution of the conventional FD-OCT method is determined by the coherence length of G(λ). In practice, FWHM (full width at half maximum) of G(λ) is used to characterize the resolution of the conventional FD-OCT method. Another consequence of the use of FFT is that the imaging depth of the conventional FD-OCT method is limited by the spectral resolution of the spectrum meter used to obtain frequency measurements [4].

The proposed method in this paper is different from the conventional FD-OCT method. Although it still uses the frequency measurements, it discretizes S(λ) in Eq. (1) on a uniform spaced wavelength grid of the OSA, which is given by (λ₁, λ₂,..., λ_N) × (z₁, z₂,..., z_M) where N denotes the dimension of frequency measurements. The depth profile is discretized on a grid given by [z₁, z₂,..., z_M] where M is the number of subsurface features of concern. Then, the depth profile of the subsurface features to be reconstructed is related to the frequency measurements by

The problem of OCT imaging is now stated as an algebraic reconstruction problem which is to obtain the unknowns in x from the measurements S(λ_k) and G(λ_k) (k = 1, 2,..., N). It uses the frequency measurements obtained in wavelength domain, and it does not need the resampling process in the reconstruction. This is the same as in [7], but the proposed method is different from the one in [7] in that solving Eq. (2) is actually performing deconvolution simultaneously. Thus, the proposed method improves the resolution regardless of FWHM of G(λ). Moreover, the conventional ART method is to solve the equation by considering an optimization problem: $\underset{x}{\text{min}}{\Vert \mathbf{x}\Vert }_2$ subject to Ax=b where || · ||₂ denotes the 2-norm of a vector.

Although the conventional ART method would give better resolution, an overdetermined equation Eq. (2) needs to solve, provided no a priori knowledge on x such as sparsity is available. This implies N ≥ M. The conventional ART method actually tries to minimize the total mean-square-error occurring to all the elements of x. However, as stated above, the OCT measurements are actually coming from a layer structure or discrete scatterers or mixture of both. Thus, x in Eq. (2) is essentially sparse. In this sense, ℓ₁-optimization of x offers an advantages of using much few number of measurements in b, because one only needs to determine those non-zero elements in x. This is feasible as shown by Restricted Isometry Property (RIP) condition in the compressed sensing theory [14]. The RIP condition specifies N ≥ log(M)μK where K is the sparsity of x and μ represents the coherence of A. The method proposed in this paper is based on the ℓ₁-optimization and it obtains the reconstruction of x as follows:

3. Performance Evaluation of the Algebraic OCT Method Based on ℓ₁-optimization

In this section, we evaluate by simulation and experiment the performance of the proposed method in terms of imaging depth, computational load, and the reconstruction error in the presence of measurement noise. The performance is benchmarked to those of the conventional FD-OCT method and the conventional ART method described in the last section.

In the subsequent simulation studies, the low coherence light is simulated with a Gaussian spectrum centered at 1550 nm with FWHM linewidth of 50 nm. The coherence length is 21 μm. This choice corresponds to the output spectrum of the light source (SLED 1550, Denselight) used in experiments reported in the later sections. In the following experiments, the frequency measurements of S(λ) and G(λ) are captured by an optical spectrum analyzer (ANDO6317B, ANDO) with a tunable spectral resolution that can go as small as 0.01 nm. Without loss of generality, only one axial scan is considered in the following simulation and experiment studies in this section.

3.1. Imaging Depth

Imaging depth is an important performance factor for frequency domain OCT. It refers to the maximum depth of the features that the OCT is capable of imaging. For the conventional FD-OCT method, the imaging depth is mainly determined by the spectral resolution of the spectrum meter [4]. Based on Shannon-Nyquist sampling theorem, the maximum imaging depth can be derived as follows [4]:

In contrast, the proposed method employs ℓ₁-optimization to reconstruct sparse data. As shown in the recent development of compressed sensing theory [?], a sparse signal can still be recovered beyond the limit imposed by Shannon-Nyquist sampling theorem. This means that the imaging depth is no longer limited to Eq. (7).

First, we examine by simulation the effect of the spectral resolution of OSA on the imaging depth. In simulation, we simulate the subsurface feature by using a mirror as the sample in Fig. 1. We set the spectral resolution of the OSA at a value, says 0.2 nm. Then the mirror is moved from z = 0 mm to z = 5 mm in a step size of Δz = 1μm. Thus, we discretize a(z) over the axial range from 0 mm to 5 mm with discretization interval of Δz = 1μm. Referring to Eq.(2), M = 5000. At the k-th step of motion of the mirror, the position of the mirror is denoted by z_k. That is, a(z_k) = 1 and a(z) = 0 for all z ≠ z_k. The frequency measurements at this step are taken at the wavelengths denoted by ${{\lambda }^{\prime }}_i={{\lambda }^{\prime }}_0+i\Delta \lambda$ for i = 1, 2,..., 771 where ${{\lambda }^{\prime }}_0=1473\hspace{0.17em}\text{nm}$. The number of 771 is the number of frequency measurements and it is returned by the OSA when its resolution is set at 0.2 nm. To use Eq. (5) to reconstruct the position of the mirror, we chose in a evenly-random manner N = 150 samples among the 771 frequency measurements. Denote the corresponding wavelengths by λ_i for (i = 1, 2,..., N = 150) for Eq. (3) and Eq. (4), without loss of generality. The reconstructed x is obtained by performing ℓ₁ optimization in Eq. (5). The optimization can be done by many well-developed algorithms such as the CVX matlab module [15]. The reconstructed position of this step is the entry position of the nonzero element in x. After all motion steps are finished, the imaging depth is identified as the maximum value of the positions which are correctly recovered. We repeat the above process for other values of the spectral resolution of the OSA. The imaging depth is plotted against the spectral resolution in Fig. 2 which also includes the imaging depth for the conventional FD-OCT method as calculated by Eq. (7). This figure shows that the imaging depth of the proposed method is no longer limited by the spectral resolution of the spectrum meter used.

 

Fig. 2 Simulated imaging range of the proposed method.

Download Full Size | PDF

Next, experiments are carried out to examine the imaging depth versus the spectral resolution of OSA. Again, a mirror is used as the sample. First, we set the spectral resolution of the OSA to 0.2 nm which gives the imaging depth for the conventional FD-OCT method as 3 mm according to Eq. (7). We put the mirror at an initial position where the optical path difference is about 0.4 mm with respect to the reference mirror. The choice of this initial position does not carry any special meaning and it is purely for convenience to see the imaging depth going through the limit at 3 mm. Then the mirror is moved away in steps to have bigger optical path difference between the mirror and the reference mirror. For each step, the mirror moves by 1 mm. In experiment, the frequency measurements are captured by the OSA at each mirror position. The axial scan reconstruction of the position of the mirror is done respectively by using the proposed method and the conventional FD-OCT method. In application of the proposed method, we discretize a(z) by 5001 points over a range of [0 mm, 5 mm] in optical path, that is M = 5001 and solve Eq. (6) with empirically determined ε = 0.5 × 10⁻³. The frequency measurements are captured in experiment at 771 wavelengths from 1473 nm to 1627 nm for the conventional FD-OCT method, out of which 150 measurements are chosen evenly randomly from 1513 nm and 1573 nm for the proposed method. The reason for choosing the center portion of the optical spectrum is that the signal there has a higher SNR. For the conventional FD-OCT method, the whole set of 771 measurements with no zero padding is used.

Figure 3 show the axial scan reconstruction obtained in experiment by the conventional FD-OCT method and the proposed method, respectively. For both figures, the magnitude of the reconstructed a(z) is plotted in logarithmic scale against the optical path difference recorded by the motion stage driving the sample mirror. Comparing the results obtained from the two methods, a clear drop trend of the magnitude is observed for both methods. This is attributed to the sensitivity roll off. However, experiment results show clearly that the mirror position reconstructed by the conventional FD-OCT method bounces back to less than 3 mm when the mirror goes actually beyond 3 mm in its third move marked by Step 3. In the Fig. 3(a), the reconstructed step profile at Step 3 has also a superimposed background or hump. This is similar to the effect observed after interpolation is performed on spectral interferogram with large time delay [6], except that in this case, the hump is attributed to the cascaded effect of interpolation error on aliased data. This indicates the imaging depth of the FD-OCT method. In contrast, the proposed method continues delivering the correct reconstruction results even when the mirror goes beyond 3 mm.

 

Fig. 3 Experiment results on the imaging depth of various OCT methods.

Download Full Size | PDF

Further studies which involve moving of the mirror further away and using different OSA resolutions are also conducted in experiment. Similar results are obtained, although they are not reported here. All these results confirm that the imaging depth of the proposed method is not limited by the spectral resolution of spectrum meter while the conventional FD-OCT method is. In practice, the achievable imaging depth depends on the reconstruction condition posed by RIP and the signal-to-noise ratio of the measurements. The analytical bound on the imaging depth of the proposed method is an interesting topic and the related work is still going on.

3.2. Computational Load

In this section, we evaluate the computational performance of the proposed method, the conventional ART method, and the conventional FD-OCT method. As evident from Eq. (2), the computational load of the proposed method is determined by the dimensions of A. In turn, the dimensions of matrix A are determined by the values of N and M. They denote respectively the number of wavelengths at which the frequency measurements are obtained, and the number of points used to discretize a(z). They basically determine the computational complexity.

Considering the RIP condition and the computational load of ℓ₁ optimization [16], we can estimate the computational load of the proposed method as O(M^3/2log²(M)μK). For the conventional ART method, SVD algorithms are usually employed for fast computational which gives the computational load of 6M³ [17]. It is easy to verify that O(M^3/2log²(M)μK) << 6M³. With this, we anticipate that the computational load of the proposed method is greatly reduced compared to the conventional ART method. In order to verify this, we perform numerical simulations to compare the computational time of the proposed method and the conventional ART method on a Pentium-4 computer with 3.0 GHz processor and 1 GB RAM, and the software used is MATLAB 7.4 under Microsoft Windows XP 32 bit.

We consider the scenario where the sample in Fig. 1 is simply a mirror. Without loss of generality, the mirror is placed at z_k = 300μm. This implies that a(z) = 0 except when z = z_k We let the spectrum G(λ) take a rectangular shape with wavelength range [600 nm, 1600 nm] and simulate S(λ_i) (i = 1, 2...2000) using Eq. (1). We chose a broader wavelength range here just in order to make sure of rank(A) = M. This is in favor of the conventional ART method in a sense that it would yield the exact and correct a(z). For the conventional ART method, we let M = 500 and find the value of corresponding N such that the condition rank(A) = M is satisfied. Then we solve the inverse problem Ax = b using the pseudo inverse function in MATLAB and record the computational time. For the proposed method, we let M = 500 and N = 20. As a reference, the computational time for the conventional FD-OCT method is also evaluated. This consists of the time needed for MATLAB to spline interpolate all the simulated G(λ) to G(κ) and the time of performing the FFT algorithm.

Table 1 summaries computational complexity of various methods and the computational time taken in experiments for various reconstruction methods. As shown, the proposed method does have a computational advantage over the conventional ART method.

Table 1. Comparison of computational time among the proposed method, the conventional ART method, and the conventional FD-OCT method

View Table

3.3. Signal to Noise Ratio (SNR)

In practice, the frequency measurements from the OSA are polluted by noises. For this case, the proposed method is implemented by making use of the constraint ||Ax – b||₂ < ε where ε is a measure of the noise power in the measurements and the value of ε is usually determined empirically. In this section, simulation studies are carried out by using Eq. (6) for the proposed method. The reconstruction error is chosen as the performance indicator, and it is examined by varying the number of measurements N under different SNR conditions.

Again considering the single mirror sample, the frequency measurements are simulated by using Eq. (1) where R_r = 2.5, [z₁ = 1μm, z₂ = 2μm, ..., z_k = 300μm, ..., z_M = 1000μm], M = 1000, and a(z = z₃00) = 1. The spectrum is simulated over the range of wavelengths from 1473 nm to 1627 nm at 771 wavelengths. The spectrum of the light source is assumed as a Gaussian over the range. Then additive Gaussian noise with zero mean and variance of ${\delta }_n^2$ is added to the spectrum generated from Eq. (2) to simulate the noisy frequency measurements. In the simulation, the variance ${\delta }_n^2$ is tuned to simulate different SNR (signal-to-noise ratio) where SNR is defined by 20 log (||x||₂ / ||n||₂) with n denoting a realization of the Gaussian noise. For each SNR value, N samples are taken from the noisy frequency measurements randomly in an even distribution. Then the ℓ₂ norm of the difference between the reconstructed x and their true values from a(z) is calculated and it is denoted as the reconstruction error. Figure 4 plots the reconstruction error against SNR for different number of frequency measurements used in Eq. (6). It is quite straightforward to see that reconstruction errors become bigger in the presence of a big noise level, given the number of frequency measurements. On the other hand, for a given SNR level, it shows that the reconstruction can be made more robust to noise by using more frequency measurements. And the dependence of the number of frequency measurements becomes lesser when a big enough value of N is chosen.

 

Fig. 4 The reconstruction error as a function of SNR.

Download Full Size | PDF

4. Onion cell imaging

In this section, the proposed method is demonstrated on an onion sample. The sample is prepared by placing a freshly sliced onion on the microscope glass slide which has 1 mm thickness. Then a thin layer of index-matching oil (Refractive index = 1.488) is applied on the top surface on the onion skin, and another microscope glass slide with thickness of 100μm is used to cover the onion sample. All glass are made from BK7. Such a packaging is to remove the reflection from the bottom surface of covering glass slide, and the top surface of the covering glass slide can be used as a reference plane which does not contribute too strong reflection.

For this application, the proposed method is implemented on a common path interferometer which is similar to the one in Fig. 1 but reference mirror M1 is removed and the collimator is replaced by a microscope objective lens (Numerical aperature=0.1). The interference is formed by the light reflected respectively from the subsurface features in the onion and the top surface of the covering glass slide. There are several reasons of using common path interferometer here. First, the reference light from the top surface of the covering glass slide is much weaker than a reference mirror, a better contrast is obtained because the non-varying component (DC) of the reflected spectrum in Eq. (1) is reduced [18]. Second, the common path interferometer provides a convenient setup for imaging. Third, a common path interferometer offer better stability in the presence of fluctuations in two interferometer arms.

In the experiment, the packaged sample is placed under the microscope lens, and the position of the sample is adjusted vertically so that the top surface of the onion is at the focal plane of the microscope lens. The whole sample sits on a nanostage (17NST101, Melles Griot) that is responsible for scanning the sample laterally under the microscope objective lens in steps of 1μm over a range of 1 mm by 1 mm. For each axial scan at a lateral position, a set of raw frequency measurements S(λ) at wavelengths over a span of 154 nm centered at 1550 nm are obtained by the OSA with its spectral resolution set at 0.4 nm. The corresponding wavelengths are denoted by ${{\lambda }^{\prime }}_1=1473\hspace{0.17em}\text{nm}$, ... ${{\lambda }^{\prime }}_{N^{\prime }}=1627\hspace{0.17em}\text{nm}$ with Δλ = 0.4 nm and N′ = 386.

To apply the proposed method to obtain a cross-section depth profile of the sample, we discretize a(z) with the discretization interval of 4μm over the axial range from 0 mm to 2 mm. That is, z₁ = 0μm, z₂ = 4μm, ..., z_M = 2000μm, and M = 501. We measure the output spectrum G(λ) of the light source at the wavelengths ${{\lambda }^{\prime }}_1,\dots ,{{\lambda }^{\prime }}_{N^{\prime }}$. Among N′ frequency measurements both for S(λ) and for G(λ), 280 frequency measurements are taken at the wavelengths which are chosen randomly in an even distribution, that is N=280. Using these measurements, we obtain A and b as in Eq. (3) and Eq. (4). Then we reconstruct each axial scan through solving the optimization problem given in Eq. (6). The image obtained by the proposed method is given in Fig. 5(a) where the vertical represents the depth but in logarithmic scale.

 

Fig. 5 Cross-section images obtained by various methods using an OSA spectral resolution of 0.4 nm (The arrow in the figures indicates the imaging depth determined by Eq. (7); the vertical axis represents the depth and the range shown is 1 mm; the horizon axis represents lateral position and the range is 1 mm; the intensity in all the images is in logarithmic scale). (a) The proposed method; (b) The conventional FD-OCT method; (c) The non-uniform Fourier transform method [19]; (d) The conventional ART method.

Download Full Size | PDF

Next we compare the reconstruction results with those of the conventional FD-OCT method, the non-uniform FFT method [19], and the conventional ART method. In these methods, we use the full set of the raw frequency measurements at N′ = 386 wavelengths given at ${{\lambda }^{\prime }}_1=1473\hspace{0.17em}\text{nm}$,..., ${{\lambda }^{\prime }}_{N^{\prime }}=1627\hspace{0.17em}\text{nm}$. For the conventional FD-OCT method and the non-uniform FFT method, the frequency measurements are first transformed from wavelength domain to the frequency domain by using spline interpolation and even sampling. For the conventional ART method, the ℓ₂ optimization is performed. The reconstructed cross-section images of using these method are shown respectively in Fig. 5(b), Fig. 5(c), and Fig. 5(d). To compare the imaging depth, we also put an arrow to depict the imaging depth of 1.5 mm as determined by Eq. (7) based on OSA resolution of 0.4 nm. But for all the images in Fig. 5, we crop the same top part without any meaningful data, and show only a part with depth range of 1 mm.

From Fig. 5(b) and Fig. 5(c), it is evident that the conventional FD-OCT method and the nonuniform FFT method suffer from aliasing and they are not able to decern layers of the onion skin when the feature is beyond the imaging depth denoted by the arrows next to the images. In contrast, the proposed method can detect the features that are missing in the Fig. 5(b). In addition, the image obtained from the proposed method also shows enhancement in resolution compared to the ones from the non-uniform FFT method and the conventional ART method.

5. Conclusions

This paper presents an OCT method that uses frequency measurements in wavelength domain to achieve tomography imaging with high resolution and longer imaging depth. Its principle is motivated by exploring the sparsity of the measurerants in OCT. The task of OCT is formulated into an algebraic reconstruction based on ℓ₁-optimization. Inheriting the merit of ART in achieving high resolution, the proposed method overcomes the bottleneck of ART method by greatly reducing its computational load. Simulation and experiment studies have demonstrated the merit performance of the proposed method. The proposed method is expected to be further develop as a fast effective method for high resolution and longer imaging depth subsurface imaging.