Development of mathematical models for quantitative OCT: A review

1.   Introduction

Optical Coherence Tomography (OCT) is a non-invasive imaging modality using non-ionizing radiation, typically light in the near infrared spectrum. Even though the first OCT imaging system was built only thirty years ago ^[23], it has attracted much interest and it is considered nowadays as a key imaging technique, especially in ophthalmology, since it can acquire images very quickly (around $50\, 000$ pixel points per second) with a high resolution of a few micrometers and a tissue penetration depth of approximately two millimetres. The OCT setup raises many challenges for a rigorous mathematical model, in particular for a quantitative analysis where we want to explicitly express the measurements in terms of the physical parameters of the sample under investigation. In addition, a concrete model needs to be implemented in three dimensions, which makes the derived algorithms computationally expensive.

In this work, we aim to give a brief overview of the different approaches to model the various OCT setups. Hereby, there are many modeling questions to be answered before formulating an appropriate model: How is the sample illuminated? Is there light absorption in the medium? Does the polarization of the light play an important role or can it be neglected? Do we take multiple scattering effects into account, or is it enough to focus on single scattering? What kind of detectors are being used? Just to name a few relevant decisions to be made.

The basis of our analysis is the mathematical description presented in the review paper ^[14] by one of the inventors of OCT, Adolf Fercher, together with ^[9] where a general formulation in this spirit was presented. We will thus describe the light as an electromagnetic wave and will be using Maxwell's macroscopic equations as our model for the light propagation, which reduces in this setting to a Helmholtz equation for the electric field. We will then start by giving an overview of some of the most prominent OCT setups and the corresponding modeling assumptions in this setting in Section 2, and show examplarily some of the resulting forward models for OCT in Section 3. In Section 4, we will present the effects of the different assumptions and simplifications to the forward simulations and compare the results to real measurement data.

2.   Different modeling choices for OCT

There exist many different variations in the experimental setup of OCT, but its basic working principle is to illuminate the sample from a fixed direction (we will use the third basis vector $e_3$ as the illumination direction) with an incident light beam $E_I$ and to measure the backreflected light $E_S$ by interferometry, that is, to superimpose $E_S$ with a copy $E_R$ of the incident field (typically produced by reflecting $E_I$ at a mirror) and to detect the resulting intensity.

The light propagation is thereby commonly described as an electromagnetic wave, which is a pair of an electric and a magnetic field evolving according to Maxwell's macroscopic equations. Since all the media the light is interacting with are typically assumed to be non-magnetic, it is sufficient to keep track of the electric field $\hat E\colon \mathbb{R}\times \mathbb{R}^3\to \mathbb{R}^3$ and the electric displacement field $\hat D\colon \mathbb{R}\times \mathbb{R}^3\to \mathbb{R}^3$, which are functions of time and space; although it is usually more convenient to do the calculations in Fourier domain, that is, for the functions

where the constant $c$ denotes the speed of light in vacuum.

2.1.   The illumination

The incident light is a solution of Maxwell's equations in the vacuum and can be seen as the electric field which would be produced in the absence of the sample. The equations can be in Fourier space reduced to the Helmholtz equation

for the electric field $E_I$ of the incident light with the additional condition that

where the differential operators $\nabla$ and $\Delta$ only act on the spatial variable $x$.

Swept Source OCT vs. Broadband Illumination: In swept source OCT, the incident wave is almost monochromatic, meaning that we have for every $x\in \mathbb{R}^3$ that

for a sufficiently small $\varepsilon > 0$, and the measurement is repeated for multiple values $k_0$ in some interval $[k_{{\rm{min}}}, k_{{\rm{max}}}]$, $0 < k_{{\rm{min}}} < k_{{\rm{max}}}$. (We remark that the function $E_I$ in Fourier space fulfills the symmetry $E_I(-k, x) = \overline{E_I(k, x)}$, since each component of the electric field $\hat E_I$ is real-valued.)

Alternatively, the illumination can be done with one illumination with a broadband laser beam with

Since we will see that we can solve the equations for the resulting electric field $E(k, \cdot)$, in particular the vector Helmholtz equation (2.8), for each value of $k$ separately, this difference is irrelevant for the mathematical model and we formulate it in the following therefore simply for a single broadband illumination.

Gaussian Beams: In the classical OCT setup, the illumination is done by a laser beam (propagating along the negative $e_3$ direction) which is focused on a certain spot $r\in \mathbb{R}^3$, which is usually chosen inside the sample. To record the full image, measurements are done for illuminations at every horizontal position $(r_1, r_2)\in \mathbb{R}^2$ (practically, on some finite grid). The depth $r_3$ of the focus is thereby normally kept constant. Thus, the incident light is typically well described by a Gaussian shaped beam, see, for example, ^[49] for a treatise of the generation of light with a laser. To obtain an expression for such a beam, we may solve the equations (2.1) and (2.2) for an arbitrary value $k$ with initial data of the form

