Quantitative Parameter Reconstruction from Optical Coherence Tomographic Data

Quantitative tissue information, like the light scattering properties, is considered as a key player in the detection of cancerous cells in medical diagnosis. A promising method to obtain these data is optical coherence tomography (OCT). In this article, we will therefore discuss the refractive index reconstruction from OCT data, employing a Gaussian beam based forward model. We consider in particular samples with a layered structure, meaning that the refractive index as a function of depth is well approximated by a piece-wise constant function. For the reconstruction, we present a layer-by-layer method where in every step the refractive index is obtained via a discretized least squares minimization. For an approximated form of the minimization problem, we present an existence and uniqueness result. The applicability of the proposed method is then verified by reconstructing refractive indices of layered media from both simulated and experimental OCT data.


Introduction
Optical coherence tomography is a non-invasive, high-precision imaging modality with micrometer resolution based on the interferometric measurement [11] of backscattered light.Since its invention in the 1990's [15], different optical coherence tomographic systems, see [5], for example, have been developed.All share the basic working principle, which can be described as follows: Laser light in the near infrared region is sent into the system.A beamsplitter then separates the light in two parts.One is directed into the sample arm and the second into the reference arm.Depending on the type of the system, the object in the sample arm is then either raster scanned, meaning that depth profiles on a lateral grid across the object are obtained, or illuminated at once.In both cases, the backscattered light by the object is then coupled into the system again and transported to the detector.There, the combined intensity of the light from the sample and the reference arm, where the light is backreflected by a perfect mirror, is detected.
Due to its (ultra-)high resolution property and its very fast acquisition rate of up to more than 50000 raster scans per second, OCT, together with its adaptions and its extensions like polarization-sensitive OCT (PS-OCT) [1,2,12], optical coherence angiography [18,19] or optical coherence elastography (OCE) [17], has become an outstanding technology in imaging of biological tissues, especially in the field of ophthalmology.
Nowadays, the structural information which is provided by such an OCT system is already successfully used for medical diagnosis.However, the obtained information is mainly of qualitative nature.For medical diagnosis it is worth to have supplementary quantitative information, like the optical (scattering) properties of the object of interest, which are considered as future key makers in the field of medical diagnosis.
The quantification of such physical properties, in our case the refractive index, from OCT data, commonly represents an ill-posed problem, which is interesting from a mathematical perspective.This is the scope of this article.In its most general form, the parameter quantification in (Fourier domain based) OCT is a severely ill-posed problem dealing with the reconstruction of the optical properties, which we describe by the refractive index which is a function of space (three dimensions) and wavenumber (one dimension), from maximally three-dimensional (two spatial dimensions and one for the wavenumber) measurement data [6].The insufficient amount of available data has driven the need for modeling possibilities of the direct problem, which are commonly expressed by assumptions on the object.A general overview of existing modeling possibilities for the direct problem in OCT is presented, for example, in [6] or [7].
The inverse problem in OCT can be seen as an inverse electromagnetic scattering problem, especially if one considers the ideal case of having direct access to the backscattered sample.It has been treated, also non-specifically related to OCT, in a theoretical context, for example, in [4] and under a number of (simplifying) assumptions, like having a weakly scattering medium, in [20].
The aim of this article is to provide a reconstruction method from experimental data obtained by a swept-source OCT system, a specific form of Fourier domain OCT, which is based on a raster scanning of the object with focused and practically monochromatic laser light centered at multiple wavelengths.Due to this raster scanning process, the (inverse) scattering problem can be restricted to the narrow illuminated region, where the refractive index is typically considered only as a depth dependent function.
One of the very first attempts on reconstructing a depth dependent refractive index from an OCT depth profile by Fourier transform has been presented in [10].The reconstruction under the, in this case, crucial assumption of a weakly scattering object (so that the Born approximation of the electromagnetic wave equation is valid) has also been treated in [23], for example, where additionally a Gaussian beam model for the incident laser light has been employed.
Another assumption on the object, which substantially simplifies the mathematical model, is that the object shows a multi-layer structure.This case typically can be identified in OCT images of the human retina [13] and in human skin imaging.For such a multi-layer structure, at least locally in the illuminated spot, the depth-dependent refractive index can be simplified to a piecewise constant function.The corresponding inverse problem has been examined in [3,8,26] and in [9] where additional inclusions have been discussed.
In this article, we consider the multi-layer structure of the sample with a Gaussian beam model, presented in [27], which resembles more efficiently the laser light illumination.Within this context, we treat the inverse problem of quantifying the refractive index by using a layer-by-layer method, where in each step the reconstruction concentrates on a pair of parameters, the refractive index and the width of the layer.
Hereby, the reconstruction is formulated as an ℓ 2 -minimization problem in the Fourier domain between the Gaussian beam prediction model and the data.That is, we match the absolute values of both in order to avoid high-frequency components which are typically present.In [8], the reconstruction within this setting for the case of a plane wave incident field only gave a unique solution when artificially adding plausible bounds on each parameter.We can omit these bounds in our analysis on the existence and uniqueness of solutions of the present minimization problem.
The analysis itself is complicated by the compact form of the Gaussian model.There the directional dependent reflection coefficients containing information about the refractive indices hinder a deeper investigation.For this reason, the analysis is divided into two parts: Firstly, we consider an approximation of the original model, where the reflection coefficient is assumed to be directional independent.We show that the modeling error caused by this approximation is small and bounded from above.For the approximation we then show the existence and the uniqueness of a solution under the strong but necessary condition that the data is in the range of the prediction.The distance in every step is obtained by the biggest overlap of the prediction model with the data once the refractive index has been calculated.
In order to support our arguments, we show the functionality of the proposed method by reconstructing the refractive indices from both simulated and experimental data.While reconstructions from simulated data, partially under very theortic assumptions, have been discussed several times, the reconstruction using experimental data has been treated fairly rare.The assumption on the object to have layered structure is a very strong assumption which turns the inverse problem to the well-posed side.Hence, regularization methods which are typically related to problems with noise, are not considered in this context.
The outline of this article is as follows: In Section 2 and Section 3 the major parts of the Gaussian beam model introduced in [27] are summarized.The section ends by giving an explicit representation of a single OCT depth profile.Based on this representation, we introduce the inverse problem which we want to treat layer-by-layer.We discuss this in Section 4. That the actual reconstruction can be formulated as a layer-by-layer method is justified in Section 5 and its reformulation as ℓ 2 -minimization problem is shown in Section 6.Before we present in Section 7 a characterization of existence and uniqueness of the solutions to an approximated version of the original problem, we show that the approximation error between the functionals is almost negligible.We conclude the section on the refractive index quantification by presenting a purely theoretical method to also retrieve uniqueness.The analysis on the minimum with respect to the width is presented in Section 8 in form of purely visual arguments by showing the behavior of the minimization functional with respect to the width.There, we consider simulated data with and without noise.Finally, in Section 9 the results of the reconstruction from simulated and experimental data of a three-layer object are presented.

