A Quantitative Model for Optical Coherence Tomography

Abstract

Optical coherence tomography (OCT) is a widely used imaging technique in the micrometer regime, which gained accelerating interest in medical imaging in the last twenty years. In up-to-date OCT literature, certain simplifying assumptions are made for the reconstructions, but for many applications, a more realistic description of the OCT imaging process is of interest. In mathematical models, for example, the incident angle of light onto the sample is usually neglected or a plane wave description for the light–sample interaction in OCT is used, which ignores almost completely the occurring effects within an OCT measurement process. In this article, we make a first step to a quantitative model by considering the measured intensity as a combination of back-scattered Gaussian beams affected by the system. In contrast to the standard plane wave simplification, the presented model includes system relevant parameters, such as the position of the focus and the spot size of the incident laser beam, which allow a precise prediction of the OCT data. The accuracy of the proposed model—after calibration of all necessary system parameters—is illustrated by simulations and validated by a comparison with experimental data obtained from a 1300 $\mathrm{n}$$\mathrm{m}$ swept-source OCT system.

1. Introduction

Optical coherence tomography (OCT) has proved to be a non-invasive, high-precision imaging technique with micrometer resolution. It emerged around 1990 for in vivo imaging of the human eye [1,2] and gained increasing interest ever since [3]. Today, extensions, such as angiography [4], polarization sensitive OCT [5] and optical coherence elastography [6], unlocked a wide range of possible applications; for example, blood vessel analysis [7] and cancer margin detection [8], while OCT endoscopes [9] are clearing the way for high-resolution imaging of internal organs and their pathologies. Multi-modal imaging techniques [10] often use OCT as morphological guidance.

While many theoretical OCT articles assume a sample geometry that is perfectly perpendicular to the OCT beam [11,12,13], commonly used OCT systems are designed for rough sample surfaces and arbitrary sample inclinations, which yield much less power at the detector. Normal incidence not only oversaturates the detector easily, especially for samples with a high refractive index and a very directed scattering profile, but can also lead to interference between the sample and optical parts inside the setup, e.g., the scan lens. To prevent imaging artifacts, normal incident is therefore usually avoided in OCT.

In addition, most works are based on a plane wave ansatz for describing the sample and reference fields [14,15,16,17]. While this approximation is valid for the immediate focus region, it is most often not true for the whole field of view of the setup or even the whole sample area. Common effects such as a focus-dependent intensity profile inside the sample cannot be described with a plane wave ansatz. While workarounds such as multiplying the spectral resolution with a sensitivity factor have been proposed [18], in this work, a Gaussian beam, similar to [19], is used as incoming wave for a more precise description of the light beams. Besides analytical models for simulations of an OCT signal, there exist also Monte Carlo methods [20,21,22], recent ones also modeling the light as Gaussian beams [23].

To make a quantitative reconstruction of the optical parameters of the sample possible, we make the simplifying assumption that the sample is not absorbing and can at least locally be described as a layered medium (with layers not necessarily perfectly perpendicular to the incident light), which is a classical assumption in this field. This simplification allows us to analytically calculate the scattered light from the sample, which is then collected by a scan lens and combined with the reference light to produce the interference pattern. We roughly model the effect of the scan lens on the scattered light by discarding plane wave components moving in wrong directions. Together with the layers and the focusing effects introduced by Gaussian beams, we derived simulations which have been in much better accordance with the experimental data obtained by a 1300 $\mathrm{n}$$\mathrm{m}$ swept-source OCT system compared to a simple plane wave approach.

The paper is structured as follows: in Section 2, we describe the OCT imaging system that was built for this work. The individual experiments performed for investigating the different effects of the system on the data are also introduced. In Section 3, we present the general problem and the governing equations. We give the forms of the sample and reference fields using the near- and far-field representation of back-scattered Gaussian waves. The section ends with the derivation of the formula for the measurement data. In Section 4, we present in two experiments the dependence of the data on the incident angle and the beam focusing and show how to use them as calibration tools to determine otherwise unknown parameters such as the beam radius in the focus. The comparison between experimental and simulated data is presented in the last section. There, we see that our model nicely predicts the behavior of the experimental data with respect to different orientations and positions of the sample relative to the focus. The proofs to all theoretical results presented in this article can be found in the Appendix A.

2. Experiments

Since there are many different variants of OCT systems around, we briefly describe the system we use for generating the data and the individual experiments.

2.1. OCT Setup and Post-Processing

All measurements where performed with a custom-built fiber-based OCT system with a central frequency of about 1300 $\mathrm{n}$$\mathrm{m}$ and 30 $\mathrm{n}$$\mathrm{m}$ bandwidth. It is built upon an akinetic swept-source from Insight Photonic Solutions, USA, which emits 60 $\mathrm{m}$$\mathrm{W}$ with a maximum repetition rate of 500 $\mathrm{k}$$\mathrm{Hz}$. This laser was selected for its high phase stability and flat power profile over its whole bandwidth. This makes it well-suited for any type of signal analysis.

As shown in Figure 1, the light is split into reference (25%) and sample arm (75%) by a fiber coupler and is released into free-space by 4 $\mathrm{m}$$\mathrm{m}$ fiber collimators. The sample arm (red) features a rotatable imaging probe with conjugated scanning and an LSM54-1310 scan lens (Thorlabs, Newton, NJ, USA) for a flat imaging plane. In the reference arm (yellow), the light is reflected by a mobile mirror. Circulators guide the reference and sample arm signal to a 50/50% fiber coupler for recombination. A dual-balance-detector (BPD-1, Insight Photonic Solutions, Lafayette, LA, USA), short DBD, records the cross-correlation term and an ATS9360 data acquisition card from Alazar Technologies, Canada, is used to digitize it. To prevent the saturation of the detector, an adjustable neutral density filter is used in the reference arm. The system achieves an axial resolution of 31 $\mathsf{\mu }\mathrm{m}$ and a lateral resolution of 24 $\mathsf{\mu }\mathrm{m}$ in air as well as a SNR of 105 dB. To ensure satisfactory oversampling, a lateral pixelsize of $9.8$ $\mathsf{\mu }\mathrm{m}$ in air is selected.

The source supplies a trigger signal used to coordinate the galvanometer movement and data acquisition. The OCT control software written in Labview and data processing is performed in MATLAB. The recorded spectrum consisting of 700 datapoints is already spaced equally in terms of wavenumber and thus can immediately be Fourier-transformed into image space. The small bandwidth of the laser makes dispersion compensation unnecessary. Zeropadding ensures the small axial pixelsize of $13.7$ $\mathsf{\mu }$$\mathrm{m}$ in air. Background removal is performed via subtraction of the average spectrum of a volume. The DBD records only the difference between both input signals, thereby removing common-mode noise and centering the signal around zero. Through the digitalization process, the signal is shifted in height, so it can be stored as unsigned integer. This shift is removed again during post-processing through the background subtraction.

2.2. Power vs. Angle

To investigate the influence of the incident angle between sample surface and OCT beam, one of the fibers entering the DBD was connected to an optical power meter (PM100C with S122C, Thorlabs, Newton, NJ, USA) instead. A mirror was fixed onto a goniometer stage with a Vernier scale (GOH-65A100RUU, OptoSigma, Santa Ana, CA, USA) with a kinematic mount. First, the goniometer was aligned to ensure that the OCT beam goes through its Pivot point. Second, the kinematic mount and a vertical stage were adjusted to put the mirror in focus and ensure normal incident of the laser beam on the mirror, using the power meter as guidance. Then, the mirror was tilted in 5 arc minute steps back and forth between $-1°$ and 1° and the power was recorded until each angular position was measured 6 times.

2.3. Power vs. Focus

For quantification of the Gaussian behavior of the focus, a motorized stage (T-LSM050A, Zaber Technologies, Vancouver, BC, Canada) was used to transport once a mirror and once a microscopy coverglass through the focus of the OCT system, with a roughly fixed tilt of about $2.75°$. The coverglass (631-0124, VWR International, Radnor, PA, USA) has a refractive index of $1.5088$ for 1300 $\mathrm{n}$$\mathrm{m}$, which needs to be taken into account during data analysis. At each position of the stage, we use 11 steps for the mirror and 7 for the coverglass, a 3D OCT volume was recorded. These 3D volumes were post-processed according to Section 2.1 and used to determine the exact incident angle and the distance of the center of the sample surface to the position where sample and reference arm would have the same length, called zero delay.

