1.

Introduction

One of the main and most promising implementations of photonic investigations of structures is optical coherence tomography (OCT). This is a noninvasive interferometric imaging technique that provides high-resolution cross-sectional information of biological tissues and materials. First generation OCT systems, called time-domain OCT (TD-OCT), was reported in 1991 by Huang et al.¹ for imaging of the human retina in vitro. In TD-OCT, the reference arm scan depth is used to measure the time-of-flight of the optical signal reflected from the sample. Fourier domain OCT (FD-OCT) was introduced by Fercher et al.² In FD-OCT, the reference arm is fixed and the internal structure of the object is encoded in the spectral response of the sample field interfering with the reference beam and can be extracted by Fourier transformation. FD-OCT improves the detection sensitivity and enables significantly higher scan speeds than earlier TD-OCT.^3–5 The short acquisition time, the elimination of moving parts, and the inherently direct access to spectral information have made FD-OCT favored over TD-OCT.

Noncontact measurement of refractive index and thickness of OCT image data are of principal importance to biomedical applications and can greatly enhance diagnosis and treatment in medicine. For example, the refractive index and thickness of the cornea are related to hydration and intraocular pressure states and thus indicate the impact of laser refractive surgery.⁶ Also, a change in thickness map of human retina layers can be related to diseases such as glaucoma,⁷ diabetes,⁸ and several neuro-ophthalmic diseases such as Alzheimer’s, Parkinson’s, and multiple sclerosis.⁹ Therefore, accurate measurement of the retina layers true thickness and refractive index profile can help for early diagnosis and prognosis of the diseases.

Various techniques have been proposed for simultaneous measurement of refractive index and thickness using an OCT system in both time domain^10–32 and Fourier domain^33–37 suitable for one-layer objects. For more than one-layer objects, the focus-tracking technique has been used to measure the refractive index and thickness of two- and three-layer systems (the cornea)^38,39 and three layers of the skin⁴⁰ but the accuracy of this technique is usually limited by the resolution of translation stages. Also, the Fourier domain approach has been used for three-layer silica–crown–dental composite–silica structure.^41,42 In the proposed method, the real spectral data are fitted to a modeled signal spectrum searching for index and thickness. This approach is purely a fitting process without measuring index and thickness.

In the search for a new approach, in this article, we introduce a theoretical method for simultaneous index and thickness measurement of multilayer systems using Fourier domain OCT without having previous information about the item under investigation. The working hypothesis of this study is that there is sufficient information in the optical spectrum acquired by currently available implementations of FD-OCT to extract both the physical structure and the refractive index profile of a stratified medium.

2.

Theory

In this section, the new theoretical framework for simultaneous measurement of refractive index and thickness of transparent multilayer systems suitable for biological applications is introduced. More detailed explanation can be found in Refs. 43 and 44. For biological samples with small refractive index contrast, the FD-OCT interference spectrum can be modeled using the well-known summation method.⁴⁵ In this modeling approximation, multiple reflections inside the sample are ignored and only single scattering events are considered. For simplicity and without losing the generality, we use the framework to extract two-layer system parameters (refractive index and thickness of the layers and the last medium index of refraction) depicted in Fig. 1. We assume that the boundaries of this system are solid and ideally flat. We also assume that the object is nonabsorbing, isotropic, and homogeneous in refractive index.

Fig. 1

Schematic representation of a two-layer system.

For such system, the normalized interference spectrum free from constant (DC) and autocorrelation (AC) terms can be described by

Once the Fresnel coefficients are found, the refractive index of the layers and the last medium index of refraction may then be extracted directly, assuming the front medium index $n_0$ is known using Eq. (2). Thicknesses can be found by dividing the OPL by the corresponding index of refraction ($d_{1,2}={\delta }_{1,2}/n_{1,2}$).