in the focus plane $\{x\in \mathbb{R}^3\mid x_3 = r_3\}$, where $f_k\colon \mathbb{R}^2\to \mathbb{C}$ is a (Gaussian like) function such that its two-dimensional Fourier transform $\hat f_k$ is supported in the disk $D_{|k|}(0)\subset \mathbb{R}^2$ with center zero and radius $|k|$ and $\eta_k\in \mathbb{R}^2\times\{0\}$ denotes a fixed polarization vector. This has a solution which is symmetric with respect to $x_3$ regarding the focus plane of the form

with

and the adapted polarization vectors

However, since the incident wave should be propagating downwards, so that before a certain time, say $t = 0$, we can ensure that $\hat E_I(t, r+x) = 0$ for all points $x$ below a certain horizontal plane where the sample is located, say in the half space $\{x\in \mathbb{R}^3\mid x_3 < 0\}$; we know for every $x\in \mathbb{R}^3$ with $x_3 < 0$ that the function $k\mapsto E_I(k, r+x) = \int_0^\infty\hat E_I(t, r+x)e^{i kct}dt$ must be holomorphically extendible to the upper half complex plane such that $\omega\mapsto\int_ \mathbb{R} |E_I(k+i\omega, r+x)|^2dk$ is bounded, the function $E_I$ should be only a combination of the terms $e^{i k\phi\cdot x}$ with $\phi\in \mathbb S^2$ and $\phi_3 < 0$ (corresponding to the plane waves $(t, x)\mapsto e^{i k(\phi\cdot x-ct)}$ moving in the direction $\phi$).

Thus, we could keep as the initial wave the downwards moving part $E_I = E_{I, -}$ of our symmetric solution, where we assume $f_k = 0$ for $|k|\notin[k_{{\rm{min}}}, k_{{\rm{max}}}]$ and impose the reality assumptions $f_{-k} = \overline{f_k}$ and $\eta_{-k} = \eta_k$. If we consider the typical case that $\hat f_k$ has a very small support around the origin, where we want to assume the Gaussian behavior

for some parameters $A_k\in \mathbb{C}$ and $a_k > 0$, we can keep only the lowest order in $\kappa$ in the integrand and find the classical form

for a Gaussian beam, depicted in Figure 1.

Plane Waves: Approximating the expression (2.6) close to the focus point, that is, for $x\approx 0$, we see that the shape of the Gaussian beam resembles the one of the plane wave

Thus, if the object is placed at the beam waist and its size is (for every $k\in[k_{{\rm{min}}}, k_{{\rm{max}}}]$) of the same order of magnitude as the Rayleigh length $2a_kk$ of the light beam, then we can approximate the illumination by a plane wave, see ^[14]. Since this assumption simplifies the mathematical analysis a lot, most of the works that followed used this form of illumination, see, for example, ^{[5,6,7,15,16,24,29,53]} and some papers from the authors of this article ^[9,10,11,13].

However, in particular, intensity changes due to the distance to the focus point are not reflected in this approximation which can have a negative impact for quantitative reconstructions.

Besides this rastering of the sample with focused illuminations, there also exist full-field OCT methods where the whole sample is illuminated at once. For these setups, the illumination is typically well modeled by plane waves, see ^[33,34], for example.

2.2.   The medium

In Maxwell's macroscopic equations, the medium is usually modeled as a non-magnetic, dielectric medium, whose optical properties are characterized by its electric permittivity $\epsilon\colon \mathbb{R}\times \mathbb{R}^3\to \mathbb{C}^{3\times 3}$, a function of the wavenumber $k$ and the position $x$, which relates the electric displacement field $D$ and the electric field $E$ in the medium via the linear relationship

In the end, we are interested in the inverse problem of reconstructing $\epsilon$ from the OCT data; however, this data seems not sufficient for a full recovery of the nine four-dimensional, complex-valued entries of $\epsilon$. It is therefore common to restrict the analysis to media with a certain structure where a reconstruction becomes possible.

There are various popular assumptions, which can in different combinations be suitable.

Isotropic Medium: In general, the inner structure of the dielectric medium can cause that the displacement field $D$ produced in the presence of an electric field $E$ is not parallel to $E$ so that the permittivity $\epsilon$ is indeed matrix-valued. To recover in this case some of the matrix entries, there exists the so-called polarization sensitive OCT method which makes illuminations with differently polarized incident waves and also detects the polarization state of the reflected light, see ^[22]. At least for orthotropic media (that is, media whose permittivity $\epsilon$ fulfill that $\epsilon(k, x)$ is symmetric and $\epsilon_{13} = 0 = \epsilon_{23}$), a functional OCT scheme has been addressed in ^[11].

Mostly, however, it is assumed that the medium is isotropic, meaning that $\epsilon$ is a multiple of the identity matrix:

with the refractive index $n\colon \mathbb{R}\times \mathbb{R}^3\to\{z\in \mathbb{C}\mid\Im z\ge0\}$.