3. Mathematical Model

Considering the workflow of the used OCT system described in the previous section, we model the parts shown in Figure 2 separately.

Firstly, in Section 3.1, we model the produced laser illumination. The laser light is split into two beams, one is sent to the sample and the other to the mirror in the reference arm. We model their scattering process in Section 3.2.1 and Section 3.2.5, respectively.

The reflected light in the sample arm is (partially) collected by a scan lens and coupled into a fiber. This aspect is the topic of Section 3.2.4.

After recombination of the scattered light beams, we model the detection via a dual-balancing detector of this superposition in Section 3.3. This in particular is discussed for the measurement related to two experiments explained in Section 2, which, in the end, are obligatory for the calibration of necessary parameters in the forward simulations.

3.1. Gaussian Beam Illumination

The shape of the light produced inside an optical resonator (we ignore at this point the finite size of the resonator and the boundary conditions) of a laser can according to [24], for example, be well described by a Gaussian beam.

We consider a Gaussian beam $\mathbf{E}:{\mathbb{R}}^3\to {\mathbb{C}}^3$ as a monochromatic solution of the electromagnetic wave equation in vacuum which reduces it to Helmholtz equation (usually it is considered as solution of the paraxial approximation of the wave equation which is not necessary here):

$$\begin{array}{ccc}\hfill \Delta \mathbf{E}\left(\mathbf{x}\right)+k_0^2\mathbf{E}\left(\mathbf{x}\right)& =0,\hfill & \hfill \mathbf{x}\in {\mathbb{R}}^3,\\ \hfill \langle \nabla ,\mathbf{E}\rangle \left(\mathbf{x}\right)& =0,\hfill & \hfill \mathbf{x}\in {\mathbb{R}}^3.\end{array}$$

(1)

It is characterized by its form

$$\mathbf{E}(x_1,x_2,r_0)=f(x_1,x_2)\mathbf{p}$$

(2)

in the focal plane $\{\mathbf{x}\in {\mathbb{R}}^3|x_3=r_0\}$ for a function $f:{\mathbb{R}}^2\to \mathbb{C}$ such that its 2D Fourier transform is compactly supported in $D_{k_0}\left(0\right)$ (the open ball with center 0 and radius $k_0$) and a polarization vector $\mathbf{p}\in {\mathbb{R}}^2\times \left\{0\right\}.$

Theorem 1.

Let $f:{\mathbb{R}}^2\to \mathbb{C}$ be a function such that its two-dimensional Fourier transform $\check{f}$ is compactly supported in $D_{k_0}\left(0\right)$ and let $\mathbf{p}\in {\mathbb{R}}^2\times \left\{0\right\}$. Then, for every $\mathbf{x}\in {\mathbb{R}}^3$ a solution of the Helmholtz problem (1) is given by

$$\begin{array}{c}\mathbf{E}\left(\mathbf{x}\right)=\frac{1}{4{\pi }^2}{\int }_{{\mathbb{R}}^2}{\check{\mathit{g}}}_{+}(k_1,k_2)e^{-i(k_1{x}_1+k_2{x}_2)}{e}^{-i\sqrt{k_0^2-(k_1^2+k_2^2)}(r_0-x_3)}d(k_1,k_2)\hfill \\ \hfill +\frac{1}{4{\pi }^2}{\int }_{{\mathbb{R}}^2}{\check{\mathit{g}}}_{-}(k_1,k_2)e^{-i(k_1{x}_1+k_2{x}_2)}{e}^{i\sqrt{k_0^2-(k_1^2+k_2^2)}(r_0-x_3)}d(k_1,k_2),\end{array}$$

(3)

with

$${\check{\mathit{g}}}_{\pm }(k_1,k_2)=\frac{1}{2}\check{f}(k_1,k_2)\left(\begin{array}{c}{p}_1\\ p_2\\ \pm \frac{p_1{k}_1+p_2{k}_2}{\sqrt{k_0^2-(k_1^2+k_2^2)}}\end{array}\right).$$

(4)

Such a wave describes well the light inside the optical resonator of the laser. Then, through one partly transparent mirror of the resonator, we obtain only the light moving in the negative $x_3$-direction of the form

$${\mathbf{E}}^{\left(0\right)}\left(\mathbf{x}\right)=\frac{1}{4{\pi }^2}{\int }_{{\mathbb{R}}^2}\check{\mathit{g}}(k_1,k_2)e^{-i(k_1{x}_1+k_2{x}_2)}{e}^{-i\sqrt{k_0^2-(k_1^2+k_2^2)}(r_0-x_3)}d(k_1,k_2),$$

(5)

with $\check{\mathit{g}}={\check{\mathit{g}}}_{+}.$ Hereby, a reasonable model for the shape of the function f is one which resembles a Gaussian function.

This laser light enters the single-mode fiber-based system. The good approximability of the laser light after exiting the single-mode fiber by a Gaussian beam shape allows us to assume that the Gaussian behavior is conserved throughout the system.

3.2. Backscattered Gaussian Fields

The laser light is split into two waves, ${\mathbf{E}}_S^{\left(0\right)}$ for the sample and $i{\mathbf{E}}_R^{\left(0\right)}$ for the reference arm, as in Equation (5), respectively, by a beam splitter and both remain in the form of a Gaussian beam, with possibly different beam parameters, for example, due an optical attenuation wheel inside the reference arm which causes a difference in light intensities between the beams.

3.2.1. The Sample Field

The beam ${\mathbf{E}}_S^{\left(0\right)}$ is now directed onto the sample and we say it is of form (5) with $f=f_S.$ Then, if the beam is sufficiently focused, meaning that the values of $|{\mathbf{E}}_S^{\left(0\right)}|$ can be neglected outside a small region, we only need to consider for the scattering process the shape of the sample inside this region. In this subregion, we denote it by $\Omega$, we assume, using the tangent plane approximation, that it can be described by a layered structure. These layers are not necessarily perpendicular to incident beam, but for simplification assumed to be parallel to each other. This is modeled by $\Omega$ being a finite union of layers:

$$\Omega =\bigcup _{j=1}^L{\Omega }_j,{\Omega }_j=\{\mathbf{x}\in \Omega |a_j\le \langle \mathbf{x},{\mathbf{\nu }}_{\Omega }\rangle <a_{j+1}\},{\left(a_j\right)}_{j=1}^L\subset \mathbb{R},$$

for some unit normal vector ${\mathbf{\nu }}_{\Omega }.$ Each of these shall be characterized by a constant refractive index $n_j\in [1,\infty ).$

Under these conditions, we model the backscattered field ${\mathbf{E}}_S$ as solution of the Helmholtz equation

$$\Delta ({\mathbf{E}}_S+{\mathbf{E}}_S^{\left(0\right)})\left(\mathbf{x}\right)+k_0^2{n}^2\left(\mathbf{x}\right)({\mathbf{E}}_S+{\mathbf{E}}_S^{\left(0\right)})\left(\mathbf{x}\right)=0,\mathbf{x}\in {\mathbb{R}}^3$$

(6)

where $n\left(\mathbf{x}\right)={\sum }_{j=1}^L{n}_j{\chi }_{{\Omega }_j}\left(\mathbf{x}\right)+{\chi }_{{\mathbb{R}}^3\backslash \Omega }\left(\mathbf{x}\right)$ and appropriate radiation conditions are assumed.

The incident field (5) is represented as a superposition of plane waves having different wave vectors. Because of the linearity of the equation it is sufficient to solve the problem for every plane wave. The result for these backscattered fields for such a sample is well known in this plane wave case, see [25,26]. For the simplest case $L=1,$ we consider an (arbitrary) plane wave as incident illumination from the top,

$${\mathbf{E}}_S^{\left(0\right),\mathrm{pl}}\left(\mathbf{x}\right)=\mathbf{\alpha }(k_1,k_2)e^{-i\langle \mathbf{k},\mathbf{x}\rangle },\mathbf{x}\in {\mathbb{R}}^3,$$

with amplitude function $\mathbf{\alpha }:{\mathbb{R}}^2\to {\mathbb{C}}^3$ and propagation vector