OCT Forward Model
The image formation in standard OCT systems nowadays is based on raster scanning the object of interest.This means that strongly focused laser light is directed by movable galvanometric mirrors to spots located on an imaginary lateral (horizontal) grid on the object's top surface.For every position then an OCT measurement, a depth profile along the vertical axis through the spot showing the microstructure inside the object, is recorded.
Latest publications in the field of OCT have pointed out that modeling such focused laser light by a single plane wave is not sufficient enough in capturing several system relevant aspects.The huge influence of the focus and the beam width, the diameter of the laser in the focus spot, for example, are a few modeling aspects which are lost when simulating the measurements of an OCT system using single plane waves.For this reason we base our analysis on a Gaussian beam model for the incident field, which is said to model laser light in an accurate way.Hereby, we adapt for our purpose the forward model which has been presented in [27].There, a model is provided which includes all relevant system parameters and allows a good approximation of an actual measurement.
In this section we summarize briefly the main parts of an OCT system, namely the light scattering, the fiber coupling and finally the measurement.Hereby, we use the fact that the object is raster scanned, which allows us to model in the following the light propagation for each raster scan independently.At the end of this section we provide an equation for a single raster scan (A-scan) data.This will be used as the data for the corresponding inverse problem.
The main parts are specifically modeled for a swept-source OCT system [5], since the data used for the numerical experiments in Section 9 has been obtained by this system type.Within such a system almost monochromatic laser light centered at multiple wavelengths in a certain spectrum is used for the illumination of the object.For each of these wavelengths the scattered light is detected by a narrow-bandwidth interferometer.This finally gives a complete measurement for each wavelength within the spectrum.
Field of incidence.We model the light propagation from the point on where the laser light has entered the sample or the reference arm respectively.The presentation restricts to the modeling of the sample arm.The reference arm, which is modeled analogously with the difference that the object is a perfectly reflecting mirror, is considered as a special case.
For each raster scan we model the incident illumination as an electromagnetic wave Ê : R × R 3 → C 3 , a function of the wavenumber k ∈ R and the spatial coordinate x ∈ R 3 , satisfying the system of Maxwell's equations in vacuum, that is (in Fourier space) where n0 = 1 is the free-space refractive index.Since we consider a swept-source OCT system, Ê is assumed to fulfill the support condition for a sufficiently small parameter ϵ0 > 0. The experiment is repeated for multiple wavelengths k0 in a spectrum S = [k1, k2].The delta-like support in wavenumber allows us to model each experiment, including the backscattering of the light by the sample, independently.This means we can solve the Helmholtz equation for the electric field for each k ∈ S separately.However, there is no difference to the mathematical formulation for a broadband illumination, an electric field which solves the Helmholtz equation for all k ∈ S simultaneously.Hence, we formulate everything simply for a broadband illumination.
In order to specify the Gaussian beam for the field of incidence, we impose additionally an initial condition at a hyperplane {x ∈ R 3 | x3 = r0}, which we locate in the focus of the beam.Hence, by r0 ∈ R we denote the focus position, which we always assume to be along the vertical line Re3, where (ei) 3  i=1 denotes the standard basis of R 3 .For a function f k : R 2 → C with compact support in D k (0) ⊂ R 2 , the ball with radius k and center zero, and a vector η ∈ S 1 × {0}, we then write where x = (x1, x2) and where denotes the two-dimensional Fourier transform with respect to the variable x.
We consider in the following illumination from the top only.Hence, we eliminate one major propagation direction from the obtained solution Ê [27] to the combined problem of Equation 1 and Equation 2 and consider only downward propagating waves, that is in direction −e3.We denote this part by E, which for any fixed wavenumber k ∈ S is then represented by where Remark 2.1 The incident wave in Equation 3 represents a weighted superposition of plane waves e iK•x , where the wave vector K represents the propagation direction of each plane wave.We note that every wave direction K in Equation 3 is implicitely given as a function of κ, that is We suppress this dependence in the following.
We complete the Gaussian beam incident field by specifying the weight function f k in the focal plane.Hereby, we use Gaussian weights only.
Assumption 1 In the following, we restrict our attention to Gaussian distributions, meaning that we consider in Equation 2 functions of the form where a > 0 is such that , where χU is the characteristic function of a set U, is almost negligible.
The Gaussian distribution in this case pronounces only those directions κ ∈ D k (0) which have a small radial deviation from zero, others are excluded from the integral in Equation 3. Hence, we can consider the quotient |κ|/k as a small parameter for which the paraxial approximation expressed by the linearization of the square-root is valid.We want to include this fact also in the representation of the polarization vector.
Assumption 2 Restricting to η = e2, we consider an approximated polarization vector with f k as in Assumption 1 in our model, which reduces the electric field to its transverse components.
We combine Equation 3 with Assumption 1 and Assumption 2 and denote the resulting incident field by E (0) = Ẽ(0) e2, with scalar complex field Light scattering.The field in Equation 6 radiates the sample, which we denote by Ω.We characterize the optical properties of Ω by the refractive index represented by a function n : The assumption that n is a real function refers to the case of a sample where the scattering dominates over the thus neglected absorption.Additionally, because of the small bandwidth of the laser, we assume that n is constant with respect to the wavenumber in the spectrum S.
The main assumption however is that Ω shows a multi-layered structure, meaning that Ω consists of a finite union of subsets Ωj, where each is characterized by a homogeneous refractive index nj ≥ 1, with nj ̸ = n l , l ∈ {j − 1, j + 1}.We claim that all layers are parallel and therefore share a single unit normal vector νΩ = (sin θΩ, 0, cos θΩ), for a small angular value of θΩ ∈ R, which we assume for simplicity to be orthogonal to the polarization vector η = e2 and which is pointing outward the object.
To summarize, we consider a sample given by for a sequence of coefficients (aj)j ⊂ R. The total number of layers J ∈ N hereby is an unknown but arbitrary finite number.We additionally assign to every layer Ωj its width, which is given by the positive real number dj = aj − aj+1.The light E : R × R 3 → C 3 in presence of the sample, Equation 7, is then modeled as a solution of the vectorial Helmholtz equation for all x ∈ R 3 : By taking the divergence of this equation, we see that the assumption of piecewise homogeneous layers implies that ∇ • E = 0 in every set Ωj, so that ∇ × (∇ × E) = −∆E and the equation reduces to the simpler Helmholtz equation (see Equation 1) for the second component Ej (the first and third components vanish by our choice of incident field) of the field E inside the layer Ωj, similar to [10,22], where k 2 n 2 0 is replaced by k 2 ñ2 .At the boundaries between the single layers, we claim that the homogeneous parts of the electric field fulfill the continuity conditions for all x ∈ ∂Ωj ∩ ∂Ωj+1 and j ∈ {1, . . ., J}.We finally denote by E (s) the second component of the backscattered field E − E (0) from the object.
Since the Fresnel equations allow us to get an explicit solution for the case of one discontinuity, we can iteratively construct a solution for a finite number of layers as a series of Fresnel solutions, which physically corresponds to multiple light reflections at the boundaries, see [3,8].Since multiply relfected light only contributes minorly, we ignore it in the model for the backreflected sample field.However, we keep in mind that these multiple reflections, at least up to a finite order, can be added to the model at any time.The assumption of a single reflection model makes it possible to use the transmitted part of the light at a certain layer directly as an incident field to the next layer.
A layer-by-layer scheme, similar to the one in [8], allows us to give an explicit representation of the backscattered field E (s) .The incident field hereby is decomposed into its plane wave parts.For each, we use the Fresnel formulas, see [16,Chapter 7], to determine the reflection coefficients (as functions of the propagation direction of the incoming plane wave) which we denote by rj : R 2 → R for every interface ∂Ωj−1 ∩ ∂Ωj (with Ω0 = ΩJ+1 = R 3 \ Ω).Finally, the reflected fields are combined again and we obtain for x, with x3 > xΩ,3, where xΩ is a point on the top surface {xΩ ∈ R 3 : xΩ • νΩ = a1}, the backscattered field with and where we denote by K = K(κ) = (κ, − k 2 − |κ| 2 ) the wave vector (introduced in Equation 4) and by Φ(K) = K − 2 (K • νΩ) νΩ the wave direction of the reflection of a plane wave with incident vector K at an interface with unit normal vector νΩ.Moreover, we introduced for j ≥ 1 the transmission coefficients and the (transmission) phase factors where the angles θ j t (κ) of transmission between the layers j and j + 1 for an incident plane wave with wave vector K(κ) can be iteratively calculated via Snell's law Remark 2.2 We remark that since |K(κ)| = k, the quotient in the definition of the angle of incidence θ 0 t (in Equation 11) is actually a function of the dimensionless parameter κ k .This implies that the same is true also for the angle of transmission θ j t and the phase factor Ψj. Hence, we will always write Ψj( κ k ) in the following.
Fiber coupling.The backscattered light is then coupled into a single-mode fiber and later transfered to the detector.This coupling is done by using a scan lens which discards all (plane wave) parts in Equation 8 of which wave direction strongly deviates from the third unit normal e3.The maximal deviation is described by the (maximal) angle of acceptance, which we denote by θ.The set of incident directions yielding an accepted wave vector is then given by This means that the area of integration in Equation 8 is reduced to the set B ⊂ R 2 .
Because of the restriction to this for typical cases small set the deviation in the directions κ is limited and small.Hence, also the reflection coefficients are only varying slighty (with κ) which motivates us to consider an approximated form of Equation 8, especially in context of the inverse problem in Section 5 and Section 6, where the reflection coefficients are directional independent.For example, if θΩ = 0, the accepted directions are given by B = D k sin θ (0), the ball centered at zero with radius k sin θ.The inclination angle of K(κ) with νΩ is maximal if |κ| = k sin θ, which means that κ is on the boundary of B. The corresponding reflection coefficient is given by (for a single boundary reflection between n and n ′ ), which we then can write as a function of the inclination angle θ 0 t .For |κ| = k sin θ, we get θ 0 t = θ and we can use the approximation .
Under the assumption of cos θ ≈ 1, which is a good approximation in our setting as θ is around 2 • , the reflection coefficient is approximated by the constant n−n ′ n+n ′ , the classical Fresnel coefficient [16].In general, the angular values of θ and θΩ are considered to be small, which allows us to assume that the variation of the reflection coefficient on B is also well approximated by a constant coefficient.Next, we show that for small angular values the difference between the reflection coefficient (as a function of the angle of incidence), see Equation 9, and its approximating constant tends to zero quadratically.
Then for small values of y we have that Proof: We determine the Taylor expansion of r locally around the point y0 = 0. Obviously, r(y) = r(−y), for all y ∈ R, which makes r an even function.Thus, for its Taylor expansion we expect only even exponents.
Estimating now the difference between r † and r for small values of y gives which proves the result.□ Applying this formula in our context, we obtain the strongest deviation from r † on the (closed) set B in the boundary point κ with maximal angular deviation θ0 = θ + 2θΩ from the unit vector −νΩ.Thus, we obtain a maximal error when approximating r by r † by the order of (θ + 2θΩ) 2 .We silently assume that the refraction at the deeper interfaces does not change the angle too much so that this approximation remains valid at all interfaces.
OCT measurement.At the detector, an interference pattern of the backscattered field from the sample combined with a reference field is measured.This reference field is obtained by the reflection of the incident light by a perfectly reflecting non-tilted mirror, which is modeled as semi-infinite object with (infinitely) large refractive index.We assume that the mirror is located around the focus of the beam where the far-field approximation [14,21] is considered to be valid [27].For the (measurement) direction s = e3 and the distance ρ ∈ R, 1 ≪ ρ, we then obtain an interference pattern where E ∞ (k, ρe3) denotes the far field approximation of the reflected incident field.Using the representation of the backscattered field, we obtain where we used the paraxial approximation, see Equation 5.The elements and ∆0 describe the position of the object with respect to the focus and the phase difference between the sample and the reference mirror respectively, see [27].The backscattered light is recorded for each wavenumber k ∈ S independently.Collecting these single measurements, leads to the measurement C as a function of the wavenumber.The interference pattern in Equation 13hereby represents the OCT measurement comprising one A-scan of the object, which is the center of interest in the OCT forward and inverse problem.
We want to assume now that all the parameters of the system are known and that the unknown is the refractive index function n of the medium.With the forward operator we can thus write the inverse problem of OCT as the determination of (nj) J+1 j=1 and (dj) J j=1 from the measurements C via the non-linear equation