Non-Absorbing Medium: For biological samples in the near-infrared spectrum (that is, with central wavelengths between $800$ and $1300$ nanometers), where OCT typically operates, we have that the scattering effects dominate the absorption effects. Thus, it is often a reasonable assumption to neglect absorption completely, which means to assume that $\epsilon(k, x)\in \mathbb{R}^{3\times3}$ for all $(k, x)\in \mathbb{R}\times \mathbb{R}^3$.

Non-Dispersive Medium: In addition, the bandwidth of the light source in a typical OCT system is within a narrow spectrum of a few tens of nanometers. Thus, it is commonly assumed that the permittivity is independent of the wavenumber $k$ over that small bandwidth of the source, see ^[6,7,37], for example. However, this way dispersion effects are excluded from the model.

Slowly Varying in Horizontal Direction: Since the laser illuminates only a small region inside the sample, we can assume that the properties of the medium within the illuminated region are almost independent of the lateral coordinates $x_1$ and $x_2$, so that the permittivity can be locally considered to depend spatially only on the depth $x_3$.

Weakly Scattering Medium: An assumption with mainly theoretical motivation, is to pretend that the medium is only weakly scattering, meaning that $\epsilon- \mathbb I_{3\times3}$ is so small that we can apply the Born approximation to linearize the forward problem and get an explicit formula for it, see ^[9,14,45,54]. Allowing for strongly scattering media in which also multiple scattering events needs to be considered, results in contrast in more complicated models ^[1,27,36,58].

Layered Medium: One of the main applications of OCT is in ophthalmology where the retina of the human eye is imaged. Another important application is the imaging of human skin. Both of them are biological tissues presenting a layered structure, where each layer is described by a distinct permittivity (usually with additional simplifying assumption such as being constant or only weakly scattering) so that we have disjoint bounded sets $\Omega_j\subset \mathbb{R}^3$, $j\in\{1, \ldots, J\}$, with a smooth boundary and local permittivities $\epsilon_j\colon \mathbb{R}\times \mathbb{R}^3\to \mathbb{C}^{3\times3}$, $j\in\{1, \ldots, J\}$, for which the permittivity has the form

where $\Omega = \bigcup_{j = 1}^J\Omega_j$ is the region of the sample and $J\in \mathbb{N}$ denotes the number of layers, see, for example, ^{[9,18,19,26,52]}.

Often one makes the additional assumption that, at least locally in the illuminated region, the domains $\Omega_j$ can be approximated by the sets

which are bounded by planes with upwards pointing normal unit vectors $\nu_{\Omega, j}\in \mathbb S^2$ and parameters $z_j\in \mathbb{R}$ (all depending on the position of the illuminated region).

Point-Like Scatterers: Similarly to the layered media, one may assume the inhomogeneities inside the object as randomly distributed point-like scatterers with constant optical properties. Then, the aim is to recover the geometry of the area where the scatterers are located ^{[35,38,46,60]}.

2.3.   The scattering

In general, Maxwell's macroscopic equations describe perfectly the wave propagation into the imaging system and the scattering theory is the main tool to derive a formula for the backscattered light from the object, see, for example, ^[8,32]. Since most OCT systems are using focused illuminations and are thus essentially performing multiple one-dimensional measurements, the electromagnetic theory can be considerably simplified, such as in ^{[3,9,15,16,44,48,51]}, to name but a few examples.

In the general case, Maxwell's macroscopic equations can be (implementing the medium via the relation (2.7)) combined into the vector Helmholtz equation

for the electric field $E\colon \mathbb{R}\times \mathbb{R}^3\to \mathbb{C}^3$, which is induced by the incident field $E_I$, that is, we request an initial condition of the form $\hat E(t, x) = \hat E_I(t, x)$ for all times $t < t_0$ with the time $t_0$ be chosen such that $\hat E_I$ is not interacting with the sample before $t_0$: $\hat E_I(t, x) = 0$ for all $t < t_0$ and $x\in\bigcup_{k\in \mathbb{R}}{\rm{supp}}\; (\epsilon(k, \cdot))$. Choosing for simplicity $t_0 = 0$, this can be again formulated as a boundedness condition on the analytic continuation of the inverse Fourier transform $E_S$ of the scattered fields $\hat E_S = \hat E-\hat E_I$ on the upper half complex plane (as indicated in the derivation of (2.6)), or more commonly expressed by the Silver–Müller radiation condition

see ^[8], for example.

The differential equation (2.8) together with this initial condition can then be rewritten as the Lippmann-Schwinger integral equation

with the integral operator $\mathcal U$ given by

where $\mathbb I$ is the identity operator and

denotes the three-dimensional outgoing Green's function of the Helmholtz operator (that is, the solution of the distributional equation $\Delta\Phi+k^2\Phi = -\delta$ fulfilling the boundedness condition for $k$ in the upper half complex plane). We note that the integration area in the operator $\mathcal U$ is restricted to the support of the functions $\epsilon(k, \cdot)- \mathbb I_{3\times3}$, which represents the sample and should thus be bounded.

