Dispersion compensation in Fourier domain optical coherence tomography Tarek A. Al-Saeed,1,* Mohamed Y. Shalaby,2 and Diaa A. Khalil2,3 1
Biomedical Engineering Department, Faculty of Engineering, Helwan University, 1 Sherif Street, Helwan, Cairo, Egypt 2
Electronics and Communications Engineering Department, Faculty of Engineering, Ain-Shams University, 11 Elsarayat Street, Abbassia, Cairo 11517, Egypt 3
Si-Ware Systems Company, 3 Khaled Ibn Al-Waleed Street, Heliopolis, Cairo, Egypt *Corresponding author: [email protected] Received 28 April 2014; revised 25 July 2014; accepted 19 August 2014; posted 3 September 2014 (Doc. ID 210914); published 2 October 2014
In this work, we propose a numerical technique to compensate for errors due to dispersion effects in Fourier domain optical coherence tomography. The proposed technique corrects for errors in depth measurements and resolution loss due to dispersion. The results show that, by using this technique, errors in thickness measurement are reduced from about 5% to less than 0.1% depending on the sample length and the amount of dispersion. Also, an improvement in the resolution from about 50 μm to less than 10 μm is demonstrated. © 2014 Optical Society of America OCIS codes: (110.4500) Optical coherence tomography; (030.1670) Coherent optical effects; (170.1650) Coherence imaging. http://dx.doi.org/10.1364/AO.53.006643
1. Introduction
Optical coherence tomography (OCT) is a competitive in vivo noninvasive imaging technique [1]. Similar to the ultrasound technique, OCT provides A-scan and B-scan images. These images are obtained using the low coherence interferometry principle. Thus, depending on the optical system, the axial and lateral resolution can be of the order of 3–10 μm, which is much better than the resolution of the ultrasound technique. On the other hand, it can provide an in-depth image to a penetration depth up to a few millimeters. In this sense, it fills the gap between ultrasound imaging and confocal microscopy. There are many applications of OCT in medicine. OCT can be used in imaging cellular and tissue microstructure, in vivo retinal imaging, imaging coronary atherosclerosis, imaging laryngeal mucosa, skin imaging, dental tissue imaging, and corneal 1559-128X/14/296643-11$15.00/0 © 2014 Optical Society of America
thickness measurement. It can also be used in monitoring the growth and development of biological tissues over time. Finally, endoscopic OCT can be used for in vivo imaging of the gastrointestinal and pulmonary tracts [1]. OCT is based on low coherence interferometry. The excellent resolution in an OCT system results from the sharp peak power received at the zero optical path difference when all the wavelengths are in phase. To achieve this resolution, a broadband optical source is used. As the spectral width of the light source increases, the axial resolution of the systems is improved. However, like any interferometric system, the existence of dispersion in one of the interferometer arms results in a spread in the white light peak point and an associated phase error [2]. Dispersion manifests itself in two ways. The first is an error in the depth measurements, i.e., in the position determination of a certain point in the sample. Thus, if the sample is composed of a layered structure, there is an error in the thickness measurements of the layers, which may be of the order of a few percent. The second 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6643
is the flattening of the envelope of the interferogram, thus worsening the resolution. In [3], the authors discussed a beating effect that results from dispersion when trying to image two closely separated surfaces, causing multiple splitting of the peaks, which may be confused with a multilayer structure. One technique to overcome sample dispersion is to introduce physically a dispersive piece in the reference arm to balance dispersion in the sample arm [4]. This technique, usually used in spectroscopy, is very hard to implement in practical OCT measurement. Other hardware dispersion compensation techniques include using grating tilt in a rapid scanning frequency domain delay line, as in [5–7]. This technique is not capable of correcting multiple dispersion orders. To compensate for higher dispersion orders, such as second and third orders, the authors of [8] inserted a delay line and introduced more length of a single-mode fiber in the sample arm in the case that there is an acousto-optic frequency modulator or an electro-optic phase modulator in the reference arm. Numerical dispersion compensation techniques include the digital resampling method in the frequency domain [9], applying phase-sensitive digital filtering [10], applying numerical correlation to the depth scan by using a depth variant kernel [11], and fitting dispersion parameters to the spectral interferogram [12]. These numerical techniques are approved to compensate for the second and third dispersion orders, but not higher. Finally, due to the fact that water is the major constituent in nearly all biological tissues [13], most of the numerical techniques assume that the dispersion in the sample is known and can be modeled by the water dispersion model. In our work, we assume also that the sample dispersion is known and that the wavelength dependence of the sample absorption is neglected. This allows proposing a numerical dispersion compensation technique that is valid for arbitrary high dispersion orders. The proposed technique is an iterative one that corrects for errors in both resolution and absolute depth. Further, this technique is applicable to multi-layered samples. The paper is composed of 6 sections. Section 1 is an introduction. Section 2 presents the basic idea of OCT. Section 3 discusses dispersion due to the wavelength dependence of refractive index. Section 4 presents our proposed dispersion compensation technique. Section 5 presents the results obtained by applying our algorithm. Finally, Section 6 concludes on the work.
where I s and I r are the light intensities coming from the sample and the reference arms, x is the optical path difference between both arms, and λ is the wavelength. jγxj is the modulus of the complex degree of coherence function of the source. The third term in Eq. (1) corresponds to the interference between the backreflected signals from the sample and reference arms. Optical sectioning in the sample is performed by recording the modulation pattern as we introduce optical delay in the reference arm. The modulation pattern has a peak at the interface between two media. The width of this peak determines the resolution in measuring the interface position. This resolution is governed by the coherence length of the source, which is given by lc
2 ln 2 λ20 ; π Δλ
(2)
where λ0 is the center wavelength, Δλ is the full width at half-maximum (FWHM) of the source power
2. Basic Theory of OCT
A typical time-domain OCT (TD-OCT) schematic, based on the Michelson interferometer, is shown in Fig. 1(a). The reading of the detector is given by p 2π x ; I I s I r 2 I s I r jγxj cos λ 6644
(1)
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
Fig. 1. Schematics of (a) the TD-OCT system and (b) the Fourierdomain OCT system.
spectral density, and the source is assumed to have a Gaussian spectrum. The Fourier domain OCT system is shown in Fig. 1(b); it is almost the same Michelson interferometer, but there is no motion in the reference arm. Instead, there is a stationary mirror and the detector is replaced by an optical spectrum analyzer. The reading of the spectrometer is called a spectral interferogram. Applying the inverse Fourier transform (IFT) on the spectral interferogram allows obtaining the depth tracking signal similar to that obtained in the TD-OCT. 3. Dispersion in an OCT System
Dispersion in OCT systems may arise from unbalanced lengths of the fiber in the sample and reference arms, in the case we are using a fiber optics interferometer. Further, dispersion may result from variation of the refractive index in the sample with wavelength. This sample dispersion is more difficult to compensate, and it is the main focus of this work. It results in the broadening of the main lobe representing the interface, thus resulting in loss in resolution. Moreover, the position of the interface might also be shifted due to the inaccurate value of the refractive index used in the extraction of the thickness, which results in error in the thickness measurement. The dispersion in a medium (sample) is characterized by the variation of the propagation constant k into the medium. The propagation constant variation with frequency can be expanded by Taylor series expansion as kω kω0 ω − ω0
addition, a nonhomogeneous medium can also be treated by dividing it into smaller layers of homogeneous refractive index as an approximation. The distance traveled by the light beam to reach the first interface of the sample is dr d0, where dr is the distance traveled by the light beam to reach the stationary mirror in the reference arm. In our technique, we assume that we do not know the thickness of any layer. Our strategy goes as follows. A. Wavelength Dependence of Refractive Index
We assume that n1 λ, n2 λ are known from Cauchy or Sellmeier forms. Then we approximate it into a constant term and a wavelength varying part as n1 λ n01 Δn1 λ;
(4)
n2 λ n02 Δn2 λ;
(5)
where
dk ω − ω0 2 d2 k 2 dω dω2
ω − ω0 3 d3 k . 6 dω3
Fig. 2. Multilayered sample.
n01 n1 λ0 ; (3)
The first term in Eq. (3), represents a constant offset and the second term represents the group velocity (first-order dispersion), and neither is related to dispersive broadening in the OCT signal. The third term represents second-order dispersion or group velocity dispersions. The fourth term represents third-order dispersion, etc. Further, in the forthcoming analysis we will not use the propagation constant k; instead, we use the wavenumber σ, which is defined as the inverse of wavelength λ, i.e., σ 1∕λ. 4. Numerical Dispersion Compensation
To describe our dispersion compensation technique we assume a stratified sample, as shown in Fig. (2). The sample is surrounded by air and is composed of two layers followed by a semi-infinite medium. The two layers are dispersive with refractive indices n1 λ and n2 λ, and thicknesses d1 and d2 , respectively. The semi-infinite medium has refractive index n3 , which is independent of wavelength; this is for simplicity only and without loss of generality. Actually, the algorithm can handle any number of layers. In
n02 n2 λ0 ;
(6)
where λ0 is the central wavelength in the range used. Further, the remaining wavelength-dependent part is fitted to a high-order polynomial in wavelength. In our practical implementation, a fifth-order polynomial in wavelength is used as Δnu λ
5 X
Cm;u λ − λ0 m ;
(7)
m1
where Cm;u are constants and u 1; 2. B. Spectral Interferogram
The spectral interferogram is the spectral values measured by the spectrometer at the output of the system. These values can be numerically expressed by Jσ Sσ
N X i0
i X ρi cos 4πσ nt dt ;
(8)
t0
where nt and dt are the refractive index and thickness of the tth layer, respectively. Moreover, σ is 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6645
J 1 σ FTI 1 x:
the wavenumber, which is the inverse of wavelength λ, and Sσ is the spectrum and is given by Sσ
1 0
σs ≤ σ ≤ σe : Otherwise
(9)
Here, σ s and σ e are the start and end wavenumbers, respectively. We assume N as the number of layers. Further, ρ is the reduced amplitude reflectivity, which accounts for the effective amplitude attenuation of the beam as it propagates through different interfaces in the nonabsorbing sample. Thus we have ρi
i−1 Y 1 − r2t ri ;
Then, we decompose the resultant spectral interferogram into one for positive and the other for negative frequencies as follows: J 1 σ
σ ≥ 0:
(11)
For negative frequencies, we have −
J σ Jjσj;
σ < 0:
J 0 σ
J σ σ ≥ 0 ; J − σ σ < 0
2
i X ρi exp −4πjσ nt dt : 2 t0 i1
(17)
J 2 σ J 1 σ exp4πjσd1 Δn1
N X ρi i1
(12)
2
exp 4πjσd1 n01 2Δn1
2
× exp 44πjσ
i X
3 nt dt 5
t0 t≠1
2
(13)
× exp 4−4πjσ
N X ρi i1
i X
2
exp −4πjσd1 n01
3 nt dt 5:
(18)
t0 t≠1
Similarly, for the negative part we have
C. Iterative Bouncing between Space and Frequency Domains
At first we obtain the IFT of the spectral interferogram as (14)
which does not suffer any dispersion at d0 and the position of it can be easily recognized. Moreover, rough approximation of d1 can be obtained as the peak of the second main lobe. Further, the region at d0 10 μm and its mirror is zeroed in I 0 x to obtain I 1 x, and this region is also extracted for further use. This step is clarified in Fig. (3). In Fig. 3(a), we consider an example of the interferogram I 0 x, where the dotted vertical line and the dashed vertical line represent the true positions of the dispersive layers. In Fig. 3(b), we plot the interferogram I 1 x for this case; dashed–dotted vertical line represents the position of d0 at which the interferogram is zeroed. Then we apply Fourier transform to the partially zeroed interferogram to get a new spectral interferogram: 6646
t0
N X
i.e., we have even spectral interferogram whose IFT is real.
I 0 x IFTJ 0 σ;
(16)
Knowing a rough value of d1, we multiply the last spectral interferogram by an exponential factor as follows:
The total spectral interferogram is given as
σ≥0 : σ ε, where ε is a certain small threshold that is equal to 0.00001, update the current value of d1 to the new value, i.e., make d1Current d1New . Then go to step 10. 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6647
15. If the condition in step 14 is not satisfied, use the new value of d1 as the exact value, i.e., d1Exact d1New . 16. From the interferogram obtained in step 11 get the following: i. an interferogram that is zero everywhere except at d1Exact 10 μm and its image. Call it I d1 x; and ii. an interferogram that is zeroed at d1Exact 10 μm and its image and assign it to the current interferogram I Current x. 17. For the current interferogram obtained in step 16(ii), apply steps 7–16 for other layers using di and Δni , where i 2 to N. 18. Finally get the corrected interferogram ΣI di x. 5. Results and Discussion A.
Input Data
To verify the validity of the proposed dispersion compensation technique, we simulated the spectral interferogram by Eq. (8). First, the distances are as follows: d0 75 μm, d1 250 μm, and d2 700 μm. Second, the refractive index of each layer is as follows. Layer (1) is silicon with refractive index that depends on wavelength as follows: n2Si λ 11.6858
bλ2 a 2 1 2; 2 λ − λ1 λ
(22)
where wavelength is in micrometers and a 9.39816 × 10−1 , b 8.10461 × 10−3 , and λ1 1.1071. Layer (2) is fused silica with refractive index that depends on the wavelength as follows: n2FS λ 1
2 X A k λ2 ; λ2 − B2k k0
(23)
where 2
3 0.6961663 A 4 0.4079426 5; 0.8974794
2
3 0.0684043 and B 4 0.1162414 5; 9.896161
in actual experimental work the measured spectral interferogram is normalized with respect to the source spectrum. B. Results
1. General Aspects of the Proposed Technique In this subsection, we demonstrate the capability of the proposed technique to compensate dispersion by its application in different cases. First, we consider the results obtained for the unknown values of layer thickness. We applied our proposed algorithm for 10 iterations at each layer and we obtained the data in Table 1. Note that, for the thickness of layer (1) d1 , stabilizes after the fourth iteration and the thickness of layer (2) d2 , stabilizes after the second iteration. Second, we consider the results for correcting the positions of the peaks of the interferogram. Figure 3(a) presents the interferogram before correction and Fig. 3(c) presents the interferogram after dispersion compensation. Note that the dashed vertical line and the dotted vertical line represent the actual positions. Third, we consider the main lobe broadening due to dispersion and the effect of our proposed algorithm on the FWHM of these lobes. First, Figs. 4(a) and 4(b), present the interferogram and its envelope at the first interface. This interface does not suffer any dispersive effects. Fourth, Figs. 5(a) and 5(b) present the interferogram and its envelope at the second interface before dispersion compensation. Figures 5(c) and 5(d) present the interferogram and its envelope at the same interface, but after correction. Figure 5(e) gives the envelope of the interferogram at the same interface before and after dispersion compensation for comparison. To quantify the effect of dispersion before and after compensation, we measured the FWHM of the main lobes at the three interfaces and present them in Table 2. Finally, we consider the time of calculations for one thickness measurement. We consider a case where d1 d2 200 μm, λ0 1.25 μm, and bandwidth 100 nm. In this case we applied FT and IFT with resolution in the wavenumber domain 2.127 cm−1, resolution in the optical path difference 0.07 μm, and number of points 67; 195. The time required to calculate the thickness d1 is equal to 2.16 s, which corresponds to only four iterations.
(24) and wavelength is in micrometers. finally, for the semi-infinite medium (3), we set the refractive index equal to 2 and assume that it is wavelength independent for simplicity. Further, we assumed that the source has a uniform spectrum centered at λ0 1.250 μm and a spectral width of 100 nm. Thus, the edges of the spectrum in Eq. (9) are given by σs
1 ; λ0 Δλ
σe
1 ; λ0 − Δλ
Table 1.
(25)
with Δλ 50 nm. Note that assuming a uniform spectrum does not result in loss of generality, since 6648
2. Artificially Increasing Dispersion Further, to ensure the power of our proposed algorithm to cancel dispersion effects, we virtually increase the dispersion effects by two different methods and show that our proposed algorithm is capable
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
d1 d2
Thickness Values Obtained by Our Proposed Algorithm.
Actual Value (μm)
Values Obtained by Our Proposed Algorithm (μm)
Number of Iterations Required
250 700
249.83 699.58
4 2
Fig. 4. (a) Interferogram at the first interface and (b) its envelope.
of compensating these effects. First, to increase dispersion effects, we modified the dispersion relation of silicon, i.e., refractive index of the first layer, in Eq. (22), by multiplying the wavelength-dependent part by a factor G, so the equation becomes n2Si λ 11.6858 G
bλ21 a : λ2 λ2 − λ21
(26)
Then, we assumed constant values of d1 d2 200 μm, while changing G and applying our proposed algorithm. To quantify the dispersion effects we measure the resolution, i.e., the FWHM of the interferogram envelope at the second interface, which is between silicon and fused silica. Moreover, we calculate the percentage error in measuring the first layer thickness d1 and define this error as Percentage Error
d1 Measured − d1 Actual d1 Actual × 100:
(27)
Figure 6(a) shows the effect of varying the constant G, in which the percentage error (without dispersion compensation) increases as expected, and it is much smaller after dispersion compensation, i.e., less than 0.12%. Also, the resolution at the second interface is much more degraded as G increases without compensation, and it is nearly constant and even smaller when the dispersion compensation is used. This is presented in Fig. 6(b). 3. Increasing Thickness Another way to demonstrate the dispersion effect is to increase the thickness of the dispersive layer (the first layer in our analysis), i.e., increasing d1. Next, we simulate the spectral interferogram and apply IFT to get the spatial interferogram. Then we calculate the percentage error in measuring the first layer and the resolution at the second interface. In Fig. 7(a), we present the percentage error as a function of the thickness of the first layer d1, with and without dispersion
compensation. Again, we also measured the FWHM at the second interface and plot it versus the layer thickness for the structure with and without dispersion compensation in Fig. 7(b). Note that, in this last analysis, we kept d2 constant at the value of 200 μm. From Fig. 7(a), we find that the percentage error in measuring d1 without dispersion compensation is nearly constant at a value of about 6%. On the other hand, with dispersion compensation, the percentage error is less than 0.09%. From Fig. 7(b), we find that the value of resolution, i.e., the FWHM of the interferogram at the second interface, increases as we increase d1 without dispersion compensation, and it is constant and smaller with dispersion compensation. A great difference is clear at large values of d1. 4. Varying the Central Wavelength Finally, we discuss the effects of varying the central wavelength λ0 on the obtained results. This is because dispersion for silicon is very high as we approach 1.1 μm, which is clear from Eq. (22). Therefore, we apply the same analysis and vary the center wavelength λ0 and get the percentage error in the first layer thickness (silicon) and the percentage error in the second layer thickness (fused silica). Moreover, we measure the resolution (FWHM) of the interferogram envelope at the second interface, between silicon and fused silica, and the third interface, between fused silica and semi-infinite layer (3). We applied this analysis for d1 1000 μm and d2 350 μm. In Fig. 8(a) we present the results obtained for the percentage error in the first layer and the second layer as a function of the center wavelength with and without dispersion compensation. Again, the error is very small after using the compensation compared to that without compensation. Further, the error in the first layer, silicon, is greater than that of the second layer, fused silica, because we are very near to 1.1 μm for silicon. In Fig. 8(b), we present the results obtained for the system resolution (FWHM) of the envelope of the interferogram at the second interface and third interface before and after dispersion compensation. From this figure, we can see that at λ0 < 1.5 μm the FWHM is much greater before compensation than that after 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6649
Fig. 5. Interferogram and its envelope at the second interface before dispersion compensation are shown in (a) and (b), respectively. The interferogram and its envelope at the same interface after dispersion compensation are shown in (c) and (d), respectively. Last, (e) shows the envelope of the interferogram at the second interface before and after dispersion compensation for comparison.
compensation for both interfaces. This is very clear since we approach the wavelength 1.1 μm for silicon. After 1.5 μm, the curves before and after compensation appear to nearly coincide since dispersion effects are much less. Finally, FWHM increases as center wavelength increases for all cases, due to the finite constant bandwidth. 6650
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
5. Error Analysis In this section we discuss errors and uncertainty in wavelength dependence of the sample refractive index and its influence on the results of dispersion compensation. This is done to check the robustness of the proposed technique. To do so, we introduce errors in
Table 2.
FWHM of Main Lobes
FWHM (μm)
First interface Second interface Third interface
Before Dispersion Compensation
After Dispersion Compensation
9.26 22.05 21.8
9.26 9.62 9.55
the parameters of Sellmeier’s equation used in the dispersion compensation of the sample under test. Thus, the simulation of the sample response is carried out using the correct parameters of Sellmeier’s equation, while the extraction of the sample response in our technique is carried out with parameters having a slight error of the order of 10% of the parameter. This error could reflect the uncertainty in the exact composition of the sample material, such as how much water it may contain or what is the exact chemical composition of a certain layer in the sample. In our current test, these errors will be introduced only in the silicon layer of the assumed sample. First we simulate the spectral interferogram for depths d1 1000 μm and d2 200 μm. The source bandwidth is equal to 100 nm, with a center
wavelength at λ0 1.35 μm. In the simulated spectral interferogram, we used Eq. (22) for silicon without introducing any errors in the parameters. Second, in the dispersion compensation, we assume that we do not know exactly the refractive index of the sample and we estimate it with an error of the order of 10% in the parameters used for approximating the refractive index of silicon by Eqs. (4), (6), and (7). The introduced errors are in the parameters of Sellmeier’s Eq. (22), specifically a, b, and λ1 . The results obtained are given in Table 3. Table 3 shows that our proposed technique is powerful in compensating for dispersion effects even if we have uncertainties in the parameters of Sellmeier’s equation of the refractive index of the sample. The error in the obtained results (thickness) is still 1 order of magnitude lower than the uncompensated results. The resolution after compensation with error 5% in λ1 is 12.7 μm. This corresponds to the worst case, which does not differ much from the case without introducing errors. Further, the error in measuring thickness after compensation is less than 0.7%. This must be compared to those values before compensation. From Table 3 we find that errors in parameter b have the least effect. The results are more sensitive to errors in λ1 . This is because the
Fig. 6. Effect of varying the factor G in Eq. (26). (a) Shows the effect on percentage error in calculating first layer thickness d1. (b) Shows the effect on resolution, i.e., FWHM at the second interface. In both figures, solid curves show the case before dispersion compensation and the dotted curves show the case after dispersion compensation.
Fig. 7. (a) Percentage error in measuring the first layer thickness as we increase d1 , with the solid curve before dispersion compensation and the dotted curve after dispersion compensation. (b) Resolution at the second interface as a function of d1, with the solid curve before dispersion compensation and the dotted curve after dispersion compensation. 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6651
Fig. 8. (a) Percentage error in thickness measurements, with solid curve and crosses for d1 before and after dispersion compensation, respectively. The dotted and dashed–dotted curves are for d2 before and after dispersion compensation, respectively. (b) Resolution, with solid curve and crosses for the second interface before and after dispersion compensation, respectively. Dotted and dashed–dotted curves are for the third interface before and after dispersion compensation.
Table 3.
Calculated Parameters Resolution before compensation (μm) Resolution after compensation (μm) Thickness error before compensation (%) Thickness error after compensation (%)
Results Obtained with Errors in Parameters of Eq. (22)
Introducing Error in λ1
−10%
10%
−10%
10%
−5%
5%
43.7 11.15 4.643 0.637
43.7 11.17 4.203 0.558
43.7 11.05 4.430 0.065
43.7 11.10 4.415 0.017
43.7 11.1 4.440 0.172
43.7 12.7 4.394 0.067
Discussion of Results
Table 1 verifies that our proposed algorithm can be used to find the unknown thicknesses if we know the dependence of refractive index on wavelength. The thickness error is less than 0.075% in not more than four iterations. Further, we corrected the positions of the interfaces and corrected the shape of the interferogram and its envelope, which is clear from Figs. 5(a) and 5(b) compared to Figs. 5(c) and 5(d). Note that we did not add a curve at the third interface because it has the same shape of that at the second interface. This is because the dispersion of silicon is much more dominant than that of fused silica at the center wavelength, i.e., 1.25 μm. Moreover, the interferograms in Figs. 5(a) and 5(b) suffer from high asymmetry about the main peak. This is due to dispersion effects, a case which can be found 6652
Introducing Error in b
43.7 11.05 4.422 0.029
working wavelength is 1.35 μm, which is near λ1 at which the refractive index has a very large value. Thus, we limited errors in this parameter to 5% because of the very high dependence of refractive value on this parameter. Moreover, error of −5% in λ1 has less effect than 5%; this is due to the resonance nature of the refractive index dependence on the wavelength. Note in Table 3 that the value of thickness error before compensation fluctuates. This is because it depends on n01 , which depends on errors in Sellmeier’s equation. C.
Introducing Error in a
Without Errors in Parameters of Sellmeier’s
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
in an ordinary Michelson interferometer when inserting a piece of glass in one arm only [14]. Finally, we have corrected the FWHM at the main lobes. As can be seen from Table 2, the erroneous values are more than twice the corrected values. Moreover, as we increased dispersion of one layer either by increasing its thickness or by artificially increasing the wavelength-dependent part by a certain factor and quantifying the power of our proposed algorithm by percentage thickness error and resolution, we found that our proposed algorithm gives marvelous results when compared to that before dispersion compensation. It should be noted that in previous work [12], the researchers compensated dispersion up to only the third order. In contrast, our algorithm compensates for the fifth order and can be extended further. The maximum order of dispersion is based on the material of the sample and the rate of variation of its refractive index with wavelength. Further, higher orders of dispersion may appear due to the operating wavelength. 6. Conclusion
We have proposed an iterative numerical dispersion compensation technique for Fourier domain OCT. In this technique, it is assumed that we know the wavelength dependence of the refractive index of the sample. Our algorithm performs well for stratified homogeneous media and corrects for arbitrary orders
of dispersion. With our algorithm we can measure distances and thicknesses with much better accuracy with no more than four iterations. Numerical results showed that errors of the order of 5% in thickness measurement without dispersion compensation are corrected to have errors of less than 0.1% with the proposed compensation technique. Further, interferogram flattening due to dispersion effects resulted in resolution of 50 μm, which is corrected by our technique to less than 10 μm. In future work, we may correct dispersion effects even if we do not know wavelength dependence of the refractive index. References 1. W. Drexler and J. G. Fujimoto, eds., Optical Coherence Tomography: Technology and Applications (Springer-Verlag, 2008). 2. T. A. Al-Saeed and D. A. Khalil, “Dispersion compensation in moving-optical-wedge Fourier transform spectrometer,” Appl. Opt. 48, 3979–3987 (2009). 3. C. K. Hitzenberger, A. Baumgartner, and A. F. Fercher, “Dispersion induced multiple signal peak splitting in partial coherence interferometry,” Opt. Commun. 154, 179–185 (1998). 4. W. Drexler, U. Morgner, F. X. Kärtner, C. Pitris, S. A. Boppart, X. D. Li, E. P. Ippen, and J. G. Fujimoto, “In vivo ultrahigh resolution optical coherence tomography,” Opt. Lett. 24, 1221–1223 (1999). 5. E. D. J. Smith, A. V. Zvyagin, and D. D. Sampson, “Real-time dispersion compensation in scanning interferometry,” Opt. Lett. 27, 1998–2000 (2002).
6. A. V. Zvyagin, E. D. J. Smith, and D. D. Sampson, “Delay and dispersion characteristics of a frequency-domain optical delay line for scanning interferometry,” J. Opt. Soc. Am. A 20, 333–341 (2003). 7. W. K. Niblack, J. O. Schenk, B. Liu, and M. E. Brezinski, “Dispersion in a grating-based optical delay line for optical coherence tomography,” Appl. Opt. 42, 4115–4118 (2003). 8. Y. Chen and X. Li, “Dispersion management up to the third order for real-time optical coherence tomography involving a phase or frequency modulator,” Opt. Express 12, 5968– 5978 (2004). 9. D. L. Marks, A. L. Oldenburg, J. J. Reynolds, and S. A. Boppart, “Digital algorithm for dispersion correction in optical coherence tomography for homogeneous and stratified media,” Appl. Opt. 42, 204–216 (2003). 10. J. F. de Boer, C. E. Saxer, and J. S. Nelson, “Stable carrier generation and phase-resolved digital data processing in optical coherence tomography,” Appl. Opt. 40, 5787–5790 (2001). 11. A. F. Fercher, C. K. Hitzenberger, M. Sticker, R. Zawadzki, B. Karamata, and T. Lasser, “Numerical dispersion compensation for partial coherence interferometry and optical coherence tomography,” Opt. Express 9, 610–615 (2001). 12. M. Wojtkowski, V. J. Srinivasan, T. H. Ko, J. G. Fujimoto, A. Kowalczyk, and J. S. Duker, “Ultrahigh-resolution, highspeed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12, 2404–2422 (2004). 13. T. R. Hillman and D. D. Sampson, “The effect of water dispersion and absorption on axial resolution in ultrahighresolution optical coherence tomography,” Opt. Express 13, 1860–1874 (2005). 14. S. Diddams and J.-C. Diels, “Dispersion measurements with white-light interferometry,” J. Opt. Soc. Am. B 13, 1120–1129 (1996).
10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6653
Biomedical Engineering Department, Faculty of Engineering, Helwan University, 1 Sherif Street, Helwan, Cairo, Egypt 2
Electronics and Communications Engineering Department, Faculty of Engineering, Ain-Shams University, 11 Elsarayat Street, Abbassia, Cairo 11517, Egypt 3
Si-Ware Systems Company, 3 Khaled Ibn Al-Waleed Street, Heliopolis, Cairo, Egypt *Corresponding author: [email protected] Received 28 April 2014; revised 25 July 2014; accepted 19 August 2014; posted 3 September 2014 (Doc. ID 210914); published 2 October 2014
In this work, we propose a numerical technique to compensate for errors due to dispersion effects in Fourier domain optical coherence tomography. The proposed technique corrects for errors in depth measurements and resolution loss due to dispersion. The results show that, by using this technique, errors in thickness measurement are reduced from about 5% to less than 0.1% depending on the sample length and the amount of dispersion. Also, an improvement in the resolution from about 50 μm to less than 10 μm is demonstrated. © 2014 Optical Society of America OCIS codes: (110.4500) Optical coherence tomography; (030.1670) Coherent optical effects; (170.1650) Coherence imaging. http://dx.doi.org/10.1364/AO.53.006643
1. Introduction
Optical coherence tomography (OCT) is a competitive in vivo noninvasive imaging technique [1]. Similar to the ultrasound technique, OCT provides A-scan and B-scan images. These images are obtained using the low coherence interferometry principle. Thus, depending on the optical system, the axial and lateral resolution can be of the order of 3–10 μm, which is much better than the resolution of the ultrasound technique. On the other hand, it can provide an in-depth image to a penetration depth up to a few millimeters. In this sense, it fills the gap between ultrasound imaging and confocal microscopy. There are many applications of OCT in medicine. OCT can be used in imaging cellular and tissue microstructure, in vivo retinal imaging, imaging coronary atherosclerosis, imaging laryngeal mucosa, skin imaging, dental tissue imaging, and corneal 1559-128X/14/296643-11$15.00/0 © 2014 Optical Society of America
thickness measurement. It can also be used in monitoring the growth and development of biological tissues over time. Finally, endoscopic OCT can be used for in vivo imaging of the gastrointestinal and pulmonary tracts [1]. OCT is based on low coherence interferometry. The excellent resolution in an OCT system results from the sharp peak power received at the zero optical path difference when all the wavelengths are in phase. To achieve this resolution, a broadband optical source is used. As the spectral width of the light source increases, the axial resolution of the systems is improved. However, like any interferometric system, the existence of dispersion in one of the interferometer arms results in a spread in the white light peak point and an associated phase error [2]. Dispersion manifests itself in two ways. The first is an error in the depth measurements, i.e., in the position determination of a certain point in the sample. Thus, if the sample is composed of a layered structure, there is an error in the thickness measurements of the layers, which may be of the order of a few percent. The second 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6643
is the flattening of the envelope of the interferogram, thus worsening the resolution. In [3], the authors discussed a beating effect that results from dispersion when trying to image two closely separated surfaces, causing multiple splitting of the peaks, which may be confused with a multilayer structure. One technique to overcome sample dispersion is to introduce physically a dispersive piece in the reference arm to balance dispersion in the sample arm [4]. This technique, usually used in spectroscopy, is very hard to implement in practical OCT measurement. Other hardware dispersion compensation techniques include using grating tilt in a rapid scanning frequency domain delay line, as in [5–7]. This technique is not capable of correcting multiple dispersion orders. To compensate for higher dispersion orders, such as second and third orders, the authors of [8] inserted a delay line and introduced more length of a single-mode fiber in the sample arm in the case that there is an acousto-optic frequency modulator or an electro-optic phase modulator in the reference arm. Numerical dispersion compensation techniques include the digital resampling method in the frequency domain [9], applying phase-sensitive digital filtering [10], applying numerical correlation to the depth scan by using a depth variant kernel [11], and fitting dispersion parameters to the spectral interferogram [12]. These numerical techniques are approved to compensate for the second and third dispersion orders, but not higher. Finally, due to the fact that water is the major constituent in nearly all biological tissues [13], most of the numerical techniques assume that the dispersion in the sample is known and can be modeled by the water dispersion model. In our work, we assume also that the sample dispersion is known and that the wavelength dependence of the sample absorption is neglected. This allows proposing a numerical dispersion compensation technique that is valid for arbitrary high dispersion orders. The proposed technique is an iterative one that corrects for errors in both resolution and absolute depth. Further, this technique is applicable to multi-layered samples. The paper is composed of 6 sections. Section 1 is an introduction. Section 2 presents the basic idea of OCT. Section 3 discusses dispersion due to the wavelength dependence of refractive index. Section 4 presents our proposed dispersion compensation technique. Section 5 presents the results obtained by applying our algorithm. Finally, Section 6 concludes on the work.
where I s and I r are the light intensities coming from the sample and the reference arms, x is the optical path difference between both arms, and λ is the wavelength. jγxj is the modulus of the complex degree of coherence function of the source. The third term in Eq. (1) corresponds to the interference between the backreflected signals from the sample and reference arms. Optical sectioning in the sample is performed by recording the modulation pattern as we introduce optical delay in the reference arm. The modulation pattern has a peak at the interface between two media. The width of this peak determines the resolution in measuring the interface position. This resolution is governed by the coherence length of the source, which is given by lc
2 ln 2 λ20 ; π Δλ
(2)
where λ0 is the center wavelength, Δλ is the full width at half-maximum (FWHM) of the source power
2. Basic Theory of OCT
A typical time-domain OCT (TD-OCT) schematic, based on the Michelson interferometer, is shown in Fig. 1(a). The reading of the detector is given by p 2π x ; I I s I r 2 I s I r jγxj cos λ 6644
(1)
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
Fig. 1. Schematics of (a) the TD-OCT system and (b) the Fourierdomain OCT system.
spectral density, and the source is assumed to have a Gaussian spectrum. The Fourier domain OCT system is shown in Fig. 1(b); it is almost the same Michelson interferometer, but there is no motion in the reference arm. Instead, there is a stationary mirror and the detector is replaced by an optical spectrum analyzer. The reading of the spectrometer is called a spectral interferogram. Applying the inverse Fourier transform (IFT) on the spectral interferogram allows obtaining the depth tracking signal similar to that obtained in the TD-OCT. 3. Dispersion in an OCT System
Dispersion in OCT systems may arise from unbalanced lengths of the fiber in the sample and reference arms, in the case we are using a fiber optics interferometer. Further, dispersion may result from variation of the refractive index in the sample with wavelength. This sample dispersion is more difficult to compensate, and it is the main focus of this work. It results in the broadening of the main lobe representing the interface, thus resulting in loss in resolution. Moreover, the position of the interface might also be shifted due to the inaccurate value of the refractive index used in the extraction of the thickness, which results in error in the thickness measurement. The dispersion in a medium (sample) is characterized by the variation of the propagation constant k into the medium. The propagation constant variation with frequency can be expanded by Taylor series expansion as kω kω0 ω − ω0
addition, a nonhomogeneous medium can also be treated by dividing it into smaller layers of homogeneous refractive index as an approximation. The distance traveled by the light beam to reach the first interface of the sample is dr d0, where dr is the distance traveled by the light beam to reach the stationary mirror in the reference arm. In our technique, we assume that we do not know the thickness of any layer. Our strategy goes as follows. A. Wavelength Dependence of Refractive Index
We assume that n1 λ, n2 λ are known from Cauchy or Sellmeier forms. Then we approximate it into a constant term and a wavelength varying part as n1 λ n01 Δn1 λ;
(4)
n2 λ n02 Δn2 λ;
(5)
where
dk ω − ω0 2 d2 k 2 dω dω2
ω − ω0 3 d3 k . 6 dω3
Fig. 2. Multilayered sample.
n01 n1 λ0 ; (3)
The first term in Eq. (3), represents a constant offset and the second term represents the group velocity (first-order dispersion), and neither is related to dispersive broadening in the OCT signal. The third term represents second-order dispersion or group velocity dispersions. The fourth term represents third-order dispersion, etc. Further, in the forthcoming analysis we will not use the propagation constant k; instead, we use the wavenumber σ, which is defined as the inverse of wavelength λ, i.e., σ 1∕λ. 4. Numerical Dispersion Compensation
To describe our dispersion compensation technique we assume a stratified sample, as shown in Fig. (2). The sample is surrounded by air and is composed of two layers followed by a semi-infinite medium. The two layers are dispersive with refractive indices n1 λ and n2 λ, and thicknesses d1 and d2 , respectively. The semi-infinite medium has refractive index n3 , which is independent of wavelength; this is for simplicity only and without loss of generality. Actually, the algorithm can handle any number of layers. In
n02 n2 λ0 ;
(6)
where λ0 is the central wavelength in the range used. Further, the remaining wavelength-dependent part is fitted to a high-order polynomial in wavelength. In our practical implementation, a fifth-order polynomial in wavelength is used as Δnu λ
5 X
Cm;u λ − λ0 m ;
(7)
m1
where Cm;u are constants and u 1; 2. B. Spectral Interferogram
The spectral interferogram is the spectral values measured by the spectrometer at the output of the system. These values can be numerically expressed by Jσ Sσ
N X i0
i X ρi cos 4πσ nt dt ;
(8)
t0
where nt and dt are the refractive index and thickness of the tth layer, respectively. Moreover, σ is 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6645
J 1 σ FTI 1 x:
the wavenumber, which is the inverse of wavelength λ, and Sσ is the spectrum and is given by Sσ
1 0
σs ≤ σ ≤ σe : Otherwise
(9)
Here, σ s and σ e are the start and end wavenumbers, respectively. We assume N as the number of layers. Further, ρ is the reduced amplitude reflectivity, which accounts for the effective amplitude attenuation of the beam as it propagates through different interfaces in the nonabsorbing sample. Thus we have ρi
i−1 Y 1 − r2t ri ;
Then, we decompose the resultant spectral interferogram into one for positive and the other for negative frequencies as follows: J 1 σ
σ ≥ 0:
(11)
For negative frequencies, we have −
J σ Jjσj;
σ < 0:
J 0 σ
J σ σ ≥ 0 ; J − σ σ < 0
2
i X ρi exp −4πjσ nt dt : 2 t0 i1
(17)
J 2 σ J 1 σ exp4πjσd1 Δn1
N X ρi i1
(12)
2
exp 4πjσd1 n01 2Δn1
2
× exp 44πjσ
i X
3 nt dt 5
t0 t≠1
2
(13)
× exp 4−4πjσ
N X ρi i1
i X
2
exp −4πjσd1 n01
3 nt dt 5:
(18)
t0 t≠1
Similarly, for the negative part we have
C. Iterative Bouncing between Space and Frequency Domains
At first we obtain the IFT of the spectral interferogram as (14)
which does not suffer any dispersion at d0 and the position of it can be easily recognized. Moreover, rough approximation of d1 can be obtained as the peak of the second main lobe. Further, the region at d0 10 μm and its mirror is zeroed in I 0 x to obtain I 1 x, and this region is also extracted for further use. This step is clarified in Fig. (3). In Fig. 3(a), we consider an example of the interferogram I 0 x, where the dotted vertical line and the dashed vertical line represent the true positions of the dispersive layers. In Fig. 3(b), we plot the interferogram I 1 x for this case; dashed–dotted vertical line represents the position of d0 at which the interferogram is zeroed. Then we apply Fourier transform to the partially zeroed interferogram to get a new spectral interferogram: 6646
t0
N X
i.e., we have even spectral interferogram whose IFT is real.
I 0 x IFTJ 0 σ;
(16)
Knowing a rough value of d1, we multiply the last spectral interferogram by an exponential factor as follows:
The total spectral interferogram is given as
σ≥0 : σ ε, where ε is a certain small threshold that is equal to 0.00001, update the current value of d1 to the new value, i.e., make d1Current d1New . Then go to step 10. 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6647
15. If the condition in step 14 is not satisfied, use the new value of d1 as the exact value, i.e., d1Exact d1New . 16. From the interferogram obtained in step 11 get the following: i. an interferogram that is zero everywhere except at d1Exact 10 μm and its image. Call it I d1 x; and ii. an interferogram that is zeroed at d1Exact 10 μm and its image and assign it to the current interferogram I Current x. 17. For the current interferogram obtained in step 16(ii), apply steps 7–16 for other layers using di and Δni , where i 2 to N. 18. Finally get the corrected interferogram ΣI di x. 5. Results and Discussion A.
Input Data
To verify the validity of the proposed dispersion compensation technique, we simulated the spectral interferogram by Eq. (8). First, the distances are as follows: d0 75 μm, d1 250 μm, and d2 700 μm. Second, the refractive index of each layer is as follows. Layer (1) is silicon with refractive index that depends on wavelength as follows: n2Si λ 11.6858
bλ2 a 2 1 2; 2 λ − λ1 λ
(22)
where wavelength is in micrometers and a 9.39816 × 10−1 , b 8.10461 × 10−3 , and λ1 1.1071. Layer (2) is fused silica with refractive index that depends on the wavelength as follows: n2FS λ 1
2 X A k λ2 ; λ2 − B2k k0
(23)
where 2
3 0.6961663 A 4 0.4079426 5; 0.8974794
2
3 0.0684043 and B 4 0.1162414 5; 9.896161
in actual experimental work the measured spectral interferogram is normalized with respect to the source spectrum. B. Results
1. General Aspects of the Proposed Technique In this subsection, we demonstrate the capability of the proposed technique to compensate dispersion by its application in different cases. First, we consider the results obtained for the unknown values of layer thickness. We applied our proposed algorithm for 10 iterations at each layer and we obtained the data in Table 1. Note that, for the thickness of layer (1) d1 , stabilizes after the fourth iteration and the thickness of layer (2) d2 , stabilizes after the second iteration. Second, we consider the results for correcting the positions of the peaks of the interferogram. Figure 3(a) presents the interferogram before correction and Fig. 3(c) presents the interferogram after dispersion compensation. Note that the dashed vertical line and the dotted vertical line represent the actual positions. Third, we consider the main lobe broadening due to dispersion and the effect of our proposed algorithm on the FWHM of these lobes. First, Figs. 4(a) and 4(b), present the interferogram and its envelope at the first interface. This interface does not suffer any dispersive effects. Fourth, Figs. 5(a) and 5(b) present the interferogram and its envelope at the second interface before dispersion compensation. Figures 5(c) and 5(d) present the interferogram and its envelope at the same interface, but after correction. Figure 5(e) gives the envelope of the interferogram at the same interface before and after dispersion compensation for comparison. To quantify the effect of dispersion before and after compensation, we measured the FWHM of the main lobes at the three interfaces and present them in Table 2. Finally, we consider the time of calculations for one thickness measurement. We consider a case where d1 d2 200 μm, λ0 1.25 μm, and bandwidth 100 nm. In this case we applied FT and IFT with resolution in the wavenumber domain 2.127 cm−1, resolution in the optical path difference 0.07 μm, and number of points 67; 195. The time required to calculate the thickness d1 is equal to 2.16 s, which corresponds to only four iterations.
(24) and wavelength is in micrometers. finally, for the semi-infinite medium (3), we set the refractive index equal to 2 and assume that it is wavelength independent for simplicity. Further, we assumed that the source has a uniform spectrum centered at λ0 1.250 μm and a spectral width of 100 nm. Thus, the edges of the spectrum in Eq. (9) are given by σs
1 ; λ0 Δλ
σe
1 ; λ0 − Δλ
Table 1.
(25)
with Δλ 50 nm. Note that assuming a uniform spectrum does not result in loss of generality, since 6648
2. Artificially Increasing Dispersion Further, to ensure the power of our proposed algorithm to cancel dispersion effects, we virtually increase the dispersion effects by two different methods and show that our proposed algorithm is capable
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
d1 d2
Thickness Values Obtained by Our Proposed Algorithm.
Actual Value (μm)
Values Obtained by Our Proposed Algorithm (μm)
Number of Iterations Required
250 700
249.83 699.58
4 2
Fig. 4. (a) Interferogram at the first interface and (b) its envelope.
of compensating these effects. First, to increase dispersion effects, we modified the dispersion relation of silicon, i.e., refractive index of the first layer, in Eq. (22), by multiplying the wavelength-dependent part by a factor G, so the equation becomes n2Si λ 11.6858 G
bλ21 a : λ2 λ2 − λ21
(26)
Then, we assumed constant values of d1 d2 200 μm, while changing G and applying our proposed algorithm. To quantify the dispersion effects we measure the resolution, i.e., the FWHM of the interferogram envelope at the second interface, which is between silicon and fused silica. Moreover, we calculate the percentage error in measuring the first layer thickness d1 and define this error as Percentage Error
d1 Measured − d1 Actual d1 Actual × 100:
(27)
Figure 6(a) shows the effect of varying the constant G, in which the percentage error (without dispersion compensation) increases as expected, and it is much smaller after dispersion compensation, i.e., less than 0.12%. Also, the resolution at the second interface is much more degraded as G increases without compensation, and it is nearly constant and even smaller when the dispersion compensation is used. This is presented in Fig. 6(b). 3. Increasing Thickness Another way to demonstrate the dispersion effect is to increase the thickness of the dispersive layer (the first layer in our analysis), i.e., increasing d1. Next, we simulate the spectral interferogram and apply IFT to get the spatial interferogram. Then we calculate the percentage error in measuring the first layer and the resolution at the second interface. In Fig. 7(a), we present the percentage error as a function of the thickness of the first layer d1, with and without dispersion
compensation. Again, we also measured the FWHM at the second interface and plot it versus the layer thickness for the structure with and without dispersion compensation in Fig. 7(b). Note that, in this last analysis, we kept d2 constant at the value of 200 μm. From Fig. 7(a), we find that the percentage error in measuring d1 without dispersion compensation is nearly constant at a value of about 6%. On the other hand, with dispersion compensation, the percentage error is less than 0.09%. From Fig. 7(b), we find that the value of resolution, i.e., the FWHM of the interferogram at the second interface, increases as we increase d1 without dispersion compensation, and it is constant and smaller with dispersion compensation. A great difference is clear at large values of d1. 4. Varying the Central Wavelength Finally, we discuss the effects of varying the central wavelength λ0 on the obtained results. This is because dispersion for silicon is very high as we approach 1.1 μm, which is clear from Eq. (22). Therefore, we apply the same analysis and vary the center wavelength λ0 and get the percentage error in the first layer thickness (silicon) and the percentage error in the second layer thickness (fused silica). Moreover, we measure the resolution (FWHM) of the interferogram envelope at the second interface, between silicon and fused silica, and the third interface, between fused silica and semi-infinite layer (3). We applied this analysis for d1 1000 μm and d2 350 μm. In Fig. 8(a) we present the results obtained for the percentage error in the first layer and the second layer as a function of the center wavelength with and without dispersion compensation. Again, the error is very small after using the compensation compared to that without compensation. Further, the error in the first layer, silicon, is greater than that of the second layer, fused silica, because we are very near to 1.1 μm for silicon. In Fig. 8(b), we present the results obtained for the system resolution (FWHM) of the envelope of the interferogram at the second interface and third interface before and after dispersion compensation. From this figure, we can see that at λ0 < 1.5 μm the FWHM is much greater before compensation than that after 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6649
Fig. 5. Interferogram and its envelope at the second interface before dispersion compensation are shown in (a) and (b), respectively. The interferogram and its envelope at the same interface after dispersion compensation are shown in (c) and (d), respectively. Last, (e) shows the envelope of the interferogram at the second interface before and after dispersion compensation for comparison.
compensation for both interfaces. This is very clear since we approach the wavelength 1.1 μm for silicon. After 1.5 μm, the curves before and after compensation appear to nearly coincide since dispersion effects are much less. Finally, FWHM increases as center wavelength increases for all cases, due to the finite constant bandwidth. 6650
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
5. Error Analysis In this section we discuss errors and uncertainty in wavelength dependence of the sample refractive index and its influence on the results of dispersion compensation. This is done to check the robustness of the proposed technique. To do so, we introduce errors in
Table 2.
FWHM of Main Lobes
FWHM (μm)
First interface Second interface Third interface
Before Dispersion Compensation
After Dispersion Compensation
9.26 22.05 21.8
9.26 9.62 9.55
the parameters of Sellmeier’s equation used in the dispersion compensation of the sample under test. Thus, the simulation of the sample response is carried out using the correct parameters of Sellmeier’s equation, while the extraction of the sample response in our technique is carried out with parameters having a slight error of the order of 10% of the parameter. This error could reflect the uncertainty in the exact composition of the sample material, such as how much water it may contain or what is the exact chemical composition of a certain layer in the sample. In our current test, these errors will be introduced only in the silicon layer of the assumed sample. First we simulate the spectral interferogram for depths d1 1000 μm and d2 200 μm. The source bandwidth is equal to 100 nm, with a center
wavelength at λ0 1.35 μm. In the simulated spectral interferogram, we used Eq. (22) for silicon without introducing any errors in the parameters. Second, in the dispersion compensation, we assume that we do not know exactly the refractive index of the sample and we estimate it with an error of the order of 10% in the parameters used for approximating the refractive index of silicon by Eqs. (4), (6), and (7). The introduced errors are in the parameters of Sellmeier’s Eq. (22), specifically a, b, and λ1 . The results obtained are given in Table 3. Table 3 shows that our proposed technique is powerful in compensating for dispersion effects even if we have uncertainties in the parameters of Sellmeier’s equation of the refractive index of the sample. The error in the obtained results (thickness) is still 1 order of magnitude lower than the uncompensated results. The resolution after compensation with error 5% in λ1 is 12.7 μm. This corresponds to the worst case, which does not differ much from the case without introducing errors. Further, the error in measuring thickness after compensation is less than 0.7%. This must be compared to those values before compensation. From Table 3 we find that errors in parameter b have the least effect. The results are more sensitive to errors in λ1 . This is because the
Fig. 6. Effect of varying the factor G in Eq. (26). (a) Shows the effect on percentage error in calculating first layer thickness d1. (b) Shows the effect on resolution, i.e., FWHM at the second interface. In both figures, solid curves show the case before dispersion compensation and the dotted curves show the case after dispersion compensation.
Fig. 7. (a) Percentage error in measuring the first layer thickness as we increase d1 , with the solid curve before dispersion compensation and the dotted curve after dispersion compensation. (b) Resolution at the second interface as a function of d1, with the solid curve before dispersion compensation and the dotted curve after dispersion compensation. 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6651
Fig. 8. (a) Percentage error in thickness measurements, with solid curve and crosses for d1 before and after dispersion compensation, respectively. The dotted and dashed–dotted curves are for d2 before and after dispersion compensation, respectively. (b) Resolution, with solid curve and crosses for the second interface before and after dispersion compensation, respectively. Dotted and dashed–dotted curves are for the third interface before and after dispersion compensation.
Table 3.
Calculated Parameters Resolution before compensation (μm) Resolution after compensation (μm) Thickness error before compensation (%) Thickness error after compensation (%)
Results Obtained with Errors in Parameters of Eq. (22)
Introducing Error in λ1
−10%
10%
−10%
10%
−5%
5%
43.7 11.15 4.643 0.637
43.7 11.17 4.203 0.558
43.7 11.05 4.430 0.065
43.7 11.10 4.415 0.017
43.7 11.1 4.440 0.172
43.7 12.7 4.394 0.067
Discussion of Results
Table 1 verifies that our proposed algorithm can be used to find the unknown thicknesses if we know the dependence of refractive index on wavelength. The thickness error is less than 0.075% in not more than four iterations. Further, we corrected the positions of the interfaces and corrected the shape of the interferogram and its envelope, which is clear from Figs. 5(a) and 5(b) compared to Figs. 5(c) and 5(d). Note that we did not add a curve at the third interface because it has the same shape of that at the second interface. This is because the dispersion of silicon is much more dominant than that of fused silica at the center wavelength, i.e., 1.25 μm. Moreover, the interferograms in Figs. 5(a) and 5(b) suffer from high asymmetry about the main peak. This is due to dispersion effects, a case which can be found 6652
Introducing Error in b
43.7 11.05 4.422 0.029
working wavelength is 1.35 μm, which is near λ1 at which the refractive index has a very large value. Thus, we limited errors in this parameter to 5% because of the very high dependence of refractive value on this parameter. Moreover, error of −5% in λ1 has less effect than 5%; this is due to the resonance nature of the refractive index dependence on the wavelength. Note in Table 3 that the value of thickness error before compensation fluctuates. This is because it depends on n01 , which depends on errors in Sellmeier’s equation. C.
Introducing Error in a
Without Errors in Parameters of Sellmeier’s
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
in an ordinary Michelson interferometer when inserting a piece of glass in one arm only [14]. Finally, we have corrected the FWHM at the main lobes. As can be seen from Table 2, the erroneous values are more than twice the corrected values. Moreover, as we increased dispersion of one layer either by increasing its thickness or by artificially increasing the wavelength-dependent part by a certain factor and quantifying the power of our proposed algorithm by percentage thickness error and resolution, we found that our proposed algorithm gives marvelous results when compared to that before dispersion compensation. It should be noted that in previous work [12], the researchers compensated dispersion up to only the third order. In contrast, our algorithm compensates for the fifth order and can be extended further. The maximum order of dispersion is based on the material of the sample and the rate of variation of its refractive index with wavelength. Further, higher orders of dispersion may appear due to the operating wavelength. 6. Conclusion
We have proposed an iterative numerical dispersion compensation technique for Fourier domain OCT. In this technique, it is assumed that we know the wavelength dependence of the refractive index of the sample. Our algorithm performs well for stratified homogeneous media and corrects for arbitrary orders
of dispersion. With our algorithm we can measure distances and thicknesses with much better accuracy with no more than four iterations. Numerical results showed that errors of the order of 5% in thickness measurement without dispersion compensation are corrected to have errors of less than 0.1% with the proposed compensation technique. Further, interferogram flattening due to dispersion effects resulted in resolution of 50 μm, which is corrected by our technique to less than 10 μm. In future work, we may correct dispersion effects even if we do not know wavelength dependence of the refractive index. References 1. W. Drexler and J. G. Fujimoto, eds., Optical Coherence Tomography: Technology and Applications (Springer-Verlag, 2008). 2. T. A. Al-Saeed and D. A. Khalil, “Dispersion compensation in moving-optical-wedge Fourier transform spectrometer,” Appl. Opt. 48, 3979–3987 (2009). 3. C. K. Hitzenberger, A. Baumgartner, and A. F. Fercher, “Dispersion induced multiple signal peak splitting in partial coherence interferometry,” Opt. Commun. 154, 179–185 (1998). 4. W. Drexler, U. Morgner, F. X. Kärtner, C. Pitris, S. A. Boppart, X. D. Li, E. P. Ippen, and J. G. Fujimoto, “In vivo ultrahigh resolution optical coherence tomography,” Opt. Lett. 24, 1221–1223 (1999). 5. E. D. J. Smith, A. V. Zvyagin, and D. D. Sampson, “Real-time dispersion compensation in scanning interferometry,” Opt. Lett. 27, 1998–2000 (2002).
6. A. V. Zvyagin, E. D. J. Smith, and D. D. Sampson, “Delay and dispersion characteristics of a frequency-domain optical delay line for scanning interferometry,” J. Opt. Soc. Am. A 20, 333–341 (2003). 7. W. K. Niblack, J. O. Schenk, B. Liu, and M. E. Brezinski, “Dispersion in a grating-based optical delay line for optical coherence tomography,” Appl. Opt. 42, 4115–4118 (2003). 8. Y. Chen and X. Li, “Dispersion management up to the third order for real-time optical coherence tomography involving a phase or frequency modulator,” Opt. Express 12, 5968– 5978 (2004). 9. D. L. Marks, A. L. Oldenburg, J. J. Reynolds, and S. A. Boppart, “Digital algorithm for dispersion correction in optical coherence tomography for homogeneous and stratified media,” Appl. Opt. 42, 204–216 (2003). 10. J. F. de Boer, C. E. Saxer, and J. S. Nelson, “Stable carrier generation and phase-resolved digital data processing in optical coherence tomography,” Appl. Opt. 40, 5787–5790 (2001). 11. A. F. Fercher, C. K. Hitzenberger, M. Sticker, R. Zawadzki, B. Karamata, and T. Lasser, “Numerical dispersion compensation for partial coherence interferometry and optical coherence tomography,” Opt. Express 9, 610–615 (2001). 12. M. Wojtkowski, V. J. Srinivasan, T. H. Ko, J. G. Fujimoto, A. Kowalczyk, and J. S. Duker, “Ultrahigh-resolution, highspeed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12, 2404–2422 (2004). 13. T. R. Hillman and D. D. Sampson, “The effect of water dispersion and absorption on axial resolution in ultrahighresolution optical coherence tomography,” Opt. Express 13, 1860–1874 (2005). 14. S. Diddams and J.-C. Diels, “Dispersion measurements with white-light interferometry,” J. Opt. Soc. Am. B 13, 1120–1129 (1996).
10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6653