Not all sets of three spectral components can be used to solve Eq. (3) since, for some sets, the determinant of the matrix ($P$) is zero and hence the equation is unsolvable. The roots of the matrix’s determinant indicate that the chosen wavenumber components do not provide linearly independent equations and are mainly related to the components that provide the same phase information from two or more interfaces. In an ideal situation in which the actual OPL is known exactly, the roots of the determinant are the only restriction on selecting wavenumber components. However, in practice, there is always uncertainty ($ϵ$) associated with the measured OPL (${\delta }^{\prime }$) so the value input into the parameter extraction algorithm may deviate from the actual OPL ($\delta$) so that

Since in practice, the matrix in Eq. (4) is constructed from the measured OPLs, there is a nonlinear frequency shift between the matrix and the signal on the right-hand side of Eq. (3). Therefore, in the presence of $ϵ$, there will be more constraints on wavenumber components and inappropriate selection of the wavenumber components leads to erroneous results. Given that this matrix equation presents a coupled set of transcendental equations containing index terms on both sides of the equal signs, the impact of $ϵ$ on index determination is not transparent so the behavior of the extracted parameters is best investigated by simulation. A framework for selecting the suitable spectral components in the presence of $ϵ$ was given in the first part of this article.^43,44 The impact of increasing $ϵ$ on the parameter extraction is investigated in this article to find out what is the maximum $ϵ$ that still allows the extraction of the parameters and how much is the accuracy loss in each parameter as $ϵ$ increases.

3.

Sample Parameters Extraction for a Range of $ϵ$

We investigate the impact of measured optical path error, $ϵ$, on the extracted parameters of a two-layer systems (i.e., $n_1$, $n_2$, $n_3$, $d_1$, $d_2$). We generate an experimental FD-OCT signal using transfer matrix method (TMM) and use the formalism to extract the parameters for a range of ${ϵ}_{1,2}$ induced to both layers measured OPL.

As an illustrative test model, we consider a low contrast test model with specification as $n_0=1.337$, $n_1=1.345$, $n_2=1.351$, $n_3=1.337$, and $d_1=10\mu \mathrm{m}$ and $d_2=30\mu \mathrm{m}$. We extract the parameters of this model for $ϵ$ ranging from ${ϵ}_{1,2}=0$ to ${ϵ}_{\mathrm{1,2}}=+10.00\mu \mathrm{m}$ for each layer with $0.1\mu \mathrm{m}$ increment and compare the extracted indices and thicknesses as per the methodology developed in the previous section. We have found that ${\mathrm{\Delta }}_{Z_R{Z}_0}$ has a neutral impact on the result and hence, for the tests in this article, we set ${\mathrm{\Delta }}_{Z_R{Z}_0}=0$.⁴³ The test result is shown in Fig. 2 where the error in the extracted parameters (vertical axis) is plotted versus $ϵ$ on the horizontal axis. The index error is shown on the left panel of Fig. 2 and the thickness error is shown on the right panel.

Fig. 2

Impact of the optical thickness measurement error, $ϵ$, on the accuracy of the extracted parameters. (a) The error of the extracted index of refraction and (b) the relative error in the extracted thicknesses.

It is seen from the right panel of Fig. 2 that even for large amounts of $ϵ$ comparable to the size of the object, the extracted refractive indices errors stay between $-0.024$ to $+0.01$ in this particular example. This is due to the fact that $ϵ$ affects only the matrix $P$ in Eq. (3) which is periodic and never diverges with $ϵ$. For higher contrast test models, the extracted index errors will be larger due to increase in multiple reflections which are ignored by the summation method. Since the physical thickness is extracted from the optical thickness divided by the corresponding index of refraction, the error of thickness on the right panel of Fig. 2 is related to $ϵ$ and error of index via

Therefore, the relative error in extracted thicknesses has steeper slope for smaller OPL.

4.

Approach to Parameter Improvement

It was shown in the previous section that in a realistic case of having no zero $ϵ$ our methodology can quickly result in significant errors. In order to enhance the parameters extraction and deal with large values of $ϵ$ which may occur in practice, an additional process is needed. We use an optimization method by constructing a modeled FD-OCT normalized signal similar to Eqs. (1) and (2) and fit the modeled signal to the actual data obtained from the digitally constructed TMM. In other words, we try to fit a summation (single scattering) spectrum to the TMM (multiple scattering) spectrum. Our modeled signal is