The Physical Parameters
So far Equation 16 together with Equation 15 has been formulated as a general inverse scattering problem for interferometric measurement data.We now want to simplify the equations by plugging in the typical values of the parameters in an OCT system and keeping only terms of significant size.In Table 1, we listed the working parameters of the OCT system used for our data acquisition.These are kept fixed for the whole modeling process and also the inverse problem.In particular, we have a rather narrow bandwidth, a small beam width, an almost normal incident angle, a small angle of acceptance, and we measure at a comparatively large distance to the object.This leads to some combinations of and relations between these parameters, which are small compared to one and which we will often simply neglect in the following analysis.

Parameter
Value Table 2. Small quantities in the experimental setup.
Assumption 3 We assume that the following quantities can be considered sufficiently small so that they can be safely neglected in the model.
(i) In general, we assume the angular values θ and θΩ to differ only slightly from zero, which allows us to consider both quantities as small parameters and to keep in the following only terms which are linear with respect to them.
(ii) We assume that the position of the detector is sufficiently far away from the object, which implies that in the relation in Equation 14the distance between the focus and the object is dominated by the distance between the object and the detector, yielding that ψ0 ≈ ρ − xΩ,3, with ρ ≫ xΩ,3.
In particular, we assume that the ratio ka ψ 0 between half the beam width k√ a (measured in multiples of the averaged wave length) and the distance ψ 0 √ a to the detector (in multiples of half the beam width) is neglible.
(iii) Similarly, we assume that the distance ∆0 between the object in the sample arm and the mirror in the reference arm is so large that the difference k2∆0 − k1∆0 in multiples of the wave lengths between measuring it with the maximal wave vector k2 and the minimal k1 is large, so that we can assume (k∆0) −1 to be small.