Born Approximation: The integral (2.9) is only an implicit equation for the scattered field. If we make the assumption of a weakly scattering medium, meaning that we consider the deviation $\epsilon(k, y)- \mathbb I_{3\times3}$ of the permittivity from the permittivity of free space to be small, we may neglect in a first order approximation the second order term $(\epsilon(k, y)- \mathbb I_{3\times3})(E(k, y)-E_I(k, y))$ in the integrand and approximate $E_S$ by its first order Born approximation, given by the explicit expression

The product of $(\epsilon(k, y)- \mathbb I_{3\times3})E_I$ at a point $y$ can be interpreted as the response of the medium there to the (unperturbed) incident illumination $E_I$, which acts as a point source that propagates according to the Green's function to an examination point $x$. The integral over the object volume then provides the total contribution of the sample.

However, the approximation does not reflect that these additional contributions also interact with the medium. This could be iteratively taken into account by using higher order Born approximations defined by

But although the usage of higher-order Born approximation in case of a strongly deviating refractive index is possible, it is due to the evaluation of many integrals not very efficient. Therefore, other possibilities such as the implementation of the radiative transfer equation and Monte Carlo simulations are often preferred, see ^{[20,27,28,36,50,58]}.

Far-Field Approximation: Maybe less an approximation, but rather attributed to the comparably large distance of the detector to the sample, it is sometimes considered that the scattered field there is only known in highest order so that barely the first term $E_S^{(\infty)}(k, \tfrac x{|x|})$ in the expansion

for $|x|\to\infty$ significantly enters into the measurements.

In particular in combination with the Born approximation, this nevertheless simplifies the expression for the scattered field as we end up with

which makes it a very common model, see, for example, ^{[9,14,15,16,17,42]}, especially, combined with plane wave illumination. When it comes to Gaussian beam illumination, the difference between far- and near-field representation becomes more noticeable. The latter counts also for system effects that are neglected when the far-field approximation is applied, see ^[5,13,59].

Linear Polarization: If we assume an isotropic medium described by a refractive index $n$, then we can write the vector Helmholtz equation (2.8) (using the vector identity $\nabla\times(\nabla\times E) = \nabla(\nabla\cdot E)-\Delta E$ together with the direct consequence $\nabla\cdot(n^2 E) = 0$ of (2.8)) in the form

from which we see that if $E\cdot\nabla n$ vanishes identically, then the vector Helmholtz equation is reduced to the scalar Helmholtz equation for every component of $E$. This is, for example, the case if $E_I(k, x) = u_I(k, x)\eta$ with $\eta\in \mathbb{R}^2\times\{0\}$ and the refractive index $n$ is assumed to depend only on $x_3$ (which may be approximated by assuming a medium which is slowly varying in horizontal direction) so that $\nabla n$ is parallel to $e_3$. Then, the electric field is of the form $E(k, x) = u(k, x)\eta$ with $u$ being a solution of

thus reducing the vector-valued problem to a scalar one, see, for example, ^{[3,14,15,16,34]}.

Fresnel Formulas and the Huygens–Fresnel Principle: We consider for simplicity a layered medium whose layers are bounded by parallel planes perpendicular to an upwards pointing normal vector $\nu_\Omega\in \mathbb S^2$ where the local permittivities are assumed to be isotropic, non-dispersive and spatially constant: $\epsilon_j(k, x) = n_j^2 \mathbb I_{3\times3}$. Then, the electric field $E$ in each layer $\Omega_j$ will be given as a superposition of plane waves of the form

with different propagation directions $\phi\in \mathbb S^2$ for appropriately chosen amplitudes $A_{j, \phi}\in \mathbb{C}$ and polarization vectors $\eta_{j, \phi}\in \mathbb S^2$ with $\phi\cdot\eta_{j, \phi} = 0$.

The amplitudes and polarizations have then to be chosen so that at every point $x_{\Omega, j}\in\partial\Omega_j\cap\partial\Omega_{j+1}$ on a discontinuity of the refractive index, the transmission conditions

are fulfilled for all sequences $(\bar y_m)_{m = 1}^\infty\subset E_j$ and $(\underline{y}_m)_{m = 1}^\infty\subset E_{j+1}$ with $\lim_{m\to\infty}\bar y_m = x_{\Omega, j}$ and $\lim_{m\to\infty}\underline{y}_m = x_{\Omega, j}$.

For only two (infinite) layers and an incident plane wave, the solution can be shown to be of the form

where the amplitudes $\tilde A, B\in \mathbb{C}$, the polarization vectors $\tilde\eta, \bar\eta\in \mathbb S^2$, and the propagation direction $\tilde\phi\in \mathbb S^2$ of the reflected and transmitted waves can be explicitly calculated via the Fresnel equations from the amplitude $A\in \mathbb{C}$, the polarization $\eta\in \mathbb S^2$ and the direction $\phi\in \mathbb S^2$ of the incident plane wave, see, for example, ^[25].