$$\mathbf{k}=\left(\begin{array}{c}{k}_1\\ k_2\\ -\sqrt{k_0^2-k_1^2-k_2^2}\end{array}\right),\left|\mathbf{k}\right|=k_0,$$

(7)

which we consider implicitly as a function of $k_1$ and $k_2.$ We obtain the reflected electric field

$${\mathbf{E}}_S^{\mathrm{pl}}\left(\mathbf{x}\right)=\mathit{\beta }(k_1,k_2)\mathbf{\alpha }(k_1,k_2)e^{-i\langle (\mathbf{k}-{\mathbf{k}}_r),{\mathbf{x}}_{\Omega }\rangle }{e}^{-i\langle {\mathbf{k}}_r,\mathbf{x}\rangle },$$

(8)

where ${\mathbf{x}}_{\Omega }$ denotes an arbitrary point of the top boundary (that is $\langle {\mathbf{x}}_{\Omega },{\mathbf{\nu }}_{\Omega }\rangle =a_1$) of the object,

$$\Phi :{\mathbb{R}}^3\to {\mathbb{R}}^3,{\mathbf{k}}_r=\Phi \left(\mathbf{k}\right)=\mathbf{k}-2\langle \mathbf{k},{\mathbf{\nu }}_{\Omega }\rangle {\mathbf{\nu }}_{\Omega }$$

(9)

the wave vector and ${\mathit{\beta }}_S$ the sum of the reflection coefficients ${\beta }_{0,\Vert },{\beta }_{0,\perp }$ of the differently polarized parts

$${\mathit{\beta }}_S(k_1,k_2)={\mathsf{\lambda }}_1{\beta }_{0,\Vert }{\mathbf{p}}_{\Vert }(k_1,k_2)+{\mathsf{\lambda }}_2{\beta }_{0,\perp }{\mathbf{p}}_{\perp }(k_1,k_2).$$

(10)

Here, we have decomposed $\mathbf{\alpha }$ into its transverse electric and magnetic polarizations, with coefficients ${\mathsf{\lambda }}_1$ and ${\mathsf{\lambda }}_2,$ respectively. Further, we use Snell’s law for the determination of the transmission angle ${\theta }_t.$

Summarizing the scattered (plane) waves for all $(k_1,k_2)$ and

$$\mathbf{\alpha }(k_1,k_2)={\check{\mathit{g}}}_S(k_1,k_2)e^{-i\sqrt{k_0^2-k_1^2-k_2^2}{r}_0}$$

then finally results in

$$\begin{array}{cc}\hfill {\mathbf{E}}_S\left(\mathbf{x}\right)& =\frac{1}{4{\pi }^2}{\int }_{{\mathbb{R}}^2}{\mathbf{E}}_S^{\mathrm{pl}}\left(\mathbf{x}\right)d(k_1,k_2)\hfill \\ \hfill & =\frac{1}{4{\pi }^2}{\int }_{{\mathbb{R}}^2}{\mathit{\beta }}_S(k_1,k_2){\check{f}}_S(k_1,k_2)e^{-i\sqrt{k_0^2-k_1^2-k_2^2}{r}_0}{e}^{-i\langle (\mathbf{k}-{\mathbf{k}}_r),{\mathbf{x}}_{\Omega }\rangle }{e}^{-i\langle {\mathbf{k}}_r,\mathbf{x}\rangle }d(k_1,k_2),\hfill \end{array}$$

(11)

with $\mathit{\beta }$ given by Equation (10).

3.2.2. Far Field Method

Since the distance between the scan lens and the sample (which is roughly 6 $\mathrm{c}$$\mathrm{m}$) is much greater than the size of of the sample itself (which is only a few millimeters), we can be tempted to simplify the integrand by using the far-field approximation.

Mathematically, this means, that we are approximating Equation (11) by its behavior at some point $r\mathbf{s},\mathbf{s}\in {\mathbb{S}}^2,$ as $r\to \infty :$

$${\mathbf{E}}_S\left(r\mathbf{s}\right)={\mathbf{E}}_{S,\infty }\left(r\mathbf{s}\right)+o(1/r)$$

To compute the dominating term ${\mathbf{E}}_{S,\infty }$, we apply the method of stationary phase, see Lemma A1, which is based on the approximation of the phase function $k_0\Psi ,$ with

$$\Psi (k_1,k_2)=\langle \frac{{\mathbf{k}}_r}{k_0},\mathbf{s}\rangle =\langle \frac{\mathbf{k}}{k_0},\mathbf{s}\rangle -2\langle \frac{\mathbf{k}}{k_0},{\nu }_{\Omega }\rangle \langle {\nu }_{\Omega },\mathbf{s}\rangle ,$$

by its Taylor series around its critical points.

Theorem 2.

Let ${\mathbf{E}}_S$ be a vector field given by Equation (11). Then, its far field approximation takes the form

$${\mathbf{E}}_{S,\infty }\left(r\mathbf{s}\right)=\frac{-i{k}_0\left|c_3\right|}{2\pi r}{\mathit{\beta }}_S(k_1,k_2){\check{f}}_S(k_1,k_2)e^{-i\sqrt{k_0^2-k_1^2-k_2^2}{r}_0}{e}^{-i\langle \mathbf{k}-{\mathbf{k}}_r,{\mathbf{x}}_{\Omega }\rangle }{e}^{i{k}_{0}sign\left(c_3\right)r},$$

(12)

where for $c_3=s_3-2\langle {\nu }_{\Omega },\mathbf{s}\rangle {\nu }_{\Omega ,3},$ the vector components $k_1$ and $k_2$ are given by

$$\left(\begin{array}{c}{k}_1\\ k_2\end{array}\right)=-k_{0}sign\left(c_3\right)\left(\begin{array}{c}{s}_1-2\langle {\nu }_{\Omega },\mathbf{s}\rangle {\nu }_{\Omega ,1}\\ s_2-2\langle {\nu }_{\Omega },\mathbf{s}\rangle {\nu }_{\Omega ,2}\end{array}\right),$$

(13)

and for $k_1,k_2$ the reflected vector ${\mathbf{k}}_r=\Phi \left(\mathbf{k}\right)$ is given by Equation (9).

3.2.3. Comparing the Near- and the Far-Fields

Comparing the representations (11) and (12) of the scattered electric field simulations, it becomes obvious that there is a difference in the visible effects provided by these methods. In this work, we are concerned with the influence of the focus in the scattered field.

In contrast to the far-field representation, the scattered field in the near field regime is heavily depending on the distance between the positions of the layer and of the focus. We neglect the vectorial quantities in Equation (11) for a moment and allow for an amplitude function

$$\check{f}(k_1,k_2)\sim e^{-a(k_1^2+k_2^2)},$$

(14)

where the parameter $a>0$ is such that the error $|\check{f}-\check{f}{\chi }_{D_{{\rho }_0}\left(0\right)}|$ is negligible. Hereby, $D_{{\rho }_0}\left(0\right)$ is a disk with small radius ${\rho }_0$ and center $0.$ Then, for a single surface (medium), we obtain the scattered field

$$\mathbf{E}\left(\mathbf{x}\right)=\frac{1}{4{\pi }^2}{\int }_{{\mathbb{R}}^2}{\beta }_0(k_1,k_2)e^{-a(k_1^2+k_2^2)}{e}^{-i\sqrt{k_0^2-k_1^2-k_2^2}{r}_0}{e}^{-i\langle (\mathbf{k}-{\mathbf{k}}_r),{\mathbf{x}}_{\Omega }\rangle }{e}^{-i\langle {\mathbf{k}}_r,\mathbf{x}\rangle }d(k_1,k_2).$$

We assume that on the small disk, the reflection coefficient ${\beta }_0$ is approximately constant and due to small deviations of $k_1,k_2$ from zero, we may approximate the root in the exponents

$$\sqrt{k_0^2-k_1^2-k_2^2}\simeq k_0-\frac{k_1^2+k_2^2}{2k_0}.$$

(15)

Further, for the sake of simplification, we restrict ${\mathbf{\nu }}_{\Omega }\in {\mathbb{R}}^2\times \left\{0\right\},{\nu }_{\Omega ,2}=0,$ fix the positions of the focus $r_0$ and the object ${\mathbf{x}}_{\Omega }=x_{\Omega ,3},x_{\Omega ,3}<0$ below the origin and evaluate at $\mathbf{x}=\mathbf{0}.$ This finally gives for the scattered field