(iv)
The ratio between the bandwidth 2k of the spectrum and its center k is assumed to be so small that it is enough to keep the zeroth order in k k−1 .
(v) Finally, we assume that the beam width is so small that the deviation k √ a sin(θ) of the tilt measured with respect to the different wave lengths in the specturm is sufficiently small to neglect all terms of higher than linear order therein.
We summarize all these small quantities and their values in our experimental setup in Table 2.

A Layer-by-layer Method for the Inverse Problem
To avoid having to solve for all the 2J + 1 parameters in Equation 16at once, we want to try to split the reconstruction, as it was also done in [8,24,25], for example, into the subproblems where Ij : R × R → L 2 (S) is the jth term in the sum of Equation 15, which corresponds to the contribution from the reflection at the boundary between the layers Ωj−1 and Ωj: Since also the coefficients (n l ) j−1 l=1 and (d l ) j−2 l=1 from the previous interfaces appear in Ij via the combined reflection coefficients r ≤j−1 and the combined phase factor Ψj of the transmissions at the previous interfaces, we want to proceed iteratively by recovering first n1 from Equation 17 for j = 1, and then obtain (nj, dj−1) from the jth problem in Equation 17 after having already recovered (n l ) j−1 l=1 and (d l ) j−2 l=1 from the previous ones.
The difficulty hereby is, however, that we a priori do not have access to the corresponding measurements C j .To get an approximation for these, we perform a Fourier transform of Equation 16with respect to the wave number k (we extend the function from S to R by zero) and obtain with k and k as in Table 1 and where we call the dual variable to the wave number k the optical distance z and write * z for the convolution with respect to z, see Figure 1.This is a combination of sinc-functions (defined by si(z) = sin(z)/z), which are centered at the frequencies of the sine under the integral in Equation 13.Since B is a small disk close to the origin, we get in a very rough approximation so that the distance between the peaks of the sinc-functions corresponds in zeroth order to twice the travel time of the light between the two interfaces.If in addition to the width of the layers, the size of the spectrum S, described by k, is sufficiently large (which is the basis of OCT), these peaks can be nicely separated from each other so that we find intervals Uj around the points ∆j( We will therefore assume that we can recast the inverse problem from Equation 16as the iterative procedure, where we start from the top and then go layer by layer deeper inside the object and recover in the jth step from the knowledge of (n l ) j−1 l=1 , (d l ) j−2 l=1 the parameters (nj, dj−1) from where Ij : R 2 → L 2 (S) is the forward operator for the jth light-layer interface interaction, defined by Equation 17and Uj is a suitably chosen interval around the peaks of the data . We note that the operator I1, corresponding to the reflection at the top boundary, only takes the refractive index as an argument and no distance.

Almost Normal Incidence
We will thus look at a specific step j of the inverse problem Equation 20.To discuss its properties, we want to approximate in Ij the refractive index rj by the expression r † j = n j−1 −n j n j−1 +n j , which corresponds to the refractive index a plane wave with normal incidence on the interface would experience.This yields the simplified forward operator Since the angle θΩ is assumed to be small, we can expand this around the value θΩ = 0, at which we have B = D k sin θ (0) and ψ1 = 0, and get an analytic expression for the leading order term of the forward operator.We will do the calculations for the operator I * 1 corresponding to the first interface, which is a bit simpler than those for deeper interfaces, since also the terms r ≤0 = 1 and Ψ1 = 0, defined in Equation 10, disappear.
By switching to polar coordinates for κ, we can explicitly calculate this integral and find the desired representation.□ Taking the Fourier transform of Equation 22 with respect to the wavenumber k, we obtain in leading order a representation of F k (I * 1 [n1]) as a sum of sinc-functions.
Lemma 5.2 Let θΩ = 0 and for k ∈ S = [k1, k2] let I * 1 be defined as in Equation 21.Then its bandlimited Fourier transform is given as the linear combination where the functions uj : R → C are defined by and we define the parameters γ and ξ as in Equation 23.
Proof: The Fourier transform of Equation 22yields ) Because of Assumption 3 we have that ak is small compared to ψ0 and we linearize the fraction We split this into the linear and the constant term in k and see that the integrals in Equation 26reduce to the two types k j e −ik(z−ϵ∆ 0 ) dk and of integrals for j ∈ {2, 3} and ϵ ∈ {−1, 1}.
• We calculate the first type by writing k j as derivative with respect to the variable z and obtain with k = 1 2 (k1 + k2) and k = 1 2 (k2 − k1) that • For the second integral, we proceed analogously which gives us together with the Fourier convolution theorem dζ.
We substitute ζ = ζ √ γ to rewrite this in the form Using that, according to Assumption 3, the term k √ γ is much smaller than one and the integral is due to the Gaussian factor e − 1 4 ζ2 restricted to a domain of order one, we approximate the sinc-function in the integral by its zeroth order in terms of ζk √ γ.The mean value theorem for Again by Assumption 3, using that k k−1 is small, we only keep terms of zeroth order and obtain Introducing the functions uj as in Equation 25, we obtain for