Accurate determination of floating parameters $b_4$ and $b_5$ is very important to achieve the best fit between the modeled and actual spectrum because they describe the frequency of the spectrum. During the search process, when the initial guess values are not close enough to the actual values, it is more likely that the algorithm is stopped at a local minimum of the residual root sum square (RSS) of difference between the modeled and actual spectrum. In this situation, the floating parameters cannot be found accurately. In order to prevent the fitting algorithm from stopping at local minima of RSS, the program randomly sweeps $b_4$ and $b_5$ within a range around the initial guess values and calculates the RSS for each set of $b_4$ and $b_5$ to search for the lowest RSS. Therefore, our fitting process is composed of three steps:

Since the algorithm must search for all possible combinations of ${\mathrm{OPL}}_1$ and ${\mathrm{OPL}}_1+{\mathrm{OPL}}_2$ within a given interval, step 2 of the algorithm given above randomly changes $b_4$ and $b_5$ and checks the RSS until the global minimum RSS is found. The algorithm in this sense is a simple implementation of a Monte Carlo search. The additional step provided in the searching algorithm forces the algorithm to jump over local minima and check for more suitable values of RSS; therefore, it is more likely to find the best estimate available. As with any Monte Carlo algorithm, the number of points searched is important to ensure finding the global minimum of the RSS. The step 3 mentioned above provides the final extracted indices $(n_1,n_2,n_3)$ and optimized OPL (${\mathrm{OPL}}_1$, ${\mathrm{OPL}}_1+{\mathrm{OPL}}_2$) from the global minimum RSS that was found in the previous step. Finally, the thicknesses are obtained from ${\mathrm{OPL}}_1/n_1$ and $[({\mathrm{OPL}}_1+{\mathrm{OPL}}_2)-{\mathrm{OPL}}_1]/n_2$.

4.1.

Optimized Results for Range of $ϵ$

We are now going to investigate the optimization method for a range of $ϵ$. We recall the test model used in Sec. 2 of this article for $ϵ$ ranging from ${ϵ}_{1,2}=0$ to ${ϵ}_{\mathrm{1,2}}=+10\mu \mathrm{m}$ for both layers with $0.4\mu \mathrm{m}$ increment. The simulation program was asked to increase the number of searches until all indices $n_1$, $n_2$, and $n_3$ are found with the absolute error $\le 0.001$. The test result is shown in Fig. 3. On the left axis of this figure is shown the extracted refractive index errors and on the right axis is shown the estimated number of searches required versus $ϵ$ on the horizontal axis. The error in this figure is the deviation of fitting parameters $\{b_1,b_2,b_3\}$ from their true values $\{n_1,n_2,n_3\}$.

Fig. 3

Impact of the number of searches on the accuracy of the extracted indices of the optimization method. On the left axis is shown the error of the extracted refractive indices and on the right axis is shown the required number of searches.

It can be seen from Fig. 3 that once we have enough numbers of searches, even for the very large amount of $ϵ$, there is a convergence to highly accurate refractive index values. For a limited number of searches, the range of $ϵ$ for which we can extract parameters with the acceptable error would be limited.

4.2.

Index Contrast

In this section, we investigate the impact of refractive index profile on two-layer system parameter errors extracted by the optimization method while the method is insensitive to thickness profile as described in Part 1, Sec. 3.2.⁴⁴ Both layers thickness is fixed at $30\mu \mathrm{m}$ and the optical paths measurement error of ${ϵ}_{1,2}=+0.1\mu \mathrm{m}$ is induced for both layers. The test models are generated with $n_1$ and $n_2$ varying from 1.337 to 1.9 in 0.011 increments (total number of 2600 test models) while $n_3=1.337$ but must still be extracted by the software. The test result is shown in the following three figures where the horizontal axis corresponds to actual $n_1$ and the vertical axis corresponds to actual $n_2$, and the color maps are the extracted parameter errors without any optimization [direct output of Eq. (3), the right panel] and after performing the optimization method (the left panel).