$$\mathbf{E}\left(\mathbf{0}\right)=\frac{{\beta }_0}{4{\pi }^2}{e}^{-i{k}_0{\psi }_0}{\int }_{\mathbb{R}}{e}^{-{\psi }_2{k}_1^2}{e}^{-i{k}_1{\psi }_1}d{k}_1{\int }_{\mathbb{R}}{e}^{-{\psi }_2{k}_2^2}d{k}_2,$$

(16)

where we defined the phase elements ${\psi }_0,{\psi }_1,{\psi }_2$ by

$${\psi }_2\left(k_0\right)=a+\frac{i}{k_0}\left({\nu }_{\Omega ,3}^2{x}_{\Omega ,3}-\frac{r_0}{2}\right),{\psi }_1=2{\nu }_{\Omega ,1}{\nu }_{\Omega ,3}{x}_{\Omega ,3},{\psi }_0=r_0-2{\nu }_{\Omega ,3}^2{x}_{\Omega ,3}.$$

(17)

Since $\Re e\left({\psi }_2\right)=a>0,$ we evaluate both integrals in Equation (16) and arrive at

$$\mathbf{E}\left(\mathbf{0}\right)=\frac{{\beta }_0}{4\pi {\psi }_2\left(k_0\right)}{e}^{-\frac{{\psi }_1^2}{4{\psi }_2\left(k_0\right)}}{e}^{-i{k}_0{\psi }_0}.$$

After complex conjugation in the exponent and taking the absolute value of the field, we find that

$$\left|\mathbf{E}\left(\mathbf{0}\right)\right|=\frac{|{\beta }_0|k_0}{4\pi \sqrt{k_0^2{a}^2+d^2}}{e}^{-\frac{k_0^2{\psi }_1^{2}a}{4(k_0^2{a}^2+d^2)}},$$

(18)

with distance $d={\nu }_{\Omega ,3}^2{x}_{\Omega ,3}-\frac{r_0}{2}.$ Considering now Equation(18) for different positions $x_{\Omega ,3}$ of the sample, and therefore for varying $d,$ corresponds to different evaluation points ${\mathbf{x}}_{\infty }=r{\mathbf{e}}_3,$ with $r=|x_{\Omega ,3}|$ in the far-field regime. Taking the absolute value of Equation (12)

$$|{\mathbf{E}}_{\infty }\left(r{\mathbf{e}}_3\right)|\simeq \frac{k_0|{\beta }_0\left|\right|c_3|}{2\pi r}{e}^{-a(k_1^2+k_2^2)}$$

we observe that opposed to the near-field representation, the far-field regime is independent of the focus position. Figure 3 provides a comparison between both for different positions of the focus and the surface.

Since the far-field approximation does not show the dependence of the electric field on the focus, we stick with the more complicated near-field representation of the scattered light ${\mathbf{E}}_S$ in Equation (11).

3.2.4. The Scan Lens

The backreflected light ${\mathbf{E}}_S$ then passes trough the scan lens and is collected by a fiber collimator. Thereby, we loose all plane wave components whose propagation directions are outside a certain angular range of the collimator. We model this by reducing the area of integration in Equation (11) to a set $\mathcal{B}$ of those scattered wave directions ${\mathbf{k}}_r$ which have an angle to the measurement direction ${\mathbf{e}}_3$ less than a certain angle of acceptance ${\theta }_{max},$ that is

$$\mathcal{B}=\left\{(k_1,k_2)\in {\mathbb{R}}^2|arccos\left(\frac{\langle {\mathbf{k}}_r,{\mathbf{e}}_3\rangle }{k_0}\right)\le {\theta }_{max}\right\},$$

(19)

which finally gives a (scattered) sample field

$$\begin{array}{c}{\mathbf{E}}_S^{\left(1\right)}\left({\mathbf{x}}_{\mathcal{D}}\right)=\frac{1}{4{\pi }^2}\times \hfill \\ \hfill {\int }_{\mathcal{B}}{\mathit{\beta }}_S(k_1,k_2){\check{f}}_S(k_1,k_2)e^{-i\sqrt{k_0^2-k_1^2-k_2^2}{r}_0}{e}^{-i\langle (\mathbf{k}-{\mathbf{k}}_r),{\mathbf{x}}_{\Omega }\rangle }{e}^{-i\langle {\mathbf{k}}_r,{\mathbf{x}}_{\mathcal{D}}\rangle }d(k_1,k_2).\end{array}$$

(20)

at the scan lens position ${\mathbf{x}}_{\mathcal{D}}=r_{\mathcal{D}}{\mathbf{e}}_3$ above the sample.

3.2.5. The Reference Field

Similarly, we model the reference field as solution to the scattering problem (6) with Gaussian incident illumination ${\mathbf{E}}_R^{\left(0\right)}$ and a medium with constant (infinitely) large refractive index. Following the same line that led to Equation (11), we obtain with $f=f_R$ a field of the form

$${\mathbf{E}}_R\left({\mathbf{x}}_{\mathcal{D}}\right)=\frac{1}{4{\pi }^2}{\int }_{D_{k_0}\left(0\right)}{\mathit{\beta }}_R(k_1,k_2){\check{f}}_R(k_1,k_2)e^{-i\sqrt{k_0^2-k_1^2-k_2^2}{r}_0}{e}^{-i\langle (\mathbf{k}-{\mathbf{k}}_r),{\mathbf{x}}_M\rangle }{e}^{-i\langle {\mathbf{k}}_r,{\mathbf{x}}_{\mathcal{D}}\rangle }d(k_1,k_2)$$

Following the experimental setup, the mirror in the reference arm is perpendicular to the incident light, so that the unit normal vector ${\mathbf{\nu }}_M={\mathbf{e}}_3,$ and positioned in the focus of the ${\mathbf{E}}_R^{\left(0\right)}$, such that, following Section 3.2.3, the far-field approximation for reference field ${\mathbf{E}}_R$ is valid. We thus have a reference field ${\mathbf{E}}_R^{\left(1\right)},$ given by

$${\mathbf{E}}_R^{\left(1\right)}\left(r_{\mathcal{D}}{\mathbf{e}}_3\right)=\frac{-i{k}_0}{2\pi r_{\mathcal{D}}}{\mathit{\beta }}_R(0,0){\check{f}}_R(0,0)e^{-i{k}_0(r_0+r_{\mathcal{D}}-2x_{M,3})}.$$

(21)

3.3. OCT Measurements

With the identities of the incident and the backscattered light in hand, we now proceed with the modeling of the detection process inside an OCT system.

In order to suppress common-mode noise and enhance the signal-to-noise ratio, a dual-balance-detector is used to record the OCT signal. After recombination of sample and reference arm light, the laser signal is split 50/50% into two different fibers ${\mathbf{F}}_1,{\mathbf{F}}_2$, each entering one of the DBDs optical inputs. The DBD then subtracts one input from the other, thereby removing everything but the cross-correlation term of the interference.

Thus, assuming that the sample and the reference fields ${\mathbf{E}}_S^{\left(1\right)}$ and $i{\mathbf{E}}_R^{\left(1\right)}$ are passing through a perfect splitter, we obtain the forms for the fields in the fibers as

$${\mathbf{F}}^1=\frac{1}{\sqrt{2}}\left({\mathbf{E}}_S^{\left(1\right)}-{\mathbf{E}}_R^{\left(1\right)}\right),{\mathbf{F}}^2=\frac{1}{\sqrt{2}}\left(i{\mathbf{E}}_S^{\left(1\right)}+i{\mathbf{E}}_R^{\left(1\right)}\right).$$

We assume, ignoring the travel paths inside the fibers, that these fields are detected at the position ${\mathbf{x}}_{\mathcal{D}}$ of the scan lens. These measurements are performed for different wavenumbers $k_0$ in a scan range $[k_{min},k_{max}].$ We therefore indicate explicitly the dependence on $k_0$ in the measurements:

$$\mathcal{M}\left(k_0\right)=\frac{1}{2}\left(|{\mathbf{F}}^1\left(x_{\mathcal{D}}\right){|}^2-{\left|{\mathbf{F}}^2\left(x_{\mathcal{D}}\right)\right|}^2\right)=-\Re e〈{\mathbf{E}}_S^{\left(1\right)},\overline{{\mathbf{E}}_R^{\left(1\right)}}〉,k_0\in [k_{min},k_{max}].$$

(22)

With the identities (20) and (21) for ${\mathbf{E}}_S^{\left(1\right)}$ and ${\mathbf{E}}_R^{\left(1\right)},$ we obtain

$$\mathcal{M}\left(k_0\right)=\frac{k_0}{8r_{\mathcal{D}}{\pi }^3}{\int }_{\mathcal{B}}\Re e\left(-i〈{\mathit{\beta }}_S(k_1,k_2){\check{f}}_S(k_1,k_2),{\mathit{\beta }}_R(0,0){\check{f}}_R(0,0)〉e^{-i\psi (k_1,k_2)}\right)d(k_1,k_2)$$

(23)

where we define the phase function

$$\psi (k_1,k_2)=\langle \mathbf{k}-{\mathbf{k}}_r,{\mathbf{x}}_{\Omega }\rangle +2k_0{x}_{M,3}+\left(\sqrt{k_0^2-k_1^2-k_2^2}-k_0\right)r_0+\langle {\mathbf{k}}_r,{\mathbf{x}}_{\mathcal{D}}\rangle -k_0{r}_{\mathcal{D}}$$

(24)

and use $\mathcal{B}$ as given in Equation (19).

4. Calibration of the Model

So far, we have investigated both the modeling of the backscattered wave from a (layered) object under Gaussian laser illumination and the measurement process of the OCT system described in Section 2.1.

However, simulations based on this explicit model and the following comparison with experimental data presupposes the knowledge of a list of system parameters. Within this list, we distinguish between parameters with values known from specifications such as the wavenumber $k_0$ and parameters which we need to calibrate from the experiment, the beam radius of the Gaussian beam and the angle of acceptance, for example. In order to be capable of extracting these parameters for the simulations, we use two calibration experiments. On the one hand, we consider an experiment which shows the behavior of the backreflected laser power for different surface tilting angles and, on the other hand, we consider the influence of varying positions of the object, with respect to the focus, on the measured data.

4.1. Angular Dependence of the Measured Power

We use a mirror as a sample and analyze the influence of the surface angle on the measured intensity of the scattered electric field. Following the measurement process in Section 2.1 and Section 2.2, the reference arm is blocked, preventing any light from the reference arm to reach the detector. Furthermore, one of the two fibers, which would normally enter the DBD, is connected to a power meter. The measured data are therefore given as the intensity of the scattered field of this mirror. Hereby, again referring to the experimental setup in Section 2.2, we model the totally reflecting mirror as a sample characterized by an infinitely large refractive index. We parameterize the unit normal vector of the mirror surface by

$${\mathbf{\nu }}_{\Omega }=\left(\begin{array}{c}sin{\theta }_{\Omega }\\ 0\\ cos{\theta }_{\Omega }\end{array}\right),$$

(25)

for small values of ${\theta }_{\Omega }\in [\underline{{\theta }_{\Omega }},\overline{{\theta }_{\Omega }}].$

We describe the measurement process for this experiment in a way that the scattered light is detected by a single scan lens point, for simplification we say ${\mathbf{x}}_{\mathcal{D}}=\mathbf{0},$ for a selected wavenumber $k_0$ in the spectrum $[k_{min},k_{max}].$ This in the end, yields a measured intensity of the form

$${\mathcal{M}}_1\left({\theta }_{\Omega }\right)={|\tau {\mathbf{E}}_S^{\left(1\right)}\left(\mathbf{0}\right)|}^2,{\theta }_{\Omega }\in [\underline{{\theta }_{\Omega }},\overline{{\theta }_{\Omega }}],$$

(26)

where ${\mathbf{E}}_S^{\left(1\right)}$ is given by Equation (20) and $\tau \in \mathbb{C}$ accounts for the traveling through the beam splitters. Additionally, we say that the function ${\check{f}}_S$ is approximately given as in Equation (14) with $a=w_0^2/4,$ where $w_0$ represents the radius of the Gaussian beam at the focus. Following the experimental setup, we fix the location $r_0<0$ (below the detector) of the focus and the mirror $x_{\Omega ,3}$ and assume that they are equal: $r_0=x_{\Omega ,3}.$ We follow the notation from Section 3.2.3, but approximate this time the exact form of the domain of integration $\mathcal{B}$ defined in Equation (19), which is an ellipse, by the rectangular domain

$$\mathcal{B}\approx [-L_1\left({\theta }_{\Omega }\right)-k_{0}sin\left(2{\theta }_{\Omega }\right),L_1\left({\theta }_{\Omega }\right)-k_{0}sin\left(2{\theta }_{\Omega }\right)]\times [-L_2,L_2]$$

with the parameters

$$L_1=k_0(-sin(2{\theta }_{\Omega }-{\theta }_{max})+sin\left(2{\theta }_{\Omega }\right)),L_2=k_{0}sin\left({\theta }_{max}\right)$$

and assume that this characterization of $\mathcal{B}$ still allows for an approximation of directions as in Equation (15). Then, using the definitions of ${\psi }_j$ for $j\in \{0,1,2\}$ in Equation (17), we obtain the intensity of the scattered field as a function of ${\theta }_{\Omega },w_0$ and ${\theta }_{max}$

$$\begin{array}{cc}\hfill {\left|\tau {\mathbf{E}}_S^{\left(1\right)}\left(\mathbf{0}\right)\right|}^2=& G({\theta }_{\Omega },w_0,{\theta }_{max}),\hfill \\ \hfill G({\theta }_{\Omega },w_0,{\theta }_{max})=& \frac{|\tau L_1\left({\theta }_{\Omega }\right){|}^2}{8|{\psi }_2{|}^2{\pi }^5}{\left|{\int }_{\mathbb{R}}{e}^{-\frac{{({\psi }_1-\zeta )}^2}{4{\psi }_2}}{e}^{i{k}_{0}sin\left(2{\theta }_{\Omega }\right)\zeta }si(L_1\left({\theta }_{\Omega }\right)\zeta )d\zeta erfi\left(L_2\sqrt{{\psi }_2}\right)\right|}^2,\hfill \end{array}$$

(27)

where $si:\mathbb{R}\to \mathbb{R}$ denotes the unnormalized sinc function given by $si\left(x\right)=\frac{sin\left(x\right)}{x}$ and $erfi:\mathbb{C}\to \mathbb{C}$ is the imaginary error function, defined by $erfi\left(z\right)=\frac{2}{\sqrt{\pi }}{\int }_0^z{e}^{{\zeta }^2}d\zeta .$ Thus, from measurements ${\mathcal{M}}_1\left({\theta }_{\Omega }\right)$ as in Equation (26), corresponding to the data provided by our power meter, for different values ${\theta }_{\Omega }\in [\underline{{\theta }_{\Omega }},\overline{{\theta }_{\Omega }}],$ we can extract the beam radius $w_0$ at the focus and the angle of acceptance ${\theta }_{max}$ as solutions of the minimization problem

$${\displaystyle (w_0,{\theta }_{max})={argmin}_{(z_1,z_2)\in {\mathbb{R}}^2}{\int }_{\underline{{\theta }_{\Omega }}}^{\overline{{\theta }_{\Omega }}}{\left|{\mathcal{M}}_1(\theta )-G(\theta ,z_1,z_2)\right|}^{2}d\theta }$$

with the function G given by Equation (27).

4.2. Reconstructing Sample Information from an OCT Experiment

In the previous section, the beam radius $w_0$ at the focus and the angle of acceptance ${\theta }_{\max }$ have been found. In order to complete the set of parameters necessary for the reconstruction from a measurement at a single point, we additionally need the normal vector ${\mathbf{\nu }}_{\Omega }$ of the tangential plane at each layer boundary.