□
By Assumption 3 we have that k∆0 ≫ 1, which holds true for the experimental data.We can then ignore for z > 0 those sinc terms in Equation 24 which are centered on the negative axis.The function , for z > 0, which reflects the minor influence of the sinc-functions uj which are centered around −∆0 at the function value close to the point z = ∆0.
(ii) In particular, we find for the norm in the highest order Proof: (i) We remark that With this, we find If we neglect herein the fourth and eighth term as they are of the order of the small quantity a k ψ 0 (which is of order 10 −3 in our setting) and plug this into our expression in Equation 27 for F k (I * 1 [n1]χS ), we end up with Equation 28.(ii) We see that the second term in Equation 28 is of lower order compared to the first one: If we thus neglect this term, we obtain Equation 29.□ In order to apply these formulas for j ≥ 2, we notice that the only difference from Equation 22is the term 2n l d l cos θ l t , where the components are calculated iteratively by using the linearization of the square-root for κ ∈ B = D k sin θ (0) making use of the smallness of the angle of acceptance θ: Thus, we can adapt Equation 29 to the contribution I * j from the jth interface by introducing the layer dependent variables for all j = 1, 2, . . ., J, where the coefficients for the first interface coincide with the system parameters: ∆1 = ∆0, ψ0,1 = ψ0, and ξ1 = ξ.Moreover, we use Lemma 2.3 to approximate the term r ≤j−1 (κ), defined in Equation 10, in the integral of I * j by With this, we obtain more generally the formula with the function for the Fourier transform of the forward operator I * j .The last error term O(θ 2 ) herein accounts for the approximation of the combined reflection coefficient r ≤j−1 by its zeroth order.
In particular, we can use this explicit expression to estimate how much the signal F k (I * j χS ) from the jth interface influences the values F k (I * l χS ) from the lth interface in an interval U l around the peak ∆ l of the second interface.We show the estimate for simplicity only for the interaction between the first and the second interface.Lemma 5.4 There exists an interval U2 such that for all z ∈ U2.
Proof: The behavior of the two Fourier transforms is mainly described by the functions U1 and U2 from Equation 33.We therefore estimate the contribution from U1 at a point z = y + ∆2 = y + ∆1 + 2n1d1 around the position ∆2, where U2 has its peak, from above and the contribution of U2 there from below.

□
The above arguments justify the method to reduce the all-at-once-approach of the inverse problem to the simpler layer-by-layer reconstruction presented in Section 4.However, the enclosed form of the integral operators Ij (and I * j respectively) still prevents us from an analytic expression for the reconstruction.Hence, we formulate our problem as a least squares minimization problem.

Least Squares Minimization for the Refractive Index
To numerically obtain a solution of the inverse problem in Equation 20 for a fixed j, we write it as a discrete least squares minimization problem of the functional for the parameters nj and dj−1.
To simplify the analysis, we will, however, replace herein the full forward operator Ij by the reduced forward operator I * j and consider thus the functional Hence, for every interface j we want to minimize J * j (nj, dj−1) with respect to nj and dj−1.Hereby, the data is provided on a discretized grid of Mj points which we denote by {zj,m} M j m=1 ⊂ Uj in the interval Uj selected in Section 4. We silently assume that the distance between the interfaces is sufficiently large so that the influence from the values F k (I l [n l , d l−1 ]χS ) for l > j in the interval Uj is negligible (as estimated in Lemma 5.4).
Finally, we declare the set of admissible values for n and d (in every step) by The reconstruction is based on the isolation of the clearly visible peaks in the data, which originates from a relatively high refractive index contrast between the single layers.If no contrast is present, meaning that nj = nj−1, the method considers the layer as a unit and goes on to the next layer.
Hence, we may exclude -for the step jnj−1, the refractive index from the previously reconstructed layer, from the admissible space Aj.
Calculating both values from this minimization functional is a slight overkill, since the width dj−1 in the jth step, can be found by the largest overlap of the sinc-function in the forward model and the peak in the data independent of their height, which itself is determined by the refractive indices of the different media.Thus, we can separate the extraction of the pair of parameters (n, d) into two parts, similar to how it is done for stepwise-gradient methods.Firstly, the width d of the layer is reconstructed, and secondly the refractive index.We delay the justification of our arguments to Section 8.There we show (visually) that the functionals (in a step j) actually allow a precise prediction of the width dj−1 without the knowledge of the refractive index nj.
We will therefore in the following consider the minimization problems with respect to nj for a known value dj−1.We thus end up with the reduced minimization problem for the functionals defined in Equation 34 and Equation 35.We hereby assume that the value of dj−1 has been reconstructed already and restrict the analysis to the refractive index argument.The second argument of the functional is therefore omitted.
To simplify the notation, we define for j and m (37) so that we can write and We note that by definition Γ * j,m is independent of the refractive index, since I * j carries the directional independent reflection coefficient r † j , see Equation 21.Before we consider the minimization problem, we want to justify that the modeling error introduced by replacing Jj with the simplified functional J * j can be controlled.
Lemma 6.1 For 2 ≤ j ≤ J + 1 let Jj and J * j be defined by Equation 34 and Equation 35, respectively.Then for a fixed value d, the following estimate holds true for any n where θ j−1 t is the (small) angle of transmission (defined in Equation 11) and where we defined the supremum norm ∥h∥ ∞, B = sup κ∈ B |h(κ)| on the reduced domain of integration B, which is defined by κ = κ k for κ ∈ B.
Proof: We first calculate the difference, which using Fj,m = Fj,m − F * j,m + F * j,m is given by By using the definition of the forward model, we find that the essential part is to find an upper bound for .
First, we need to estimate the integral After a change of coordinates κ = kκ, where κ is in the reduced set B, we can replace B by the disk D sin(2θ Ω +θ) in the integral and find an upper bound For any value of n it holds that the supremum norm The difference between the actual function and its approximated form then is estimated by Up.
We combine all the above estimates to find Finally, we find an upper bound for Thus, we get We summarize all coefficients of B in a constant C and obtain the desired result.□ The estimate and its proof for j = 1 is similar