The error distribution corresponding to the extracted $n_1$ is shown in Fig. 4. It can be seen from Fig. 4(b) of that $n_1$ can be extracted by the optimization method with absolute error $\le 0.001$, suitable for biological applications, in the entire $n_1:n_2$ space tested. Without optimization, Fig. 4(a), we observe very poor extraction of $n_1$ where error varies between $1.92\times {10}^{-5}$ and 2.13. The figure shows that the optimization method improves greatly the usefulness of this solution.

Fig. 4

Error distribution of the extracted $n_1$ with ${ϵ}_{1,2}=+0.1\mu \mathrm{m}$. The horizontal axis corresponds to varying $n_1$ and the vertical axis corresponds to varying $n_2$. (a) The error map of initial (guess) values corresponding to $n_1$ and (b) the error map of extracted $n_1$ by the optimization method.

The error distribution corresponding to the extracted $n_2$ is shown in Fig. 5. From Fig. 5(b), $n_2$ can be extracted by the optimization method with absolute error $\le 0.001$ in area limited by $n_1<1.60$ and $n_2<1.80$. Without optimization, Fig. 5(a), error in $n_2$ varies between $4.57\times {10}^{-5}$ and 3.18.

Fig. 5

Error distribution of the extracted $n_2$ with ${ϵ}_{\mathrm{1,2}}=+0.1\mu \mathrm{m}$. The horizontal axis corresponds to varying $n_1$ and the vertical axis corresponds to varying $n_2$. (a) The error map of initial (guess) values corresponding to $n_2$ and (b) the error map of extracted $n_2$ by the optimization method.

The error distribution corresponding to the extracted $n_3$ is shown in Fig. 6. From Fig. 6(b), $n_3$ can be extracted by the optimization method with absolute error $\le 0.001$ in area limited by $n_1<1.55$ and $n_2<1.60$. Without optimization, error in $n_3$ varies between $1.90\times {10}^{-2}$ and 1.78 and never reaches the targeted accuracy.

Fig. 6

Error distribution of the extracted $n_3$ with ${ϵ}_{1,2}=+0.1\mu \mathrm{m}$. The horizontal axis corresponds to varying $n_1$ and the vertical axis corresponds to varying $n_2$. (a) The error map of initial (guess) values corresponding to $n_3$ and (b) the error map of extracted $n_3$ by the optimization method.

5.

Conclusion

In conclusion, we introduced in this article (together with its first part)⁴⁴ a framework for simultaneous extraction of refractive index and thickness of transparent multilayer systems without any previous information about the item under investigation from its FD-OCT data. No additional arrangement, reference reflector, or mechanical scanning is needed, which is the main advantage of this approach.

It was shown that the measured OPL uncertainty has a significant impact on parameters extraction, but this factor can be overcome through the introduced optimization method even for the relatively large amount of uncertainty in measured OPL comparable to the size of the object. The optimization approach enhances the final results superlatively and the simulation shows that the potential accuracy of the method is high enough for biological applications.

The accuracy of the method depends mainly on the refractive index contrast of the object while it is independent of the sample’s physical thickness. This is because in the theoretical model, only single scattering events are taken into account and this approximation is not efficient for large index contrast objects. It was also shown that using the optimization method, all parameters of a two-layer system and the last medium index of refraction can be obtained with the absolute error $\le 0.001$ when $n_{1,2}<1.55$. This implies that the maximum tolerable index contrast of the optimization method is in the range of $|1.55-1.337|=0.213$. Therefore, this method can be suitable for biological applications in which the refractive index rarely exceeds 1.45.

Some practical implementation limitations need to be investigated more deeply to better understand the real limits of the solution. The two most important are likely speckle noise and the impact of the focusing optical beam (nonuniform sample illumination) which are left for future work.