In swept-source OCT, one in-depth profile of the sample, that is a measurement of the form of $\mathcal{M}$ in Equation (22) (in this case centered at $x_1=x_2=0$), called an A-scan, is acquired during one wavenumber sweep of the laser. To get 3D information, raster scanning in $x_1$ and $x_2$ direction over a certain field of view is performed. The data used in the following are considered as a B-scan, a line of A-scans where only $x_1$ varies at a fixed position $x_2$. Since we assume our layer boundaries to be planes with a certain normal vector ${\mathbf{\nu }}_{\Omega },$ the surface points fulfill an equation of the form $\langle {\mathbf{x}}_{\Omega },{\mathbf{\nu }}_{\Omega }\rangle =c.$ If we can therefore determine at every raster position $x_{\Omega ,1},x_{\Omega ,2}$ the third component $x_{\Omega ,3},$ this determines the normal direction ${\mathbf{\nu }}_{\Omega }.$

Since the single A-scans along those lines are performed independently, we treat these A-scan as single measurements. We shift the coordinate system always so that the incident beam is located at $x_1=x_2=0$ and therefore have ${\mathbf{x}}_{\Omega }=x_{\Omega ,3}{\mathbf{e}}_3.$ We recall that the mirror in the reference arm is modeled as a medium described by an infinitely large refractive index with unit normal vector ${\mathbf{\nu }}_M={\mathbf{e}}_3$ and fixed position at ${\mathbf{x}}_M=x_{M,3}.$

Under these assumptions, we rewrite Equations (23) and (24) as

$${\mathcal{M}}_2\left(k_0\right)=-\frac{k_0}{8r_{\mathcal{D}}{\pi }^3}{\int }_{\mathcal{B}}{\beta }_{\Omega }(k_1,k_2){\check{f}}_S(k_1,k_2){\check{f}}_R(0,0)sin\left(\psi (k_1,k_2)\right)d(k_1,k_2)$$

(28)

with ${\beta }_{\Omega }=\langle {\mathit{\beta }}_S(k_1,k_2),{\mathit{\beta }}_R(0,0)\rangle$ and

$$\psi (k_1,k_2)=(k_3-k_{r,3})x_{\Omega ,3}+2k_0{x}_{M,3}+\left(\sqrt{k_0^2-k_1^2-k_2^2}-k_0\right)r_0+(k_{r,3}-k_0)r_{\mathcal{D}}.$$

For fixed mirror position ${\mathbf{x}}_M$, focus $r_0$ and detector $r_{\mathcal{D}}$, $\psi$ only varies with respect to different depth positions $x_{\Omega ,3}$ of the the sample. Thus, if we can determine the function $\psi$ from the measurements ${\mathcal{M}}_2$ for different A-scans, we also obtain the depth information about the sample.

Under the simplifying assumption that the far-field approximation of the scattered sample field, using $\mathbf{s}={\mathbf{e}}_3$ in Theorem 2, is a reasonable approximation in this case, we rewrite Equation (28) as

$${\mathcal{M}}_2\left(k_0\right)={\left(\frac{k_0}{2r_{\mathcal{D}}\pi }\right)}^2\left|c_3\right|{\beta }_{\Omega }(k_1,k_2){\check{f}}_S(k_1,k_2){\check{f}}_R(0,0)cos\left(k_0\frac{\psi (k_1,k_2)}{k_0}\right),$$

(29)

where the point $(k_1,k_2)$ is defined by Equation (13). Since $\psi$ depends linearly on $k_0$, the measurements are then given as a harmonic oscillation with respect to $k_0$ and with frequency $\psi /k_0.$ To solve for this frequency, we want to Fourier transform with respect to $k_0$, which we define by

$$\mathcal{F}\left(m\right)(\kappa )=\frac{1}{\sqrt{2\pi }}{\int }_{\mathbb{R}}m\left(k_0\right)e^{-i{k}_0\kappa }d{k}_0.$$

However, since we only have band-limited data, we will study the function

$$I:\mathbb{R}\to {\mathbb{R}}_{+},\kappa \mapsto {\left|\frac{1}{\sqrt{2\pi }}{\int }_{k_{\min }}^{k_{\max }}{\mathcal{M}}_2\left(k_0\right)e^{-i{k}_0\kappa }d{k}_0\right|}^2.$$

(30)

We will show in the following that $\psi /k_0$ is determined as the argument where the maximum is located, that is, $\psi /k_0={argmax}_{\kappa }I(\kappa )$. (The absolute value is used to avoid real- or imaginary parts with higher frequent oscillations in order to stably calculate a maximal point.)

We assume that ${\check{f}}_S$ and ${\check{f}}_R$ in Equation (29) are of exponential form as in Equation (14) and define the measurement function

$$M({\Theta }_0;k_0)={\mathcal{M}}_2\left(k_0\right)=K{k}_0^2{e}^{-k_0^2{\sigma }^2}cos\left(k_0{\Theta }_0\right),$$

with the parameters

$${\sigma }^2=\frac{w_0^2}{4}{sin}^2\left(2{\theta }_{\Omega }\right),K=\frac{|c_3|}{{(2\pi r_{\mathcal{D}})}^2}{\beta }_{\Omega }(k_1,k_2),{\Theta }_0=\frac{\psi }{k_0}.$$

By rewriting

$$M({\Theta }_0;k_0)=-K{\partial }_{{\sigma }^2}{e}^{-k_0^2{\sigma }^2}cos\left(k_0{\Theta }_0\right)$$

and interchanging the integral and differentiation in the Fourier transform $\mathcal{F}\left({\partial }_{{\sigma }^2}{e}^{-k_0^2{\sigma }^2}\right),$ we find a form for the Fourier integral in Equation (30) as a convolution

$$\begin{array}{c}\mathcal{F}\left(M\right)({\Theta }_0;\kappa )=\frac{K\delta }{2\pi }{\int }_{\mathbb{R}}-{\partial }_{{\sigma }^2}\frac{1}{\sqrt{2{\sigma }^2}}{e}^{-\frac{1}{4{\sigma }^2}{\zeta }^2}(si(\delta (\kappa -{\Theta }_0-\zeta ))e^{-i\overline{k}(\kappa -{\Theta }_0-\zeta )}\hfill \\ \hfill +si(\delta (\kappa +{\Theta }_0-\zeta ))e^{-i\overline{k}(\kappa +{\Theta }_0-\zeta )})d\zeta ,\end{array}$$

(31)

for $\overline{k}=\frac{k_{\max }+k_{\min }}{2}$ and $\delta =\frac{k_{\max }-k_{\min }}{2}.$ To simplify this expression, we introduce the values

$${\sigma }_{\overline{k}}=\frac{1}{\overline{k}},{\sigma }_{\delta }=\frac{1}{\delta },$$

(32)

and observe from

Table 1. List of parameters: The central wavelength and the wavenumbers $k_{\min }$ and $k_{\max }$ are determined through the specifications of the used swept-source. The beam radius and the angle of acceptance were found through calibration (see Section 5.1). The numerical aperture NA was calculated from the angle of acceptance. The parameter $\sigma$ is given for a typical tilting angle ${\theta }_{\Omega }$ in our experiments.

Parameter	Value	Unit
${\mathsf{\lambda }}_0$ central wavelength	1300	$\mathrm{n}\mathrm{m}$
$w_0$ beam radius at focus	$14.15$	$\mathsf{\mu }\mathrm{m}$
${\theta }_{max}$ angle of acceptance	1.5709	${}^{°}$
$NA$	0.037
$k_{\min }$	$4.7835$	$\mathsf{\mu }{\mathrm{m}}^{-1}$
$k_{\max }$	$4.8973$	$\mathsf{\mu }{\mathrm{m}}^{-1}$
$\overline{k}=\frac{1}{2}(k_{\max }+k_{\min })$	$4.8404$	$\mathsf{\mu }{\mathrm{m}}^{-1}$
$\delta =\frac{1}{2}(k_{\max }-k_{\min })$	$0.056858$	$\mathsf{\mu }{\mathrm{m}}^{-1}$
$\sigma =\frac{w_0}{2}\left|sin\left(2{\theta }_{\Omega }\right)\right|$, ${\sigma }_{\overline{k}}={\overline{k}}^{-1}$, ${\sigma }_{\delta }={\delta }^{-1}$	$0.2753$, $0.2066$, $17.588$	$\mathsf{\mu }\mathrm{m}$

