Quantitative OCT reconstructions for dispersive media

TL;DR

This paper develops a method for reconstructing the optical properties of dispersive media using OCT data, addressing phase retrieval and parameter identification in layered, inhomogeneous media.

Contribution

It introduces an iterative algorithm for frequency-dependent parameter reconstruction in dispersive media, including phase retrieval and absorption considerations.

Findings

Effective reconstruction of dispersive media properties from OCT data.

Algorithm handles phase retrieval and absorption effects.

Applicable to layered inhomogeneous media.

Abstract

We consider the problem of reconstructing the position and the time-dependent optical properties of a linear dispersive medium from OCT measurements. The medium is multi-layered described by a piece-wise inhomogeneous refractive index. The measurement data are from a frequency-domain OCT system and we address also the phase retrieval problem. The parameter identification problem can be formulated as an one-dimensional inverse problem. Initially, we deal with a non-dispersive medium and we derive an iterative scheme that is the core of the algorithm for the frequency-dependent parameter. The case of absorbing medium is also addressed.

Tables

Table 1. Table 1: Reconstructed values using Algorithm 1 for a three-layer medium.

length (mm)	${\ell }_0$	${\ell }_1$	${\ell }_2$	${\ell }_3$
exact	0.50000	0.20000	0.30000	0.10000
reconstructed (noise free)	0.50000	0.20004	0.29990	0.10006
reconstructed ($5\%$ noise)	0.49960	0.20021	0.30030	0.09976
refractive index	$n_1$	$n_2$	$n_3$	$n_4$
exact	1.55000	1.41000	1.48000	1.00000
reconstructed (noise free)	1.55107	1.41070	1.48087	1.00014
reconstructed ($5\%$ noise)	1.55272	1.40740	1.48170	0.99851

Table 2. Table 2: Reconstructed values using Algorithm 2 for a three-layer medium.

length (mm)	${\ell }_0$	${\ell }_1$	${\ell }_2$	${\ell }_3$
exact	0.70000	0.15000	0.40000	0.13000
reconstructed (noise free)	0.69885	0.15210	0.39628	0.13451
reconstructed ($5\%$ noise)	0.70340	0.15387	0.39337	0.13669
refractive index	$n_1$	$n_2$	$n_3$	$n_4$
exact	1.55000	1.40500	1.55000	1.00000
reconstructed (noise free)	1.55107	1.40568	1.55107	0.99782
reconstructed ($5\%$ noise)	1.55164	1.40599	1.55139	0.99695

Table 3. Table 3: Reconstructed values using Algorithm 3 for a three-layer medium. The values of the refractive indices are given at ω ∗ = 2.351 × 10 15 . superscript 𝜔 ∗ 2.351 superscript 10 15 \omega^{\ast}=2.351\times 10^{15}.

length (mm)	${\ell }_0$	${\ell }_1$	${\ell }_2$	${\ell }_3$
exact	0.70000	0.20000	0.30000	0.10000
reconstructed ($2\%$ noise)	0.69790	0.20139	0.29569	0.10329
refractive index	$n_1\left({\omega }^{\ast }\right)$	$n_2\left({\omega }^{\ast }\right)$	$n_3\left({\omega }^{\ast }\right)$	$n_4$
exact	1.54488	1.41000	1.49675	1.00000
reconstructed ($2\%$ noise)	1.54499	1.41094	1.49884	1.00484

Equations

$$\nabla \times E(t,\mathrm{x})+\frac{1}{c}\frac{\partial H}{\partial t}(t,\mathrm{x})=0,\nabla \times H(t,\mathrm{x})-\frac{1}{c}\frac{\partial D}{\partial t}(t,\mathrm{x})=0,$$

$$D(t,\mathrm{x})=E(t,\mathrm{x})+4\pi {\int }_{\text{mathbbm}R}\chi (\tau ,\mathrm{x})E(t-\tau ,\mathrm{x})d\tau .$$

$$E(t,\mathrm{x})=E_0(t,\mathrm{x}),H(t,\mathrm{x})=H_0(t,\mathrm{x}),t<0,\mathrm{x}\in \text{mathbbm}{R}^3.$$

$${\int }_{\text{mathbbm}R}|(E-E_0)(t,\mathrm{x})+(E_r-E_0)(t,\mathrm{x}){|}^{2}dt,r\in \text{mathbbm}R,\mathrm{x}\in \mathcal{D}.$$

$$(\hat{E}-{\hat{E}}_0)(\omega ,\mathrm{x}),\omega \in \text{mathbbm}R,\mathrm{x}\in \mathcal{D},$$

$$\mathcal{F}(f)(\omega )={\int }_{\text{mathbbm}R}f(t)e^{i\omega t}dt.$$

$$\hat{m}(\omega ,\mathrm{x})=|(\hat{E}-{\hat{E}}_0)(\omega ,\mathrm{x})+({\hat{E}}_r-{\hat{E}}_0)(\omega ,\mathrm{x}){|}^2,\omega \in \text{mathbbm}R,\mathrm{x}\in \mathcal{D}.$$

$${\hat{m}}_s(\omega ,\mathrm{x})=|(\hat{E}-{\hat{E}}_0)(\omega ,\mathrm{x}){|}^2,\omega \in \text{mathbbm}R,\mathrm{x}\in \mathcal{D},$$

$$E(t,\mathrm{x})=\left(\begin{array}{c}0\\ u(t,z)\\ 0\end{array}\right),H(t,\mathrm{x})=\left(\begin{array}{c}v(t,z)\\ 0\\ w(t,z)\end{array}\right).$$

$$\Delta u(t,z)-\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}{\int }_{\text{mathbbm}R}ϵ(\tau ,z)u(t-\tau ,z)d\tau =0,$$

$$u(t,z)=u_0(t,z),t<0,z\in \text{mathbbm}R.$$