Now, similarly to the higher Born approximation, we can iteratively calculate the solution of the general scattering problem by—starting with the incident wave, decomposed as a sum of plane waves—successively reflecting and scattering the different plane waves at the discontinuities of the layers. Since each reflection step goes along with a loss in amplitude, one can in practice only keep those plane wave components which were produced from a low number of reflections. For a more detailed explanation of the construction, we refer to ^[7,13,43].

We want to mention at this point that even if we take in this construction only contributions from waves into account which have been reflected at most once, the so obtained single scattering solution would be different from the corresponding Born solution, since it would still depend non-linearly on the refractive index (as it appears in the exponent of every plane wave component).

We refer to ^[6,7,21,47] for works using similar modeling. The layered structure can be also combined with different differential equations, like the radiative transfer equation ^[20,55].

Instead of calculating the solution inside the layers completely, there also exists the approach of applying the extended Huygens-Fresnel principle, see ^[4,31,61], to propagate the values of the electric field from one layer boundary directly to the next layer boundary and in that way also back to the detector, see ^[3,44,51].

2.4.   The measurements

Since the time resolution is usually not sufficient to determine from measurements of the rapidly oscillating electromagnetic waves $\hat E_S$ the Fourier transform $E_S$, the way to get information about the phase is to split off a part of the incident wave via a beam splitter and reflect it at a mirror, which produces a well-known electromagnetic field $E_R$ (which can be calculated from the Fresnel equations), and recombine it with the scattered field $E_S$ to measure the superposition $E_S+E_R$. There are slight variations on how this is exactly done.

Time-Domain vs. Spectral-Domain OCT: In time-domain OCT, the mirror to produce the field $\hat E_R$ is movable and for every shifted position $\mu\in[-\mu_0, \mu_0]$ of the mirror, we get for an induced scattered field $\hat E_S$ at a detector point $x\in \mathbb{R}^3$ the averaged intensity measurement

The dominant first term is simply the intensity of the incident light beam, which is known, and the last term can be obtained by an intensity measurement of the scattered field without superimposing it with the reference field $\hat E_R$ (but often it is simply ignored as it is of second order in the field $\hat E_S$, which is usually much weaker than $\hat E_R$). Therefore, we can assume to know by this the second, cross-correlation term.

Using the Parseval–Plancherel identity, we can write it in the form

so that we can calculate from this by doing a Fourier transform with respect to $\mu$ (provided the interval of $\mu$ is sufficiently large—in general, a convolution operator will additionally appear, see ^[16]) the product

The expression $\Re(M)$ can alternatively be directly measured by sending the sum $\hat E_S+\hat E_R$ through a spectrometer giving the values $|E_S(k, x)+E_R(k, x)|^2$ of which the values $|E_S(k, x)|^2$ and $|E_R(k, x)|^2$ can be again obtained separately leaving only $\Re(M)$. This is the way spectral-domain OCT is working, which has the experimental advantage, that there is no more necessity of moving the mirror which speeds up the measurements.

We want to remark that if not a broadband signal is used, but a swept source method, then for each illumination with an almost monochromatic incident wave at wavenumber $k$, only one value $\Re(M(k, x))$ is recovered with this spectral method, and the full signal is obtained by repeating the procedure for every value $k\in[k_{{\rm{min}}}, k_{{\rm{max}}}]$.

Polarization Sensitive OCT: For a polarization sensitive setup, every measurement is repeated with an incident wave where just the polarization vectors in every Fourier component are changed. Additionally, the superposition $E_S+E_R$ is splitted and sent to different polarization filters so that we can measure the single components $\Re(E_{S, j}(k, x)\overline{E_{R, j}(k, x)})$, $j\in\{1, 2, 3\}$, of the cross-correlations.

Focused Illumination vs. Full-field OCT: In full-field OCT, the whole sample is, as we have seen, illuminated at once and we measure the resulting scattered field (in one of the previously described ways) at all points $x$ on a detector plane $\mathcal D = \{x\in \mathbb{R}^3\mid x_3 = d\}$ above the sample.

In the classical setup, we generate for each point $(r_1, r_2, d)\in\mathcal D$ an incident wave $(k, x)\mapsto E_I(k, r+x)$ with focus in $r = (r_1, r_2, r_3)\in \mathbb{R}^3$ for some depth focus $r_3$ (inside the sample) and detect the hereby induced scattered field $E_S$ at the detector position $(r_1, r_2, d)$ only. Hereby, we produce at every point on the detector a measurement which contains information of the sample in a small neighbourhood of the vertical line $\{(r_1, r_2, x_3)\mid x_3\in \mathbb{R}\}$, which is typically called an "A-scan".