, that ${\sigma }_{\overline{k}}$ and $\sigma$ are of the same order and ${\sigma }_{\delta }$ is considerably larger compared to both of them, meaning that

$${\sigma }_{\overline{k}}=Q\sigma ,\sigma \ll {\sigma }_{\delta },$$

(33)

for some $Q\in \mathbb{R}$ which is close to one. Writing the functions under the integral (31) in terms of these values gives us with

$$\begin{array}{cc}\hfill u_{{\sigma }_{\overline{k}},\sigma }(\zeta )& =\frac{1}{{(\sqrt{2}\sigma )}^3}{e}^{-\frac{1}{Q^2}}\left(1-\frac{{\zeta }^2}{2{\sigma }^2}\right)e^{-{\left(\frac{\zeta }{2\sigma }-i\frac{1}{Q}\right)}^2},\hfill \\ \hfill g_{{\sigma }_{\delta },\pm }(\zeta )& =si\left(\frac{\zeta }{{\sigma }_{\delta }}\right)e^{-i\frac{1}{{\sigma }_{\overline{k}}}(\kappa \pm {\Theta }_0)},\hfill \end{array}$$

(34)

the expression

$$\mathcal{F}\left(M\right)({\Theta }_0;\kappa )=\frac{K\delta }{2\pi }{\int }_{\mathbb{R}}{u}_{{\sigma }_{\overline{k}},\sigma }(\zeta )\left(g_{{\sigma }_{\delta },-}(\kappa -{\Theta }_0-\zeta )+g_{{\sigma }_{\delta },+}(\kappa +{\Theta }_0-\zeta )\right)d\zeta .$$

(35)

Considering Equation (33), we will expand this around $\sigma =0$.

Lemma 1.

Let ${\sigma }_{\overline{k}},{\sigma }_{\delta },\sigma$ as in Equation (32) satisfying Equation (33). Further, let $f_{{\sigma }_{\overline{k}},\sigma }$ be defined as in Equation (34). Then, we have for small values of σ the approximation

$${\int }_{\mathbb{R}}{u}_{{\sigma }_{\overline{k}},\sigma }(\zeta )g_{{\sigma }_{\delta },\pm }(\zeta )d\zeta \simeq \frac{1}{2\sqrt{2}}{e}^{-\frac{1}{Q^2}}\frac{4\sqrt{\pi }}{{\sigma }^2{Q}^2}{g}_{{\sigma }_{\delta },\pm }\left(0\right).$$

(36)

Thus, by applying Lemma 1 to Equation (35), we obtain after changing back to the original system of coordinates

$$\mathcal{F}\left(M\right)({\Theta }_0;\kappa )\simeq \frac{K\delta }{\sqrt{2\pi }}{e}^{-{\overline{k}}^2{\sigma }^2}{\overline{k}}^2\left(si(\delta (\kappa -{\Theta }_0))e^{-i\overline{k}(\kappa +{\Theta }_0)}+si(\delta (\kappa +{\Theta }_0))e^{-i\overline{k}(\kappa -{\Theta }_0)}\right),$$

resulting in

$$\begin{array}{c}{\left|\mathcal{F}\left(M\right)({\Theta }_0;\kappa )\right|}^2\simeq \frac{K^2{\delta }^2}{2\pi }{e}^{-2{\overline{k}}^2{\sigma }^2}{\overline{k}}^4(si{(\delta (\kappa -{\Theta }_0))}^2\hfill \\ \hfill +si{(\delta (\kappa +{\Theta }_0))}^2+2si(\delta (\kappa -{\Theta }_0))si(\delta (\kappa +{\Theta }_0))cos\left(2\overline{k}{\Theta }_0\right)).\end{array}$$

(37)

Note that the dominant sinc terms are centered symmetrically with respect to the origin. In order to derive an explicit expression for the maximum of Equation (37), we want to assume that ${\Theta }_0$ is far away from the origin (which can be accomplished experimentally by tuning the position of the sample), then, these sinc functions do not influence each other strongly. We shift one of them to the origin by setting ${\kappa }^{\prime }=\kappa -{\Theta }_0$ and obtain

$$\begin{array}{c}F({\Theta }_0;{\kappa }^{\prime })\simeq {\left|\mathcal{F}\left(M\right)({\Theta }_0;{\kappa }^{\prime }+{\Theta }_0)\right|}^2\hfill \\ \hfill \simeq si{(\delta {\kappa }^{\prime })}^2+si{(\delta ({\kappa }^{\prime }+2{\Theta }_0))}^2+2si(\delta {\kappa }^{\prime })si(\delta ({\kappa }^{\prime }+2{\Theta }_0))cos\left(2\overline{k}{\Theta }_0\right).\end{array}$$

(38)

Lemma 2.

Let F be defined by Equation (38). Then, for ${\Theta }_0\to \infty ,$ the function F attains a local maximum at ${\kappa }^{\prime }=0.$

Shifting back to the original coordinates and using Lemma 2 yields that Equation (37) attains a maximum at $\kappa ={\Theta }_0,$ that is ${\Theta }_0={argmax}_{\kappa }I(\kappa ),$ which finally gives a representation of Equation (30) as

$$I\left({\Theta }_0\right)\simeq \frac{K^2{\delta }^2}{2\pi }{e}^{-2{\overline{k}}^2{\sigma }^2}{\overline{k}}^4\left(1+si{(2\delta {\Theta }_0)}^2+2si\left(2\delta {\Theta }_0\right)cos\left(2\overline{k}{\Theta }_0\right)\right).$$

(39)

Thus, from the definition of ${\Theta }_0$ we can uniquely determine $\psi .$

We use this information for the reconstruction of the surface angle ${\theta }_{\Omega }.$ For two different, but known lateral positions $x_{\Omega ,1}^j,j\in \{1,2\},$ we consider A-scans leading to measurements ${\mathcal{M}}_2^j$ of the form (28), for different depth positions $x_{\Omega ,3}^j,$ for $j\in \{1,2\},$ respectively. By using the above analysis (under the assumption that the far-field approximation is valid), we determine from the Fourier transform of these two the phase contribution $\psi$ in dependence of $x_{\Omega ,3}^1$ and $x_{\Omega ,3}^2.$ Under the assumption that ${\theta }_{\Omega }$ is considered small, the subtraction of these two then leads to

$$\psi \left(x_{\Omega ,3}^2\right)-\psi \left(x_{\Omega ,3}^1\right)\simeq 2(x_{\Omega ,3}^2-x_{\Omega ,3}^1),$$

which gives the difference in depth $(x_{\Omega ,3}^2-x_{\Omega ,3}^1).$ Together with known lateral information and using that the unit normal vector on the surface satisfies $\langle {\mathbf{\nu }}_{\Omega },({\mathbf{x}}_{\Omega }^2-{\mathbf{x}}_{\Omega }^1)\rangle =0,$ we determine ${\theta }_{\Omega }$ as

$${\theta }_{\Omega }=arctan\left(\frac{x_{\Omega ,1}^2-x_{\Omega ,1}^1}{x_{\Omega ,3}^1-x_{\Omega ,3}^2}\right).$$

5. Results

Finally, we want to validate our model by comparing the simulations with experimental data. We focus on the two previously addressed experiments (see Section 4.1 and Section 4.2, respectively, Section 2.2 and Section 2.3). This quantitative approach shows the dependence of the data on the surface tilting and the focus position.

First we will use the calibration measurement to calculate the beam radius at the focus and the angle of acceptance, then, we look (using the just calibrated parameters) at the data from Section 2.3. The main parameters are presented in Table 1.

5.1. Power vs. Angle—Experiment

Following Section 4.1, we first calibrate the beam radius and the angle of acceptance from the experimental data for different values of the surface angle. This procedure is presented in Algorithm 1. The algorithm is based on the approximated form (27) shortening the computation time by evaluating a one-dimensional integral, instead of the two-dimensional integration presented in Equation (20). As discussed in Section 2.2, the sample arm power arriving at one of the DBD entrances was measured $M=6$ times at $J=25$ angular positions. While the laser power behaves very stable, due to the expected error from the rotational stage and the Gaussian dependence from the angle, some error is observed in the power vs. angle data. Thus, in the following, the data will be plotted with errorbars, representing the standard deviation.

Algorithm 1: Extraction method for the beam radius $w_0$ and the angle ${\theta }_{max}$ of acceptance.

