|
1.IntroductionOptical coherence tomography (OCT)1 is a subsurface imaging technique that has gained popularity in many applications.2–4 Since it was invented, time-domain OCT (TD-OCT) and Fourier domain OCT (FD-OCT) methods have been developed, among which the FD-OCT method has been recognized as the superior method due to its better sensitivity and higher speed.5,6 In the FD-OCT method, the spectral interferogram measurements are regarded as the Fourier transform of the depth profile. However, the mutual interference of the optical field scattered from all the sample features also contributes to a part of the spectral interferogram. This gives rise to autocorrelation noise (ACN) in the FD-OCT reconstruction results. As the resultant ACN overlaps with the actual subsurface features in OCT images, the presence of ACN adds ambiguities when interpreting the OCT image.7 This is especially so in the case of full range detection or a sample consisting of several strong scattering layers with separations such that the ACN appears at the same location as the features of the depth profile. Thus, ACN is a major concern to the OCT community. In order to eliminate the ACN, several approaches have been proposed. The simplest approach is to move the reference mirror away from the sample surface in terms of optical path. This approach reduces the influence of ACN but does not eliminate it. The second approach is to remove the mutual interference by using two or more phase-shifted spectral interferograms.8,9 This approach eliminates ACN but requires multiple data acquisition. Furthermore, additional phase-shifted errors would complicate the matter.10 The third approach6 is to avoid ACN by reconstructing the optical fields scattered from the sample subsurface features and the reference mirror. Since the optical fields and the depth profile are direct Fourier transform pairs, the depth profile is reconstructed by Fourier transform without ACN. To uniquely reconstruct the optical fields, however, the minimum phase condition is required. This limits the class of samples in practice. The forth kind of approache is based on the cepstrum signal processing.7 The cepstrum of the optical fields is reconstructed first, from which the depth profile is obtained by Fourier transform. This cepstrum approach is applicable for a limited class of samples, as the approach requires the non-overlap of the cepstrums. In this paper, the removal of ACN is carried out by exploiting the sparsity of the depth profile. The optical fields scattered from the subsurface features are regarded as the responses of a sparse finite impulse response (FIR) filter. Then, the OCT problem of reconstructing the depth profile from the spectral interferogram is formulated into one of identifying the parameters of a sparse FIR filter subject to its known spectrum. This problem can be solved by -norm optimization by virtue of the sparsity. The modeling of scattered optical fields by a sparse FIR filter helps to remove the ACN. This is because the OCT reconstruction is carried out on the basis of a precise expression of spectral interferogram in terms of the depth profile and the mutual interference has already been taken into account in the reconstruction. But the challenge is that, in order to identify the sparse parameters of the FIR filter, -norm optimization should be done subject to a quadratic constraint that relates the spectral interferogram to the depth profile. To overcome the challenge, an iterative soft-thresholding algorithm of -norm optimization is developed. The performance of the proposed reconstruction method is evaluated by experiments, and the results are benchmarked against those of the cepstrum method and the minimum phase method. It is shown that the proposed method suppresses the ACN effectively. 2.PrincipleThe principle of the proposed method is illustrated by using the FD-OCT setup as shown in Fig. 1. The setup uses a fiber-coupled, low-coherence light source. The light beam goes through the fiber circulator and collimator, and it is split into the reference beam and the sample beam by the beam splitter. The sample beam penetrates the sample along the axial direction and it is reflected by the subsurface features at different depths. The reference beam is reflected at the reference mirror. Eventually, the light reflected from all the sample features and the reference mirror are coupled back into the fiber by the collimator. The coupled beams are then transmitted into an optical spectrum analyzer (OSA) that is usually designed with a fiber-coupled input. At OSA, the sample and light beams interfere and the raw spectral interferogram in the lambda domain is obtained. The interferogram would be used for constructing a single depth profile of the sample. Finally, the OCT image of the sample is constructed by assembling all the A-scan depth profiles with the assistance of the motion stage used to scan the sample laterally over the axial beam. Without loss of generality, the following discussion only considers recovering a single A-scan depth profile of OCT. In Fig. 1, we set up a -coordinate with the origin set as the position given by the zero optical path difference relative to the reference mirror. The coordinate system is used to indicate the positions of the sample features relative to the reference mirror. Furthermore, we introduce the following denotations. denotes the optical field of low coherence light with denoting time instant, denotes the optical field reflected from the reference mirror, the optical field from the sample, the output optical spectrum of the low-coherence light source with denoting the light wavelength, the reflectivity of the reference mirror, and the depth profile to reconstruct. From the system theory point of view, we can model the optical field of the low-coherence light as the output of a filter with spectrum given by , that is, , where is the transfer function of the filter, denotes a unit-time delay operator ( i.e., ), and is an identically independent random sequence with unit power spectrum. Hence, the spectrum of the filter is given by , where denotes the light speed. Similarly, the scattered optical field from the sample can be modeled as the output of a system with its input as the low-coherence light. Specifically, . Assume that absorption is weak or negligible. If we discretize the depth profile into where is an integer and is the axial sampling rate, the scattered optical field can be considered as the sum of the input white light reflected from the scatters at the positions given by with the reflectance of . In this sense, the unit-time delay can be considered as and the transfer function is given by with . Then, the spectral interferogram is related to the sample features by the following expression: where denotes the expectation operator and the spectrum is obtained by substituting into the transfer functions and . It is worth pointing out that Eq. (1) actually tallies with the traditional expression of spectral interferogram in FD-OCT.11 Expanding the equation, the ACN term refers to , which represents the mutual interference of all the optical fields scattered by the features of the sample.12 However, the traditional expression of is in wavenumber domain in FD-OCT.13 Here we express the spectral interferogram in wavelength domain for the reason that the spectral measurements are usually obtained in terms of wavelength. This has the additional advantage of avoiding the interpolation operation commonly needed in FD-OCT.With the above models, the task of OCT reconstruction can be formulated as one of identifying the parameters of filter based on the available spectral measurements given by . In order to eliminate the ACN, the parameter identification should be done subject to the precise expression of in Eq. (1). In practice, the subsurface features observed by OCT are of sparse scatters. In view of this sparsity, we can consider a FIR filter of , where denotes the order of the filter.14 The FIR representation means that there is only a finite number of optical path differences between the scattered optical field from the sample and the optical field reflected from the reference mirror. We can vary the order of the filter to consider different sparsity of the subsurface features. Assume the spectral measurements are recorded by an OSA at wavelengths denoted by . The spectral measurements available at an OSA are thus expressed as where with .By considering the sparsity, the OCT reconstruction of the depth profile can be solved by the following -norm optimization subject to quadratic constraint: where a is a vector with its th element and its -norm is defined by .In practice, the spectral interferogram obtained in the experiment is always corrupted by noise. This leads us to consider and solve the proposed method in the following mixed -optimization:15 where denotes the noise level.The Lagrangian formulation of the optimization problem of Eq. (5) can be formulated as follows: where , and μ is a regularization parameter that is required to control how stringent the quadratic constraint is μ performs the same role as does in Eq. (5).Equation (6) can be solved by an iterative soft thresholding method.16 For this method, each iteration performs three main procedures: synthesizing, updating, and thresholding. In the procedure of synthesizing, an estimated value of is substituted into to obtain the synthesized spectral interferogram denoted by at th iteration. This procedure is denoted by applying a synthesis operator on , that is, where is a vector with entries given by . In the next procedure, the estimate of is updated by applying an adjunct operator on the innovation given by . The updated estimate is given by where denotes a vector with entries given by , and is the gain used to control the updating procedure.In the last procedure, the elements of are thresholded by using a soft thresholding function, , given by where is the function to be soft thresholded and μ is the threshold level. The choice of μ is experimentally determined, but it can also be determined by various criteria like the Akaike information criterion17 or the Bayesian information criterion.18It can be seen from the above schematic of soft thresholding method that the key challenge in this method is to construct the adjunct operator in the updating procedure. By taking advantage of the FIR representation introduced in this paper, we propose one realization of as follows: where , , and is the th element of .To summarize, the proposed OCT reconstruction method characterized in Eq. (4) or Eq. (5) can be implemented by the Algorithm 1. Algorithm 1OCT reconstruction with ACN removal.
3.Performance of the Proposed MethodIn this section, the performance of the proposed method is evaluated in terms of ACN removal capability, resolution, and computing demand. The performance is also compared to those of the FD-OCT method, the minimum phase method, and the cepstrum method. For a fair comparison, all the methods perform the reconstruction on the same set of the A-scan data, and the reconstruction is carried out on an Intel Core(TM)2 Quad CPU computer with 2.66 GHz processor and 4 GB RAM, and the software used is MATLAB 7.4 under 32 bits Microsoft Windows 7 operating system. We begin with the evaluation of the effectiveness of the proposed method in removing the ACN. This is demonstrated and benchmarked by applying the method in two experiments to obtain OCT images. The first experiment uses a well-controlled sample formed by a two-mirror system. The second experiment applies the OCT method to image an onion sample. In all the experiments, the OCT is implemented on an interferometer setup as shown in Fig. 1. But for the second experiment, an additional microscope objective lens () is inserted between the collimator and the sample. This is to help collect the light reflected by the sample better. The light source used in experiments is a broadband Superluminescent diode (SLED 1550, Denselight) with the center wavelength at 1550 nm and the coherence length of 21 μm in air. The optical spectrum analyzer (ANDO6317B, ANDO) records the spectral interferogram with resolution 0.1 nm over a span of 154 nm centered at 1550 nm. It returns spectral measurements consisting of 1541 readings. In the first experiment, a two-mirror system given in Fig. 2 is used as the sample in Fig. 1. The two mirrors labeled by and simulate the two interfaces of the sample. Mirror is carried by a motorized nanostage (Melles Griot 17NST101) and remains stationary. In this case, the depth profile consists of two spikes. Different depth profiles can be simulated by moving mirror . To carry out the experiment, the positions of mirrors , , and the reference mirrors are calibrated. The calibration of and is done with the reference mirror blocked and the OSA replaced with an optical power meter. The position of mirror is tuned until a maximum interference intensity is read at the optical power meter. The calibration of the reference mirror involves three steps. First is to block the mirror and unblock the reference mirror. Then, the position of the reference mirror is tuned until a maximum interference intensity is read at the optical power meter. Third, the position of the reference mirror is shifted until it is “in front” of both and by 120 μm. This gives the usual arrangement in OCT setup.13 After the calibration, all the mirrors are unblocked and the optical power meter is replaced with the OSA. The OCT imaging is first performed to reconstruct an initial depth profile characterized by a position of mirror . Then, the OCT imaging is performed on different depth profiles simulated by moving mirror in steps. Finally, all the reconstructed depth profiles are assembled to give one OCT image. In each step, mirror is moved by 1 μm, and it is moved a total of 20 steps. At each step, the spectral resolution of the OSA is set to 0.1 nm, and a raw spectral interferogram is recorded over a span of 154 nm centered at wavelength 1550 nm. After the spectral interferogram is obtained, the depth profiles are reconstructed by different reconstruction methods including the proposed method, the Fourier transform–based FD-OCT method, the minimum phase method, and the cepstrum method. For the proposed method, the raw spectral interferogram is input into the Algorithm 1 with setting of , , , and . The reconstructed OCT image is shown in Fig. 3(d). For the remaining three methods, the raw spectral interferogram obtained at each mirror step is first spline interpolated to obtain a spectral interferogram that is uniformly sampled on wavenumber. For the Fourier transform–based FD-OCT method, each resampled interferogram is Fourier transformed using FFT, and all the depth profiles are assembled to give the depth profile image in Fig. 3(a). For the minimum phase method, the depth profile at each mirror step is reconstructed by using 1000 iterations of the algorithm proposed in the paper by Ozcan6 and the obtained depth profile image is displayed in Fig. 3(b). Figure 3(c) gives the corresponding result obtained by the cepstrum method.6 For all the images in Fig. 3, the vertical axis denotes the dimension measurement over a range of 180 μm and the horizontal axis denotes the different positions of mirror . Comparing the images, the first observation is that the image obtained by the proposed method appears sharper than the others. Second, the ACN appears clearly in the top part of the image obtained by the Fourier transform–based FD-OCT method. They are suppressed in the image obtained by the minimum phase method and the cepstrum method, and they are totally eliminated in the image obtained by our proposed method. To see more closely, Fig. 4 shows a single depth profile obtained, respectively by the four methods. It is evident that an autocorrelation peak appears at μm position for the FD-OCT method, the minimum phase method, and the cepstrum method. But the ACN is eliminated by the proposed method. Summarizing the experimental results, it is confirmed that the proposed method can effectively eliminate the ACN. It is also noted that one of the feature peaks reconstructed by the proposed method is slightly lower than the one obtained by the FD-OCT method. This is attributed to the soft thresholding algorithm, which tends to underestimate the actual value. In the second experiment, the OCT imaging of a piece of onion skin is considered. Onion skin consists of hexagonal onion cell at the epithelium layer.19 The onion cell membrane forms the scatters which act as features in the depth profile. Thus, onion samples are often used to verify the OCT methods in the literature.20 In this paper, the onion sample is prepared by placing a freshly sliced onion on the microscope glass slide having 1 mm thickness. Then, a thin layer of index-matching oil () is applied on the top surface of the onion skin, and another microscope glass slide with thickness of 100 μm is used to cover the onion sample. All glasses 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. In the experiment, the packaged sample sitting on a nanostage (17NST101, Melles Griot) 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 sample is scanned laterally over the microscope objective lens in steps of 1 μm over a range of 1 mm. At each lateral position, a set of A-scan spectral interferogram is measured over wavelengths with a span of 154 nm centered at 1550 nm. The measurements are obtained with the spectral resolution of OSA set at 0.1 nm. After the A-scan spectral interferograms at all the lateral positions are obtained, we proceed to reconstruct the images by the proposed method, the FD-OCT method, the minimum phase method,6 and the cepstrum method,7 respectively. The reconstruction is done line by line, each line denoting one A-scan of depth profile. Then, the reconstructed cross-section OCT image resulting from the lateral scan is obtained by assembling all the A-scan depth profiles according to their lateral position. For each reconstruction method, the settings are the same for individual A-scans in the cross-section image. For the proposed method, the settings used in the Algorithm 1 are chosen as , , , and . The cross-section image is shown in Fig. 5(a). For the other three methods, the spectral interferograms are first spline interpolated line by line to obtain another set of spectral interferograms, each of which is uniformly sampled on wavenumber. Then, Fourier transform is applied for the FD-OCT method, the algorithm described in the paper6 with 1000 iterations for the minimum phase method, and the procedure described in the paper7 is applied for the cepstrum method, respectively. The cross-section images reconstructed by the three methods are given, respectively in Fig. 5(b), 5(c), and 5(d). Comparing the images in Fig. 5, it can be seen that there is a large area of artifacts on the top part of reconstructed image obtained from the Fourier transform–based FD-OCT. These artifacts are suppressed in the reconstructed cross-section images resulting from the proposed method, the minimum phase method, and the cepstrum method. However, since the artifacts are due to not only ACN but also other noises in practice, some left-over artifacts are still visible in the top parts of Fig. 5(a), 5(c), and 5(d). The second observation is the vertical artifacts appearing in Fig. 5(c) and 5(d). This may be due to two reasons. One is that the condition of non-overlap of the cepstrums is not satisfied for the depth profile at the corresponding positions. The second is the spurious artifacts as mentioned in the literature.7 But the vertical artifacts are totally removed in Fig. 5(a), which is obtained by our proposed method. Another observation when comparing the four methods is that the reconstructed image obtained by using our proposed method exhibits better clarity. This is because the method actually is based on an algebra reconstruction scheme which factored the point spread function in the reconstruction. In the next part of this section, the resolution and computing demand of the various methods are compared. To evaluate the resolution, a quantitative measure called resolution gain (RG) is deployed here. This measure is commonly used in ultrasound imaging.21 It is calculated as follows. First, we calculate the autocorrelation function of one A-scan with respect to the depth position. Then, the autocorrelation function is normalized, with its maximum normalized to 1. The number of positions where the autocorrelation function is greater than 0.75 were counted to give the autocorrelation value. Second, we calculate the autocorrelation values, respectively on the A-scans of the proposed method, FD-OCT method, the minimum phase method, and the cepstrum method. Next, we use the one of the FD-OCT method as a reference, and calculate the ratio between this reference value and the autocorrelation values of other methods. This ratio is called the resolution gain. Without lost of generality, we consider the scenario of a mirror surface as the sample. The spectrum measurement is first simulated according to Eq. (1). Then, we apply respective method to reconstruct one A-scan, based on which we calculate the resolution gain for the method. The obtained resolution gains are summarized in Table 1. It is clear that the proposed method achieves the highest resolution gain. This implies that proposed method achieves highest resolution among various method. We attribute this to the inherent deconvolution properties of solving Eq. (5). Table 1Different performance measures.
The computing demand of various methods is evaluated by the computational time that is required by various methods to finish the reconstruction from one set of spectrum measurement. For each reconstruction method, we use the timer function in MATLAB to record the time needed. For the FD-OCT method, the minimum phase method, and the cepstrum method, the computational time also includes the time for spline interpolation of all the spectrum measurement and the time for FFT transform. The third column of Table 1 lists the computational time needed by the four methods to obtain the onion images in Fig. 5. It shows that the proposed method requires more time to perform the reconstruction. This is understandable, as the proposed method is based on an iterative reconstruction algorithm. This may mean that the proposed method achieves the high resolution at expense of computational time. The overall performance of the proposed method and its comparison to other methods is summarized in Table 1. 4.ConclusionsA reconstruction method is proposed for the removal of the ACN in the reconstruction of depth profile by using OCT. The method is based on the modeling of the scattered optical fields in OCT by a sparse FIR filter. The depth profile reconstruction is formulated as one parameter identification problem of a sparse FIR filter given the spectrum response. The parameter identification is carried out via an -norm optimization, but with a quadratic constraint. To overcome the challenge, an iterative soft thresholding scheme is developed in the optimization to identify the filter parameters. Experiments and comparison studies show that the proposed method is able to remove the ACN from the OCT images effectively. ReferencesD. Huanget al.,
“Optical coherence tomography,”
Science, 254
(5035), 1178
–1181
(1991). http://dx.doi.org/10.1126/science.1957169 SCIEAS 0036-8075 Google Scholar
O. P. BrunoJ. Chaubell,
“Inverse scattering problem for optical coherence tomography,”
Opt. Lett., 28
(21), 2049
–2051
(2003). http://dx.doi.org/10.1364/OL.28.002049 OPLEDP 0146-9592 Google Scholar
P. E. Andersenet al.,
“Advanced modelling of optical coherence tomography systems,”
Phys. Med. Biol., 49
(7), 1307
–1327
(2004). http://dx.doi.org/10.1088/0031-9155/49/7/017 PHMBA7 0031-9155 Google Scholar
D. Stifter,
“Beyond biomedicine: a review of alternative applications and developments for optical coherence tomography,”
Appl. Phys. B, 88
(3), 337
–357
(2007). Google Scholar
R. LeitgebC. K. HitzenbergerA. F. Fercher,
“Performance of fourier domain vs. time domain optical coherence tomography,”
Opt. Express, 11
(8), 889
–894
(2003). http://dx.doi.org/10.1364/OE.11.000889 OPEXFF 1094-4087 Google Scholar
A. OzcanM. J. F. DigonnetG. S. Kino,
“Minimum-phase-function-based processing in frequency-domain optical coherence tomography systems,”
J. Opt. Soc. Am. A, 23
(7), 1669
–1677
(2006). http://dx.doi.org/10.1364/JOSAA.23.001669 JOAOD6 0740-3232 Google Scholar
C. S. Seelamantulaet al.,
“Exact and efficient signal reconstruction in frequency-domain optical-coherence tomography,”
J. Opt. Soc. Am. A, 25
(7), 1762
–1771
(2008). http://dx.doi.org/10.1364/JOSAA.25.001762 JOAOD6 0740-3232 Google Scholar
M. Wojtkowskiet al.,
“Full range complex spectral optical coherence tomography technique in eye imaging,”
Opt. Lett., 27
(16), 1415
–1417
(2002). http://dx.doi.org/10.1364/OL.27.001415 OPLEDP 0146-9592 Google Scholar
A. F. Fercheret al.,
“Complex spectral interferometry OCT,”
Proc. SPIE, 3564 173
–178
(1999). https://doi.org/http://link.aip.org/link/doi/10.1117/12.339152 PSISDG 0277-786X Google Scholar
P. S. HuangS. Zhang,
“Fast three-step phase-shifting algorithm,”
Appl. Opt., 45
(21), 5086
–5091
(2006). http://dx.doi.org/10.1364/AO.45.005086 APOPAI 0003-6935 Google Scholar
A. F. Fercheret al.,
“Optical coherence tomography—principles and applications,”
Rep. Prog. Phys., 66
(2), 239
–303
(2003). http://dx.doi.org/10.1088/0034-4885/66/2/204 RPPHAG 0034-4885 Google Scholar
V. M. Gelikonovet al.,
“Coherent noise compensation in spectral-domain optical coherence tomography,”
Optic. Spectros., 106
(6), 895
–900
(2009). Google Scholar
G. HäuslerM. W. Lindner,
“Coherence radar and spectral radar—new tools for dermatological diagnosis,”
J. Biomed. Opt., 3
(1), 21
–31
(1998). http://dx.doi.org/10.1117/1.429899 JBOPFO 1083-3668 Google Scholar
F. J. Taylor,
“,”
Principles of Signals and Systems, Mcgraw Hill, Inc., New York
(1994). Google Scholar
M. ZibulevskyM. Elad,
“L1-L2 optimization in signal and image processing,”
IEEE Signal Process. Mag., 27
(3), 76
–88
(2010). ISPRE6 1053-5888 Google Scholar
S. Mallet,
“,”
A Wavelet Tour of Signal Processing: The Sparse Way, Academic Press, Burlington, MA
(2009). Google Scholar
H. AkaikeB. N. PetroxF. Caski,
in Proc. of 2nd International Symposium on Information Theory,
267
(1973). Google Scholar
G. Schwarz,
“Estimating the dimension of a model,”
Ann. Statist., 6
(2), 461
–464
(1978). http://dx.doi.org/10.1214/aos/1176344136 ASTSC7 0090-5364 Google Scholar
M. D. Kulkarniet al.,
“Image enhancement in optical coherence tomography using deconvolution,”
Electron. Lett., 33
(16), 1365
–1367
(1997). http://dx.doi.org/10.1049/el:19970913 ELLEAK 0013-5194 Google Scholar
H. L. SeckY. ZhangY. C. Soh,
“Optical coherence tomography by using frequency measurements in wavelength domain,”
Opt. Express, 19
(2), 1324
–1334
(2011). http://dx.doi.org/10.1364/OE.19.001324 OPEXFF 1094-4087 Google Scholar
C. YuC. ZhangL. Xie,
“An envelope signal based deconvolution algorithm for ultrasound imaging,”
Signal Process., 92
(3), 793
–800
(2012). http://dx.doi.org/10.1016/j.sigpro.2011.09.024 SPRODR 0165-1684 Google Scholar
|