The reconstruction of the medium is thus essentially reduced to a (one-dimensional) reconstruction of the sample on each of these vertical lines from its corresponding A-scan. Thus, the A-scans are treated independently and the image quality reduces far from the focal plane. Interferometric synthetic aperture microscopy (ISAM) on the other hand uses a more sophisticated model and takes into account scattering and diffractions effects in order to produce images with resolution independent of the focus position. ISAM may be implemented by an existing OCT system with the addition of a spectral detector. We refer the interested reader to ^{[34,38,40,41]} and references therein.

3.   Examples of forward models

To give a brief overview, what assumptions are typically made and how they are usually combined, we tried to collect and classify some of the theoretical works on OCT in Table 1.

Illumination

p.w. = plane wave, G.b. = Gaussian beam

Medium

w.s. = weakly scattering, (n.-)abs. = (non-)absorbing, (n.-)disp. = (non-)dispersive, lay. = layered medium, p.sca. = point scatterers

Scattering

B.A. = Born approximation, F.A. = far-field approximation, Huy. = Huygens principle l.p. = linear polarization, FF. = Fresnel formulas

Although we cannot present all the works in detail here, we would like to give a brief summary of the main ideas. The first set of articles (see the first row of Table 1) assumes plane wave illumination for a weakly scattering medium. The works of Fercher and colleagues ^{[14,15,16,17]} are the ones that every researcher interested in the fundamentals of OCT should start with. The basics of reflection tomography and backscattering are presented but also how an OCT system is implemented. All different types of OCT are addressed and how other factors (like noise) affect the signal. In ^[9] we see one of the first attempts by mathematicians to describe how OCT operates. In addition, iterative methods are presented for reconstructing the refractive index from OCT data even for anisotropic media. The case of a medium consisting of discrete reflectors is addressed in ^[24] and the formula of the measured intensity is derived for both TD- and FD-OCT. Here, signal processing techniques required in practice (like phase shifting) are also analysed. In ^[26] a model is presented that deals with the non-linear sampling of the measured signal both theoretically and numerically. The novelty of the work ^[34] has to do with the implementation of a method that produces images with spatially independent resolution (see subsection 2.4 where we discuss about ISAM).

The second set of papers refer to mathematical models based on plane wave illumination and assume layered structure for the medium. In ^[6,7] the authors use an one-dimensional model to reconstruct the refractive indices of a multi-layered medium. The work ^[12] can be seen as an extension to more dimensions with the addition of point-like inhomogeneities inside the layers. A similar inverse problem to reconstruct both the thickness and the refractive index of every layer is studied in ^[37]. The application of Fresnel's equations was discussed in ^[52,53]. A theoretical approach for various OCT systems is presented in ^[52] whereas in ^[53] a matrix formulation is used and the solution of a minimization problem produces accurate reconstructions.

The next set of papers consider the scattering of a Gaussian beam by a medium in a FD-OCT system. The works ^[38,39] (just to name a few) present the advantages of ISAM for a medium consisting of multiple points scatterers. The derived solution of the inverse problem is independent of the focus position. In ^[57] the authors of this review article collaborated with researchers from the Medical University of Vienna to validate the accuracy of the presented mathematical model for layered media with experimental data.

The last row in Table 1 refers to works where the analytical model describing OCT is based on the extended Huygens-Fresnel (EHF) principle and how multiple scattering events affect the signal. The main advantages is that it is independent of the sample arm geometry and it can be combined with Monte Carlo simulations to extract the scattering properties of the sample. We refer only to the early works since there have been tens of later articles using this model.

In the following, we consider two particular cases and how these assumptions are combined to a complete forward model.

3.1.   Single scattering plane wave model

One of the classical models, used for example in ^[14], assumes (locally) a plane wave illumination and considers a medium which is isotropic, non-absorbing, and non-dispersive. Then, we have seen that we can describe the medium by a refractive index $n\colon \mathbb{R}^3\to \mathbb{C}$, being a function of $x$.

Assuming further that the medium is weakly scattering and applying the Born and far-field approximations, we obtain for the far-field pattern $E_S^{(1, \infty)}$ of the scattered field $E_S$ induced by an incident plane wave $E_I(k, x) = Ae^{-ikx_3}\eta$ with $A\in \mathbb{C}$ and $\eta\in \mathbb{R}^2\times\{0\}$ in the direction $\vartheta\in \mathbb S^2$ with $\vartheta_3 > 0$ the expression

As expected, the measurements thus are proportional to the Fourier transform $\mathcal F[n^2-1](k(\vartheta+e_3))$ of the reflectivity $n^2-1$. We stress here once more that OCT is a band-pass technique and we actually detect only high frequency Fourier components because only backscattered light is measured within a small frequency range.