Result: $w_0$ and ${\theta }_{\max }$

Input: wavenumber $k_0=\frac{2\pi }{{\mathsf{\lambda }}_0}$ with the central wavelength ${\mathsf{\lambda }}_0=1300 \mathrm{n}\mathrm{m}$, ${\mathcal{M}}_1\left({\theta }_j\right)$,

for $j=1,\cdots ,J$;

$(w_0,{\theta }_{max})={argmin}_{(z_1,z_2)}\frac{1}{J}{\sum }_{j=1}^J{\left({\mathcal{M}}_1\left({\theta }_j\right)-G({\theta }_j,z_1,z_2)\right)}^2;$

The simulated data ${\mathcal{M}}_1,$ see Equation (26), is given for ${\theta }_{\Omega }={\theta }_j,j=1,\dots ,J$. The match between experimental and simulated data is presented in Figure 4. We remark that both data and simulation follow a Gaussian behavior and attain the maximum at normal incidence, as expected.

5.2. Power vs. Focus—Experiment

Using the calibrated spot size and angle of acceptance from the previous experiment, we compare the simulations with experimental data for multiple B-scans of a mirror and a coverglass as samples of interest.

Following Section 4.2, we first determine the surface angle ${\theta }_{\Omega }$ and adapt the integration area in Equation (19). The coverglass, which is described by a medium with constant refractive index $n_1=1.5088$ and perfectly parallel surfaces, has a thickness d, which is also determined from the experimental data, see Algorithm 2.

The experimental data are measured at a series of different wavelengths ${\mathsf{\lambda }}_j\in (1282.86 \mathrm{n}\mathrm{m},$ $1313.71 \mathrm{n}\mathrm{m})$, $j=1,\cdots ,J,J=700,$ equally spaced in wavenumber $k_{0,j}=\frac{2\pi }{{\mathsf{\lambda }}_j}.$ As described in Section 2.3, the sample was imaged at different positions $x_{\Omega ,3}=x_n,n=1,\cdots ,N,$ along the depth axis.

We ignore polarization effects in the following and use the form (20) with ${\beta }_S=1$ for the simulations of the scattered field of the mirror data. However, for the coverglass experiment, we extend the form to a layer model with two parallel surfaces

$$\begin{array}{c}\hfill {\mathbf{E}}_S^{\left(1\right)}\left(\mathbf{x}\right)=\frac{1}{4{\pi }^2}{\int }_{\mathcal{B}}{\beta }_S(k_1,k_2){\check{f}}_S(k_1,k_2)e^{-i\sqrt{k_0^2-k_1^2-k_2^2}{r}_0}{e}^{-i\langle (\mathbf{k}-{\mathbf{k}}_r),{\mathbf{x}}_{\Omega }\rangle }{e}^{-i\langle {\mathbf{k}}_r,\mathbf{x}\rangle }d(k_1,k_2),\end{array}$$

where we have given the reflection coefficient

$${\beta }_S(k_1,k_2)={\beta }_0-{\beta }_0(1-{\beta }_0^2)e^{-i{k}_{0}2n_{1}dcos{\theta }_t}.$$

We assume again that ${\check{f}}_S$ can be well approximated by Equation (14) and define $\mathcal{B}$ as in Equation (19).

Algorithm 2: Extraction scheme for the tilting angle ${\theta }_{\Omega }$ and the thickness d from the power vs. focus experiment.

Although the dual-balance detection already lowers the noise level in the data, we need to minimize the effects of the residual noise for a meaningful comparison with the simulations.

Thus, for every individual B-scan, we use the mean value of all maximum intensities of it’s A-scans, i.e., we consider the map $\frac{1}{\tilde{N}}{\sum }_{n=1}^{\tilde{N}}{max}_{\kappa \in \mathbb{R}}{I}_n(\kappa ),$ where we use $I_n=I,$ as defined in Equation (30), with $n\in \{1,\cdots ,\tilde{N}\}$ accounting for the number of A-scans in every B-scan and for the different position $x_{\Omega ,3}^n$ (corresponding to this A-scan), as a reference for our simulations. The errorbars in Figure 5 represent the standard deviation of these variations over each B-scan. The highest value of the experimental data is matched to the highest value of the simulation for comparison.

Figure 5 shows, that the Gaussian behavior of the data for the mirror experiment follows the theory.

In Figure 6, we see the comparison between experimental data and simulations for both boundaries of the coverglass. Unfortunately, the calibrated parameters $(w_0,{\theta }_{max})$ from the previous subsection do not yield optimal results, see Figure 6. Similar to Figure 5, a sufficiently strong decrease of the averaged maximum values away from the focus position can be identified in the simulations for the coverglass experiment as well.

At this point, we remark by comparing the experimental datasets, see Figure 7, that the Gaussian curve for the coverglass experiment shows a slightly stretched behavior. This is explained by the measurement of the returning laser light in a diffusive regime originating from the reflection at a slightly rough coverglass boundary surface.

However, updated parameters can be found using an experiment similar to the calibration of $w_0$ and ${\theta }_{max}$, see Figure 8. In contrast to the calibration, the power measurement for the coverglass includes information of both boundaries and, therefore, the measured field intensity is provided as a sum

$${\left|{\mathbf{E}}_S^{\left(1\right)}\right|}^2={\left|{\mathbf{E}}_{S,1}^{\left(1\right)}+{\mathbf{E}}_{S,2}^{\left(1\right)}\right|}^2,$$

where we consider first order reflections only. Due to additional scattering events inside the coverglass material, the background information ${\mathbf{E}}_{S,2}^{\left(1\right)}$ is smaller than ${\mathbf{E}}_{S,1}^{\left(1\right)}$ and therefore neglected for the calibration.

A comparison for the coverglass experiment—after updating the system parameters—shows that almost all simulated data points lie inside the estimated range for both boundaries, see Figure 9.

For a better visualization, in Figure 10, a direct comparison between experimental data and the simulations (using recalibrated system parameters) for a single A-scan is provided. In order to obtain a reasonable comparison, both experimental and simulated datasets are normalized to one.

6. Discussion

We have considered samples with a very distinct scattering profile, but still a slight difference in the angular reflectivity profile could be observed for the mirror and the coverglass. For diffusely scattering samples, the proposed description can easily be generalized, especially with the mentioned automatic angular power measurement. An automated measurement would also reduce the error, so that a fast angle scanning procedure could be implemented prior to OCT imaging.

The error in the power vs. focus experiment (see Figure 5 and Figure 6) is caused by power variations inside single B-scans. These are explained by a combination of reasons. The tilt of the sample with respect to the OCT illumination generates a slight continuous change in distance to the focus for each A-scan inside a B-scan. In addition, the scan lens induces a certain curvature of field, resulting in a change of illumination depending on the raster scanning position $x_1$.

Although the presented results in this work show suitable correspondence with the provided OCT data, we note that the algorithm for solving the minimal-error-solution problem in extracting the beam radius and the angle of acceptance from the power-angle experiment suffers from the fact that the function in Equation (27) is highly oscillating and it is therefore difficult to find an “optimal” set of parameters.

Nevertheless, for a rather large range of values of these parameters, we observe a reasonable match with the experimental data, at least, for simple, layered examples. The model should, however, work nicely also for more complicated samples with different geometries (that is, with more and potentially curved layers).

We expect that this model can be used as a forward model, which simulates the OCT data from the knowledge of the refractive index of the medium and can be used for the corresponding inverse problem of reconstructing the refractive index in the layers.

7. Conclusions

We presented a method to model the image formation in OCT based on a real-life 1300 $\mathrm{n}$$\mathrm{m}$ swept-source setup. In contrast to publications based on plane wave models for the OCT system, the proposed model includes the effect of additional system relevant parameters such as the focus and the beam radius of the incident laser light and the angle of acceptance. We also suggested a way how to determine these (not necessarily a priori known) parameters, either from the OCT data themselves and from calibration measurements.

A comparison between simulation and experiment shows that the presented model together with the derived system and sample parameters produces a quantitatively correct prediction of the OCT data. We therefore expect that this model can serve as a forward model for an inversion algorithm to quantitatively reconstruct the optical material properties of a sample from OCT data, in particular its refractive index.