Existence and Uniqueness Results
We finally turn to the question if the minimization problem ( * ) has a unique minimizer.Actually, under the assumption that the data is an element of the range of the forward model, that is for any m we have we can show that there exists a unique solution to the minimization problem ( * ) corresponding to the approximated functional in Equation 35.
For j = 1, we have that Λ1,m = 0, for all m, which reduces the functional J * 1 to This functional, even in the best case, where the data is in the range of the forward model, attains two global maxima.In order to exclude one of the possible solutions, we use the fact that the object is surrounded by air and that the refractive index of the first layer thus has to satisfy n1 > 1.
For j ≥ 2, we obtain the uniqueness from the definition of the functional, which allows -in contrast to j = 1 -the calculation of the reflection coefficient from a second order polynomial equation.For this reason, we discuss the cases j = 1 and j ≥ 2 separately.Theorem 7.1 For j = 1 let J * 1 be defined as in Equation 35and let y1 satisfy the range condition Equation 38 for some ñ1 > 1.Then there exists a unique solution to the minimization problem ( * ).
Proof: We take the derivative with respect to n and find Since we have ∂nr Since we have excluded n = 1 (as refractive index from the previous layer) from the possible solutions, we find that r † 1 (n) = ±r † 1 (ñ1), which further yields .
Since sgn(r † 1 (ñ1)) = −1, due to the assumption that ñ1 > 1, we can also exclude n− from the set of possible solutions.Otherwise it would hold that n < 1.
The non-negativity of the functional J * 1 and the fact that J * 1 (n+) = 0 implies that n+ is actually the global minimum.□ Lemma 7.2 Let n be a minimum of J * 1 , then we have Proof: Let n satisfy ∂nJ * 1 (n) = 0. Then according to the proof of Theorem 7.1 the second derivative of J * 1 with respect to n is given by Thus, together with r † 1 ∂nr † 1 = 1−n (1+n) 3 ̸ = 0, we immediately find that We formulate the uniqueness for the deeper interfaces in the case j = 2.The generalization to general j is straightforward.Proposition 7.3 Let the functional J * 2 be defined as in Equation 35 and let the data y2 satisfy the range condition Equation 38 for some r2 : Then, for a suitable choice of discretization points z2,m, m ∈ {1, . . ., M2} with M2 > 1, the value ñ2 is the unique minimizer of J * 2 : J * 2 (ñ2) = 0 and J * 2 (n) > 0 for all n ̸ = ñ2.
Proof: Using the abbreviations from Equation 37, the non-negative functional J * 2 can be written as Since all terms in this sum are non-negative, it can only vanish if all of them are zero, that is, if we have for all m ∈ {1, . . ., M2} where Γ * 2,m ̸ = 0. Inserting the definitions of Λ2,m and Γ * 2,m , this becomes the condition that the function is constant, at least restricted to the values z2,m.However, we see from the asymptotic expression Equation 28 for F k (I * 1 [n1]χS ) and the correspondingly adapted expression for the second interface (as in Equation 32) that this function (in the leading order terms) is not constant.Therefore, for a suitable choice of discretization points z2,m, we can guarantee that there are at least two values m1, m2 for which , thus excluding any zero values of J * 2 other than r2.□ Remark 7.4 So far we have considered the case where the data is in the range of the forward model without any noise.Noisy data is treated similarly.We note that since the positive functional J * j (also with noise) remains a fourth order polynomial in r † j , the existence of the local minima with respect to n is guaranteed.One is located on the half (1, nj−1), the other one on (nj−1, ∞).The functional J * j (as a continuous function of n) changes also continuously depending on how the data is disturbed.This means, that up to a certain (small) level of noise, the position of the global minimum remains in the correct position.However, cases where the noise is such that the position of the global minimum might jump to the other side cannot be excluded.
Sign determination.A difficulty, which we are facing in the reconstruction of the first layer in Theorem 7.1 (and which has been identified also in the reconstruction method shown in [8]) is that the minimum of the functional can only be determined up to its absolute value.In [8], we had in each step two local (and also global) minima, where we excluded one by using a priori information on the refractive indices of the layers.
Additional noise on the data may lead to the case where the functional (for the refractive reconstruction from layers deeper inside the object) attains equal values for the two local minima.In this case, a distinction of the correct solution to the minimization problem ( * ) is no longer possible.Now we want to show from a very theoretical point of view how one may exclude one of the possible solutions.
In detail, let j ≥ 2 and let the refractive indices and the widths n l , d l−1 , 1 ≤ l ≤ j − 1, be determined and assume that the data in the jth step is within the range of the forward operator I * j for some Hence, rj (or ñj) is the objective to be reconstructed in this step.What we have seen so far is that the minima with respect to n are distributed such that one is located on each side of the point nj−1, the refractive index from the previous layer.
Then, the position of ñj with respect to nj−1, in detail the difference nj−1 − ñj, determines Thus, the correct solution ñj of the minimization problem -in step j -is completely connected to the sign of the reflection coefficient.
If one is able to determine the sign of the reflection coefficient rj from the data, one already the determines the position of the exact solution with respect to nj−1.Therefore, while searching for the minimum of the functional, we may restrict to a certain half, (1, nj−1) or (nj−1, ∞), of the admissible values.
Basically the goal is to read off the sign from the measurement data by avoiding the highly frequent parts originating from the exponential factor e −i k ., see Equation 20.To do so, we multiply the equation with the exponential factor with opposite sign and obtain the sign of rj in the neighbourhood of z = ∆j.We present this method for small values θΩ, where we by using Lemma 5.3 can calculate the Fourier transform of the integral in Equation 21 explicitly in the leading order.For simplicity, we consider in the following the single reflection from the first interface, meaning we consider the case j = 1.
Lemma 7.5 Let the asymptotic behavior as in Lemma 5.3 hold true, meaning that the Fourier transform of the data, F k C 1 χS , is given by Equation 28 for a r † 1 .Further, let ∆1 and ξ1 be defined as in Equation 31.Then we obtain for z + ∆1 > 0. Furthermore, in a sufficiently small neighbourhood of z = 0 we have that Proof: Using the real part of Equation 28, we see that is of lower order compared to − si(kz) + e − k2 γ cos( kξ1) si(k(z + ξ1)).
Hence, we have that Then, in a sufficiently small neighbourhood of z = 0, the sign of the dominant term is determined by the sign of the first sinc-function, meaning that sgn si (kz) − e − k2 γ cos( kξ1) si(k(z + ξ1)) = sgn si (kz) = 1.