In Figure 2 we see the few accessible Fourier data (in red) for two different OCT systems compared to the available Fourier data (gray upper cone): In standard OCT with focused illumination, we will only detect the directly back-scattered values, that is, $\vartheta = e_3$ (at each horizontal position), whereas in full-field OCT we would (ideally) obtain the contributions from all $\vartheta\in \mathbb S^2$ with $\vartheta_3 > 0$.

This result is the main reason for the few mathematical-oriented papers, since for many years the inverse problem in OCT, to recover the reflectivity $n^2-1$ from the OCT data, was just accomplished by the application of the Fourier transform. However, the need to improve image quality and enhance quantitative analysis forced different mathematical groups to relax some of the assumptions and consider different cases to handle the more general scheme.

3.2.   A Gaussian illumination model

In ^[57], we used a Gaussian incident wave $E_I$ for analyzing a layered medium with constant local refractive indices and whose layers $\{x\in \mathbb{R}^3\mid z_{j+1} < x\cdot\nu_\Omega\le z_j\}$ are bounded by parallel planes with normal vector $\nu_\Omega$.

Using (2.3) for the data in the focal plane, this leads (we only consider the A-scan at the position $(r_1, r_2) = (0, 0)$ to simplify the notation) to

where the following notation was introduced:

● $\beta_j\colon \mathbb{R}\times \mathbb{R}^2\to \mathbb{R}^3$ denotes the dispersive reflection coefficient, accounting for the differently polarized parts of the electric field;

● for a $\kappa\in \mathbb{R}^2$ with $K = (\kappa, -\sqrt{k^2-|\kappa|^2})\in \mathbb{R}^3,$ we define the scattered vector

● $\theta_{t, j}$ denotes the angle of transmission from the layer $j$ to the layer $(j+1)$;

● $x_\Omega = (0, 0, z_1)$ denotes the point on the top surface of the medium.

A part which often has been neglected in articles on OCT, is the transport of the reflected light to the detector. In practice, a common setup is that the back-reflected light is coupled into a single-mode fiber and later transferred to the detector. Since both fibers connecting the sample and the reference arm with the detector are assumed to have equal length, the additional phase factor, which commonly represents the way of light through a fiber, cancels out when being combined in (2.10).

While the way through the fiber only takes a minor role in an OCT measurement, the coupling into the fiber has a more prominent position in here ^[38,57]. This coupling is done by using a scan lens which discards all (plane wave) parts in equation (3.1) whose propagation directions strongly deviate from the illumination direction $e_3$ ^[30]. The maximal deviation is described by the (maximal) angle of acceptance, which we denote by $\theta.$ The set of incident directions yielding an accepted wave vector is then given by

Thus, we change the area of integration in (3.1) to $\mathcal B_k\subset \mathbb{R}^2$, which gives us (after multiplying with the reflected field) our forward model.

To simplify the expression, the integral may be replaced by its far-field approximation: For a direction $\vartheta\in \mathbb S^2_+$ and (large) $\rho\in \mathbb{R},$ the far-field approximation of the reflected field by a single interface, that is $j = 1$ in (3.1), is given by ^[57]

where $\kappa_0$ is such that $\psi(K(\kappa_0)) = \vartheta$ and $C$ is constant depending on $\kappa_0.$

It has been shown in ^[57] that the far-field approximation is meaningful only if the medium is located in the focus of the laser beam. However, there are still cases where this approximation fails as we demonstrate in Section 4.1 with a numerical example.

Although we have for this model no closed form for the reconstruction of the positions $z_j$ of the layer boundaries and the refractive indices $n_j$ of the layers, the parameters can numerically be calculated to optimally fit the data which gives for suitable samples satisfying the modeling assumptions good reconstructions.

4.   Numerical experiments

In this section we present numerical experiments using both simulated and real OCT data. The goal is to show that some mathematical models even if they are commonly used might not be appropriate for some specific OCT setups and some objects.

4.1.   Near- vs. Far-Field representation

We address a special case where near- and far-field representations of the sample field coincide at a specific point but deviate in general.

Let the imaging object be a totally reflecting mirror that is not tilted, meaning $\beta\approx 1$ and $\nu_\Omega = e_3$. We assume further that it is placed below the examination point $x = 0.$ Then, after linearizing the root in (3.1), we obtain the following form for the reflected field

where we used the definition of $\hat f_k$ in equation (2.5). It is straightforward to show that $\mathcal B_k,$ in this simple case, is given by the ball with radius $k \sin(\theta).$ Thus, we parametrize the integral to get

For $z_1 = r_3,$ we observe that the absolute value of (3.3) is only a good approximation of $|E_S(k, 0)|$ if $\theta$ increases to the point that almost all back-reflected waves are collected. The difference between a small (left) and a large (right) value of $\theta$ is shown in the bottom row of Figure 3.

4.2.   Plane wave vs. Gaussian shaped illumination

We present an example with real data generated from a swept-source OCT system. We show that the Gaussian beam seems advantageous. In addition, we observe that even multiple reflections (included in the model) do not improve the simulations. We compare the real data with the single-reflection model for the plane and the Gaussian shaped illumination. Simulations with second-order reflections are also presented.

