[t]11em$^1$Faculty of Mathematics

University of Vienna

Oskar-Morgenstern-Platz 1

A-1090 Vienna, Austria

[t]11em$^2$ Project-team Makutu

Inria Bordeaux Sud-Ouest,

E2S–UPPA, UMR CNRS 5142,

[t]11em$^3$Institute of Mathematics

Technical University Berlin

Straße des 17. Juni 136

D-10623 Berlin, Germany

[t]15em$^4$Johann Radon Institute for Computa-

tional and Applied Mathematics

Altenbergerstraße 69

A-4040 Linz, Austria

In this paper, we study the mathematical imaging problem
of diffraction tomography (DT), which is an inverse scattering
technique used to find material properties of an object by
illuminating it with probing waves and recording the scattered waves.
Conventional DT relies on the Fourier diffraction theorem, which is
applicable under the condition of weak scattering.
However, if the object has high contrasts or is too large compared
to the wavelength, it tends to produce multiple scattering, which
complicates the reconstruction.
We give a survey on diffraction tomography and compare
the reconstruction of low and high contrast objects.
We also implement and compare the reconstruction using the full
waveform inversion method which, contrary to the Born and Rytov
approximations, works with the total field and is more robust to multiple scattering.

§ INTRODUCTION

Diffraction tomography (DT) is a technique for reconstructing the scattering potential of an object from measurements of waves scattered by that object.
DT can be understood as an alternative to, or extension of, classical computerized tomography. In computerized tomography a crucial assumption is that the radiation, X-rays for instance, essentially propagates along straight lines through the object. The attenuated rays are recorded and can be related to material properties $f$ of the object by means of the Radon, or X-ray, transform. A central result for the inversion of this relation is the Fourier slice theorem. Roughly speaking, it says that the Fourier transformed measurements are equal to the Fourier transform of $f$ evaluated along slices through the origin, <cit.>.

The straight ray assumption of computerized tomography can be considered valid as long as the wavelength of the incident field is much smaller than the size of the relevant details in the object. As soon as the wavelength is similar to or greater than those details, for instance in situations where X-rays are replaced by visible light, diffraction effects are no longer negligible.
As an example of a medical application, an optical diffraction experiment in <cit.> utilized a red laser of wavelength 633nm to illuminate human cells of diameter around 10, which include smaller subcellular organelles.
One way to achieve better reconstruction quality in such cases is to drop the straight ray assumption and adopt a propagation model based on the wave equation instead.

The theoretical groundwork for DT was laid more than half a century ago <cit.>. The central result derived there, sometimes called the Fourier diffraction theorem, says that the Fourier transformed measurements of the scattered wave are equal to the Fourier transform of the scattering potential evaluated along a hemisphere. This result relies on a series of assumptions: (i) the object is immersed in a homogeneous background, (ii) the incident field is a monochromatic plane wave, (iii) the scattered wave is measured on a plane in $\R^3$, and (iv) the first Born approximation of the scattered field is valid.

On the one hand the Born approximation greatly simplifies the relationship between scattered wave and scattering potential. On the other hand, however, it generally requires the object to be weakly scattering, thus limiting the applicability of the Fourier diffraction theorem. An alternative is to assume validity of the first Rytov approximation instead <cit.>. While mathematically this amounts to essentially the same reconstruction problem, the underlying physical assumptions are not identical to those of the Born approximation, leading to a different range of applicability in general <cit.>. Nevertheless, the restriction to weakly scattering objects remains.

Full waveform inversion (FWI) is a different approach that can overcome
some of the limitations of the first-order methods, typically at the cost of being
computationally more demanding.
It relies on the iterative minimization of a cost functional which penalizes the
misfit between measurements and forward simulations of the total field,
cf. <cit.>.
Here, the forward model consists of the solution of the full wave equation,
without simplification of first-order approximations.
It results in a nonlinear minimization problem to be solved, typically with
Newton-type methods <cit.>.

In practical experiments, there are sometimes only measurements of the intensity, i.e., the absolute value of the complex-valued wave, available.
Different phase retrieval methods were investigated, e.g., in <cit.>.
For this paper, we assume that both the phase and amplitude information is present, which can be achieved by interferometry, cf. <cit.>.

§.§ Contribution and outline

In this paper, we present a numerical comparison of three reconstruction approaches for diffraction tomography on simulated data, based on (i) the Born approximation, (ii) the Rytov approximation, and (iii) FWI. The setting we use for this comparison is 2D transmission imaging in a homogeneous background with (approximate) plane wave irradiation. The object is assumed to make a full turn during the experiment, providing measurements for a uniform set of incidence angles. In addition we investigate how providing additional data by varying the wavelength affects the reconstruction. The scattering potentials considered here are test phantoms of varying sizes, shapes, and contrasts. Moreover, for data generation purposes we compare several forward models.

For numerical reconstruction under the Born and Rytov approximations a well-known method is the backpropagation algorithm <cit.>, which is widely used in practice, cf. <cit.> and also <cit.>.
Our algorithms rely on the nonuniform discrete Fourier transform (NDFT), which was used in 3D Fourier diffraction tomography yielding better results than discrete backpropagation <cit.>.
Our FWI-based reconstruction uses an iterative Newton-type method on
an $L^2$ distance between data and simulations. Here, the discretization
of the partial differential equations associated with the wave propagation
uses the hybridizable discontinuous Galerkin method (HDG), <cit.>.
It is implemented, together with the inverse procedure, in the open-source parallel software [<https://ffaucher.gitlab.io/hawen-website/> ],

The outline of this paper is as follows. The conceptual experiment is detailed in <ref>. Forward models are presented in <ref>, and their numerical performance is compared in <ref>. The Fourier diffraction theorem is formulated and discussed in <ref>. Further, <ref> covers the reconstruction algorithms used for the numerical experiments, which are presented in <ref>. A concluding discussion of our findings is given in <ref>.

§ EXPERIMENTAL SETUP

We consider the tomographic reconstruction of a two-dimensional object taking into account diffraction of the incident field. The object is assumed to be embedded in a homogeneous background and illuminated or insonified by a monochromatic plane wave.
In fact, for the computational experiments, we implement and compare
several approaches to approximate the incident plane wave,
see <ref> and <ref>.
We restrict ourselves to transmission imaging, where the incident field propagates in direction $\be_2=(0,1)^\top$ and the resulting field is measured on the line $x_2 = \rM$. The distance between the measurement line and the origin, $\rM>0$, is sufficiently large so that it does not intersect the object.
From the measurements, we aim to reconstruct the object's scattering properties. In order to improve the reconstruction quality, we generate additional data by rotating the object or changing the incident field's wavelength.
See <ref> for an illustration of the experimental setup.

[->] (-5,0)–(5,0) node[right] $x_1$;
[->] (0,-5.5)–(0,4.5) node[right] $x_2$;

[red,opacity=0.3] (0,0) circle (2cm and 2cm);
[red,opacity=0.8] (-.8,.8) circle (.4);
[red,opacity=0.3,rotate=45] (0,0) circle (1.2cm and .3cm);
[red,opacity=1] (0.5,-0.5) circle (.15);
[red,opacity=1] (1,-0.5) circle (.15);
[red,opacity=1] (0.5,-1) circle (.15);
[red,opacity=1] (1.0,-1) circle (.15);
at (3,-2) object;
[->]  (0.43, 2.46) arc (80:10:2.5);

[dashed] (-5,-3.5) – (5,-3.5);
[dashed] (-5,-4) – (5,-4);
[dashed] (-5,-4.5) – (5,-4.5);
at (-8,-4) incident field;
[->] (5.5,-4.5) – (5.5,-3.5);
at (6,-4) $\be_2$;

[line width = 3pt] (5,3.5) – (-5,3.5) node [left] measurement line;
at (7,3.5) $x_2 = \rM$;

Experimental setup.

We now introduce the physical quantities needed subsequently. Let $\lambda>0$ denote the wavelength of the incident wave and $k_0={2\pi}/{\lambda}$ the background wave number. Furthermore, let $n(\bx)$ denote the refractive index at position $\bx \in \R^2$ and $n_0$ the constant refractive index of the background.
From these quantities, we define the wave number

\begin{equation*} %\label{eq:refractive}
k(\bx)\coloneqq k_0 \frac{n(\bx)}{n_0}.
\end{equation*}

Furthermore, the wave number $k$ can be equivalently expressed in terms of the
angular frequency $\omega$ and the wave speed $c$ such that

\begin{equation} \label{eq:refractive:omega-c}
k(\bx) \,=\, \dfrac{\omega}{c(\bx)} \, \qquad \text{and} \qquad
k_0 \,=\, \dfrac{\omega}{c_0} \,,
\end{equation}

where $c_0$ is the constant wave speed in the background.
The scattering potential $f$ is obtained by
subtracting the background wave number $k_0$:

\begin{equation}\label{eq:f}
f(\bx)\coloneqq k^2(\bx)-k_0^2
= k_0^2 \left( \frac{n(\bx)^2}{n_0^2}-1\right) \, .
\end{equation}

Note that, for all practical purposes, $f$ can be assumed to be
bounded and compactly supported in the disk
$\mathcal B_{\rM} = \{\bx\in\R^2: \norm{\bx}<\rM\}$.

In our subsequent reconstructions with Born and Rytov approximations,
$f$ is the quantity to be reconstructed from the measured data and $k_0$ is
known. On the other hand, with FWI, we reconstruct $c$, see <ref>.
These two quantities can be related to each other via

\begin{equation} \label{eq:fwi:c_and_f}
c(\bx) \,=\, \sqrt{\dfrac{\omega^2}{k_0^2 \,+\, f(\bx)}} \, .
\end{equation}

§ FORWARD MODELS

In this section, we propose several forward models for the
experiment presented above. For all of them, the starting point
is the system of equations

\begin{equation} \label{eq:start}
\left. %\lbrace
\begin{gathered}
\big(-\Delta \,-\, k(\bx)^2 \big) \, \utot(\bx) \,=\, g(\bx), \\
\big(-\Delta \,-\, k_0^2 \big) \, \ui(\bx) \,=\, g(\bx), \\
\utot(\bx)\,=\,\ui(\bx)\,\,+\,\,\us(\bx),
\end{gathered}
\right\} \quad \bx \in \R^2.
\end{equation}

Here, $\ui$ is the given incident field, the total field $\utot$ is what is recorded on the measurement line $\{\bx\in\R^2: x_2=\rM\}$, and the difference between the two constitutes the scattered field $\us$. We describe different sources $g$ in the following subsections. The scattered field $\us$ is assumed to satisfy the Sommerfeld radiation condition which requires that

\begin{equation*}
\lim_{\norm{\bx} \to \infty} \sqrt{\norm{ \bx} } \left( \frac{\partial\us}{\partial \norm{ \bx} } - \i k_0 \us \right) = 0
\end{equation*}