□
In the best case possible, where the measurement data (on basis of a plane wave model) is known for all wavenumbers k ∈ R, the position ∆1 is perfectly determined by the duality of the Fourier transform However, for the Gaussian model, see Equation 28, the situation is different.In order to determine ∆1, the position of the peak, from the measurement, the absolute value of Equation 28 is used.There, the sinc-term centered around z = ∆1 − ξ1 overlaps (and interacts) with the one centered around z = ∆1 and therefore eventually shifts the position of the global maximum slightly to ∆1 + µ0 for a µ0 ≪ 1.Even for a small change, say kµ0 > π/2, we then get from Equation 42sgn An increasing bandwidth 2k of wavenumbers, which yields that increases the precision which is needed for determining ∆1.However, for an actual OCT setup, the precision is limited and provides only non-sufficient reconstructions.

Behavior of the Functional with Respect to the Width
In the previous sections, we showed that in every step the minimization problem ( * ) for the functional J * j in Equation 35 attains a unique solution.Hereby, we have assumed that the determination of the width dj−1 (in step j) can be carried out independently of n and therefore can be assumed to be reconstructed before the jth step, j > 1.
In this section we want to justify our argument, that this reconstruction may be separated, by presenting the behavior of the corresponding functionals for different values with respect to d.That is, we want to argue that the functional J ( * ) j (n, d) in the jth step (almost) independently of the refractive index n yields a unique minimum with respect to d.
We have seen that the actual data is nicely predicted by the forward model showing also the sinc-function structure.To this end let us assume that the data y be given and in the range of the forward model, meaning that for j and m let yj,m = Λj,m + F k I * j [ ñj, dj−1 ]χS (zj,m) for some (ñj, dj−1) ∈ Aj and let rj = r † j (ñj) denote the corresponding directional independent reflection coefficient.Clearly, we obtain the best match, meaning that the according minimization functional is zero, if the amplitude, which is determined by the refractive index n, and the width d satisfy (n, d) = (ñj, dj−1).
In the actual situation the refractive index in the jth step is reconstructed after the width of the layer is determined.This means that the reflection coefficient corresponding to the forward model I * j in Equation 35, which further determines the absolute value after Fourier transform, is certainly not matched.We want to show that even though the amplitude is not matched perfectly, the functional (and the corresponding minimization problem with respect to d) still yields a unique minimum at the correct position.
We cannot explicitly calculate the extremal values of the functional J ( * ) j .Instead, we first provide, in Figure 2, for a better understanding a comparison between the simulated data y, for a two layer model, and the forward model.This is shown for different values of d and r † 2 , and the L 2 -error between the data and the model is calculated, see Table 3.
The influence of r † 2 , especially for a large value, on the L 2 -error has a more severe impact compared to the cases where the distance in the model is chosen incorrectly.In order to retrieve the width correctly, we have to use a test functional for suitably chosen refractive index, since the correct is not determined at this stage.In order to minimize the L 2 -error, the test model with respect to r † 2 (and n) should satisfy |r † 2 | < |r2|, which guarantees that the second peak of the model is below the one of the data.
In Figure 3, we show the (approximated) functional for the same two layer medium on a grid of values d.From left to right, the number of grid points M2 for the integration used in the functional in Equation 35 increases.In order to point out the independence of the existence of a (global) minimum on the refractive index, we plotted for each M2,i, i = 1, 2, 3, the functional for different values of r † 2 , satisfying |r † 2 | ≤ |r2|, where |r2| is the upper bound for the peak height in the data.Clearly, for the actual value r † 2 = r2, the functional is zero for all M2,i.By increasing the number of grid points we additionally obtain a sharpening effect which pronounces the minimum even better.Further, we see that for r † 2 with the wrong sign and where the simulated peak is ultimately on the level of the sideloops of the sinc-functions, the brown curve in Figure 3, the minimum is shifted slightly for larger values of M2,i.The location of the global minimum of the functional for each M2,i and every r † 2,l , l ∈ {1, . . ., 4}, is presented in the first three lines of Table 4. There, we see that for an intermediate number of grid points, the localization of the minimum position is not shifted for resonable values of r † 2 .Finally, we consider the case where the data is disturbed  and number of grid points M 2,i .The first three lines correspond to the case where the data satisfies the range condition and which are plotted in Figure 3.The two bottom lines show the minimas of the functional, where an additional term was added to the data.
for some Φ representing possible noise or a reflection from deeper inside the sample.We use the forward model |Λj,m + r † j (n)Γ * j,m | 2 for the functional J * j in Equation 35, and we still obtain a unique global minimum.For the simulation of the data, we used for Φ the forward model (for an additional layer reflection) with r3 = −r2 and distance d2 = d1.In Figure 4, we plot the functional for different values of r † 2 and M2,i ∈ {150, 300}.In Table 4 the two bottom rows show the position of the global minimum for the disturbed data.Here, again, an intermediate value of M2 is the most effective one.Aside from the global minimum, which we obtain in every case above, there exist also local minima.These occur if the peak overlaps with the sideloops of the data-sinc-function and therefore yield a smaller L 2 -error.Hence, for the minimization (process) of the functionals, the initial guess is crucial.Otherwise the solution is stuck in the neighbourhood of one of these local minima.