The OCT experiment was performed by a fiber-based OCT system built at the Medical University of Vienna. It uses a swept-source that emits 60 mW at a central frequency of around $1300$nm and bandwidth from approximatively $1282.86$nm to $1313.76$nm. A fiber coupler splits the light into the reference ($25\%$) and sample arm ($75\%$). A dual-balance-detector measures the cross-correlation term and it is digitized by a data acquisition card. For more details on the OCT system see ^[57].

The medium was a three-layered sample consisting of glass, water and glass. The sample was raster-scanned $1024 \times 1024$ times in lateral directions. For each raster scan we record a spectrum consisting of $1498$ data points, equally spaced with respect to wavenumber. Due to the small bandwidth of the laser, dispersion compensation is not needed. In order to adjust all necessary parameters for the Gaussian incident field the object must be shifted multiple times along its axial direction. For each shift, the measurement process is repeated.

The optical properties of the sample are known. The used coverglass and the water have, at the wavelength $1300$ nm, refractive indices $1.5088$ and $1.3225$, respectively. The thickness of each layer was calculated manually using the optical distance in the image space.

In Figure 4 we present a comparison between the Fourier transform (with respect to wavenumber) of a single raster scan and the simulations for a plane wave and a Gaussian beam illumination.

Regarding the major peaks, the maxima of the sinc functions, we observe that the Gaussian model performs clearly better. This point is crucial for the inverse problem where the quality of the reconstruction of the refractive index is highly dependent on the height difference between real data and forward simulations ^[7,13]. However, the plane wave model seems to catch better the minor peaks appearing far from the centres of the sinc functions and this is reflected in Table 2 where we see that this model leads to smaller root-mean-square error (RMSE).

That this lack of information cannot be compensated by adding multiple reflections in the simulations is presented in Figure 5, which shows a comparison, again in Fourier domain, between the real data and the plane wave model with first and multiple reflections. The multiple reflections were calculated up to second order. Hereby, we consider only the multiple reflections that travel at most the time needed for light to travel from the front to the back surface. The effect of the multiple reflections is minor as expected from the analytic representation. This is also true for higher order reflections.

4.2.1.   Focusing effect

The effect showing responsible for this inconsistency between the two presented simulations in Figure 4 is the focus of the laser beam. It is a strongly affecting factor that does not appear in the plane wave illumination but it is dominant in the Gaussian beam case. Thus, it is essential to consider this effect when we want to solve the inverse problem of reconstructing the material parameters from the OCT data.

The difference between the plane wave and the Gaussian beam model is even more highlighted when the object of interest is moved out of the focus. Figure 6 shows the absolute value of the Fourier transform of the cross-correlation term. Here, one can easily identify a strong decrease of the intensity as we move further from the focus position. For the simulations, we assumed that the laser light is reflected by a non-tilted totally reflecting mirror. The mirror is shifted to different positions $x^j_\Omega, \, j = 1, \dots, 9,$ with respect to the focus position $r_3.$ For each position the cross-correlation term (2.10) is recorded on the exact same spectrum as in Section 4.2.

The lacking intensity decrease in the plane wave model as the object moves out of the focus point can be artificially compensated for by multiplying with an axial point spread function (PSF) which imitates the exponential behavior of the actual Gaussian beam with respect to the focus position ^[56]. For a perfectly reflecting object, the axial PSF takes the form

where $d$ accounts for the distance between the object and the focus position and $z_r$ is the Rayleigh length, see [56,equation (4)].

In order to show the similarity of the two models, we compare the intensity in the sample arm as a function of the distance between the object and the focus. In Figure 7, we present the simulations for both models at $41$ points. The Gaussian model is plotted for multiple values of the maximum angle of acceptance. Due to the strong focus of the OCT setup, which is accompanied by a small beam width in the center, the Gaussian model converges, as the angle of acceptance increases, to a state of saturation. This happens for angular values larger than $4^\circ$ and then the Gaussian model coincides with the PSF correction.

5.   Conclusions

In this work we presented different possibilities for modeling an Optical Coherence Tomography system. The performed analysis was based on different mathematical models originated from Maxwell's equations or their simplification, the Helmholtz equation. We tried to cover as many as possible cases considering various illumination techniques, medium properties, scattering behaviours and measurement ways. The addressed models were also tested/compared using simulated and real OCT data.

Acknowledgments

We would like to thank Wolfgang Drexler and Lisa Krainz both from the Medical University of Vienna for providing us with the experimental data presented in Section 4.2. P. Elbau and L. Veselka acknowledge support from the Austrian Science Fund (FWF) within the project F6804–N36 of the Special Research Programme SFB F68 "Tomography Across the Scales". L. Mindrinos was supported by the Christos Papakyriakopoulos bequest.

Conflict of interest

All authors declare no conflicts of interest in this paper.