uniformly for all directions $\bx/\norm{\bx}$. It guarantees that $\us$ is an outgoing wave. Further details concerning derivation and analytical properties of problems like <ref> can be found, for instance, in <cit.>.

The models considered below are based on the following specifications of <ref>. Their numbers agree with the corresponding subsection numbers, where the models will be explained in more detail.

1. Plane wave No source ($g=0$) and $\ui$ is an ideal plane wave, see <ref>.

1.1 Born model No source, $\ui$ is an ideal plane wave and $\us$ satisfies the Born approximation.

1.2 Rytov model No source, $\ui$ is an ideal plane wave and $\us$ satisfies the Rytov approximation.

In addition we propose the following two models with sources.

2.1 Point source $g$ represents a point source located far from the object.

2.2 Line source $g$ represents simultaneous point sources positioned along
a straight line. We refer to this configuration as a `line source'.

<ref> contains a numerical comparison of these forward models.

The proposed selection of equations is motivated in part by the availability of methods for their numerical inversion. While the Born and Rytov models can be inverted using nonuniform Fourier methods, the point and line source models are well-suited for FWI.

§.§ Incident plane wave

Monochromatic plane waves are basic solutions $u$ of the homogeneous Helmholtz equation

\begin{equation*}
\big(-\Delta \,-\, k_0^2 \big) \, u \,=\, 0.
\end{equation*}

They take the form $u(\bx) = \e^{\i k_0 \bx \cdot \bs}$, where the unit vector $\bs$ specifies the direction of propagation of $u$. Plane waves are widely studied in imaging applications and theory and we refer to <cit.> for further information.

In the first model we consider the incident field is a monochromatic plane wave propagating in direction $\be_2$

\begin{equation}\label{eq:planewave}
\ui(\bx) = \e^{\i k_0 x_2}.
\end{equation}

In this case, we obtain from <ref> the following equation for the scattered field

\begin{equation} \label{eq:uscat-rhs-uinc}
\left(-\Delta \,-\, k(\bx)^2\right) \, \us(\bx) \,=\, f(\bx) \, \e^{\i k_0 x_2} \, .
\end{equation}

§.§.§ The Born approximation

<ref> can be written as

\begin{equation*}
\left(-\Delta \,-\, k_0^2\right) \, \us(\bx) \,=\, f(\bx) \, \left( \e^{\i k_0 x_2} + \us(\bx) \right).
\end{equation*}

If the scattered field $\us$ is negligible compared to the incident field $\e^{\i k_0 x_2}$,
we can ignore $\us$ on the right-hand side and obtain

\begin{equation} \label{eq:uborn-rhs-uinc}
(-\Delta \,-\, k_0^2) \, \uborn(\bx) \,=\, f(\bx) \, \e^{\i k_0 x_2} \,,
\end{equation}

where $\uborn$ is the (first-order) Born approximation to the scattered field.
Supplementing this equation with the Sommerfeld radiation condition we have a unique solution
corresponding to an outgoing wave <cit.>. It can be written as a convolution

\begin{equation}\label{eq:conv-2d}
\uborn(\bx) =
%[G * f\ui] (\bx) \coloneqq
\int_{\R^2} G(\bx-\by) \, f(\by)\,\e^{\i k_0 y_2} \dd\by,
\end{equation}

where $G$ is the outgoing Green's function for the Helmholtz equation.
In $\R^2$, it is given by

\begin{equation} \label{eq:G2}
= \frac{\i}{4} H^{(1)}_0(k_0 \norm{\bx})
,\qquad \bx\in\R^2\setminus\{\bo\},
\end{equation}

where $H^{(1)}_0$ denotes the zeroth order Hankel function of the first kind, see <cit.>. We note that, in spite of a singularity at the origin, $G$ is locally integrable in $\R^2$.

The second-order Born approximation can be obtained by replacing the plane wave $\e^{\i k_0 y_2}$ in <ref> by the sum $\e^{\i k_0 y_2} + \uborn(\by)$. Iterating this procedure yields Born approximations of arbitrary order. For more details we refer to <cit.> and <cit.>.

§.§.§ The Rytov approximation

In this subsection, we derive an alternative approximation for the scattered field.

Introducing the complex phases $\phii$, $\phitot$ and $\phis$ according to

\begin{equation} \label{eq:exp}
\utot= \e^{\phitot}, \quad \ui = \e^{\phii}, \quad  \phitot= \phii + \phis,
\end{equation}

one can derive from <ref>, with $g=0$, the following relation

\begin{equation} \label{eq:totphase}
(-\Delta -k_0^2)(\ui \phis)=\left(f + \left(\nabla\phis\right)^2\right)\ui,
\end{equation}

where $\left(\nabla\phis\right)^2 = \left( \partial \phis / \partial x_1 \right)^2 + \left( \partial \phis / \partial x_2 \right)^2$.
The details of this derivation can be found, for instance, in <cit.>.
Neglecting $\left(\nabla\phis\right)^2$ in <ref> we obtain

\begin{equation} \label{eq:Rytov}
(-\Delta -k_0^2)(\ui \phir)=f\ui,
\end{equation}

where $\phir$ is the Rytov approximation to $\phis$. Note that we still assume $\ui$ to be a monochromatic plane wave, as given in <ref>. Thus the product $\ui \phir$ solves the same equation as $\uborn$.
If we define the Rytov approximation to the scattered field, $\urytov$, in analogy to <ref> via

\begin{equation*}
\urytov = \e^{\phir + \phii} - \ui,
\end{equation*}

and replace $\phir$ by $\frac{\uborn}{\ui}$, we obtain a relation between the two approximate scattered fields that can be expressed as

\begin{equation}\label{eq:Born-Rytov}
\uborn
%= \ui\phir
%= \ui (\phitot-\phii)
%= \ui\, \log \frac{\utot}{\ui}
= \ui\, \log \left( \frac{\urytov}{\ui}+1 \right).
\end{equation}

The relation between Born and Rytov in <ref>
is not unique because of the multiple branches of the complex logarithm.
In practical computations, this is addressed by a phase unwrapping as we will see in <ref>.

There have been many investigations about the validity of the Born and Rytov approximations, see, e.g., <cit.> or <cit.>.
The Born approximation is reasonable only for a relatively (to the wavelength) small object.
In particular, for a homogeneous cylinder of radius $a$, the Born approximation is valid if $a (n-n_0) < \lambda/4$, where $\lambda$ is the wavelength of the incident wave and $n$ is the constant refractive index inside the object.
In contrast, the Rytov approximation only requires that
$ n-n_0 > (\nabla \phis)^2 / k_0^2$,
i.e., the phase change of the scattered phase $\phis$, see <ref>, is small over one wavelength, but it has no direct requirements on the object size and is therefore applicable for a larger class of objects.
The latter is also observed in numerical simulations in <cit.>.
However, for objects that are small and have a low contrast $n-n_0$, the Born and Rytov approximation produce approximately the same results.

§.§ Modeling the total field using line and point sources

As an alternative to ideal incident plane waves we consider models with one or several point sources.
If arranged the right way, the resulting incident field can resemble a monochromatic plane wave. We refer to <ref> for a numerical comparison of the different models presented here.

§.§.§ Point source far from object

In this case the right-hand side in <ref> is a Dirac delta function so that we obtain

\begin{equation} \label{eq:fwi:helmholtz-point}
\left\lbrace \begin{aligned}
\big(-\Delta \,-\, k(\bx)^2 \big) \, \utot_{\mathrm{P}}(\bx) \,=\, \delta(\bx -\bx_0)\, , \\
\big(-\Delta \,-\, k_0^2 \big) \, \ui_{\mathrm{P}}(\bx) \,=\, \delta(\bx -\bx_0) \, .
\end{aligned} \right. \end{equation}

If the position of the source is given by $\bx_0 = -r_0\be_2$ with $r_0>0$ sufficiently large, then, after appropriate rescaling, $\ui_{\mathrm{P}}$ approximates a plane wave with wave number $k_0$ and propagation direction $\be_2$ in a neighborhood of $\bo$.

§.§.§ Line source

Alternatively, we let $g$ be a sum of Dirac functions and consider

\begin{equation} \label{eq:fwi:helmholtz-full}
\left\lbrace \begin{aligned}
\big(-\Delta \,-\, k(\bx)^2 \big) \, \utot_{\mathrm{L}}(\bx) \,=\, \sum_{j=1}^{N_{\mathrm{sim}}}
\delta(\bx-\bx_j)\, , \\
\big(-\Delta \,-\, k_0^2 \big) \, \ui_{\mathrm{L}}(\bx) \,=\, \sum_{j=1}^{N_{\mathrm{sim}}}
\delta(\bx-\bx_j) \, ,
\end{aligned} \right. \end{equation}

where the number $N_{\mathrm{sim}}$ of simultaneous point sources should be sufficiently large. Moreover, the positions $\bx_j$ should be arranged uniformly along a line perpendicular to the propagation direction $\be_2$ of the plane wave.
This is illustrated in <ref>.

§ NUMERICAL COMPARISON OF FORWARD MODELS

In this section, we numerically compare
the forward models presented above.
For the discretization of partial differential
equations, several approaches exist, we mention for
instance the finite differences that approximates
the problem on a nodal grid (e.g., <cit.>),
or methods that use the variational formulation, such as
finite elements <cit.> or discontinuous Galerkin
methods, <cit.>.
In our work, we use the hybridizable discontinuous Galerkin
method (HDG) for the discretization and refer to
<cit.> for more details.
The implementation precisely follows the steps prescribed
in <cit.>, and it is carried out in the
open-source parallel software ,
see <cit.> and <ref>.
While the propagation is assumed on infinite space,
the numerical simulations are performed on a finite
discretization domain $\Omega\subset\R^2$, with absorbing boundary
conditions <cit.> implemented to simulate
It corresponds to the following Robin-type condition applied
on the boundary $\Gamma$ of the discretization domain $\Omega$:

\begin{equation} \label{eq:abc}
-\mathrm{i} \, k(\bx) \, u(\bx) \,+\ \partial_n u(\bx) \,=\,0, \qquad
\text{for $x$ on $\Gamma$.}
\end{equation}

where $\partial_n u$ denotes the normal derivative of $u$. The test sample used below is a homogeneous medium encompassing a circular object of radius $4.5$ around the origin with contrast $f=1$.
This corresponds to the characteristic function

\begin{equation}\label{eq:f0}
\mathbf 1^{\mathrm{disk}}_a(\bx) := % \mathbf 1_{\mathcal{B}_{a}}(\bx),
%  \qquad \bx\in\R^2,
\begin{cases}
1, & \bx\in\mathcal B_a,\\
0, & \bx\in\R^2\setminus\mathcal B_a,
\end{cases}
\end{equation}

of the disk $\mathcal B_a$ with radius $a>0$.
The incident plane wave has wave number $k_0 = 2\pi$.

§.§ Modeling the scattered field assuming incident plane waves

We consider the solutions $\us$ of <ref> and $\uborn$ of <ref>,
both satisfying boundary condition <ref>, and simulated following the HDG discretization
indicated above.
As an alternative for computing the Born approximation $\uborn$, we discretize the convolution <ref>
with the Green's function $G$ given in <ref>.
In particular, applying an $N\times N$ quadrature on the uniform grid $\mathcal R_N = \{ -\rM, -\rM + 2\rM/N, \dots,\rM-2\rM/N \}^2$ to <ref>, we obtain

\begin{equation} \label{eq:uborn-conv}
\uborn(\bx)
\approx
u^{\mathrm{Born}}_{\mathrm{conv},N} (\bx) \coloneqq
\left(\frac{2\rM}{N}\right)^2
\sum_{\by\in\mathcal R_N} G(\bx-\by) \, f(\by)\,\e^{\i k_0 y_2}
\end{equation}

In <ref>, we illustrate the solutions
obtained with the different formulations.
We observe that the transmission waves contain the most energy, that is,
waves that “cross” the object. On the other hand,
the solution has very low amplitude elsewhere, including the reflected
In <ref>, we see that the Born
approximation leads to an incorrect amplitude of the solution, in
particular the imaginary part.
In addition, the imaginary part of $u^{\mathrm{Born}}_{\mathrm{conv},200}$
does not match the one of $u^{\mathrm{Born}}$.

-1   1
Perturbation model $f=\mathbf 1_{4.5}^{\mathrm{disk}}$.

-0.50  0.50

Real part of $\us$.
-0.50  0.50

Real part of $\uborn$.

1  1

enlarge y limits=true,
enlarge x limits=false,
scale only axis,
ylabel style = yshift =0mm, xshift=0mm,
xlabel style = yshift =0mm, xshift=0mm,
clip mode=individual,
legend columns=1,
label style=font=,
tick label style=font=,
legend style=font=,
legend pos=south east,

[color=black,line width=2]
table[x expr = xrcv,
y expr = ,

[color=blue!80!black,line width=2]
table[x expr = xrcv,
y expr = ]

[color=red,line width=1.50,dotted]
table[x expr = xrcv,
y expr = *]

[color=green!70!black,line width=1.50,dotted]
table[x expr = xrcv,
y expr = *]

[color=purple,line width=1.5,dotted]
table[x index=0,y index=1]images/fwi/modeling_2d_Uscat/data_rcv/f1-1000mHz-u.txt; $\mathrm{Re}(\uborn_{\mathrm{conv},200})$

[color=cyan!80!black,line width=1.5,dotted]
table[x index=0,y index=2]images/fwi/modeling_2d_Uscat/data_rcv/f1-1000mHz-u.txt; $\mathrm{Im}(\uborn_{\mathrm{conv},200})$

Comparison of the solutions at fixed height $z=$10.

Comparison of the scattered wave $\us$ and the Born approximation $\uborn$.
The computations are performed on a
domain $[-25,25]\times[-25,25]$
with boundary conditions given in
Further, we display $\uborn_{\mathrm{conv},200}$ which is the Born approximation obtained by the convolution <ref>.

§.§ Modeling the total field using line and point sources

The objective is to evaluate how considering line and point sources
differs from using incident plane waves, and if the data
obtained with both approaches are comparable.
To compare the scattered fields obtained from <ref> and <ref> with the solution $\us$ of
<ref>, one needs to rescale according to

\begin{align*} %\label{eq:scaling_uscat}
%\ushat \,=\, \tilde{\alpha} \big( \,\uttilda - \uitilda \,\big) \,=\, \tilde{\alpha}\, \ustilda \,.
\us_{\mathrm{P}}		&=	\alpha_{\mathrm{P}} \left( \utot_{\mathrm{P}} - \ui_{\mathrm{P}} \right), \\
\us_{\mathrm{L}}		&=	\alpha_{\mathrm{L}} \left( \utot_{\mathrm{L}} - \ui_{\mathrm{L}} \right),
\end{align*}

where $\alpha_{\mathrm{P}}$ and $\alpha_{\mathrm{L}}$ are constants depending only on the number and positions of the point sources $\bx_j$.
We illustrate in <ref>, where
the line source is positioned at fixed height $x_2=-15$ and composed
of 441 points between $x_1=-22$ and $x_1=22$.
For the case of a point source, we have to consider a very wide
domain, namely $[-500,500]\times[-500,500]$, and the point source is
positioned in $(x_1=0,x_2=-480)$.
In <ref>, we plot the
corresponding solutions on a line at height $x_2=10$, i.e.,
for measurements of transmission waves.

-1   1

Computational domain $[-25,25]^2$
for line source with perturbation model $f_{4.5}^{\mathrm{disk}}$.

-0.01  0.01

Real part of $\utot_{\mathrm{L}}$.
-0.005  0.005

Real part of $\us_{\mathrm{L}}$.

-1   1

Computational domain $[-500,500]^2$
for point source with perturbation model
-0.03  0.03

Real part of $\utot_{\mathrm{P}}$ near origin.
-0.01  0.01

Real part of $\us_{\mathrm{P}}$ near origin.

Comparison of the solutions at fixed height $x_2=$10.

Total field $\utot_{\mathrm{L}}$ and scattered field $\us_{\mathrm{L}}$ (line source) and
total field $\utot_{\mathrm{P}}$ and scattered field $\us_{\mathrm{P}}$ (point source).
The line source is composed of $N_{\mathrm{sim}} = 441$
points at fixed height $x_2=-15$.
The computational domain for the point source is
very large such that the perturbation is barely
visible and the source is positioned in $\bx_0 = (0,-480)^\top$.

We see that the simulation using the line source is very
close to the original solution $\us$, in fact, it is a more
accurate representation than the Born approximation pictured
in <ref>.
The simulation using a point source positioned far away is
also accurate, except for the middle area of the imaginary
Furthermore, the major drawback of using a single point source
is that it necessitates a very big domain, hence largely
increasing the computational cost.

§ FOURIER DIFFRACTION THEOREM

In this section we discuss the inverse problem of recovering the scattering potential from measurements of the scattered wave under the Born or Rytov approximations.
Before stating the fundamental result in this context, see <ref> below, we have to introduce further notation.

We denote by $\ktran$ the two-dimensional Fourier transform and by $\ktran_1$ the partial Fourier transform with respect to the first coordinate,

\begin{align*}
\mathcal{F} \phi(\bk)	&=  (2\pi)^{-1} \int_{\R^2} \phi(\bx) \e^{-\i \bx \cdot \bk} \dd \bx, \qquad \bk \in \R^2, \\
\mathcal{F}_{1} \phi(k_1,x_2)	&=  (2\pi)^{-\frac{1}{2}} \int_{\R} \phi(\bx) \e^{-\i k_1 x_1} \dd x_1, \qquad (k_1,x_2)^\top \in \R^2.
\end{align*}

For $k_1\in[-k_0,k_0]$, we define

\begin{align*}
\kappa(k_1) \coloneqq  \sqrt{k_0^2-k_1^2}.% \qquad \text{and} \qquad \bh^\pm (k_1) \coloneqq \begin{pmatrix} k_1 \\ \pm \kappa \end{pmatrix}.
\end{align*}

We can now formulate the Fourier diffraction theorem, see for instance <cit.>, <cit.> or <cit.>.

Let $f$ be bounded with $\supp f\subset \mathcal{B}_{\rM}$. For $k_1 \in (-k_0,k_0)$, we have

\begin{equation} \label{eq:recon}
\ktran_{1}\uborn (k_1, \rM) = \sqrt{2 \pi} \frac{\i \e^{\i \kappa \rM}}{2 \kappa} \ktran f\left(k_1, \kappa-k_0 \right).
\end{equation}

According to the Fourier diffraction theorem, <ref>, the measurements of the scattered wave $\uborn$ can be related to the scattering potential $f$ on a semicircle in k-space. Below we discuss how this so-called k-space coverage of the experiment is affected by (i) rotating the object and (ii) varying the wave number $k_0$ of the incident field $\ui$.

§.§ Rotating the object

Suppose the object rotates around the origin during the experiment. Then the resulting orientation-dependent scattering potential can be written as

\begin{equation*}
f^\alpha (\bx) = f(R_\alpha \bx), \quad \bx \in \R^2,
\end{equation*}

where $\alpha$ ranges over a (continuous or discrete) set of angles $A\subset [0,2\pi]$ and

\begin{equation*} %\label{eq:R-mat}
R_\alpha = \begin{pmatrix} \cos \alpha & -\sin \alpha \\ \sin \alpha & \phantom{-}\cos \alpha \end{pmatrix}
\end{equation*}

is a rotation matrix. If we let $u^\alpha$ be the Born approximation to the wave scattered by $f^\alpha$, the collected measurements are given by

\begin{equation*}
u^\alpha(x_1,\rM), \quad x_1 \in \R, \quad \alpha \in A.
\end{equation*}

Exploiting the rotation property of the Fourier transform, $\ktran (f \circ R_\alpha) = (\ktran f) \circ R_\alpha$, we obtain
from <ref>

\begin{equation*} %\label{eq:recon_n}
\ktran_{1}u^\alpha(k_1,\rM)= \sqrt{2\pi} \frac{\i \e^{\i \kappa \rM}}{2 \kappa} \ktran f \left(R_\alpha (k_1,\kappa-k_0)^\top \right).
\end{equation*}

Thus the k-space coverage, that is, the set of all spatial frequencies $\by \in \R^2$ at which $\ktran f$ is accessible via the Fourier diffraction theorem, is given by

\begin{equation*}
\mathcal{Y} = \left\{ \by = R_\alpha (k_1,\kappa-k_0)^\top \in \R^2 : \abs{k_1} < k_0, \, \alpha \in A \right\}.
\end{equation*}

It consists of rotated versions (around the origin) of the semicircle $(k_1,\kappa-k_0)^\top$, $\abs{k_1}<k_0$, see <ref>.

[red,opacity=.3] (-1.5,1.5) arc (135:-135:2.12) – (-1.5,-1.5) arc (180:90:1.5) – (0,0) arc (270:180:1.5);

[thick, red, dotted] (-1.5,1.5) arc (180:360:1.5);
[red] (-1.5,-1.5) arc (180:0:1.5);
[red, dashed] (1.5,1.5) arc (90:270:1.5);

[->]   (-3,0) – (3,0) node[below,black] $y_1$;
[->]   (0,-2.5) – (0,2.8) node[below left,black] $y_2$;

(2.12,-0.1) – (2.12,0.1);
[above] at (1.65,0) $\sqrt{2}k_0$;
(-0.1,1.5) – (0.1,1.5);
[right] at (0.1,1.5) $k_0$;
(-0.1,-1.5) – (0.1,-1.5);
[right] at (0.1,-1.5) $-k_0$;

[red,opacity=.3] (0,0) circle (2.12);

[thick, red, dotted] (-1.5,1.5) arc (180:360:1.5);
[red] (-1.5,-1.5) arc (180:0:1.5);
[red, dashed] (1.5,1.5) arc (90:270:1.5);
[red, dashdotted] (-1.5,1.5) arc (90:-90:1.5);

[->]   (-3,0) – (3,0) node[below,black] $y_1$;
[->]   (0,-2.5) – (0,2.8) node[below left,black] $y_2$;

(2.12,-0.1) – (2.12,0.1);
[above] at (1.65,0) $\sqrt{2}k_0$;

k-space coverage for a rotating object.
Left: half turn, $A = [0,\pi]$. Right: full turn, $A = [0,2\pi]$. The k-space coverage (light red) is a union of infinitely many semicircles, each corresponding to a different orientation of the object. Some of the semicircles are depicted in red: solid arc ($\alpha = 0$), dashed ($\alpha = \pi/2$), dotted ($\alpha = \pi$), dash-dotted ($\alpha = 3 \pi/2$).

§.§ Varying wave number

Now suppose the object is illuminated or insonified by plane waves with wave numbers ranging over a set $K\subset (0,+\infty)$. Recall the definition of the scattering potential $f_{k_0} = k_0^2 (n^2/n_0^2 -1)$ from <ref>, but note that we have now added a subscript to indicate the dependence of $f$ on $k_0$. If the variation of the object's refractive index $n$ with $k_0 \in K$ is negligible, we can write

\begin{equation} \label{eq:f-multifreq}
f_{k_0}(\bx) = k_0^2 f_1(\bx), \quad \bx \in \R^2.
\end{equation}

If no confusion arises, we may write $f = f_1$. Denoting by $u_{k_0}$ the Born approximation to the wave scattered by $f_{k_0}$, the resulting collection of measurements is

\begin{equation*}
u_{k_0}(x_1,\rM), \quad x_1 \in \R, \quad k_0 \in K.
\end{equation*}

Then, according to <ref>, we have

\begin{equation*}
\ktran_{1}u_{k_0}(k_1,\rM)= \sqrt{2\pi} \frac{\i \e^{\i \kappa \rM}}{2 \kappa} k_0^2 \ktran f_1 \left( k_1, \kappa - k_0 \right).
\end{equation*}

Notice that now $\kappa$ also varies with $k_0$. The resulting k-space coverage

\begin{equation*}
\mathcal{Y} = \left\{ \by = \left( k_1, \kappa - k_0 \right)^\top \in \R^2 : k_0 \in K, \, \abs{k_1} < k_0 \right\}
\end{equation*}

is a union of semicircles scaled and shifted (in direction of $-\be_2$) according to $k_0 \in K$.

Consider, for example,

$K = [k_{\mathrm{min}},k_{\mathrm{max}}]$. Then the k-space coverage consists of all points $\left(y_1, y_2 \right)^\top\in\R^2$ such that
\abs{y_1} \le k_{\mathrm{max}}$

\sqrt{k_{\mathrm{max}}^2-y_1^2} - k_{\mathrm{max}}
\ge y_2
\ge
\begin{cases} -\abs{y_1},& \abs{y_1}\ge k_{\mathrm{min}},\\
\sqrt{k_{\mathrm{min}}^2-y_1^2} - k_{\mathrm{min}}, & \text{otherwise},
\end{cases}

see <ref>. Note that, in contrast to the scenarios of <ref>, there are large missing parts near the origin.

[red,opacity=.3] (-1.5,-1.5) arc (180:0:1.5) – (.75,-.75) arc (0:180:.75);

[red] (-1.5,-1.5) arc (180:0:1.5);
[red] (-.75,-.75) arc (180:0:.75);

[->]   (-2.5,0) – (2.5,0) node[below,black] $y_1$;
[->]   (0,-2.5) – (0,1) node[right,black] $y_2$;

[black] (-0.12,-.75) – (0.12,-.75);
[left] at (0,-0.75) $-k_{\mathrm{min}}$;
[black] (-0.12,-1.5) – (0.12,-1.5);
[left] at (0,-1.5) $-k_{\mathrm{max}}$;

k-space coverage for $k_0$ covering the interval $[k_{\mathrm{min}},k_{\mathrm{max}}]$.

§.§ Rotating the object with multiple wave numbers

We combine the two previous observations by picking a finite set of wave numbers $K \subset(0,\infty)$ and performing a full rotation of the object for each $k_0\in K$. Let $u_{k_0}^\alpha$ be the Born approximation to the wave scattered by $f^\alpha_{k_0} = k_0^2 f_1 \circ R_\alpha $. Then the full set of measurements is given by

\begin{equation} \label{eq:u_a_k}
u_{k_0}^\alpha(x_1,\rM), \quad x_1 \in \R, \quad \alpha \in [0,2\pi], \quad k_0\in K,
\end{equation}

and the Fourier diffraction theorem yields

\begin{equation} \label{eq:fdt_final}
\ktran_{1}u_{k_0}^\alpha (k_1,\rM)= \sqrt{2\pi} \frac{\i \e^{\i \kappa \rM}}{2 \kappa} k_0^2 \ktran f_1 \left(R_\alpha (k_1,\kappa - k_0)^\top \right).
\end{equation}

We deduce from <ref> and <ref> that the resulting k-space coverage $\mathcal Y$ is the union of disks with radii $\sqrt2 k_0$, all centered at the origin. Hence, $\mathcal Y$ is just the largest disk, that is, the one corresponding to the largest wave number $\max K$. However, smaller disks in k-space are covered more often, which might improve reliability of the reconstruction for noisy data.

§ RECONSTRUCTION METHODS

In the following, we assume data generated by line-sources according to the
setup described in <ref>.
We simulate the total fields solutions to <ref>,
which are the synthetic data used for the reconstruction.

§.§ Reconstruction using Full Waveform Inversion

For the identification of the physical properties of the medium,
the Full Waveform Inversion (FWI) relies on an iterative minimization
of a misfit functional which evaluates a distance between numerical
simulation and measurements of the total field.
The Full Waveform Inversion method arises in the context of
seismic inversion for sub-surface Earth imaging,
cf. <cit.>,
where the measured seismograms are compared to simulated waves.

With FWI, we invert with respect to the
wave speed $c$, from which the wave number
is defined according to <ref>.
It further connects with the model perturbation $f$ according to
In our experiment, $c_0$ is used as an initial guess
(i.e., we start from constant background) and then
$c$ is inverted rather than $f$, as discussed in <ref>.
Given some measurements $\mathbf{d}$ of the
total field, the quantitative reconstruction of
the wave speed $c$ is performed following the
minimization of the misfit functional $\misfit$ such that

\begin{equation} \label{eq:fwi:main-minimization}
\min_c \, \misfit(c) \, , \qquad
\misfit  \,=\, \mathrm{dist}\Big(\restrict \utot , \, \mathbf{d}\Big) \, ,
\qquad \text{where $u$ solves \autoref{eq:fwi:helmholtz-full}} \, .
\end{equation}

Here, $\mathrm{dist}(\cdot)$ is a distance function to evaluate
the difference between the measurements and the simulations,
and $\restrict$ is a linear operator to restrict the solution
to the positions of the receivers.
For simplicity, we do not encode a regularization term
in <ref> and refer the readers
to, e.g., <cit.> and the
references therein.

Several formulations of the distance function have been
studied for FWI (in particular for seismic applications),
such as a logarithmic criterion, <cit.>, the
use of the signal's phase or amplitude, <cit.>,
the use of the envelope of the signal, <cit.>,
criteria based upon cross-correlation,
or optimal transport distance, <cit.>.
Here, we rely on a least-squares approach where the
misfit functional is defined as the $L^2$ distance
between the data and simulations:

\begin{equation} \label{eq:fwi:misfit-l2}
\misfit(c) \,:=\,  \dfrac{1}{2} \,
\sum_{\omega\in c_0 K} \, \sum_{\alpha\in A} \Vert \, \restrict
\utot(c, \omega,\, \alpha)
\,-\, {\mathbf{d}(\omega,\alpha)} \, \Vert_{L^2(-\lM,\lM)}^2   \, ,
\end{equation}

where $\mathbf{d}(\omega,\alpha)$ refers to the measurement data of the total field at the measurement plane
with respect to the object rotated with angle $\alpha$,
and $\utot(c, \omega,\, \alpha)$ is the solution of <ref> with $k(\bx) = \omega / c(R_\alpha^\top\bx)$.
The last term $R_\alpha^\top \bx$ encodes the rotation of the object.
We note that a rotation of the object is equivalent of the rotation of both the measurement line and the direction of the incident field.
We have encoded
a sum over the frequencies $\omega$, which are
chosen in accordance with the frequency content
available in measurements.
In the computational experiments,
we further investigate uni- and multi-frequency

The minimization of the misfit functional <ref>
follows an iterative Newton-type method as depicted in <ref>.
Due to the computational cost, we use first-order information
and avoid the Hessian computation, <cit.>: Namely, we
rely on the nonlinear conjugate gradient method for the model update,
cf. <cit.>.
Furthermore, to avoid the formation of the dense Jacobian, the
gradient of the misfit functional is computed using the adjoint-state
method, cf. <cit.>.
In <ref>, we further implement a progression
in the frequency content, which is common to mitigate the
ill-posedness of the nonlinear inverse problem, <cit.>.
We further invert each frequency independently, from low
to high, as advocated by <cit.>.
For the implementation details using the HDG discretization,
we refer to <cit.>.

Initial wave speed model $c_0$.
Initiate global iteration number $\ell:= 1$;

frequency $\omega \in c_0 K$
iteration $j = 1, \, \ldots, \, n_\text{iter}$
Compute the solution to the wave equation using
current wave speed model $c_\ell$ and frequency
$\omega$, that is, the solution to <ref>
with $k(\bx)=\omega/c_\ell(\bx)$
Evaluate the misfit functional $\misfit$ in <ref>
Compute the gradient of the misfit functional using the adjoint-state method
Update the wave speed model using nonlinear conjugate gradient method to obtain $c_{\ell+1}$
Update global iteration number $\ell\gets \ell+1$

Approximate wave speed $c$, from which the scattering potential $f$ can be computed via Equations <ref> and <ref>.
Iterative reconstruction of the wave speed model
following the minimization of the misfit functional.
At each iteration, the total field solution to
<ref> is computed and the gradient
of the misfit functional is used to update the wave speed
The algorithm stops when the prescribed number of
iterations is performed for all of the frequencies
of interest.

In the computational experiments, the reconstruction
with FWI assumes the availability of the total fields
which are solutions to <ref>,
and we invert with respect to the (frequency
independent) wave speed $c$ defined in
We could instead use the representation with
relation $k^2=k_0^2 + f$, and invert with respect
to the perturbation $f$, imposing the (known) smooth
background $c_0$.
Inverting with respect to $c$ rather than $f$ is
mainly motivated by consistency with existing
literature in FWI <cit.>, in which the
background model ($c_0$) is usually unknown.
Nonetheless, reformulating the minimization with respect
to $f$ and imposing $c_0$ could improve the efficiency of
FWI, as advocated by the data-space reflectivity inversion
of <cit.>.

§.§ Reconstruction based on the Born and Rytov approximations

In this section, we present numerical methods for the computation of the Born and Rytov approximations from <ref> and <ref>, respectively, as well as the reconstruction of the scattering potential.
We concentrate on the case of full rotations of the object using incident waves with different wave numbers $k_0\in K$, see <ref>.
The tomographic reconstruction is based on the Fourier diffraction theorem, <ref>,
and the nonuniform discrete Fourier transform.
Nonuniform Fourier methods have also been applied in computerized tomography <cit.>, magnetic resonance imaging <cit.>, spherical tomography <cit.> or surface wave tomography <cit.>.

In the following, we describe the discretization steps we apply.
For $N \in 2 \mathbb N$, let

$$\mathcal I_N \coloneqq \left\{-\frac{N}{2} + j: j=0,\ldots,N-1\right\}.$$

We sample the scattering potential $f$ on the uniform grid $\mathcal R_N\coloneqq \frac{2\rs}{N} \mathcal I_N^2$ in the square $[-\rs,\rs]^2$ for some $\rs>0$.
We assume that we are given measurements of the Born approximation

u^{\alpha}_{k_0}(x_1,\rM),\qquad x_1\in [-\lM, \lM],

for $\alpha\in A\subset[0,2\pi]$ and $k_0\in K$, cf. <ref>.
We want to reconstruct the scattering potential $f=f_1$, recall <ref>, utilizing <ref>.
We adapt the reconstruction approach of <cit.>, which is written for the 3D case.
First, we need to approximate the partial Fourier transform

\begin{equation} \label{eq:Fuborn}
\mathcal F_{1} u^{\alpha}_{k_0}(k_1,\rM)
= \frac{1}{\sqrt{2\pi}} \int_{\R} u^{\alpha}_{k_0}(x_1,\rM)\, \e^{-\i x_1k_1} \dd x_1
\end{equation}

The discrete Fourier transform (DFT) of $u(\cdot,\rM)$ on $m$ equispaced points $x_1\in({2\lM}/{m}) \mathcal I_m$ can be defined by

\begin{equation} \label{eq:dft}
\mathbf F_{1,m} u (k_1,\rM)
\frac{1}{\sqrt{2\pi}}
\frac{2\lM}{m}
\sum_{x_1\in \frac{2\lM}{m} \mathcal I_m} u(x_1,\rM)\, \e^{-\i x_1k_1}
,\qquad k_1 \in \frac{\pi}{\lM} \mathcal I_m,
\end{equation}

which gives an approximation of <ref>.
Then <ref> yields

\begin{equation}\label{eq:fdt2}
k_0^2 \ktran f (R_\alpha(k_1,\kappa-k_0)^\top)
-\i \sqrt{\frac{2}{\pi}} \kappa \e^{-\i\kappa \rM}
\ktran_{1} u^{\alpha}_{k_0}(k_1,\rM)
\end{equation}

for $\abs{k_1}\le k_0$.
Considering that we sample the angle $\alpha$ on the equispaced, discrete grid $A=({2\pi}/{n_A}) \{0,1,\dots,n_A-1\}$
and some finite set $K\subset(0,\infty)$,
<ref> provides an approximation of $\ktran f$ on the non-uniform grid

\begin{multline*}
\mathcal{Y}_{m,n_A} := \Big\{ R_{\alpha}(k_1,\kappa-k_0)^\top :
\\ k_1 \in \frac{\pi}{\lM} \mathcal I_m, \, \abs{k_1}\le k_0, \, \alpha \in \frac{2\pi}{n_A} \{0,1,\dots,n_A-1\},\, k_0\in K \Big\}
\end{multline*}

in k-space,
from which we want to reconstruct the scattering potential $f$.

Let $M$ be the cardinality of $\mathcal{Y}_{m,n_A}$. The two-dimensional nonuniform discrete Fourier transform (NDFT) is the linear operator
$\mathbf F_N\colon \R^{N^2} \rightarrow \R^M$
defined for the vector
$\mathbf f_N
\coloneqq
\left( f ( \bx ) \right)_{\bx \in \mathcal R_N}
elementwise by

\begin{equation}\label{eq:ndft}
\mathbf F_N \mathbf f_N (\by)
\coloneqq \frac{1}{2\pi} \frac{(2\rs)^2}{N^2} \sum_{\bx \in \mathcal R_N} f(\bx) \e^{-\i \bx \cdot \by}
, \quad \by \in \mathcal Y_{m,n_A},
\end{equation}

see <cit.>.
It provides an approximation of the Fourier transform

\begin{equation}\label{eq:ndft-approx}
\ktran f(\by) \approx \mathbf F_N \mathbf f_N(\by), \qquad \by\in\mathcal Y_{m,n_A}.
\end{equation}

Solving an equation $\mathbf F_N \mathbf f_N (\by) = \bf b$ for $\mathbf f_N$ amounts to applying an inverse NDFT, which usually utilizes an iterative method such as the conjugate gradient method on the normal equations (CGNE), see <cit.> and <cit.>.
One should be aware that the notation regarding conjugate gradient algorithms varies in the literature:
the algorithm called CGNE in <cit.>, is known as CGNR in <cit.>.
Conversely, the algorithm CGME in <cit.> is known as CGNE in <cit.>.

In conclusion, our method for computing $f$ given the Born approximation $\uborn$ is
summarized in <ref>.

Measurement data

$$u^\alpha_{k_0}(x_1,\rM),\quad x_1\in \frac{2\lM}{m} \mathcal I_m,\ \alpha\in A=\frac{2\pi}{n_A} \{0,\dots,n_A-1\},\; k_0\in K.$$

$k_0\in K$
$\alpha\in A$
Compute $-\i \sqrt{\frac{2}{\pi}} \kappa \e^{-\i\kappa \rM}
\mathbf F_{1,m} u^\alpha_{k_0}(k_1,\rM),$ $k_1\in \frac{\pi}{\lM}\mathcal I_m,$ with a DFT in <ref>

Solve <ref> with <ref> for $\mathbf f_N$ using the conjugate gradient method
Approximate scattering potential $\mathbf f_N\approx (f(\bx))_{\bx\in\mathcal R_N}$.

Iterative reconstruction of the scattering potential $f$ based on the Born approximation using an inverse NDFT.

The Rytov approximation $\urytov$, see <ref>, is closely related to the Born approximation,
but it has a different physical interpretation.
Assuming that the measurements arise from the Rytov approximation, we apply <ref> to obtain $\uborn$ from which we can proceed to recover $f$ as shown above.
We note that the actual implementation of <ref> requires a phase unwrapping because the complex logarithm is unique only up to adding $2\pi\i$, cf. <cit.>.
In particular, we use in the two-dimensional case

\begin{equation} \label{eq:unwrap}
\uborn
= \ui \left( \i \operatorname{unwrap\left(\arg\left(\frac{\urytov}{\ui}+1\right)\right)} + \ln\abs{\frac{\urytov}{\ui}+1}  \right),
\end{equation}

where $\arg$ denotes the principle argument of a complex number and $\operatorname{unwrap}$ denotes a standard unwrapping algorithm.
For the reconstruction with the Rytov approximation, we can use <ref> as well, but we have to preprocess the data $u$ by <ref>.

§ NUMERICAL EXPERIMENTS

In this section, we carry out numerical experiments
comparing the reconstruction obtained with FWI
(<ref>) and
Born and Rytov approximations (<ref>),
using single and multi-frequency data-sets.
We consider different media with varying shapes and
amplitude for the embedded objects.
Our experiments use synthetic data with added noise:
Firstly, synthetic simulations are carried out for the known
wave speeds using the software  <cit.>.
The discretization relies on a fine mesh
(usually a few hundred thousands cells in the discretized domain)
and polynomials of order $5$ to ensure accuracy.
Then, white Gaussian noise is incorporated in the synthetic
data, with a signal-to-noise ratio of 50.
The reconstruction with FWI also relies on software ,
but uses different discretization setups to foster the computational
time: the discretization mesh is coarser (usually a few tens of thousand
cells) and the polynomial order varies with the cells, depending on the
(local to the cell) wavelength, in order to remain as small as
possible, as detailed in <cit.>. The computational
cost of FWI is further discussed in <ref>.

We perform a full rotation of the object for a single or for multiple
frequencies $\omega$ and thus wave numbers $k_0$, cf. <ref>.
For different frequencies, the scattering potential is scaled according to <ref>.
We always reconstruct the rescaled scattering potential $f_1$, which we will simply denote by $f$ in the following.
In all numerical experiments, we rely on forward data generated with the forward model of line sources in <ref>.

We compare the reconstruction quality based on
the peak signal-to-noise ratio (PSNR) of the
reconstruction $\mathbf g$ with respect to the ground
truth $\mathbf f$ determined by

\begin{equation*}
\operatorname{PSNR}(\mathbf f,\mathbf g)
\coloneqq 10 \log_{10} \frac{\max_{\bx\in\mathcal R_N} \abs{\mathbf f(\bx)}^2}{N^{-2} \sum_{\bx\in\mathcal R_N} \abs{\mathbf f(\bx) - \mathbf g(\bx)}^2},
\end{equation*}

where higher values indicate a better reconstruction quality.

§.§ Reconstruction of circular contrast with various amplitudes and sizes

For the initial reconstruction experiments,
we consider a circular object in a homogeneous
background, namely the scattering potential $f$ of <ref>.
We investigate different sizes and contrasts for
the object, as shown in <ref>.
The data are generated for $n_A=\num{40}$ angles of incidence
equally partitioned between 0 to 351,
every 9,
and the measurement line is sampled on the 200 point uniform grid $ {10}^{-1}\, \mathcal I_{200} \subset [-\lM,\lM]$
with $\lM=10$.
Let us note that in this
context of a circularly symmetric object,
the data of each angle are similar and correspond to that of
<ref> for $f=\mathbf 1_{4.5}^{\mathrm{disk}}$.

-20 20
-20 20

-5  5
Perturbation $f$ for radius 4.5 and amplitude 5: model $f=\fdiskbb$.

-5  5
Perturbation $f$ for radius 4.5 and amplitude -5: model $f=\fdiskbc$.

-1  1
Perturbation $f$ for radius 2 and amplitude 1: model $f=\fdiskaa$.

Different perturbation models $f$
used for the computational experiments,
given for frequency $\omega/(2\pi)=\num{1}$,
with the relation to the wave speed given in
Both the size and contrast
vary: we consider two radii (4.5 and 2)
and three contrasts ($\num{1}$, $\num{5}$ and $\num{-5}$
with corresponding wave speeds $c=\num{0.9876}$, $c=\num{0.9421}$ and $c=\num{1.0701}$,
respectively), for a total of six configurations.
The computations are carried out on the domain $[-50,50]\times[-50,50]$,
i.e., a slightly larger setup than <ref>,
and we only picture the area near the origin for clearer visualization.

§.§.§ Reconstruction using FWI with single-frequency data-sets

We first only use data at frequency $\omega/(2\pi)=\num{1}$,
that is, wave number $k_0=2\pi$
for the reconstruction of the different perturbations illustrated
in <ref>.
With the background $k_0=2\pi$ (i.e., wave speed of 1),
it means that we
only rely on waves with wavelength 1
when propagating in the (homogeneous) background.
Then all measurements in $\bx$ are in multiples of the wavelength.
In the case of a single frequency, only the inner loop remains
in <ref>, and we perform 50 iterations.
In <ref>, we picture the
reconstruction obtained for the six different perturbations $f$.
We observe that the reconstructions of the smaller object
of radius 2 (Figures <ref>,
<ref> and
are more accurate, both in terms of the circular shape,
and in terms of amplitude.
In the case of the larger object,
the mild amplitude (<ref>)
is accurately recovered, while the stronger contrasts
(Figures <ref> and
are only partially retrieved.
Here, the outer part of the disk appears,
but the amplitude is incorrect with a ring effect
and incorrect values in the inner area.
Therefore, the reconstruction using single-frequency data
is limited and its success depends on two factors: the size
of the object and its contrast.

-7 7
-7 7

-1  1
Reconstruction for model $f=\fdiskaa$ (PSNR 23.50).

-5  5
Reconstruction for model $f=\fdiskab$
(PSNR 24.43).

-5  5
Reconstruction for model $f=\fdiskac$
(PSNR 24.25).

-1  1
Reconstruction for model $f=\fdiskba$
(PSNR 14.79).

-5  5
Reconstruction for model $f=\fdiskbb$
(PSNR 15.71).

-5  5
Reconstruction for model $f=\fdiskbc$
(PSNR 10).

Reconstruction using iterative minimization using data
of frequency $\omega/(2\pi)=\num{1}$ only.
In each cases, 50 iterations are performed and the initial model consists
in a constant background where $k_0=2\pi$. The data consist of $n_A=\num{40}$ different
angles of incidence from 0 to 351.

§.§.§ Reconstruction using FWI with multiple frequency data-sets

The difficulty of recovering a large object with a
strong contrast can be mitigated by the use of multi-frequency
data-sets, allowing a multi-scale reconstruction <cit.>.
We carry out the iterative reconstruction using
increasing frequencies, starting with $\omega/(2\pi)=\num{0.2}$ and
up to $\omega/(2\pi)=\num{1}$.
Following <cit.>, we use a sequential
progression, that is, every frequency is inverted
The reconstructions for the object of radius 4.5
and contrast $f=\pm\num{5}$ are pictured
in <ref>.
Contrary to the case of a single frequency
(see <ref>), the
reconstruction is now accurate and stable: the
amplitude is accurately retrieved and
the circular shape is clear, avoiding the circular
artifacts observed in Figures
and <ref>.

-7 7
-7 7

-5  5
Reconstruction for model $f=\fdiskbb$
(PSNR 19.43).

-5  5
Reconstruction for model $f=\fdiskbc$
(PSNR 19.02).

Reconstruction using multi-frequency data from
$\omega/(2\pi)=\num{0.2}$ to $\omega/(2\pi)=\num{1}$.
The initial model consists in a constant wave speed $c_0=1$.
The data consist of $n_A=\num{40}$ different angles of incidence from
0 to 351.

It is possible to recover the model with a single frequency,
that needs to be carefully chosen depending on the size of
the object and the amplitude of the contrast.
We have seen in <ref> that
the frequency $\omega/(2\pi)=\num{1}$ is sufficient for the
object of radius 2 but for the radius
4.5, we need a lower frequency (i.e., larger wavelength)
to uncover the larger object.
We illustrate in <ref>
the reconstruction using data at only $\omega/(2\pi)={0.7}$, where
we see that the shape and contrast are retrieved accurately.
Nonetheless, it is hard to predict this frequency a-priori and
we believe it remains more natural to use multiple frequencies
(when available in the data), to ensure the robustness of the

-7 7
-7 7

-5  5
Reconstruction for model $f=\fdiskbb$
(PSNR 15.33).

-5  5
Reconstruction for model $f=\fdiskbc$
(PSNR 15.03).

Reconstruction using frequency $\omega/(2\pi)=\num{0.7}$.
The initial model consists in a constant wave speed $c_0=1$.
The data consist of $n_A=\num{40}$ different angles of incidence from
0 to 351.

§.§.§ Reconstruction using Born and Rytov approximations

For the reconstruction with <ref>, which relies on the Born or Rytov approximation,
we use the same data as in the above experiment.
We use a grid with $N=240$ and $\rs = {N}/({8\sqrt2})\approx10$.
The numerical results indicate that $\rs$ should not be smaller than $\lM$.
Since we have $k_1^2\le k_0^2$ and the distance between two grid points of $k_1$ is ${\pi}/{\lM}$, only around ${2k_0 \lM}/{\pi}\approx40$ of them contribute to the data of the inverse NDFT.

In the following reconstructions, we use a fixed number of 20 iteration steps in the conjugate gradient method.
Initially, we use the frequency $\omega/(2\pi)=1$ of the incident wave, therefore $k_0=2\pi$.
Reconstructions of the circular model $f=\mathbf 1^{\mathrm{disk}}_a$
are shown in <ref>.
We note that all reconstructions are reasonably good,
where the Rytov reconstruction looks slightly better inside the object.

-7.0711 7.0711
-7.0711 7.0711
-7.0711 7.0711
-7.0711 7.0711
-1  1

Model $f=\mathbf 1^{\mathrm{disk}}_{4.5}$

Born reconstruction (PSNR 19.35)

Rytov reconstruction (PSNR 19.28)

Model $f=\mathbf 1^{\mathrm{disk}}_{2}$

Born reconstruction (PSNR 24.13)

Rytov reconstruction (PSNR 24.01)

Reconstructions with with the Born and Rytov approximation, where the data $u(\cdot,\rM)$ is generated with the line source model.
The incident field has the frequency $\omega/(2\pi)=1$.
Visible is only the cut out center, where we compute the PSNR.

For a higher amplitude of the model function $f$, the limitations of the linear models become apparent.
In <ref>, we see that the Born reconstruction of the larger object fails,
and for the smaller object only the Rytov approximation yields a good reconstruction,
which is consistent with <ref>.
With the Born approximation, we recognize the object's shape but not its amplitude, which is consistent with the observations in <cit.>.
However, as we see in <ref>, even the FWI reconstruction makes a considerable error in the object's interior, and we cannot expect the linear models to be better.

-7.0711 7.0711
-7.0711 7.0711
-7.0711 7.0711
-7.0711 7.0711
-5  5

Model $f=5\cdot\mathbf 1^{\mathrm{disk}}_{4.5}$

Born reconstruction (PSNR 4.68)

Rytov reconstruction (PSNR 11.91)

Model $f=5\cdot\mathbf 1^{\mathrm{disk}}_{2}$

Born reconstruction (PSNR 19.39)

Rytov reconstruction (PSNR 21.31)

Same setting as in <ref>, but with a higher amplitude of $5$.

We see that the FWI and the Born/Rytov reconstructions contain different kinds of artifacts.
Therefore, a comparison of the visual image quality perception does not necessarily yield the same conclusions as for the computed PSNR values.
Furthermore, the size of the object has a considerable effect on the PSNR,
e.g., the images in <ref> (c) and (f) show a comparable visual quality, but the latter's PSNR is considerably better because of the lower error in the background farther away from the object,
see also <cit.> for a study on the PSNR.

In <ref>, we use the same setup as before, but with the
frequency $\omega/(2\pi)=0.7$ instead of $\omega/(2\pi)=1$ and thus the wave number $k_0=\omega$.
Apparently, the reconstruction becomes worse with lower frequency, because it provides a smaller k-space coverage.

-7.0711 7.0711
-7.0711 7.0711
-7.0711 7.0711
-7.0711 7.0711
-1  1

Born reconstruction (PSNR 18.05)

Rytov reconstruction (PSNR 16.52)

Reconstructions of $\mathbf 1^{\mathrm{disk}}_{4.5}$, where the incident field has the frequency $\omega/(2\pi)=0.7$ instead of $1$.
The rest of the setting is from <ref>.

§.§ Reconstruction of embedded shapes: Phantom 1

We consider a more challenging scenario with
shapes embedded in the background medium.
In <ref>, we picture
the perturbation $f$ consisting of a disk and heart
included in an ellipse,
with $f$ varying from 0 to 0.5.
The computational domain corresponds to
$[-20,20]\times[-20,20]$, with line-sources positioned at
a distance $R=\num{10}$ and receivers in $\rM=\num{6}$ to
capture the data.
The data are generated using $n_A=\num{100}$ incidence angles $\alpha$ equispaced on $[0,2\pi]$,
following the steps described in <ref>.
This is illustrated in <ref>.

-0.1  0.5
Perturbation model at frequency $\omega/(2\pi)=1$,
the wave speed is equal to 1 in the background.
The positions of the source and the receivers recording
transmission data are pictured in white.

-25   25

Real part of the global solution to <ref>
at frequency $\omega/(2\pi)=\num{1}$, the source is discretized by
$N_{\mathrm{sim}}=\num{1361}$ simultaneous excitations at fixed height $x_2=-10$.
-5   5

Real part of the scattering solution at frequency $\omega/(2\pi)=\num{1}$.

-25   25

Zoom near origin of figure panel b)

-5   5

Zoom near origin of figure panel c)

Scattered solution measured at the 201
receivers positioned at fixed height $x_2=\num{6}$.

Illustration of the acquisition setup and
generated data. The computations are carried
out on the domain $[-20,20]\times[-20,20]$.
While FWI uses the total field,
the reconstruction based upon Born and Rytov
approximations use the scattered solutions,
obtained after removing a reference solution
corresponding to a propagation in an homogeneous
medium, cf. <ref>.

§.§.§ Reconstruction using FWI

We carry out the reconstruction following <ref>, and
the results are pictured in <ref>,
where we compare the use of single and multi-frequency data.
In this example, we see that with relatively low frequency data
(i.e., relatively large wavelength), such as
for frequency $\omega/(2\pi)=\num{0.7}$ and $\omega/(2\pi)=\num{1}$, the reconstruction is smooth,
see <ref> and <ref>,
and one needs to use smaller wavelengths to obtain a better reconstruction,
see <ref> and <ref>.
In <ref>, we see that multi-frequency data
gives the best reconstruction, it is also the most robust as one does not need
to anticipate the appropriate wavelength before carrying out the reconstruction.
Here both the shapes and contrast in amplitude are accurately obtained.
We notice some oscillatory noise in the reconstructed models, which
could certainly be reduced by incorporating a regularization criterion
in the minimization <cit.>.

-0.1  0.5

Using single-frequency,
(PSNR 22.38).

Using single-frequency,
(PSNR 22.91).

Using single-frequency,
(PSNR 23.10).

Using single-frequency, $\omega/(2\pi)=\num{1.4}$.
(PSNR 23.31).

Using multi-frequency, $\omega/(2\pi)\in\{\num{0.7}, \num{1}, \num{1.2}, \num{1.4} \}$.
(PSNR 27.28).

Reconstruction of the model
with FWI and starting from a homogeneous
background with $f=0$. The models are given
at frequency $\omega/(2\pi)=1$ and the wave speed
is equal to 1 in the background.

In <ref>, we conduct a
similar computational experiment, but increasing the contrast in the
included heart shape where $f$ has now a value of
2, see <ref>.
We provide single and multi-frequency reconstructions,
and observe that large wavelengths still provide
a smooth reconstruction. The high contrast in the
heart is well recovered.

-0.1   2

True model.

Using single-frequency,
(PSNR 25.04).

Using single-frequency,
(PSNR 25.91).

Using single-frequency,
(PSNR 26.19).

Using single-frequency,
(PSNR 26.75).

Using multi-frequency,
$\omega/(2\pi)\in\{\num{0.7}, \num{1}, \num{1.2}, \num{1.4} \}$
(PSNR 28.76).

Reconstruction with FWI starting
from a homogeneous background
with $f=0$. The models are given
at frequency $\omega/(2\pi)=1$ and the wave speed
is equal to 1 in the background.

§.§.§ Reconstruction using Born and Rytov approximations

We perform the reconstruction with <ref> of the test model $f$ from <ref>.

In the following tests, we discretize $f$ on a finer grid of resolution $N=720$, which covers the radius $\rs = 15/\sqrt2$.
The PSNR is computed only on the central part of the grid that is visible in the image.
Since we know that the $f$ is real-valued, we take only the real part of the reconstruction.
For simplicity, we use a constant number of $12$ iterations in the conjugate gradient method of the inverse NDFT.

The Born and Rytov reconstructions are shown in <ref>,
where the data $u$ is the same as in <ref>.
The reconstruction with a higher frequency of the incident wave is more accurate, since it provides a larger k-space coverage, which is the disk of radius
$\sqrt{2} k_0 = \sqrt{2} \omega$,
see <ref>.
Moreover, the multi-frequency reconstruction is shown in <ref> (g) and (h).
Even though the multi-frequency setup covers the same disk in k-space, it still seems superior because we have more data points of the Fourier transform $\mathcal Ff$.

-0.1   0.5

True model $f$

Born reconstruction at frequency $\omega/(2\pi)={1} $. (PSNR 24.69)

Rytov reconstruction at frequency $\omega/(2\pi)={1} $. (PSNR 24.66)

Rytov reconstruction at frequency $\omega/(2\pi)={0.7} $. (PSNR 22.72)

Rytov reconstruction at frequency $\omega/(2\pi)={1.2} $. (PSNR 25.32)

Rytov reconstruction at frequency $\omega/(2\pi)={1.4} $. (PSNR 26.14)

Born reconstruction using multi-frequency, $\omega/(2\pi)\in\{0.7,1,1.2,1.4\}$. (PSNR 26.37)

Rytov reconstruction using multi-frequency, $\omega/(2\pi)\in\{0.7,1,1.2,1.4\} $. (PSNR 26.37)

Reconstructions for different frequencies of the incident wave.
The PSNR is computed on the visible part of the grid for the real part of the reconstruction, since we know that $f$ must be real.

For the similar model from <ref> with a higher contrast,
the reconstructions with Born and Rytov approximation differ more from the FWI reconstruction because of the more severe scattering, see <ref>.
In general, we can expect the FWI reconstruction to be better since it is a numerical approximation of the wave equation,
of which the Born or Rytov approximations are just linearizations.

-0.1   2

True model $f$

Born reconstruction at frequency $\omega/(2\pi)={1} $. (PSNR 23.53)

Rytov reconstruction at frequency $\omega/(2\pi)={1} $. (PSNR 24.47)

Rytov reconstruction at frequency $\omega/(2\pi)={0.7} $. (PSNR 21.78)

Rytov reconstruction at frequency $\omega/(2\pi)={1.2} $. (PSNR 25.55)

Rytov reconstruction at frequency $\omega/(2\pi)={1.4} $. (PSNR 26.40)

Born reconstruction using multi-frequency, $\omega/(2\pi)\in\{0.7,1,1.2,1.4\} $. (PSNR 23.77)

Rytov reconstruction using multi-frequency, $\omega/(2\pi)\in\{0.7,1,1.2,1.4\} $. (PSNR 25.92)

Reconstructions with a higher contrast, where the rest of the setting is the same as in <ref>.

§.§ Reconstruction of embedded shapes: Phantom 2

We now consider the case with combinations of smaller convex and
non-convex shapes included in the background medium.

§.§.§ Reconstruction using FWI

In <ref> and
<ref>, we show
the model perturbation, which consist in small
objects buried in the background.
FWI is carried out with single and multiple
frequencies, while we investigate a mild
contrast in <ref>
(where $f$ is at most 0.5) and a
stronger contrast in <ref>
(where $f$ is at most 2).
We see that the model can be recovered using a single frequency, which has to be
selected depending on the contrast. As an alternative, multi-frequency data appears
to be a robust candidate and always provides a good reconstruction, for both the
object's shape and amplitude.
The reconstruction quality with high contrast in <ref> seems to be of a similar level as with low contrast.

-0.1   0.5

True contrast.

Using single-frequency,
(PSNR 18.91).

Using single-frequency,
(PSNR 19.50).

Using single-frequency,
(PSNR 19.90).

Using single-frequency,
(PSNR 20.10).

Using multi-frequency, $\omega/(2\pi)\in\{\num{0.7}, \num{1}, \num{1.2}, \num{1.4} \}$
(PSNR 23.04).

Reconstruction with FWI starting
from a homogeneous background
with $f=0$. The models are given
at frequency $\omega/(2\pi)=1$ and the wave speed
is equal to 1 in the background.

-0.1   2

True contrast.

Using single-frequency,
(PSNR 20.75).

Using single-frequency,
(PSNR 21.50).

Using single-frequency,
(PSNR 22.27).

Using single-frequency,
(PSNR 22.72).

Using multi-frequency,
$\omega/(2\pi)\in\{\num{0.7}, \num{1}, \num{1.2}, \num{1.4} \}$
(PSNR 22.82).

Reconstruction with FWI starting
from a homogeneous background
with $f=0$. The models are given
at frequency $\omega/(2\pi)=1$ and the wave speed
is equal to 1 in the background.

§.§.§ Reconstruction using Born and Rytov approximations

In <ref>, we show the reconstruction
using Born and Rytov approximations.
Here, we can clearly see that we need a higher frequency in
order to resolve small features of the object.
However, the reconstructions are still inferior to the FWI.

-0.1   0.5

True model $f$

Born reconstruction at frequency $\omega/(2\pi)={1} $. (PSNR 20.14)

Rytov reconstruction at frequency $\omega/(2\pi)={1} $. (PSNR 20.10)

Rytov reconstruction at frequency $\omega/(2\pi)={0.7} $. (PSNR 18.00)

Rytov reconstruction at frequency $\omega/(2\pi)={1.2} $. (PSNR 21.29)

Rytov reconstruction at frequency $\omega/(2\pi)={1.4} $. (PSNR 22.17)

Born reconstruction using multi-frequency, $\omega/(2\pi)\in\{0.7,1,1.2,1.4\} $. (PSNR 21.91)

Rytov reconstruction using multi-frequency, $\omega/(2\pi)\in\{0.7,1,1.2,1.4\} $. (PSNR 21.93)

Reconstructions of the more complicated shapes.
The models are given at frequency $\omega/(2\pi)=1$.

With a higher contrast, the reconstruction with the Rytov approximation is considerably better than the one with the Born approximation, see <ref>.
This is consistent with <ref>.
Interestingly, the shapes reconstruction in high contrast barely profits from taking frequencies higher than 1, even though the k-space coverage is larger.
In this situation, the Rytov reconstruction is almost comparable to the one with lower contrast, but still worse than the FWI.

-0.1   2

True model $f$

Born reconstruction at frequency $\omega/(2\pi)={1} $. (PSNR 17.16)

Rytov reconstruction at frequency $\omega/(2\pi)={1} $. (PSNR 18.73)

Rytov reconstruction at frequency $\omega/(2\pi)={0.7} $. (PSNR 16.10)

Rytov reconstruction at frequency $\omega/(2\pi)={1.2} $. (PSNR 20.14)

Rytov reconstruction at frequency $\omega/(2\pi)={1.4} $. (PSNR 21.15)

Born reconstruction using multi-frequency, $\omega/(2\pi)\in\{0.7,1,1.2,1.4\} $. (PSNR 16.92)

Rytov reconstruction using multi-frequency, $\omega/(2\pi)\in\{0.7,1,1.2,1.4\} $. (PSNR 20.33)

Reconstructions of the more complicated shapes with a higher contrast than in <ref>.
The models are given at frequency $\omega/(2\pi)=1$.

§.§ Computational costs

Computational cost of FWI.
The computational cost of FWI comes from the discretization
and resolution of the wave problem <ref>
for each of the sources in the acquisition, coupled with the
iterative procedure of <ref>.
In our numerical experiments, we use the software  for the
iterative inversion, <cit.>, <ref>,
which relies on the Hybridizable discontinuous Galerkin
discretization, <cit.>.
The number of degrees of freedom for the discretization
depends on the number of cells in the mesh, and the
polynomial order.
In the inversion experiments, we use a fixed mesh for
all iterations, with about fifty thousands cells. On
the other hand, the polynomial order is selected depending
on the wavelength on each cell. That is, each of the cells
in the mesh is allowed to have a different order (here
between $3$ to $7$), see <cit.>.
Then, when the frequency changes, while the mesh remains
the same, the order of the polynomial evolves accordingly
to the change of wavelength.
Once the wave equation, <ref>, is
discretized, we obtain a linear system which size is the
number of degrees of freedom, that must be solved for the
different sources (i.e., the different incident angles).
We rely on the direct solver MUMPS, <cit.>, such
that once the matrix factorization is computed, the numerical
cost of having several sources (i.e., multiple right-hand sides
in the linear system) is drastically mitigated, motivating
the use of a direct solver instead of an iterative one.

Our numerical experiments have been carried out on the
Vienna Scientific Cluster [<https://vsc.ac.at/>],
using 48 cores.
For the reconstructions of <ref>,
<ref>, <ref>
and <ref>, the size of the computational
domain is 40$\times$40, with about 350000 degrees of
Using single-frequency data, 50 iterations are performed
in <ref> and the total computational time is of
about 40min.
In the case of multiple frequencies, we have a total of 120
iterations and the computational time is of about 1 45min.

Computational cost of Born and Rytov approximations.
The conjugate gradient method used in the inverse NDFT requires in each iteration step the evaluation of an NDFT <ref> and its adjoint.
We utilize the nonequispaced fast Fourier transform (NFFT) algorithm <cit.>, implemented in the open-source software library  <cit.>, which can compute an NDFT in $\mathcal{O}(N^2 \log N + M)$ arithmetic operations, which is considerably less than the $\mathcal{O}(N^2M)$ operations of a straightforward implementation of <ref>.

The numerical simulations of <ref> have been carried out on a 4-core Intel Core i5-6500 processor.
We used 12 iterations of the conjugate gradient method and noted that the reconstruction quality hardly benefits from a higher number of iterations.
The reconstruction of an image took about one second, with a grid size $720\times720$ of $f$ and $200\times100$ data points of $u$ for each frequency.
Therefore, the numerical computation of the Born and Rytov approximations is much faster than the FWI.

§ CONCLUSION

We study the imaging problem for diffraction tomography, where
wave measurements are used to quantitatively reconstruct the
physical properties, i.e., the refractive index.
The forward operator that describes the wave propagation corresponds
with the Helmholtz equation, which, under the assumption of small
background perturbations, can be represented via the Born and Rytov

Firstly, we have compared different forward models in terms of the resulting
measured data $u$. It highlights that, even in the case of a small circular
object, the Born approximation is not entirely accurate to represent the
total wave field given by the Helmholtz equation.
In addition, the source that initiates the phenomenon (e.g., a point source
located very far from the object, or simultaneous point source along a line),
also plays an important role as it changes the resulting signals,
hence leading to systematic differences depending on the choice of forward model.
We found that the line source model approximates the plane wave pretty well.

Secondly, we have carried out the reconstruction using data from the total field $\utot$,
and compared the efficiency of the Full Waveform Inversion method (FWI) with that
of the Born and Rytov approximations.
FWI works directly with the Helmholtz problem, <ref>,
hence giving a robust approach that can be implemented in all configurations,
however at the cost of possibly intensive computations.
On the other hand, the Born and Rytov are computationally cheap, but lack
accuracy when the object is too large or when the contrast is too strong.
We have also noted that the Rytov approximation gives better results than the Born one.
Furthermore, for all reconstruction methods, we have shown that
using data of multiple frequencies allows to improve the robustness
of the reconstruction by providing information on multiple wavelengths.

§.§ Acknowledgments

This work is supported by the Austrian Science Fund (FWF)
within SFB F68 (“Tomography across the Scales”), Projects
F68-06 and F68-07.
FF is funded by the Austrian Science Fund (FWF) under the Lise
Meitner fellowship M 2791-N.
Funding by the DFG under Germany's Excellence Strategy – The Berlin Mathematics Research
Center MATH+ (EXC-2046/1, Projektnummer: 390685689) is gratefully acknowledged.
For the numerical experiments, we acknowledge the use of the
Vienna Scientific Cluster  (<https://vsc.ac.at/>).

§ REFERENCES

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

ACM New York, NY, USA

journaltitleACM Transactions on Mathematical Software (TOMS)
titlePerformance and scalability of the block low-rank multifrontal factorization on multicore architectures

journaltitleJournal of Applied Mathematics and Optimisation
titleAbout the stability of the inverse problem in the 1-d wave equation

IOP Publishing

journaltitleInverse Problems
titleA priori estimates of attraction basins for nonlinear least squares, with application to Helmholtz seismic inverse problem

Blackwell Publishing Ltd

journaltitleGeophysical Prospecting
titleComparison of waveform inversion, part 2: phase approach

titleTotal variation-based phase retrieval for diffraction tomography

titleMultiscale seismic waveform inversion

journaltitleApplied Optics
shortjournalAppl. Opt.
titleValidity of diffraction tomography based on the first Born and the first Rytov approximations

Society of Exploration Geophysicists

titleMigration-based traveltime waveform inversion of 2-D simple structures: A synthetic example

journaltitleSIAM Journal on Numerical Analysis
titleUnified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems

seriesApplied Mathematical Sciences
titleInverse Acoustic and Electromagnetic Scattering Theory

journaltitleUltrasonic Imaging
shortjournalUltrason. Imaging
titleA filtered backpropagation algorithm for diffraction tomography

Cambridge University Press

titleMathematical Foundations of Imaging, Tomography and Wavefield Inversion

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleAbsorbing boundary conditions for numerical simulation of waves

booktitleIEEE Photonics Conference (IPC)
titleAn iterative reconstruction algorithm for optical diffraction tomography

Université de Pau et Pays de l'Ardour

titleContributions to seismic full waveform inversion for time harmonic wave equations: Stability estimates, convergence analysis, numerical experiments involving large scale optimization algorithms

journaltitleJournal of Open Source Software
title: time-harmonic wave modeling and inversion using hybridizable discontinuous Galerkin discretization

Society of Exploration Geophysicists

titleFull Reciprocity-Gap Waveform Inversion, enabling sparse-source acquisition

titleA priori estimates of attraction basins for velocity model reconstruction by time-harmonic Full Waveform Inversion and Data-Space Reflectivity formulation

Society of Exploration Geophysicists

titleReciprocity-gap misfit functional for Distributed Acoustic Sensing, combining data from passive and active sources

journaltitleComputer Methods in Applied Mechanics and Engineering
titleAdjoint-state method for Hybridizable Discontinuous Galerkin discretization, application to the inverse acoustic wave problem

journaltitleGeophysical Journal International
titleEigenvector Models for Solving the Seismic Inverse Problem for the Helmholtz Equation

Blackwell Publishing Ltd Oxford, UK

journaltitleGeophysical Journal International
titleTheoretical background for continental-and global-scale full-waveform inversion in the time–frequency domain

journaltitleJ. Opt. Soc. Am. A
titleHybrid diffraction tomography without phase information

Longman Scientific & Technical

seriesPitman Research Notes in Mathematics Series
titleConjugate Gradient Type Methods for Ill-Posed Problems

Springer Science & Business Media

titleNodal discontinuous Galerkin methods: algorithms, analysis, and applications

booktitleNew Trends in Parameter Identification for Mathematical Models
seriesTrends in Mathematics
titleAn SVD in Spherical Surface Wave Tomography

journaltitleInverse Problems
titleOptimal Mollifiers for Spherical Deconvolution

journaltitleInverse Probl. Imaging
titleReconstructing a Function on the Sphere from Its Means Along Vertical Slices

titleDiffraction tomography with Fourier ptychography

journaltitleTelecommunication Systems
titleThe accuracy of PSNR in predicting video quality for different video scenes and frame rates

IOP Publishing

journaltitleJapanese Journal of Applied Physics
titleCalculation of Refractive Index Distribution from Interferograms Using the Born and Rytov's Approximation

Philadelphia, PA

Society for Industrial
Applied Mathematics (SIAM)

noteReprint of the 1988 original
seriesClassics in Applied Mathematics
titlePrinciples of Computerized Tomographic Imaging

journaltitleSIAM Journal on Optimization
titleMinimization Based Formulations of Inverse Problems and Their Regularization

howpublished<https://www.tu-chemnitz.de/ potts/nfft>
noteContributors: F. Bartel, M. Fenn, T. Görner, M. Kircheis, T. Knopp, M. Quellmalz, M. Schmischke, T. Volkmer, A. Vollrath
titleNFFT 3.5, C subroutine library

journaltitleACM Trans. Math. Software
titleUsing NFFT3 - a Software Library for Various Nonequispaced Fast Fourier Transforms
pagesArticle 19130

journaltitleJournal of Scientific Computing
titleTo CG or to HDG: a comparative study

journaltitleInverse Problems
shortjournalInverse Probl.
titleFourier reconstruction for diffraction tomography of an object rotated into arbitrary orientations

journaltitleInt. J. Biomed. Imag.
noteID 24727
titleA note on the iterative MRI reconstruction from nonuniform k-space data

journaltitleSIAM J. Sci. Comput.
titleStability results for scattered data interpolation by trigonometric polynomials

Society for Industrial
Applied Mathematics

booktitleConference on Inverse Scattering: Theory and Application
titleThe seismic inverse problem as a sequence of before stack migrations

Society of Exploration Geophysicists

titleWave-equation traveltime inversion

journaltitleJ. Opt. Soc. Am. A
titlePhase-retrieval and intensity-only reconstruction algorithms for optical diffraction tomography

The Royal Astronomical Society

journaltitleGeophysical Journal International
titleMeasuring the misfit between seismograms using an optimal transport distance: application to full waveform inversion

Oxford University Press

titleFinite element methods for Maxwell's equations

titleThe Theory of Diffraction Tomography

journaltitleBMC Bioinformatics
titleODTbrain: a Python library for full-view, dense diffraction tomography

B. G. Teubner

titleThe mathematics of computerized tomography

Philadelphia, PA

seriesMonographs on Mathematical Modeling and Computation
titleMathematical Methods in Image Reconstruction

Springer Series in Operations Research

titleNumerical Optimization

journaltitleGeophysical Journal International
titleA review of the adjoint-state method for computing the gradient of a functional with geophysical applications

seriesApplied and Numerical Harmonic Analysis
titleNumerical Fourier Analysis

journaltitleIMA J. Numer. Anal.
titleA new linogram algorithm for computerized tomography

journaltitleGeophysical Journal International
titleGauss–Newton and full Newton methods in frequency–space seismic waveform inversion

Blackwell Publishing Ltd

journaltitleGeophysical Prospecting
titleComparison of waveform inversion, part 3: amplitude approach

Blackwell Publishing Ltd

journaltitleGeophysical Prospecting
titleComparison of waveform inversion, part 1: conventional wavefield vs logarithmic wavefield

journaltitleIEEE Transactions on Microwave Theory and Techniques
shortjournalIEEE Trans. Microw. Theory Tech.
titleLimitations of Imaging with First-Order Diffraction Tomography

Optical Society of America

journaltitleOptics express
titleOptical diffraction tomography for high resolution live cell imaging

titleInversion of seismic reflection data in the acoustic approximation

Blackwell Publishing Ltd Oxford, UK

journaltitleGeophysical Journal International
titleA correlation-based misfit criterion for wave-equation traveltime tomography

Society of Exploration Geophysicists

titleSH-wave propagation in heterogeneous media: velocity-stress finite-difference method

Society of Exploration Geophysicists

titleAn overview of full-waveform inversion in exploration geophysics

journaltitlePure Appl. Opt.
titleComparison of phase retrieval methods for optical diffraction tomography

journaltitleOptics Communications
shortjournalOpt. Commun.
titleThree-dimensional structure determination of semi-transparent objects from holographic data