$$n(t,z)=\{\begin{array}{ll}{n}_0,& z\in \text{mathbbm}R\setminus \overline{L},\\ n_j(t),& z\in L_j.\end{array}$$

$$n(z)=\{\begin{array}{ll}{n}_0,& z\in \text{mathbbm}R\setminus \overline{L},\\ n_j,& z\in L_j,\end{array}$$

$${\partial }_{tt}u(t,z)=\frac{c^2}{n^2(z)}\Delta u(t,z),t\in \text{mathbbm}R,z\in \text{mathbbm}R,$$

$$u_0(t,z)=f_0(z-c_{0}t),$$

$$n(z)=\{\begin{array}{ll}{n}_0,& z<z_1,\\ n_1,& z>z_1.\end{array}$$

$${\partial }_{tt}u(t,z)$$

$$u(t,z)$$

$$u(t,z)=\{\begin{array}{ll}{u}_{-}(t,z)=f_0(z-c_{0}t)+g_0(z+c_{0}t),& t\ge 0,z<z_1,\\ u_{+}(t,z)=f_1(z-c_{1}t),& t\ge 0,z>z_1.\end{array}$$

$$\underset{z\uparrow z_1}{\lim }{u}_{-}(t,z)=\underset{z\downarrow z_1}{\lim }{u}_{+}(t,z),\underset{z\uparrow z_1}{\lim }{\partial }_z{u}_{-}(t,z)=\underset{z\downarrow z_1}{\lim }{\partial }_z{u}_{+}(t,z),$$

$$R:f_0\mapsto g_0,R[f_0](z+c_{0}t)=\frac{c_1-c_0}{c_1+c_0}{f}_0\left(2z_1-(z+c_{0}t)\right),$$

$$T:f_0\mapsto f_1,T[f_0](z-c_{1}t)=\frac{2c_1}{c_1+c_0}{f}_0\left(z_1+\frac{c_0}{c_1}((z-c_{1}t)-z_1)\right).$$

$${\partial }_{tt}u(t,z)$$

$$u(t,z)$$

$$u(t,z)=\{\begin{array}{ll}{u}_{-}(t,z)=g_0(z+c_{0}t),& t\ge 0,z<z_1,\\ u_{+}(t,z)=f_1(z-c_{1}t)+g_1(z+c_{1}t),& t\ge 0,z>z_1.\end{array}$$

$$R_{-}:g_1\mapsto f_1,R_{-}[g_1](z-c_{1}t)=\frac{c_0-c_1}{c_1+c_0}{g}_1\left(2z_1-(z-c_{1}t)\right),$$

$$T_{-}:g_1\mapsto g_0,T_{-}[g_1](z+c_{0}t)=\frac{2c_0}{c_1+c_0}{g}_1\left(z_1+\frac{c_1}{c_0}((z+c_{0}t)-z_1)\right).$$

$$U_3[f_1](t,z)=f_1(z-c_{1}t)+g_1(z+c_{1}t).$$

$$n(z)=\{\begin{array}{ll}{n}_0,& z<z_1,\\ n_1,& z\in (z_1,z_2),\\ n_2,& z>z_2.\end{array}$$

$$u(t,z)$$

$$+\text{mathds}{1}_{(y,\infty )}\sum _{j=0}^{\infty }{U}_3\left[(R_{-}{R}_{+}{)}^{j}T{f}_0\right](t,z),t\ge 0,y\in (z_1,z_2),$$

$$U_1[f](t,z)$$

Full text

1Faculty of Mathematics

University of Vienna

Oskar-Morgenstern-Platz 1

A-1090 Vienna, Austria

2Johann Radon Institute for Computational

 and Applied Mathematics (RICAM)

Altenbergerstrasse 69

A-4040 Linz, Austria

1. Introduction

Optical Coherence Tomography (OCT) is nowadays considered as a well-established imaging modality producing high-resolution images of biological tissues. Since it first appeared in the beginning of 1990s [11, 17, 30], OCT has gained increasing acceptance because of its non-invasive nature and the use of non-harmful radiation. Main applications remain tissue diagnostics and ophthalmology. It operates at the visible and near-infrared spectrum and the measurements consist mainly of the backscattered light from the sample. OCT is analogous to Ultrasound Tomography where acoustic waves are used and differs from Computed Tomography (where electromagnetic waves are also used) because of its limited penetration depth (few millimeters) due to the lower energy radiation. As OCT data we consider the measured intensity of the backscattered light at some detector area usually far from the medium.

However, the intensity of light, undergoing few scattering events, is not measured directly, but the OCT setup is based on low coherence interferometry. The incoming broadband and continuous wave light passes through a beam-splitter and it is split into two identical beams. One part travels in a reference path and is totally back-reflected by a mirror and the second part is incident on the sample. The backscattered from the sample and the back-reflected light are recombined and their superposition is then measured at a detector. The maximum observed intensity refers to constructive interference, and this happens when the two beams travel equal lengths. For a detailed explanation of the experimental setup we refer to [10, 34] and to the book [4].

The way the measurements are performed characterizes and differentiates an OCT system. We summarize here the different setups considered in this work:

Time-domain OCT:

The reference mirror is moving and for each position a measurement is performed. By scanning the reference arm, different depth information from the sample is obtained.

Frequency-domain OCT:

The mirror is placed at a fixed position and the detector is replaced by a spectrometer, which captures the whole spectrum of the interference pattern.

State-of-the-art OCT:

The incoming light is focused, through objective lenses, to a specific region at a certain depth in the sample. The backscattered light is measured at a point detector.

Standard OCT:

The vector nature of light is ignored and the electromagnetic wave is treated as a scalar quantity. Then, only the total intensity is measured.

Time- and Frequency-domain OCT provide almost equivalent measurements that are connected through a Fourier transform. The advantage of the later is that no mechanical movement of the mirror is required, improving the acquisition time. The last two cases simplify the following mathematical analysis since we can consider scalar quantities and depth-dependent optical parameters. For an overview of the different mathematical models that can be used in OCT we refer to the book chapter [6].

We consider Maxwell’s equations to model the light propagation in the sample, which is assumed to be a linear isotropic dielectric medium. We deal with dispersive and non-dispersive media. Firstly, using a general representation for the initial illumination, we present the direct problem of computing the OCT data, given the optical properties of the sample. Then, we derive reconstruction methods for solving the inverse problem of recovering the refractive index, real or complex valued. Motivated by the layer stripping algorithms [28, 31], we present a layer-by-layer reconstruction method that alternates between time and frequency domain and holds for dispersive media.

Without loss of generality, the OCT system can be simplified by placing the beam-splitter and the detector at the same position. The medium is contained in a bounded domain $\Omega\subset\mathbbm{R}^{3},$ such that $\mbox{supp}\chi(t,\cdot)\subset\Omega,$ for all $t\in\mathbbm{R},$ where $\chi$ is the electric susceptibility, a scalar quantity describing the optical properties of a linear dielectric medium. We set $\chi$ to zero for negative times. Also, the medium for $t\leq 0$ is assumed to be in a stationary state with zero stationary fields. Then, the electric field $E\in C^{\infty}(\mathbbm{R}\times\mathbbm{R}^{3};\mathbbm{R}^{3})$ and the magnetic field $H\in C^{\infty}(\mathbbm{R}\times\mathbbm{R}^{3};\mathbbm{R}^{3}),$ in the absence of charges and currents, satisfy the Maxwell’s equations

[TABLE]

where $c$ is the speed of light and $D$ is the electric displacement, given by

[TABLE]

This relation models a linear dielectric, dispersive medium with inhomogeneous, isotropic and non-stationary parameter.

The two identical laser pulses, one incident on the sample and the other on the mirror, are described initially, before the time $t=0,$ as vacuum solutions of the Maxwell’s equations, meaning (1) for $D\equiv E,$ defined by $E_{0},\,H_{0}\in C^{\infty}(\mathbbm{R}\times\mathbbm{R}^{3};\mathbbm{R}^{3}).$ In practice, the medium is illuminated by a Gaussian light, however at the scale of the sample the laser pulse can be approximated by a linearly polarized plane wave [9]. We assume that the incident wave does not interact with the medium until $t=0,$ resulting to the condition

[TABLE]

The mirror is modeled as a medium with (infinitely) large constant electric susceptibility, with surface given by the hyperplane placed at $r\in\mathbbm{R}$ distance from the source. Given the form of the incident wave, the reference field (back-reflected field), denoted by $E_{r},$ can be explicitly calculated.

The sample wave (backscattered wave) is given as a solution of the system (1) – (3). Then, the two backward traveling waves are recombined at the beam splitter, assumed to be at the detector position. In time-domain OCT, the sum of these two fields, integrated over all times, is measured at each point of the two-dimensional detector array $\mathcal{D}\subset\mathbbm{R}^{2}$. Thus, as observed quantity we consider

[TABLE]

Under some assumptions on the incident field [6], we may recover from the above measurements, the quantity

[TABLE]

where $\hat{f}=\mathcal{F}(f)$ denotes the Fourier transform of $f$ with respect to time

[TABLE]

In frequency-domain OCT, the detecting scheme is different. The mirror is not moving ($r$ is fixed) and the detector is replaced by a spectrometer. Then, the intensity of the sum of the Fourier transformed fields at every available frequency (corresponding to different pixels at the CCD camera) is measured

[TABLE]

In practice, we obtain data only for few frequencies restricted by the limited bandwidth of the spectrometer. The OCT system allows also for measurements of the intensities of the two fields independently, by blocking one arm at a time. Thus, we assume that the quantity

[TABLE]

is also available. The main difference between the two setups is that (5) provide us with the full information of the backscattered field, amplitude and phase, which is not the case in (7), where we get phase-less data. We address later the problem of phase retrieval, meaning how to obtain (5) from (7).

Up to now, what we have modeled is known as full-field OCT where the whole sample is illuminated by an extended field. The main problem is that we want to reconstruct a $(1+3)$-dimensional function $\chi$ from OCT data, either (4) or (6), which are $(1+2)$-dimensional. Thus, we have to impose additional assumptions in order to compensate for the lack of dimension. To solve this problem, we consider a medium which admits a multi-layer structure. This assumption is not far from reality since OCT is mainly used in ophthalmology (imaging the retina) and human skin imaging. In both cases the imaging object consists of multiple layers with varying properties and thicknesses [14, 20].

If the medium is non-dispersive, meaning that the optical parameter is stationary, the function $\chi$ can be modeled as a $\delta-$distribution in time, so that its Fourier transform (temporal) does not depend on frequency. Then, even if we have enough information (theoretically), in OCT, as in any tomographic imaging technique, we deal with the problem of inverting partial and limited-angle data. This is the result of measuring only the back-scattered light for a limited frequency spectrum. In OCT, a narrow beam is used, resulting to an almost monochromatic illumination centered around a frequency.

In the following, we focus on data provided from a state-of-the-art and standard OCT system, where point-like illumination is used. In this case, only a small region inside the object is illuminated so that the function $\chi$ can be assumed depth-dependent and constant in the other two directions. Again we assume that locally the illumination is still properly described by a plane wave.

Let $\mathrm{x}=(x,y,z),$ where the $z-$direction denotes the depth direction. We model the light as a transverse electric polarized electromagnetic wave of the form

[TABLE]

Then, the Maxwell’s equations (1) together with (2) are simplified to

[TABLE]

for the scalar valued function $u,$ where $\Delta=\partial^{2}/\partial z^{2}.$ Here, we define the time-dependent electric permittivity $\epsilon(t,z)=\delta(t)+4\pi\chi(t,z),$ which varies also with respect to depth. The condition (3) is replaced by

[TABLE]

The medium admits a multi-layered structure with $N$ layers orthogonal to the $z-$direction, having spatial-independent but time-dependent refractive indices $n=\sqrt{\epsilon},$ and varying lengths. We define $L=\cup_{j=1}^{N}L_{j}$ and we set

[TABLE]

This setup is commonly used for modeling the problem of parameter identification from OCT data. Even the volumetric OCT data consist of multiple A-scans, which is a one-dimensional cross-sectional of the medium across the $z-$direction. Under the assumption of a layered medium, the multiple A-scans are averaged over the $x-$ and $y-$directions producing a profile of the measured intensity with respect to frequency or depth (post-processed image).

In Fig. 1, we see the experimental data for a three-layer medium with total length $0.7$mm, having two layers (top and bottom) filled with Noa61 $(n_{1}=n_{3}\approx 1.55)$ and a middle one filled with DragonSkin $(n_{2}\approx 1.405)$. The spectrometer uses a grating with central wavelength $840$nm, going from $700$ to $960$nm. On top right, we see an A-scan of the $``$raw$"$ data (depth information), meaning the intensity of the combined sample and reference fields at a given point on the surface plane. The left picture is the post-processed B-scan (two-dimensional cross-sectional of the volumetric data) where we see clearly the four reflections from the boundaries. The bottom right picture presents the averaged (over lateral dimension) post-processed version of the data on the left. All $x-$axes are in pixel units.

We refer to [2, 25, 32, 33] for recent works using similar setup and assumptions. Our work differs from previous methods in that we consider dispersive medium. We deal also with absorbing medium, a property that is usually neglected. We address three different cases for the layered medium:

• •

$n_{j}(t)\equiv n_{j},\quad j=1,...,N$ (non-dispersive),

• •

$\hat{n}_{j}(\omega)\in\mathbbm{R},\quad j=1,...,N$ (dispersive),

• •

$\hat{n}_{j}(\omega)\in\mathbbm{C},\quad j=1,...,N$ (dispersive with absorption).

The paper is organized as follows. In Sec. 2 we present the forward problem, meaning given the medium (location and properties) find the measurement data. We derive formulas that are also needed for the corresponding inverse problem, which we address in Sec. 3. Iterative schemes are presented for dispersive media and a mathematical model is given for the case of absorbing media. In Sec. 4 we give numerical results for simulated data, and we show that the parameter identification problem can be solved under few assumptions.

2. The forward problem

We derive mathematical models for the direct problem in OCT for multi-layer media with piece-wise inhomogeneous refractive index. We start with a single-layer medium and then we generalize to more layers. The multiple reflections are also taken into account. Most of the formulas presented in this section, like the solutions of the initial value problems or the reflection and transmission operators (analogue to the Fresnel equations) can be found in classic books on partial differential equations [8, 29] and optics [1, 3, 13], respectively. However, we summarize them here, on one hand because we want to derive a rigorous mathematical model in both time and frequency domains and on the other hand because they are needed for the corresponding inverse problems. The easier but essential time-independent case is treated first. Then, we consider the time-dependent case by moving to the frequency-domain for real and complex valued parameters.

2.1. Non-dispersive medium

Here, we simplify (10), and we consider the following form for the refractive index

[TABLE]

for $j=1,...,N.$ We describe the light propagation using (8) together with (9). Under the above assumption, we obtain

[TABLE]

the one-dimensional wave equation. In the following, we use $c_{j}=c/n_{j},\,j=0,...,N.$ Let us assume that the initial field is given by the form

[TABLE]

together with the assumption that $\mbox{supp}f_{0}\subset(-\infty,z_{1})$, where $z_{1}$ represents the surface (first boundary point) of the medium $L$. This assumption on the support of the function reflects the condition that the laser beam does not interact with the probe until time $t=0$.

We model the single-layer medium as $L=(z_{1},\,z_{2}),$ for $z_{1}<z_{2},$ but initially we consider the case

[TABLE]

Then, we obtain the system

[TABLE]

The above system of equations describes a wave traveling from the left incident on the interface at $z_{1}\in\mathbbm{R}$, see the left picture of Fig. 2. It is easy to derive the solution, which is given by

[TABLE]

Here, the function $g_{0}$ and $f_{1}$ describes the reflected and transmitted field, respectively. Given the continuity condition at $z=z_{1}$, meaning

[TABLE]

we find a representation of $g_{0}$ and $f_{1}$ via operators. We denote the reflection operator by $R$ and the transmission operator by $T$, defined by

[TABLE]

and

[TABLE]

The fact that $\mbox{supp}f_{0}\subset(-\infty,z_{1})$ implies that for every $t<0$, $R[f_{0}]=0,$ in $(-\infty,z_{1}),$ and $T[f_{0}]=0,$ in $(z_{1},\infty)$. This is true, since neither the reflected nor the transmitted wave exists before the interaction of the initial wave with the boundary. Finally, we define the operator $U_{1}:f_{0}\mapsto u,$ mapping the initial function $f_{0}$ to the solution $u,$ given by (16), of the problem (15).

Now we consider the following problem

[TABLE]

for an initial wave $g_{1}$ with $\mbox{supp}g_{1}\subset(z_{1},\infty),$ for $n$ as in (14).

This problem refers to the case of a wave incident from the right on the boundary $z=z_{1}$, see the right picture in Fig. 2. Again, we obtain a reflected and a transmitted part of the wave. The solution of this problem is given by

[TABLE]

As previously, we find a representation of the reflected and transmitted waves using operators acting on the initial wave. Here, we denote the reflection operator by $R_{-}$ and the transmission operator by $T_{-},$ having the forms

[TABLE]

and

[TABLE]

We define the solution operator $U_{2}:g_{1}\mapsto u$, mapping the initial $g_{1}$ to the solution of the problem (18).

It is trivial to model an operator $U_{3}:f_{1}\mapsto u,$ where $f_{1}$ satisfies $\operatorname{supp}f_{1}\subset(-\infty,z_{2}),$ for an interface at $z=z_{2},$ with $n(z)=n_{2},$ for $z>z_{2}.$ This setup models how a transmitted, from the boundary at $z=z_{1},$ wave propagates for $t\geq 0.$ We know that on $(-\infty,z_{2})$, $U_{3}[f_{1}]$ is of the form

[TABLE]

We define in addition the operator $R_{+}:f_{1}\mapsto g_{1}$. These three different cases are combined to produce the following result.

Proposition \theproposition.

Let $L=(z_{1},\,z_{2})$ be a single-layer medium, and let the refractive index be given by

[TABLE]

If the initial wave $f_{0},$ see (13), satisfies $\operatorname{supp}f_{0}\subset(-\infty,z_{1}),$ then the solution of (12), together with $u(t<0,z)=f_{0},$ is given by

[TABLE]

with $U_{1},\,U_{2},$ and $U_{3}$ defined as before.

Proof**:**

The function $u,$ given by (19), is by construction a solution to the wave equation problem in both $(-\infty,y)$ and $(y,\infty)$. Thus, we have only to check if both parts coincide in the interval $(z_{1},z_{2})$. To do so, we recall the definitions

[TABLE]

for $t\geq 0,$ and $z\in(z_{1},\,z_{2}).$ By plugging these formulas in (19) we get for $(-\infty,y)$ the term

[TABLE]

and for $(y,\infty)$ the term

[TABLE]

Under the assumption that both series are convergent, we see that the two terms coincide. The last thing to show is that (19) also satisfies the initial condition $u(t,z)=f_{0}(z-c_{0}t),$ for all $t<0$ and $z\in\mathbbm{R}$. This can be seen by considering the supports of the operators $U_{1},\,U_{2}$ and $U_{3},$ as they are defined previously.

The formula (19), consists of three terms and accounts also for the multiple reflections occurring in the single-layer medium. Each term neglects either the boundary at $z_{1},$ or the one at $z_{2}$. The operator $U_{1}$ maps the initial wave $f_{0}$ to the solution $u$, in the real line. Then, the transmitted wave $T[f_{0}]$ is traveling back and forth between $z_{1}$ and $z_{2}$, describing the multiple reflections, given by the field

[TABLE]

The operator $U_{2}$ now uses for every $j\in\mathbbm{N}$ the reflection of (20) at $z_{2}$ as initial and gives back a solution of the sub-problem (18). The last term models the wave interacting with the boundary at $z=z_{2}$, by application of $U_{3},$ which uses (20) as an initial function.

In the following example, we present the forms of the single- and double-reflected wave from the boundary at $z=z_{2},$ measured in $(-\infty,z_{1}).$

Example \theexample.

We know already that the reflected wave from the boundary $z=z_{1},$ is given by (17). We present now the reflected waves $\varrho_{\eta_{r}},\,\eta_{r}=1,2$ in the interval $(-\infty,z_{1}),$ where $\eta_{r}$ counts for the numbers of the undergoing reflections, meaning

[TABLE]

For $\eta_{r}=1,$ using the definition of the operator $T$ applied to $f_{0}$ we get

[TABLE]

The argument of $f_{0}$ is now a function of $z-c_{1}t$, and we can apply $R_{+},$ resulting to

[TABLE]

Thus, we have

[TABLE]

a function of $z+c_{1}t,$ where the operator $T_{-}$ can be to applied to give

[TABLE]

Following the same procedure for $\eta_{r}=2$ now, we end up with the following form

[TABLE]

Now, we move to the case of a multi-layer medium. The solution will be derived using the above formulas and consider the problem layer-by-layer. We define $L_{j}=(z_{j},\,z_{j+1}),$ for $j=1,...,N.$ The refractive index is given by (11) and we define

[TABLE]

The parameter $\tilde{n}_{+}$ represents the refractive index in the remaining $N-1$ layers.

Next we want to find a solution of (12) for

[TABLE]

and $y\in(z_{1},\,z_{2}).$ The first case is exactly the same as the problem already discussed for the single-layer case, meaning that the application of the operators $U_{1}$ and $U_{2}$ is still valid. For the later case, we consider an initial wave $f$ with $\mbox{supp}f\subset(-\infty,z_{2})$ and by $U_{3}[f]$ we denote the solution of this sub-problem. We know that in $(-\infty,z_{2})$, $U_{3}[f]$ gives

[TABLE]

We define an operator $R_{+}$ with $R_{+}[f]=g,$ corresponding to the multi-reflected light from the boundaries $z_{2},\dots,z_{N+1}$, if we illuminate illumination by $f$. Then, we get the following result.

Proposition \theproposition.

Let $n$ be defined by (11), and $n_{-}$, $n_{+}$ as in (21). Then the solution of (12), together with $u(t<0,z)=f_{0},$ is given by

[TABLE]

We know that the solution of (12) in $(-\infty,z_{1})$ admits the form

[TABLE]

and we define an operator $\tilde{R},$ through $\tilde{R}[f_{0}]=g$.

Proposition \theproposition.

Let the operators $R,R_{-}$ and $T$ be defined as previously and $\tilde{R}[f_{0}]$ be given. Then, the following holds

[TABLE]

where

[TABLE]

Proof**:**

From Proposition 2.1 we know that

[TABLE]

describes a solution. It also holds that

[TABLE]

Then, using $If_{0}(t,x)=f_{0}(z-c_{0}t)$, we get

[TABLE]

which admits the equivalent form

[TABLE]

This results to (23) for $\tilde{f}$ and $\tilde{g},$ as in (24).

Remark \theremark:

*The function $\tilde{g}$ describes the total amount of light which travels back from the remaining $N-1$ layers, meaning it considers all multiple reflections. *

Lemma 2.1**.**

Let $L$ be a multi-layer medium consisting of $N\in\mathbbm{N}$ layers, and let the refractive index be given by (11). Then, the solution of (12), together with $u(t<0,z)=f_{0},$ can be computed layer-by-layer.

Proof**:**

Starting with the first layer, we use $n$ defined in (11) and (21) and we apply Proposition 2.1. We thus obtain $\tilde{f}$ and $\tilde{g},$ presented in Proposition 2.1. Then, the function $\tilde{f}$ is the initial wave for the corresponding problem with parameter now given by

[TABLE]

where $\tilde{n}$ represents the refractive index of the next $N-2$ layers. Repeating the same argument, we use (22), with $f_{0}$ replaced by $\tilde{f},$ for the updated operators. We continue this procedure for the new parameters and operators and we end up with the solution for $n$ given by (11).

After some lengthy but straightforward calculations, we can generalize the formulas of Example 2.1 for the $k$-th layer of the medium and $\eta_{r}\in\mathbbm{N},$ resulting to the field

[TABLE]

valid for an initial function $f_{0},$ with $\mbox{supp}f_{0}\subset(-\infty,z_{1}).$

2.2. Dispersive medium

In this section, we consider the form (10) for the refractive index and we set $\hat{n}_{0}(\omega)=\hat{n}_{0}>0.$ We assume $\operatorname{\Re e}\{\hat{n}(\omega)\}>0,$ and $\operatorname{\Im m}\{\hat{n}(\omega)\}\geq 0,$ for all $\omega\in\mathbbm{R}.$ We recall (8). Unfortunately, an explicit solution, as in Sec. 2.1, cannot be derived here for a time dependent parameter. However, applying the Fourier transform with respect to time, we get the Helmholtz equation

[TABLE]

For $y\in(z_{1},\,z_{2}),$ we define

[TABLE]

with

[TABLE]

The refractive index $\tilde{n}$ accounts for the parameter of the remaining $N-1$ layers.

Initially, we consider the problem of a right-going incident wave of the form

[TABLE]

incident at the interface $z=z_{1}.$ Then, the corresponding problem reads

[TABLE]

for an artificial boundary point $z_{1}^{+}>z_{1}$. The boundary radiation condition is such that there is no left-going wave at the region $(z_{1},+\infty).$

The solution admits the form

[TABLE]

where we define the reflection and transmission operators $R$ and $T,$ respectively, by

[TABLE]

The solution operator is then given by $V_{1}:\hat{u}_{0}\mapsto\hat{u}$. The next sub-problem is described by

[TABLE]

for an incident left-going wave of the form $\hat{u}_{0}(\omega,z)=\beta_{1}(\omega)e^{-i\frac{\omega}{c}\hat{n}_{1}(\omega)z},$ and an artificial boundary point at $z_{1}^{-}<z_{1}$. The boundary radiation condition is that the left-going wave in $(-\infty,z_{1})$ is zero. The solution now is given by

[TABLE]

where $R_{-}$ and $T_{-}$ are defined by

[TABLE]

Let again $V_{2}:\hat{u}_{0}\mapsto\hat{u},$ denote the corresponding solution operator.

The final sub-problem, deals with the scattering of a right-going wave of the form $\hat{u}_{0}(\omega,z)=\alpha_{1}(\omega)e^{i\frac{\omega}{c}\hat{n}_{1}(\omega)z},$ by a medium supported in $(z_{2},\,+\infty)$ with refractive index $\tilde{n}.$ The governing equations are

[TABLE]

The radiation condition now ensures that in $(z_{2},\,+\infty)$ exist only right-going waves. The solution is given by

[TABLE]

We define $R_{+}:\alpha_{1}(\omega)\mapsto\beta_{1}(\omega),$ and the relevant operator $V_{3}:\hat{u}_{0}\mapsto\hat{u},$ mapping the incident field to the solution of this specific problem. We remark that the operator $R_{+}$ cannot be computed explicitly because it contains also the information from the remaining $N-1$ layers.

Proposition \theproposition.

Let the incident wave be of the form (27). We define

[TABLE]

Then, the field

[TABLE]

for $y\in(z_{1},z_{2}),$ is the solution of the Helmholtz equation (26), for the refractive index defined as above.

Proof**:**

By construction, $\hat{u}$ fulfills (26) in $(-\infty,y)$ and $(y,\infty)$. Thus, it remains to show that the two parts coincide in $(z_{1},z_{2}).$ Recalling the definitions of $V_{1},\,V_{2},$ and $V_{3},$ restricted in $(z_{1},z_{2})$, we get

[TABLE]

and

[TABLE]

We reorder the terms and we observe that they coincide in $(z_{1},z_{2}).$

Remark \theremark:

If $L=(z_{1},z_{2})$ denotes a single-layer medium, with material parameter $\hat{n}$, given by

[TABLE]

*then we can compute $R_{+}$ explicitly, and also the operator $V_{3}$. *

Example \theexample.

The amplitude of the $j$-th reflection in a certain layer is described by the term

[TABLE]

The single reflected wave from the most left boundary of $L$ is given by (28). For the $k$-th layer of the medium, we obtain the back-reflected field

[TABLE]

where $\alpha_{0}$ is the amplitude of the incident wave $\hat{u}_{0}$.

The solution of (26) in $(-\infty,z_{1})$ admits the form

[TABLE]

and we define $\tilde{R}(\omega):\alpha_{0}(\omega)\mapsto\tilde{\beta}(\omega)$.

Lemma 2.2**.**

Let the incident wave be of the form (27), and let $\tilde{R}[\alpha_{0}]$ be known. Then, the following relation holds

[TABLE]

for

[TABLE]

calculated from the previously defined operators $R_{-}$, $T$ and $T_{-}$.

Proof**:**

We know that in $(-\infty,z_{1})$,

[TABLE]

for $\hat{u}^{j}_{0,-},$ defined as in Proposition 2.2. Using the definitions of $V_{1}$ and $V_{2},$ we get

[TABLE]

This results to

[TABLE]

which is equivalent to

[TABLE]

This completes the proof.

The amplitudes $\tilde{\alpha}_{0}$ and $\tilde{\beta}_{0}$, defined in (31), correspond to the amplitudes of the Fourier transforms of $\tilde{f}$ and $\tilde{g},$ given in Proposition 2.1. Furthermore, one can derive an analogue of Lemma 2.1 also for a dispersive medium.

3. The inverse problem

We address the inverse problem of recovering the position, the size and the optical properties of a multi-layer medium with piece-wise inhomogeneous refractive index. We identify the position by the distance from the detector to the most left boundary of the medium, and the size by reconstructing the constant refractive index $n_{0}$ of the background medium. Initially, we discuss the problem of phase retrieval and possible directions to overcome it and then we present reconstruction methods for non-dispersive and dispersive media. We end this section by giving a method, which with the use of the Kramers-Kronig relations, makes the reconstruction of a complex-valued refractive index (absorbing medium) possible. Let $z=z_{d}$ denote the position of the point detector.

3.1. Phase retrieval and OCT

The phase retrieval problem, meaning the reconstruction of a function from the magnitude of its Fourier transform, has attracted much attention in the optical imaging community, see [27] for an overview. When dealing with experimental data, additional problems arise, like different types of noise and incomplete data. Mathematically speaking, the problem corresponds to a least squares minimization problem for a non-convex functional. In our case, where we are given one-dimensional data of the form (6) or (7), unique reconstruction of the phase is not possible [15]. However, there exist convergent algorithms that produce satisfactory results under some assumptions on the signal, like bounded support and non-negativity constraints. These algorithms are alternating between time and frequency domains, using usually less coefficients than samples, which makes the exact recovery almost impossible.

In OCT, this problem has been also well studied, see for example [21, 23, 26]. The main idea is either to consider a phase-shifting device in the reference arm or to combine OCT with holographic techniques. The first case, the one we consider here, produces different measurements by placing the mirror at different positions, meaning by changing the path-length difference between the two arms.

As already discussed in Sec. 1, we have measurements of the form

[TABLE]

for $r$ fixed, where $u_{r}$ denotes the $y-$component of the reference field $E_{r},$ and we also acquire the data

[TABLE]

Since we know the incident field $\hat{u}_{0}$ explicitly, we can also compute the reference field $\hat{u}_{r}-\hat{u}_{0}$ at the point detector. Then, the problem of phase retrieval we address here is to recover $\hat{u}-\hat{u}_{0}$ from the knowledge of $\hat{m}$ and $\hat{m}_{s}$ for all $\omega\in\mathbbm{R}.$ We know, from [18, 19], that if $u_{r}-u_{0}$ is compactly supported, then there exist at most two solutions $u-u_{0}$. See the left picture of Fig. 3, where we visualize graphically the two solutions by plotting in the complex plane the two above equations at specific frequency for the setup of the third example presented later in Sec. 4.

If in addition, there exist a constant $\gamma\in[0,1),$ such that

[TABLE]

then, there exists at most one solution in $L^{2}(\mathbbm{R})$ with compact support in $[0,\infty).$ However, it is hard to verify that the reference field fulfills this condition and we observed that, in all numerical examples, this inequality does not hold, for an incident plane wave. Thus, in order to decide which solution of the two is correct, extra information is needed. Motivated by the phase-shifting procedure, we consider data for two different positions of the mirror, let us say $r=r_{1},r_{2}.$ Then, we get the data

[TABLE]

Using $\hat{m}(r_{1};\cdot)$ and $\hat{m}_{s},$ we get two possible solutions, and from $\hat{m}(r_{2};\cdot)$ and $\hat{m}_{s},$ other two. But since $\hat{m}_{s}$ is the same in both cases, we find the unique solution as the common solution of the two pairs. This is illustrated at the right picture of Fig. 3, where we plot the three above relations at two different frequencies. This way, we get unique solutions at every available frequency.

Thus, having measurements for two different positions of the reference mirror, we may consider that frequency-domain OCT provide us with the quantity

[TABLE]

the equivalent measurements of a time-domain OCT system, see (5).

3.2. Reconstructions in time domain

We consider initially a non-dispersive medium. Then, the refractive index is given by (11) and the time-dependent OCT data admit the form

[TABLE]

for $T>0.$ Here, we assume that the initial wave (known explicitly) does not contribute to the measurements. The following presented algorithms are based on a layer-by-layer procedure. At the first step, we do reconstruct the parameters for a given layer and then we update the data, to be used for the next layer. Thus, we assume that we have already recovered the boundary point $z_{k-1}$ and the coefficient $c_{k-1}$ of the layer $L_{k-1}.$ We denote by $m^{(k)}(t,z_{d}),$ the data corresponding to a $(N-k+1)-$layer medium, with the most left layer being the $L_{k}.$

As we see in Fig. 1, the (time-dependent) data consist mainly of $N+1$ major $``$peaks$"$, for a $N-$layer medium, and some minor $``$peaks$"$, related to the light undergoing multiple-reflections in the medium. The experimental data are, of course, also noisy and may show some small $``$peaks$"$ because of the OCT system. The first two $``$peaks$"$ correspond, for sure, to the single-reflected light from the first two boundaries. We propose a scheme to neglect minor $``$peaks$"$ due to multiple reflected light. Then, the first major $``$peak$"$ in $m^{(k)},$ corresponds to the back-reflected light from the interface at $z_{k}$.

Step 1:

We isolate the first $``$peak$"$ by cutting off around a certain time interval $[T_{1},T_{2}],$ meaning we consider

[TABLE]

The time interval can be fixed for all layers and depends on the time support of the initial wave. On the other hand, this wave can be described by (25), if we use $\eta_{r}=1$ and $k:=k-1.$ Therefore, we obtain the equation

[TABLE]

with $\tilde{f}_{0}$ given by

[TABLE]

The supremum of (34), using its shift-invariance property, gives the value of $c_{k}$. The position of the interface $z_{k}$ can then be recovered from (34), by solving the minimization problem

[TABLE]

Since both parameters are time-independent, there exist also other variants for solving this overdetermined problem.

Step 2:

Before moving to the layer $L_{k+1}$, we have to update the data function. We could just remove the contribution of the current layer, meaning $\tilde{m}^{(k)}.$ However, the first $``$peak$"$ might not correspond to the reflection from $z_{k},$ for $k>1,$ but to contributions of multiple reflected wave from previous layers, arriving at the detector before the major wave. Since, we have recovered the properties of $L_{k},$ we can compute all future multiple reflections from this layer using (25), let us call them $\mathcal{R}[c_{k},\,z_{k}].$

Then, we update the data as

[TABLE]

Repeating these steps, we end up with the following result.

Lemma 3.1**.**

Let $L$ be a multi-layer medium, with $N\in\mathbbm{N}$ layers, characterized by $n,$ given by (11). Then from the knowledge of $n_{0}$, the initial wave $f_{0},$ and the measurement data (33), following the above iterative scheme, we can uniquely reconstruct $n_{j}$ and $z_{j}$ for $j=1,...,N+1.$

The above scheme can be written in an operator form, by the application of Propositions 2.1 and 2.1. For the sake of presentation, we consider the case of the first layer in order to avoid redefining all operators.

Step 1:

Recall the definition of the operator $\tilde{R}$, applied to the initial wave $f_{0},$ which describes the total amount of the reflected light. Considering the data (33), we get

[TABLE]

Following the same procedure, using (17), we can recover $c_{1}$ and $z_{1},$ from the reduced data equation

[TABLE]

Step 2:

We update all operators and from Proposition 2.1, we obtain $\tilde{f}$ and $\tilde{g}$, given by (24). From the definition of $\tilde{g}$, we see that we obtain the function

[TABLE]

describing the updated data. The advantage here is that we do not need to subtract the multiple reflections term, since they are already included in the updated version of $\tilde{g}.$

3.3. Reconstructions in frequency domain

We consider the case of a dispersive medium, with a piece-wise inhomogeneous refractive index $\hat{n}(\omega),$ for $\omega\in\mathbbm{R}.$ Initially, we assume $\hat{n}\in\mathbbm{R}.$ The derived iterative scheme can be applied also to the simpler non-dispersive case, giving a reconstruction method in the frequency domain for a time-independent parameter.

3.3.1 Dispersive medium

The refractive index is given by (10) and initially, we restrict ourselves to the case of real-valued $\hat{n}.$ The data $\hat{m}$ are given by (32). As before, we present a layer-by-layer scheme. We assume that the boundary point $z_{k-1}$ and the coefficient $\hat{n}_{k-1}$ of the layer $L_{k-1}$ are already recovered. We denote by $\hat{m}^{(k)}(\omega,z_{d})$, the data corresponding to the $(N-k+1)-$layer medium. Here, we need to choose the time interval $\mathds{1}_{[T_{1},T_{2}]}$ wisely, compared to Sec. 3.2, because of the presence of dispersion. However, this results only to slightly longer time interval, since in the wavelength range, where OCT operates, scattering dominates absorption.

Step 1:

As we have seen in Fig. 1, for example, from the data $\hat{m}^{(k)}$ we cannot distinguish the different $``$peaks$".$ Thus, we have to switch back to the time domain in order to isolate the first $``$peak$".$ We apply

[TABLE]

We use (30) for $k:=k-1$ and $\eta_{r}=1$ to get

[TABLE]

with

[TABLE]

describing an already known quantity. The absolute value of (35) gives

[TABLE]

which together with a suitable condition on $\hat{n}_{k}$ allow us to recover the refractive index.

To reconstruct the position $z_{k}$ we define the function

[TABLE]

and we observe that the absolute value of its derivative, together with (35), results to

[TABLE]

Together with a non-negativity constrain, we obtain $z_{k}.$

Step 2:

As in the time domain case, we update the data by subtracting (35) and the terms representing the multiple reflections from the already recovered layers, called $\mathcal{R}[\hat{n}_{k},z_{k}].$ We define

[TABLE]

Repeating the steps, a reconstruction of the properties and the lengths of all the remaining layers is obtained.

Lemma 3.2**.**

Let $L$ be a multi-layer medium, with $N\in\mathbbm{N}$ layers, characterized by $\hat{n}$, given by the Fourier transform of (10). If we restrict ourselves to the case $\hat{n}(\omega)\in\mathbbm{R},$ for all $\omega\in\mathbbm{R}$, then the above iterative scheme, allows us to uniquely reconstruct $\hat{n}_{j}$ and $z_{j}$ for $j=1,...,N+1,$ given $\hat{n}_{0}$, the incident wave $\hat{u}_{0}$ and the measurement data (32).

3.3.2 Absorbing medium

Here, we consider the case of a complex-valued material parameter, $\hat{n}(\omega)\in\mathbbm{C},$ for every $\omega\in\mathbbm{R}.$ The real part describes how the medium reflects the light and the imaginary part (wavelength dependent) determines how the light absorbs in the medium. In [5, 7] we considered the multi-modal PAT/OCT system, meaning that we had additional internal data from PAT, in order to recover both parts of the refractive index. Here, the multi-layer structure allows us to derive an iterative method that requires only OCT data. We decompose $\hat{n}$ as

[TABLE]

The measurements are again given by (32). The above presented iterative scheme, fails in this case. Indeed, recall the formula (35), which we considered for recovering the parameters of the $k-$th layer. Taking the absolute value, for $\hat{n}_{k-1}(\omega),\,\hat{n}_{k}(\omega)\in\mathbbm{C},$ gives

[TABLE]

The last term in the above expression, describes how the amplitude of the wave decreases, i.e. attenuation, and prevents us from a step-by-step solution, since $z_{k}$ still appears. Thus, we propose a different scheme that takes into account also the relations between the parts $\nu$ and $\kappa,$ meaning the Kramers-Kronig relations. We stress here that $\hat{n}$ is holomorphic in the upper complex plane, satisfying $\hat{n}(\omega)=\hat{n}(-\omega)^{\ast}.$ The parts of the complex-valued refractive index are connected through

[TABLE]

In addition, defining the reflection coefficient,

[TABLE]

and using its expression in polar coordinates $\rho_{k}=|\rho_{k}|e^{i\theta_{k}},$ we obtain

[TABLE]

a function that diverges logarithmically as $\omega\rightarrow\infty,$ and is not square-integrable [22]. However, the following relation holds for the phase of the complex-valued reflectivity

[TABLE]

We define the operator

[TABLE]

and we get the relations in compact form

[TABLE]

Of course, when working with the Kramers-Kronig relations (36) and (38), one has to deal with the problem that they assume that information is available for the whole spectrum, something that is not true for experimental data. Another problem could be the existence of zero’s of the reflection coefficient in the half plane. There exist generalizations of there formulas that can overcome these problems, like the subtractive relations which require few additional data. We refer to [12, 16, 24, 35] for works dealing with the applicability and variants of the Kramers-Kronig relations. This practical problem is out of the scope of this paper and will be considered in future work, where we will examine numerically, with simulated and real data, the validity of the proposed scheme.

Step 1:

At first, we consider the reconstruction of the interface $z_{k}.$ We apply once the logarithm to the absolute value of (35) and then we take imaginary part of the logarithm of (35). Using the definition (37), we obtain the system of equations

[TABLE]

We define the data functions

[TABLE]

and the system (39) takes the form

[TABLE]

to be solved for $z_{k}.$ We rewrite the last equation using the formulas (36) and (38) and the first equation, to obtain

[TABLE]

The last equation is solved at given frequency $\omega^{\ast}\neq 0,$ in order to obtain the location of the interface

[TABLE]

We can now recover $\hat{n}_{k}$ from (35) which admits the from

[TABLE]

by equating the real and imaginary parts.

Step 2:

We update the data as in the non-absorbing case.

4. Numerical Implementation

We solve the direct problem considering two different schemes depending on the properties of the refractive index, see (10) and (11). First, we consider the time-independent refractive index and phase-less data from a time-domain OCT system.

4.1. Reconstructions from phase-less data in time domain

We model the incident field as a gaussian wave centered around a frequency $\omega_{0}$ moving in the $z$-direction of the form

[TABLE]

with width $\sigma,$ where $z_{0}$ denotes the source position and $c\approx 3\times 10^{8}$m/s is the speed of light.

The simulated data are created by solving (12) using a finite difference scheme. We restrict $z\in[0,1.5]$mm and we set $T>0,$ the final time. We consider absorbing boundary conditions at the end points and we set $u(0,z)=u_{0}(0,z)$ and $\partial_{t}u(0,x)=0$ as initial conditions. The left-going wave is ignored. We consider equidistant grid points with step size $\Delta z=\lambda_{0}/100,$ where $\lambda_{0}=2\pi c/\omega_{0},$ is the central wavelength, and time step $\Delta t$ such that the CFL condition is satisfied.

The measurement data are given by

[TABLE]

where $z_{d}$ denotes the position of the point detector. We add noise with respect to the $L^{2}$ norm

[TABLE]

where $v$ is a vector with components normally distributed random variables and $\delta$ denotes the noise level. We have to stress that the total time is such that the data contain also information from the multiple reflections inside the medium. We define the length of the $k-$th layer

[TABLE]

and we set $\ell_{-1}=|z_{d}-z_{0}|$, and $\ell_{0}=z_{0}-z_{d}.$ We denote by

[TABLE]

the reflection coefficient at the interface $z=z_{k}.$

Here, we assume that we know only $n_{0}=1,$ and the positions of the source and the detector. Thus, we aim for recovering the position, the size and the optical properties of the medium.

The proposed iterative scheme for a $N$-layer medium is presented in Algorithm 1, where the output is the reconstructed refractive indices and lengths of the layers. First, we order the observed $``$peaks$"$ at the image with respect to time, producing the set of data $(t_{l},\,p_{l}),$ for $l=1,2,...,\Lambda.$ The number $\Lambda\geq N$ describes the number of single and multiple reflections arrived at the detector before the final time $T$. In order to obtain a physically compatible solution we impose some bounds $[\underline{n},\,\overline{n}]$ on the refractive index. This condition is not necessary for data with phase information. We update the data by neglecting the multiple reflections. To do so, once we have recovered the length and the refractive index of a layer, we neglect the $``$peaks$"$ appearing later referring to multiple reflections inside this layer. Of course, because of numerical error and noisy data, we give a tolerance depending on the time duration of the wave.

We define the error function

[TABLE]

where $(n_{k},\,\ell_{k})$ and $(\tilde{n}_{k},\,\tilde{\ell}_{k}),$ for $k=1,...,N$ are the exact and the reconstructed values, respectively.

In the first example, we consider a three-layer medium positioned at $\ell_{0}=0.5$mm, with $\ell_{-1}=0.$ We set $(n_{1},\,n_{2},\,n_{3})=(1.55,\,1.41,\,1.48)$ and lengths $(\ell_{1},\,\ell_{2},\,\ell_{3})=(0.2,\,0.3,\,0.1)$mm. The obtained data (41) for this example are given at the left picture in Fig. 4.

The results are presented in Table 1 for $[\underline{n},\,\overline{n}]=[1.345,2]$ and $\mbox{tol}=0.1$ps. We obtain accurate and stable reconstructions, with $\epsilon=3.34\times 10^{-7}.$ This algorithm can be easily applied to multi-layer media and it is presented here since it will be the core of the more complicated algorithm in the Fourier domain.

4.2. Reconstructions having phase information in frequency domain

We aim to reconstruct the time-dependent refractive index, meaning its frequency-dependent Fourier transform. As discussed already in Sec. 3.1, it is possible from the phase-less OCT data to recover the full information, implying that we consider as measurement data the function

[TABLE]

Here, the frequency interval $[\underline{\omega},\overline{\omega}],$ with $\overline{\omega}>\underline{\omega}\gg 0,$ models the OCT data, recorded by a CCD camera placed after a spectrometer with wavelength range $[2\pi c/\overline{\omega},2\pi c/\underline{\omega}],$ in a frequency-domain OCT system.

4.2.1 Non-dispersive medium

In order to construct the data (42), we consider the time-dependent back-reflected field derived in the previous section, we add noise and we take its Fourier transform with respect to time. Then, we truncate the signal at the interval $[\underline{\omega},\overline{\omega}],$ see the right picture in Fig. 4.

In Algorithm 2 we present the main steps of the iterative scheme as described in Sec. 3.3. In Step 1, we take advantage of the causality property of the time-dependent signal and we zero-pad $\hat{m}(\omega),$ for all $\omega\in\mathbbm{R}\setminus[\underline{\omega},\overline{\omega}],$ and then we recover the signal as two times the real part of the inverse Fourier transformed field. In the second step, we initially approximated the derivative with respect to frequency using finite differences but it did not produce nice reconstructions due to the highly oscillating signal. We replace the derivative with a high-order differentiator filter taking into account the sampling rate of the signal. In Fig. 6, we plot the function $\phi(\omega),\,\omega\in[\underline{\omega},\overline{\omega}],$ and we see that it is constant in a central interval (called trusted) and oscillates close to the end points. Thus, we denote by $\omega^{\ast}$ either a chosen frequency in the trusted interval or the mean of the frequencies in this trusted interval. We address here that one could also average over the whole spectrum and still get reasonable results.

In the second example, we consider again a three-layer medium with parameters $(n_{1},\,n_{2},\,n_{3})=(1.55,\,1.405,\,$ $1.55)$ and lengths $(\ell_{1},\,\ell_{2},\,\ell_{3})=(0.15,\,0.5,$ $0.13)$mm. Here, $\ell_{0}=0.7$mm and $\ell_{-1}=0.2$mm. The central frequency is given by $\omega_{0}=2\pi c/\lambda_{0},$ with $\lambda_{0}=800$nm. The sampling rate is $f_{s}=100c/\lambda_{0}.$ The recovered parameters for noise-free and noisy data are presented in Table 2. The relative error is $\epsilon=2.164\times 10^{-5}.$ In Fig. 6, we see how the data change as the Algorithm 2 progresses. The picture on the right shows the data in frequency domain with respect to the $``$peaks$"$ presented in the left picture where we see the time-domain data. The yellow curve (right) represents the full data, the red curve (right) the data if we neglect the red $``$peak$"$ (left), the blue curve shows the data if we neglect also the blue $``$peak$"$, and so on. The green curve (right) represents the signal from the multiple reflections. We observe that the OCT signal maintains the Gaussian form of the incident wave, centered around the central frequency, and the different reflections result to the oscillations of the field.

4.2.2 Dispersive medium

The incident field in the frequency domain takes the form

[TABLE]

which is the Fourier transform with respect to time of $u_{0},$ given by (40), restricted to positive frequencies. This field describes a plane wave moving in the $z-$direction having a Gaussian profile perpendicular to the incident direction, centered around $\omega_{0}.$ We generate the data considering the formula (30) and then we add noise. The Algorithm 3 summarizes the steps of the iterative scheme, which for a frequency-independent refractive index simplifies to Algorithm 2.

We model the wavelength-dependent refractive index of the medium using the standard formula, known as Cauchy’s equation,

[TABLE]

for some fitting coefficients $\beta_{j},\,j=1,2,3.$ In Fig. 8, we see the exact refractive index of the first (left) and the third (right) layer for the medium used in the third example. Afterwards, we consider the refractive index as a function of frequency.

As already discussed, the calculations close to the end points were not stable and since here we are interested in reconstructing the frequency-dependent refractive index, we restrict the computational domain and then we extrapolate the recovered functions in order to update the data.

The medium lengths are given by $(\ell_{1},\,\ell_{2},\,\ell_{3})=(0.2,\,0.3,\,0.1)$mm, and we set $\ell_{0}=0.7$mm and $\ell_{-1}=0.$ The second layer has constant refractive index given by $n_{2}(\omega)=1.41.$ The reconstructions of $n_{1}(\omega)$ and $n_{3}(\omega)$ are presented in Fig. 8 for data with $2\%$ noise. In Table 3, we see the recovered lengths and the refractive indices at specific frequencies.

5. Conclusions

In this work we addressed the inverse problem of recovering the optical properties of a multi-layer medium from simulated data modelling a frequency-domain OCT system. We considered the cases of non-dispersive, dispersive and absorbing media. We proposed reconstruction methods and we presented numerical examples justifying the feasibility of the derived schemes. Stable reconstruction with respect to noise were presented. The methods are based on standard equations, equivalent to the Fresnel equations, and to ideas from stripping algorithms. The originality of this work lies in the combination of them into a new iterative method that addresses also the frequency-dependent case, which needs special treatment. As a future work, we plan to examine the applicability of the iterative schemes for experimental data and test numerically the method for absorbing media.

acknowledgement

The work of PE and LV was supported by the Austrian Science Fund (FWF) in the project F6804–N36 (Quantitative Coupled Physics Imaging) within the Special Research Programme SFB F68: $``$Tomography Across the Scales$".$