Numerical Experiments
OCT setup and phantom description.In the following examples we illustrate the applicability of the proposed method for both simulated and experimental data.The OCT experimental data, shown Figure 5, was generated by a swept-source OCT system with a laser light source centered at the wavelength λ0 = 1300nm.The bandwidth around the central wavelength ranges from approximatively 1282.86nm to 1313.76nm.During a single experiment the object was raster scanned on a grid of 1024 × 1024 points in horizontal direction and for each raster scan a spectrum consisting of 1498 data points, equally spaced with respect to the wavenumber, was recorded.Due to the relatively small bandwidth of the laser, the approximation that n does not dependent on the wave number is valid.For the reconstruction, a single depth profile was chosen from the data set and no averaging has been done.Additionally, for the calibration of all necessary system parameters, introduced in [27], the object was shifted multiple times along the axial direction, where for each position a full measurement was recorded.
In order to examine the correctness of the recovered values, the object needs to consist of materials where the ground truths of the refractive indices, on the used spectrum, are available.The sample is a three-layer medium where two coverglasses enclose water.The refractive indices of both materials, at the central wavelength, are given by ng ≈ 1.5088 and nw ≈ 1.3225.The thickness of every layer is estimated between 0.14mm and 0.19mm.
The simulated data were generated for the same three-layer sample with refractive indices n0 = 1, n1 = 1.5088, n2 = 1.3225, n3 = 1.5088 and widths d = [0.174,0.186, 0.173]mm by using the forward model presented in Section 2. Given the list of system and object parameters described in Assumption 3, the formula for the reflected field Equation 8 was implemented where the domain of accepted wave directions, see Equation 12, was replaced by its discretized version.
In addition we use a sample showing characteristica in lateral directions, which are then visible in the outcome for different raster scans across the object.On a grid of 20 × 20, we define the object for every (l, m) ∈ {1, . . ., 20} 2 as a two layer object with n(l, m) = [1, 1.5088, n(l, m)], where the distribution n is shown in Figure 12.For every grid element, we collected the simulated measurement data.For the reconstruction 5% uniformly distributed noise was added to the data.
Multiple reflections, due to their minor contribution to the final measurement, have been omitted.
Inverse problem and minimization.For the inverse problem, the functionals Jj and J * j , based on the direct model, defined in Equation 34 and Equation 35 respectively, were implemented.The Simulated vs. Real Data In Section 7, we discussed the minimization problem for the data yj, 1 ≤ j ≤ J + 1, which so far, has been considered as an element of the range of forward operator (plus a certain noise level δ).However, experimental data fails to satisfy this relation because of the unknown intensity in the focal plane: In Assumption 1 we have set the amplitude of the Gaussian distribution in our forward model equal to one, which further yields that the maximal intensity of the incident field in Equation 3, that is E (0) 2 in the focal spot x3 = r0, is also one.For the experimental data the maximal intensity represents the power of the laser light in the focal spot and is in general an unknown quantity.We consider this quantity as real positive number Q0.Even for a single layer reflection, the factor Q0 causes incorrect reconstructions if it is not included in the model.Let the data be given by ỹ = Q0y, for y = |r † (ñ1)Γ * 1,m 0 | 2 , for a single grid point z1,m 0 .Then from Theorem 7.1, we deduce that |r † (n)Γ * 1,m 0 | 2 = Q0|r † (ñ1)Γ * 1,m 0 | 2 , which then yields |r † (n)| = √ Q0|r † (ñ1)|.We correct the mismatch between the model and the experimental data by the quantifying Q0 from the first air-glass reflection.Hereby, we use the coverglass plate as a calibrational layer.We recover instead of the refractive index n1, which we assume in this case to be determined perfectly before, the quantity Q0 as the quotient between the data and the forward model for j = 1 in Equation 34.

Minimization in Two
Steps: In Section 4 (which is carried out in Section 7) we proposed a layerby-layer method where in each step a single pair of parameters (nj, dj−1), 1 ≤ j ≤ J + 1, is recovered.Each step itself is then again split into two steps.After the minimum of the functional J * j is determined, the output is used as an input in form of an initial guess for the actual minimization problem corresponding to the functional Jj.In Figure 6 and Figure 7 a comparison of the two functionals (for the minimum value of d) for every reconstruction step is provided for simulated and experimental data, respectively.
Results.The outcome of both minimization problems, for simulated and experimental data, are presented in this section.In Figure 8 we plot the functionals Jj, j = 1, 2, 3, on a (n, d)-grid on the half of the admissible values Aj where the global minimum is located.
For the experimental data, the functionals evaluated over a (n, d)-grid are shown in Figure 9. Additionally, a contour plot of the close neighbourhoods of each pair of minima are shown in Figure 10.The reconstructed values are plotted in Figure 11.The laterally reconstructed refractive index distribution is shown in Figure 12.Comparison between the functionals J j (blue line) and J * j (red circles).The x-axis shows possible values of the refractive index n j in each step j = 1 to j = 3 for the reconstruction for simulated data.The order is from left to right.

Conclusion
The inverse problem in quantitative optical coherence tomography for layered media has been discussed.Imaging", Lisa Krainz and Wolfgang Drexler have been supported via the subproject F6803-N36 "Multi-Modal Imaging".

Figure 1 .
Figure 1.A single A-scan in wavenumber domain showing high-frequent components (left).The absolute value of the dataset after Fourier transform avoiding these high-frequency components in the data (right).

Figure 5 .
Figure 5.The OCT experimental data for a glass-water-glass sample produced by a 1300nm system.A cross-sectional image, a B-scan, (left) showing 150 depth profiles (A-scans) of the object.A single depth scan (right), which is taken along the red line of the B-scan, shows the single layers of the object.The intensity in the A-scan is normalized such that the first peak has magnitude one.

Figure 6 .
Figure 6.Comparison between the functionals J j (blue line) and J * j (red circles).The x-axis shows possible values of the refractive index n j in each step j = 1 to j = 3 for the reconstruction for simulated data.The order is from left to right.

Figure 7 .Figure 8 .
Figure 7.Comparison between the functionals J j (blue line) and J * j (red circles).The x-axis shows possible values of the refractive index n j in each step j = 1 to j = 3 for the reconstruction for experimental data.The order is from left to right.
Based on a Gaussian beam forward model a discrete least squares minimization problem for the reconstruction of the refractive index and the thickness of each layer is formulated.Existence and uniqueness of solutions for the minimization problem are discussed.The method is validated by numerical examples considering the refractive index reconstruction from a three-layer object from both simulated and experimental data.

Figure 9 .
Figure 9.The functionals J j (for experimental data) for each reconstruction step evaluated on a grid of 150 × 100 values for the refractive index n and the width d.

Figure 10 .
Figure 10.A contour plot of the functionals J j (for experimental data) for each reconstruction step evaluated within a close neighbourhood of each pair of minima (n, d).

Figure 11 .
Figure 11.A comparison of the reconstructed values (red) with the ground truth (blue).The reconstruction for the simulated data without noise (left) and the experimental data (right) match nicely the ground truth.

Figure 12 .
Figure 12.The ground truth refractive index distribution n on a 20 × 20 grid, showing the discrete versions of two circles nC 1 = 1.45, nC 2 = 1.4 and an ellipse nE = 1.1 surrounded by a homogeneous medium with refractive index nB = 1.37 (left).The reconstructed distribution from simulated data with 5% (random) noise (right).

Table 1 .
Parameters of the experimental setup of the considered OCT system.

Table 3 .
The different values used for the plot in Figure2and the calculated L 2 -error.The first line in the table showing the actual values used for the simulated of y.
The simulated data y for two reflections in blue is compared with the model for different values of d and r † 2 : for a relatively large value r † 2 and small distance d (brown), the actual distance d and small r † 2 (red) and where both values are off (black).

Table 4 .
The minimum values of the functional J * 2 for different reflection coefficients r † 2