When topological derivatives met regularized Gauss-Newton iterations in holographic 3D imaging

TL;DR

This paper introduces an automatic 3D electromagnetic imaging algorithm combining topological derivatives and regularized Gauss-Newton iterations, enabling accurate nanoscale reconstructions from holograms with fewer detectors.

Contribution

It develops a novel method integrating topological derivatives with Gauss-Newton iterations for improved 3D holographic imaging, allowing automatic object number updates during optimization.

Findings

Achieves fast, accurate 3D shape reconstructions with nanoscale resolution.

Handles multiple components and non-convex geometries effectively.

Requires fewer detectors for reliable imaging.

Abstract

We propose an automatic algorithm for 3D inverse electromagnetic scattering based on the combination of topological derivatives and regularized Gauss-Newton iterations. The algorithm is adapted to decoding digital holograms. A hologram is a two-dimensional light interference pattern that encodes information about three-dimensional shapes and their optical properties. The formation of the hologram is modeled using Maxwell theory for light scattering by particles. We then seek shapes optimizing error functionals which measure the deviation from the recorded holograms. Their topological derivatives provide initial guesses of the objects. Next, we correct these predictions by regularized Gauss-Newton techniques. In contrast to standard Gauss-Newton methods, in our implementation the number of objects can be automatically updated during the iterative procedure by new topological derivative computations. We show that the combined use of topological derivative based optimization and iteratively regularized Gauss-Newton methods produces fast and accurate descriptions of the geometry of objects formed by multiple components with nanoscale resolution, even for a small number of detectors and non convex components aligned in the incidence direction. The method could be applied in general imaging set-ups involving other waves (microwave imaging, elastography...) provided closed-form expressions for the topological and Frechet derivatives are determined.

Equations

$$\displaystyle\begin{array}[]{rcl}\mathbf{curl}\,({1\over\mu_{e}}\mathbf{curl}\,\mathbf{E})-{\kappa_{e}^{2}\over\mu_{e}}\mathbf{E}=0&\mbox{in}&\mathbb{R}^{3}\setminus\overline{\Omega},\\[4.30554pt] \mathbf{curl}\,({1\over\mu_{i}}\mathbf{curl}\,\mathbf{E})-{\kappa_{i}^{2}\over\mu_{i}}\mathbf{E}=0&\mbox{in}&\Omega,\\[4.30554pt] \hat{\mathbf{n}}\times\mathbf{E}^{-}=\hat{\mathbf{n}}\times\mathbf{E}^{+},&\mbox{on}&\partial\Omega,\\[4.30554pt] {1\over\mu_{i}}\hat{\mathbf{n}}\times\mathbf{curl}\,\mathbf{E}^{-}={1\over\mu_{e}}\hat{\mathbf{n}}\times\mathbf{curl}\,\mathbf{E}^{+},&\mbox{on}&\partial\Omega,\\[4.30554pt] {\rm lim}_{|\mathbf{x}|\rightarrow\infty}|\mathbf{x}|\big{|}\mathbf{curl}\,(\mathbf{E}-\mathbf{E}_{\rm inc})\times{\mathbf{x}\over|\mathbf{x}|}-\imath\kappa_{e}(\mathbf{E}-\mathbf{E}_{\rm inc})\big{|}=0,&&\end{array}$$

$$\displaystyle\begin{array}[]{rcl}\mathbf{curl}_{\tilde{\mathbf{x}}}\,(\mathbf{curl}_{\tilde{\mathbf{x}}}\,\tilde{\mathbf{E}})-k_{e}^{2}\tilde{\mathbf{E}}=0&\mbox{in}&\mathbb{R}^{3}\setminus\overline{\tilde{\Omega}},\\ \mathbf{curl}_{\tilde{\mathbf{x}}}\,(\mathbf{curl}_{\tilde{\mathbf{x}}}\,\tilde{\mathbf{E}})-k_{i}^{2}\tilde{\mathbf{E}}=0&\mbox{in}&\tilde{\Omega},\\ \tilde{\mathbf{n}}\times\tilde{\mathbf{E}}^{-}=\tilde{\mathbf{n}}\times\tilde{\mathbf{E}}^{+},&\mbox{on}&\partial\tilde{\Omega},\\ \beta\tilde{\mathbf{n}}\times\mathbf{curl}_{\tilde{\mathbf{x}}}\,\tilde{\mathbf{E}}^{-}=\tilde{\mathbf{n}}\times\mathbf{curl}_{\tilde{\mathbf{x}}}\,\tilde{\mathbf{E}}^{+},&\mbox{on}&\partial\tilde{\Omega},\\ {\rm lim}_{|\tilde{\mathbf{x}}|\rightarrow\infty}|\tilde{\mathbf{x}}|\big{|}\mathbf{curl}_{\tilde{\mathbf{x}}}\,(\tilde{\mathbf{E}}-\tilde{\mathbf{E}}_{\rm inc})\times{\tilde{\mathbf{x}}\over|\tilde{\mathbf{x}}|}-\imath k_{e}(\tilde{\mathbf{E}}-\tilde{\mathbf{E}}_{\rm inc})\big{|}=0,&&\end{array}$$

$$I_{\mathrm{meas}}({\mathbf{x}}_j)=|{\mathbf{E}}_{\Omega ,k_i}({\mathbf{x}}_j){|}^2,j=1,\dots ,N,$$

$$I_{\mathrm{meas}}({\mathbf{x}}_j)=|{\mathbf{E}}_{\Omega }({\mathbf{x}}_j){|}^2,j=1,\dots ,N,$$

$$J({\mathbb{R}}^3\setminus \overline{\Omega })=\frac{1}{2}\sum _{j=1}^N|I_{\Omega }({\mathbf{x}}_j)-I_{\mathrm{meas}}({\mathbf{x}}_j){|}^2.$$

$$\displaystyle\begin{array}[]{l}\|I_{\rm meas}-I_{\rm meas}^{\eta}\|_{2}=\big{(}\sum_{j=1}^{N}|I_{\rm meas}(\mathbf{x}_{j})-I_{\rm meas}^{\eta}(\mathbf{x}_{j})|^{2}\big{)}^{1/2}\leq\eta\,\|I_{\rm meas}\|_{2},\end{array}$$

$$J(\mathcal{R}\setminus \overline{B_{\epsilon }(\mathbf{x})})=J(\mathcal{R})+\frac{4}{3}\pi {\epsilon }^3{D}_T(\mathbf{x},\mathcal{R})+o({\epsilon }^3),\epsilon \to 0,$$

$${\Omega }_{\mathrm{ap}}:=\{\mathbf{x}\in \mathcal{R}|D_T(\mathbf{x},\mathcal{R})<-C\}.$$

$${\Omega }_{\mathrm{ap}}:=\{\mathbf{x}\in {\mathcal{R}}_{\mathrm{obs}}|D_T(\mathbf{x},{\mathbb{R}}^3)<(1-C_0)\underset{\mathbf{y}\in {\mathcal{R}}_{\mathrm{obs}}}{\min }{D}_T(\mathbf{y},{\mathbb{R}}^3)\},$$

$$D_T(\mathbf{x},{\mathbb{R}}^3)=3\mathrm{Re}\left[\frac{k_e^2(k_e^2-k_i^2\beta )}{(k_i^2\beta +2k_e^2)}\mathbf{E}(\mathbf{x})\cdot \overline{\mathbf{P}}(\mathbf{x})-\frac{1-\beta }{1+2\beta }\mathbf{curl}\mathbf{E}(\mathbf{x})\cdot \mathbf{curl}\overline{\mathbf{P}}(\mathbf{x})\right]$$

$$\displaystyle\begin{array}[]{r}\mathbf{curl}\,(\mathbf{curl}\,\mathbf{E})-k_{e}^{2}\mathbf{E}=0\quad\mbox{in $\mathbb{R}^{3}$},\\ {\rm lim}_{|\mathbf{x}|\rightarrow\infty}|\mathbf{x}|\big{|}\mathbf{curl}\,(\mathbf{E}-\mathbf{E}_{\rm inc})\times\hat{\mathbf{x}}-\imath k_{e}(\mathbf{E}-\mathbf{E}_{\rm inc})\big{|}=0,\end{array}$$

$$\displaystyle\begin{array}[]{r}\mathbf{curl}\,(\mathbf{curl}\,\mathbf{P})-k_{e}^{2}\mathbf{P}=2\sum_{j=1}^{N}(I_{\rm meas}-|\mathbf{E}|^{2})\mathbf{E}\,\delta_{\mathbf{x}_{j}}\quad\mbox{in $\mathbb{R}^{3}$},\\[4.30554pt] {\rm lim}_{|\mathbf{x}|\rightarrow\infty}|\mathbf{x}|\big{|}\mathbf{curl}\,\mathbf{P}\times\hat{\mathbf{x}}+\imath k_{e}\mathbf{P}\big{|}=0,\end{array}$$

$$\overline{\mathbf{P}}(\mathbf{x})=\frac{1}{k_e^2}\sum _{j=1}^N{\mathbf{curl}}_{\mathbf{x}}{\mathbf{curl}}_{\mathbf{x}}\left(G_{k_e}(\mathbf{x}-{\mathbf{x}}_j)2(I_{\mathrm{meas}}({\mathbf{x}}_j)-|{\mathbf{E}}_{\mathrm{inc}}({\mathbf{x}}_j){|}^2)\overline{{\mathbf{E}}_{\mathrm{inc}}({\mathbf{x}}_j)}\right),$$

$${\Omega }_{\mathrm{new}}:=\{\mathbf{x}\in {\Omega }_{\mathrm{ap}}|D_T(\mathbf{x},{\mathbb{R}}^3\setminus {\overline{\Omega }}_{\mathrm{ap}})<(1-c_1)\underset{\mathbf{y}\in {\Omega }_{\mathrm{ap}}}{\max }{D}_T(\mathbf{y},{\mathbb{R}}^3\setminus {\overline{\Omega }}_{\mathrm{ap}})\}$$

$$\cup \{\mathbf{x}\in {\mathcal{R}}_{\mathrm{obs}}\setminus {\overline{\Omega }}_{\mathrm{ap}}|D_T(\mathbf{x},{\mathbb{R}}^3\setminus {\overline{\Omega }}_{\mathrm{ap}})<(1-C_1)\underset{\mathbf{y}\in {\mathcal{R}}_{\mathrm{obs}}\setminus {\overline{\Omega }}_{\mathrm{ap}}}{\min }{D}_T(\mathbf{y},{\mathbb{R}}^3\setminus {\overline{\Omega }}_{\mathrm{ap}})\}.$$

$$\begin{array}[]{ll}J(\mathbb{R}^{3}\setminus(\overline{\Omega\setminus\overline{B_{\varepsilon}(\mathbf{x})}}))&=J((\mathbb{R}^{3}\cup B_{\varepsilon}(\mathbf{x}))\setminus\overline{\Omega})\\[4.30554pt] &=J(\mathbb{R}^{3}\setminus\overline{\Omega})-{4\over 3}\pi\varepsilon^{3}D_{T}(\mathbf{x},\mathbb{R}^{3}\setminus\overline{\Omega})+o(\varepsilon^{3})\end{array}$$

$$\displaystyle D_{T}(\mathbf{x},\mathbb{R}^{3}\setminus\overline{\Omega})=\left\{\begin{array}[]{l}3\,{\rm Re}\big{[}{k_{e}^{2}(k_{e}^{2}-k_{i}^{2})\over(k_{i}^{2}+2k_{e}^{2})}\mathbf{E}(\mathbf{x})\cdot\overline{\mathbf{P}}(\mathbf{x})\big{]},\quad\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Omega},\\ 3\,{\rm Re}\big{[}{k_{i}^{2}(k_{e}^{2}-k_{i}^{2})\over(k_{e}^{2}+2k_{i}^{2})}\mathbf{E}(\mathbf{x})\cdot\overline{\mathbf{P}}(\mathbf{x})\big{]},\quad\mathbf{x}\in\Omega,\end{array}\right.$$

$$\displaystyle\begin{array}[]{rcl}\mathbf{curl}\,(\mathbf{curl}\,\mathbf{E})-k_{e}^{2}\mathbf{E}=0,&\mbox{in}&\mathbb{R}^{3}\setminus\overline{\Omega},\\ \mathbf{curl}\,(\mathbf{curl}\,\mathbf{E})-k_{i}^{2}\mathbf{E}=0,&\mbox{in}&\Omega,\\ \hat{\mathbf{n}}\times\mathbf{E}^{-}=\hat{\mathbf{n}}\times\mathbf{E}^{+},&\mbox{on}&\partial\Omega,\\ \beta\hat{\mathbf{n}}\times\mathbf{curl}\,\mathbf{E}^{-}=\hat{\mathbf{n}}\times\mathbf{curl}\,\mathbf{E}^{+},&\mbox{on}&\partial\Omega,\\ {\rm lim}_{|\mathbf{x}|\rightarrow\infty}|\mathbf{x}|\big{|}\mathbf{curl}\,(\mathbf{E}-\mathbf{E}_{\rm inc})\times\hat{\mathbf{x}}-\imath k_{e}(\mathbf{E}-\mathbf{E}_{\rm inc})\big{|}=0,&&\end{array}$$

$$\displaystyle\begin{array}[]{rcl}\mathbf{curl}\,(\mathbf{curl}\,\overline{\mathbf{P}})-k_{e}^{2}\overline{\mathbf{P}}=2\sum_{j=1}^{N}(I_{\rm meas}-|\mathbf{E}|^{2})\overline{\mathbf{E}}\,\delta_{\mathbf{x}_{j}}&\mbox{in}&\mathbb{R}^{3}\setminus\overline{\Omega},\\ \mathbf{curl}\,(\mathbf{curl}\,\overline{\mathbf{P}})-k_{i}^{2}\overline{\mathbf{P}}=0,&\mbox{in}&\Omega,\\ \hat{\mathbf{n}}\times\overline{\mathbf{P}}^{-}=\hat{\mathbf{n}}\times\overline{\mathbf{P}}^{+},&\mbox{on}&\partial\Omega,\\ \beta\hat{\mathbf{n}}\times\mathbf{curl}\,\overline{\mathbf{P}}^{-}=\hat{\mathbf{n}}\times\mathbf{curl}\,\overline{\mathbf{P}}^{+},&\mbox{on}&\partial\Omega,\\ {\rm lim}_{|\mathbf{x}|\rightarrow\infty}|\mathbf{x}|\big{|}\mathbf{curl}\,\overline{\mathbf{P}}\times\hat{\mathbf{x}}-\imath k_{e}\overline{\mathbf{P}}\big{|}=0,\end{array}$$

$$\displaystyle\begin{array}[]{l}\xi_{k+1}={\rm Argmin}_{\xi\in X}\|{\cal F}(x_{k}^{\delta})+{\cal F}^{\prime}(x_{k}^{\delta})\xi-y^{\delta}\|_{Y}^{2}+\alpha_{k}\|x_{k}^{\delta}+\xi-x_{0}\|_{X}^{2},\end{array}$$

$$\displaystyle\begin{array}[]{l}\xi_{k+1}=-({\cal F}^{\prime}(x_{k}^{\delta})^{*}{\cal F}^{\prime}(x_{k}^{\delta})+\alpha_{k}I)^{-1}[{\cal F}^{\prime}(x_{k}^{\delta})^{*}({\cal F}(x_{k}^{\delta})-y^{\delta})+\alpha_{k}(x_{k}^{\delta}-x_{0})],\end{array}$$

$$\displaystyle\begin{array}[]{ll}\mathbf{q}_{\rm ap}=\mathbf{c}_{\rm ap}+\boldsymbol{\xi}_{\rm ap},&\quad\mathbf{c}_{\rm ap}\in\mathbb{R}^{3},\;\boldsymbol{\xi}_{\rm ap}(\hat{\mathbf{x}})=r_{\rm ap}(\hat{\mathbf{x}})\hat{\mathbf{x}},\\[4.30554pt] &\quad r_{\rm ap}(\hat{\mathbf{x}})=\sum_{n=0}^{n_{max}}\sum_{m=-n}^{n}\gamma_{n,m}^{\rm ap}Y_{n,m}(\hat{\mathbf{x}}),\end{array}$$

$$\mathcal{I}({\mathbf{q}}_k)+{\mathcal{I}}^{\prime }({\mathbf{q}}_k)\mathbf{\xi }\approx I_{\mathrm{meas}}^{\eta }$$

$$\displaystyle\begin{array}[]{l}\boldsymbol{\xi}_{k+1}=\underset{\boldsymbol{\xi}}{\operatorname{Argmin}}\big{\|}I^{\eta}_{\rm meas}-{\cal I}(\mathbf{q}_{k})-{\cal I}^{\prime}(\mathbf{q}_{k})\boldsymbol{\xi}\big{\|}^{2}_{\rm 2}+\alpha_{k}\big{\|}\mathbf{q}_{k}+\boldsymbol{\xi}-\mathbf{q}_{\rm ap}\big{\|}^{2}_{H^{s}(\mathbb{S}^{2})},\end{array}$$

$$\displaystyle\begin{array}[]{l}\big{(}{\cal I}^{\prime}(\mathbf{q}_{k})^{*}{\cal I}^{\prime}(\mathbf{q}_{k})+\alpha_{k}I\big{)}\boldsymbol{\xi}_{k+1}={\cal I}^{\prime}(\mathbf{q}_{k})^{*}(I^{\eta}_{\rm meas}-{\cal I}(\mathbf{q}_{k}))+\alpha_{k}(\mathbf{q}_{\rm ap}-\mathbf{q}_{k}),\end{array}$$

$$\displaystyle\begin{array}[]{l}\mathbf{q}(\boldsymbol{\xi})=\mathbf{q}_{\rm ap}+\boldsymbol{\xi},\quad\boldsymbol{\xi}(\hat{\mathbf{x}})=r(\hat{\mathbf{x}})\hat{\mathbf{x}},\end{array}$$

$${\mathcal{I}}^{\prime }({\mathbf{q}}_{\mathrm{ap}})\mathbf{\xi }=\left(\begin{array}{c}2\mathrm{Re}\left[\overline{{\mathbf{E}}_{{\Omega }_{\mathrm{ap}}}({\mathbf{x}}_1)}\cdot \dot{\mathbf{E}}({\mathbf{x}}_1)\right]\\ ⋮\\ 2\mathrm{Re}\left[\overline{{\mathbf{E}}_{{\Omega }_{\mathrm{ap}}}({\mathbf{x}}_N)}\cdot \dot{\mathbf{E}}({\mathbf{x}}_N)\right]\end{array}\right),$$

$$\displaystyle\begin{array}[]{ll}\mathbf{curl}\,(\mathbf{curl}\,\dot{\mathbf{E}})-k_{e}^{2}\dot{\mathbf{E}}=0&\quad\mbox{in $\mathbb{R}^{3}\setminus\overline{\Omega}$},\\ \mathbf{curl}\,(\mathbf{curl}\,\dot{\mathbf{E}})-k_{i}^{2}\dot{\mathbf{E}}=0&\quad\mbox{in $\Omega$},\\ \hat{\mathbf{n}}\times\dot{\mathbf{E}}^{+}-\hat{\mathbf{n}}\times\dot{\mathbf{E}}^{-}=g_{D}&\quad\mbox{on $\partial\Omega$},\\ \hat{\mathbf{n}}\times\mathbf{curl}\,\dot{\mathbf{E}}^{+}-\beta\,\hat{\mathbf{n}}\times\mathbf{curl}\,\dot{\mathbf{E}}^{-}=g_{N}&\quad\mbox{on $\partial\Omega$},\\ {\rm lim}_{|\mathbf{x}|\rightarrow\infty}|\mathbf{x}|\big{|}\mathbf{curl}\,\dot{\mathbf{E}}\times{\mathbf{x}\over|\mathbf{x}|}-\imath k_{e}\dot{\mathbf{E}}\big{|}=0,&\end{array}$$

$$g_D$$

$$g_N$$

$${\mathbf{\xi }}^{\ast }(\hat{\mathbf{x}})=r^{\ast }(\hat{\mathbf{x}})\hat{\mathbf{x}}={\mathcal{I}}^{\prime }({\mathbf{q}}_{\mathrm{ap}}{)}_{|L^2}^{\ast }\mathbf{h},$$

$$\displaystyle\begin{aligned} r^{*}=\hat{\mathbf{x}}\cdot{\cal I}^{\prime}(\mathbf{q}_{\rm ap})_{|L^{2}}^{*}\boldsymbol{h}=&\;r^{2}_{\rm ap}{\rm Re}\left[-(1-\beta^{-1})\hat{\mathbf{n}}\times\mathbf{curl}_{\mathbf{x}}\,\overline{\mathbf{E}^{+}_{\boldsymbol{h}}}\cdot\hat{\mathbf{n}}\times\mathbf{curl}_{\mathbf{x}}\,\overline{\mathbf{E}^{+}}\right.\\ &\left.\hskip 28.45274pt+(\frac{1}{k_{e}^{2}}-\frac{1}{\beta k_{i}^{2}}){\rm div}_{\partial\Omega}\big{(}\hat{\mathbf{n}}\times\mathbf{curl}\,\overline{\mathbf{E}^{+}_{\boldsymbol{h}}}\big{)}\,{\rm div}_{\partial\Omega}\big{(}\hat{\mathbf{n}}\times\mathbf{curl}\,\overline{\mathbf{E}^{+}}\big{)}\right.\\ &\left.\hskip 28.45274pt-({k_{e}^{2}}-{\beta k_{i}^{2}})\hat{\mathbf{n}}\times\overline{\mathbf{E}^{+}_{\boldsymbol{h}}}\cdot\hat{\mathbf{n}}\times\,\overline{\mathbf{E}^{+}}\right.\\ &\left.\hskip 28.45274pt+(1-\beta){\rm div}_{\partial\Omega}\big{(}\hat{\mathbf{n}}\times\overline{\mathbf{E}^{+}_{\boldsymbol{h}}}\big{)}\,{\rm div}_{\partial\Omega}\big{(}\hat{\mathbf{n}}\times\,\overline{\mathbf{E}^{+}}\big{)}\right]\circ\mathbf{q}_{\rm ap},\end{aligned}$$

Full text

When topological derivatives met regularized Gauss-Newton

iterations in holographic 3D imaging

A. Carpio 111Universidad Complutense de Madrid, Madrid 28040, Spain, and Courant Institute, New York University, NY 10012, USA, T.G. Dimiduk 222Tesla, Palo Alto, CA 94304, USA, F. Le Louër 333Université de Technologie de Compiègne, 60203 Compiègne, France, M.L. Rapún 444Universidad Politécnica de Madrid, Madrid 28040, Spain

Abstract. We propose an automatic algorithm for 3D inverse electromagnetic scattering based on the combination of topological derivatives and regularized Gauss-Newton iterations. The algorithm is adapted to decoding digital holograms. A hologram is a two-dimensional light interference pattern that encodes information about three-dimensional shapes and their optical properties. The formation of the hologram is modeled using Maxwell theory for light scattering by particles. We then seek shapes optimizing error functionals which measure the deviation from the recorded holograms. Their topological derivatives provide initial guesses of the objects. Next, we correct these predictions by regularized Gauss-Newton techniques. In contrast to standard Gauss-Newton methods, in our implementation the number of objects can be automatically updated during the iterative procedure by new topological derivative computations. We show that the combined use of topological derivative based optimization and iteratively regularized Gauss-Newton methods produces fast and accurate descriptions of the geometry of objects formed by multiple components with nanoscale resolution, even for a small number of detectors and non convex components aligned in the incidence direction. The method could be applied in general imaging set-ups involving other waves (microwave imaging, elastography…) provided closed-form expressions for the topological and Fréchet derivatives are determined.

1 Introduction

Digital in-line holography is a promising tool for high speed three dimensional (3D) imaging of live cells and soft matter [21, 44]. It can achieve high temporal (microseconds) and spatial (nanometers) resolution while avoiding the usage of toxic stains and fluorescent markers. Holograms are two-dimensional (2D) light interference patterns that contain information about the 3D positions and optical properties of an object or set of objects [56]. Figure 1 illustrates the formation of an in-line hologram from the interference of the light field scattered from a sample and the undiffracted beam [34]. In tracking experiments, one expects to infer the position of objects as a function of time analyzing a time-series of holograms. Instead, in characterization experiments, one aims to extract the size, shape and refractive index of the particles under study. Traditional optical reconstructions shine light back through the hologram to produce a 3D image, though this process may introduce a number of artifacts for sizes comparable to the employed light wavelength [54]. In contrast, digital holography aims to achieve numerical reconstructions, facing the inherent difficulty of computationally recovering 3D geometries from the 2D holograms they generate. From the mathematical point of view, finding objects and their optical properties from holograms is an ill-posed inverse scattering problem.

Recent work has relied on scattering theory to analyze holograms. As demonstrated in [52, 53], in-line holograms can be predicted combining Lorenz-Mie scattering theory with a model for the propagation and interference of the light fields in the microscope. Later, spherical colloidal particles were successfully tracked and characterized by fitting scattering models based on the Lorenz-Mie theory to holograms [34]. More general scattering approximations allow to treat non-spherical particles like rods as well as clusters of spherical particles [21, 59] fitting to the data a forward model, that is, a model that can evaluate a hologram based on a theory of scattering and propagation inside the microscope, using solvers such as discrete dipole approximations (DDA), see for instance [58, 59]. These methods proceed by ‘least squares fitting’: a few parameters representing radius, orientation, position, refractive index and so on are adjusted iteratively by a Levenberg-Marquardt algorithm to minimize the error when comparing the synthetic holograms generated by the approximate objects as predicted by the selected forward model and the true measured hologram. Alternative bayesian and machine learning approaches are discussed in [15, 57]. Successful reconstructions typically require significant a priori knowledge about the objects being imaged, such as their approximate positions in the field of view and a simple parametrization (sphere, cylinder…), for example.

To track and characterize objects without a priori knowledge (other than the optical properties of the ambient medium and the incident light) we may formulate more comprehensive optimization problems. The idea is to optimize the error functionals with respect to arbitrary shapes and arbitrary functions representing the imaged objects and their optical properties. When illuminating simple object configurations with time harmonic and polarized light, the problem somewhat simplifies, since we can use scalar Helmholtz approximations to simulate the polarized component of light [5, 6]. Within this simplified framework, Ref. [6] succeeds in producing first guesses of objects and of their optical properties from the holograms they generate by combining topological derivative and gradient based optimization procedures, provided the size of the objects is of the same order or smaller than the employed light wavelength. Helmholtz equations for the forward problems are solved in [6] by coupled boundary element (BEM)/finite element (FEM) methods [48, 51], whereas shapes are constructed by means of blobby molecule coverings and signed distance functions as in [5], without imposing any specific parametrization.

Instead of relying on particular scattering theories, here we implement topological derivative based optimization using the 3D vector Maxwell equations, which requires new closed-form formulas adapted to the holographic setting. We obtain them extending ideas developed in [6] to deal with holograms $I=|{\mathbf{E}}|^{2}$ under scalar approximations for polarized light to derivations for error functionals using complex measured data $\bf E$ and vector Maxwell constraints [47, 8, 36, 37]. As it happens in other imaging problems (acoustics, elastography…), topological derivatives allow us to generate first guesses of holographied objects and to correct the number of boundaries by creating, merging or destroying components. However, topological derivative based iterations may get stuck without providing a precise description of the shapes unless data distributed over a wide enough angle for a wide enough range of frequencies or incoming incident directions are available [1, 7, 8, 26]. We find that this is the case in holography. Due to the way microscopes are built, only one incident direction and one frequency are usually available and data are measured on a limited screen behind the object [5, 34, 59]. Moreover, our tests show that the results do not really improve switching to optimization procedures that deform the approximate object contours along vector fields, such as shape derivative based deformations or level set techniques [9, 17, 28, 46]. In this paper, we overcome this difficulty by combining topological derivative based optimization with iteratively regularized Gauss-Newton methods, see the videos in the supplemental material. For the true objects we will only require ${\cal C}^{2}$ regularity. However, those objects will be approximated by star-shaped components during the iterative procedure. Star-shaped objects are defined by rays emerging from one fixed position $\mathbf{c}$. Their boundary is located at distances of that point which vary with the angles. Therefore, they can be described using a spherical coordinate system, as it happens for ellipsoids and smoothed polyhedra, for instance. A general object formed by several star-shaped components is represented as $\Omega=\cup_{\ell=1}^{\cal L}\Omega_{\ell},$ $\partial\Omega_{\ell}=\mathbf{c}_{\ell}+\mathbf{r}_{\ell}(\mathbb{S}^{2}),$ where $\mathbf{r}_{\ell}:\mathbb{S}^{2}\rightarrow\mathbb{R}^{3}$ is a combination of spherical harmonics. Our method identifies automatically the number of defects ${\cal L}$ (i.e., we do not assume ${\cal L}$ to be known), and provides the centers $\mathbf{c}_{\ell}$ and the radii functions $\mathbf{r}_{\ell}$ defining our approximation of the true configuration, which does not need to consist in star-shaped objects. Furthermore, the method could be extended to a broader class of parameterizations dropping the star-shaped constraint, see [27]. To simplify, we assume that the optical properties are known and constant inside each component.

The starting point of our algorithm are star-shaped parametrizations fitted to initial guesses of the holographied objects provided by topological derivative optimization of the error functional comparing the synthetic hologram generated by any object with the true hologram. Next, we linearize the synthetic hologram contribution in this quadratic error functional about the current parametrization to obtain a nonlinear least squares problem for the next star-shaped parametrization, which is regularized including a Tikhonov quadratic term. Gauss-Newton iterations are a very powerful tool for finding the best descent direction for minimizing nonlinear least squares problems [2, 31, 32, 27]. There is work on coupling topological gradients and Gauss-Newton iterations [19] to solve a 2D inverse problem in elasticity, applied to parameter reconstruction instead of shape reconstruction. That study does not consider the iteratively regularized Gauss-Newton method whose regularizing Tikhonov parameter is exponentially decreasing to ensure the convergence of the algorithm [31, 32] when the stoppping rule is given by Morozov’s discrepancy principle, as we do here. Moreover, the classical application of the iteratively regularized Gauss-Newton method requires an initial guess formed by the true number of components, while our hybrid method does not: it automatically updates this number by topological derivative computations. This is the key starting point for constructing automatic/smart inverse algorithms. The hybrid inverse algorithm we are proposing is therefore new, as it is its application to holography.

Many other approaches have been developed for inverse scattering problems involving time harmonic electromagnetic waves of wavelengths larger than those of light. We may mention qualitative techniques such as linear sampling, factorization and MUSIC methods [4] and a variety of methods tracking permittivity variations [41]. In such cases, the complex amplitude (modulus and phase) of the scattered field is measured. When working with light, phases are not measurable due to its high frequency, only intensities in the form of interference patterns such as holograms are measurable. Ref. [5] showed that in case the complex amplitude was known the performance of topological methods might improve. Ref. [6] proposes a procedure based on gradient optimization and gaussian filtering to numerically approximate the modulus and the phase from the measured hologram. However, the applicability of methods exploiting this numerical approximation of the complex field is currently limited by the numerical error and by the constrained use of a single incident wave.

The paper is organized as follows. Section 2 recalls the formulation of the inverse holography problem. Section 3 adapts topological techniques to generate initial approximations of the holographied objects in the absence of a priori information, other than the emitted light and the refractive index of the ambient medium. Section 4 formulates an iteratively regularized Gauss-Newton method employing star-shaped parametrizations. This allows to use closed-form expressions of the Fréchet derivatives with respect to the parametrization in the linearized operators. It also permits the use of fast spectral solvers for the Maxwell equations. Section 5 proposes a hybrid algorithm combining topological derivative based iterations to create, merge, or destroy objects with regularized Gauss-Newton iterations to sharpen their shape. We illustrate numerically the performance of this scheme for configurations containing multiple and non necessarily convex neither star-shaped components. The numerical solution of the auxiliary Maxwell systems and details on the computation of Fréchet derivatives and their adjoints are discussed in Section 6. Finally, Section 7 summarizes our conclusions and comments on possible further developments. A few appendices contain complementary information and technical details. Appendix A briefly compares with scalar approximations based on Helmholtz equations. Appendix B establishes formulas for the Fréchet, shape and topological derivatives relevant in holography. Appendix C collects some background on the selected spherical harmonics and related Mie expansions for ease of the reader.

2 The inverse holography problem

The general form of an inverse scattering problem is the following. An incident wave ${\mathbfcal E}_{\rm inc}$ (electromagnetic, elastic, acoustic, thermal…) interacts with objects contained in a medium. The resulting wave field ${\mathbfcal E}$ is somehow measured at a set of detectors. Knowing the emitted waves, the data measured at the detectors, and the properties of the background medium, we aim to reconstruct the geometry of the objects and their relevant material properties with regard to the employed wave (permittivities and permeabilities, elastic constants, thermal diffusivities…). In other words, we seek objects $\Omega$ with material parameters $\kappa$ such that, when the emitted waves ${\mathbfcal E}_{\rm inc}$ interact with such objects, the resulting wave field ${\mathbfcal E}_{\Omega,\kappa}$ generates at the detectors data ${\cal D}_{\Omega,\kappa}$ which agree with the measured data ${\cal D}_{\rm meas}.$ In practice, one never knows the exact value of the true field at the detectors due to errors and noise. Therefore, it is often enough to find objects $\Omega$ with material parameters $\kappa$ such that the error $J(\Omega,\kappa)$ when comparing ${\cal D}_{\Omega,\kappa}$ and ${\cal D}_{\rm meas}$ at the detectors is as small as possible, that is, we seek objects $\Omega\subset{\mathbb{R}}^{3}$ and parameter functions $\kappa:\Omega\rightarrow{\mathbb{R}}$ minimizing the selected cost functional $J(\Omega,\kappa)$.

In holography, the emitted waves are light beams, the measured data take the form of a hologram ${\mathbfcal I}_{\rm meas}$ recorded at a screen behind the object, and the material parameters of interest are usually its permittivity or its refractive index. The resulting wave field is governed by the linear time dependent Maxwell equations. When the emitted light beams are time harmonic, that is, ${\mathbfcal E}_{\rm inc}(\mathbf{x},t)={\rm Re}[e^{-\imath\omega t}{\mathbf{E}}_{\rm inc}(\mathbf{x})]$, the resulting wave fields also happen to be time harmonic ${\mathbfcal E}_{\Omega,\kappa}(\mathbf{x},t)={\rm Re}[e^{-\imath\omega t}{\mathbf{E}}_{\Omega,\kappa}(\mathbf{x})]$ and the complex amplitude $\mathbf{E}_{\Omega,\kappa}(\mathbf{x})$ satisfies a stationary version of the time dependent Maxwell equations, the so-called forward problem:

[TABLE]

where $\mu_{i},\varepsilon_{i},\kappa_{i}$ and $\mu_{e},\varepsilon_{e},\kappa_{e}$ are the permeabilities, permittivities and wavenumbers $\kappa^{2}=\omega^{2}\varepsilon\mu$ of the objects and the ambient medium, respectively [3, 11, 30, 49, 51]. In biological media, $\mu_{i}\sim\mu_{e}\sim\mu_{0}$, $\mu_{0}$ being the vacuum permeability [40]. The signs $+$ and $-$ denote the values from outside and inside $\Omega$, respectively. The vector $\hat{\mathbf{n}}$ represents the outer unit normal vector. We have imposed transmission conditions at the interface $\partial\Omega$ between the objects and the ambient medium, together with the Silver-Müller radiation condition at infinity. We will consider incident plane waves polarized in a direction $\hat{\mathbf{p}}$ orthogonal to the direction of propagation $\hat{\mathbf{d}}$, that is, with $\mathbf{E}_{\rm inc}(\mathbf{x})=E_{0}\hat{\mathbf{p}}e^{\imath\kappa_{e}\hat{\mathbf{d}}\cdot\mathbf{x}}$ where $E_{0}$ represents the units and magnitude of the incident field. As sketched in Fig. 1(b), in practice the direction of propagation is the $z$ axis, orthogonal to the hologram recording screen.

It is well known that system (6) has a unique solution for any real positive $\kappa_{e}>0$ [51]. Elliptic regularity implies that this solution belongs to the Sobolev space $H^{2,0}(\Omega^{\prime})=\{{\mathbf{E}}\in H^{2}(\Omega^{\prime}),\;{\rm div}\,{\mathbf{E}}=0\}$ for any smooth bounded domain $\Omega^{\prime}\subset\mathbb{R}^{3}\setminus\overline{\Omega}$ [23, 25]. Sobolev’s embeddings ensure then continuity in $\Omega^{\prime}$.

The solution of system (6) is the forward field ${\mathbf{E}}_{\Omega,\kappa_{i}}$. The total field ${\mathbf{E}}_{\Omega,\kappa_{i}}$ is the sum of the scattered field ${\mathbf{E}}_{\rm sc,\Omega,\kappa_{i}}$ and the incident field ${\mathbf{E}}_{\rm inc}$ outside the objects, but becomes the transmitted field ${\mathbf{E}}_{\rm tr,\Omega,\kappa_{i}}$ inside. The hologram is $I_{\Omega,\kappa_{i}}=|{\mathbf{E}}_{\rm inc}+{\mathbf{E}}_{\rm sc,\Omega,\kappa_{i}}|^{2}=|{\mathbf{E}}_{\Omega,\kappa_{i}}|^{2}$ evaluated at detectors placed at the screen.

Strategies to solve numerically (6) are discussed in Section 6. For numerical purposes, it is convenient to nondimensionalize the problem. We set $\mathbf{\tilde{x}}:=\mathbf{x}/L$, $\tilde{\Omega}:=\Omega/L$, $\tilde{\mathbf{E}}(\mathbf{\tilde{x}})=\mathbf{E}(\mathbf{x}/L)/E_{0}$, $\tilde{I}(\mathbf{\tilde{x}})=I(\mathbf{x}/L)/E_{0}^{2}$, $\tilde{\mathbf{E}}_{inc}(\mathbf{\tilde{x}})=\mathbf{E}_{inc}(\mathbf{x}/L)/E_{0},$ where $L$ is a reference length unit. Choosing $L$ equal to a typical diameter of the object $\Omega$, for instance, we obtain the dimensionless forward system:

[TABLE]

$k_{i}=\kappa_{i}L$ and $k_{e}=\kappa_{e}L$ being the ‘dimensionless wavenumbers’ (size parameters in the terminology of [3]) inside and outside the object, respectively. This assumes that the permeabilities are constant in the ambient medium and in the connected components of $\Omega$, so that they can be scaled out. We will also assume $\beta={\mu_{e}\over\mu_{i}}\sim 1.$ The standard refractive indexes become $n_{i}={\kappa_{i}\lambda\over 2\pi}={k_{i}\lambda\over 2\pi L}$ and $n_{e}={\kappa_{e}\lambda\over 2\pi}={k_{e}\lambda\over 2\pi L}$, where $\lambda$ is the employed light wavelength. For ease of notation, in the sequel we will drop the  $\tilde{}$ symbol. In what follows, all the magnitudes are dimensionless. We wish to image objects using light, ranging from sizes of a few nanometers (viruses, colloidal particles, cell structures) to a few microns (prokaryotic cells). The reference length will be set to $L=1\mu$m. We will work with laser lights of wavelengths varying from $405$ nm (violet light) to $660$ nm (red light). For refractive indexes typical of cellular structures we find $k_{i}$ and $k_{e}$ in the ranges $12-15$ and $20-25$, respectively, for instance.

In the general framework described above and assuming the ambient refractive index known, the inverse holography problem consists in finding objects $\Omega=\cup_{\ell=1}^{\cal L}\Omega_{\ell}$ and functions $k_{i}:\Omega\rightarrow\mathbb{R}^{+}$ satisfying the equation

[TABLE]

where $\mathbf{E}_{\Omega,k_{i}}={\mathbf{E}}_{\rm inc}+{\mathbf{E}}_{\rm sc,\Omega,k_{i}}$ is the solution of the forward problem (12) with object $\Omega$ and dimensionless wavenumber $k_{i}$ (i.e. refractive index $n_{i}$) whereas $I_{\rm meas}$ represents the hologram measured at the screen points $\mathbf{x}_{j},$ $j=1,\ldots,N.$

To simplify, in this paper we consider that the refractive indexes are constant and known inside each component of the objects. Therefore, the problem becomes finding objects $\Omega$ such that the equation:

[TABLE]

is satisfied. Alternatively, we can reformulate this equation as a constrained optimization problem: Find the global minimum $\Omega$ of

[TABLE]

Here, $I_{\Omega}=|{\mathbf{E}}_{\Omega}|^{2}$ and $\mathbf{E}_{\Omega}$ is the solution of the dimensionless forward system (12). The object $\Omega$ is the design variable. The stationary Maxwell system (12) is the constraint. The true objects are a global minimum at which the cost functional vanishes.

In practice, we never know the true hologram $I_{\rm meas}$ but a hologram $I_{\rm meas}^{\eta}=I_{\rm meas}+I_{\eta}$ affected by noise $I_{\eta}$ of magnitude $\eta$, meaning that

[TABLE]

where $\delta=\|I_{\rm meas}-I_{\rm meas}^{\eta}\|_{2}$ is the noise level. Thus, $I_{\rm meas}$ is replaced in (14) and (15) by $I_{\rm meas}^{\eta}$. In the next section we explain how to approximate the number, size and location of holographied objects $\Omega$ from a hologram using the topological derivative of the cost functional (15).

3 Topological derivative based imaging

In the absence of any information on the holographied objects, other than the measured hologram and the nature of the ambient medium, some key features can be captured by a topological sensitivity analysis of the cost functional (15).

3.1 First guesses

The topological derivative of functional (15) measures its sensitivity to including and removing points in reconstructions of an object. Given a region ${\cal R}$ and a point $\mathbf{x}\in{\cal R}$, we have the expansion

[TABLE]

for any ball $B_{\varepsilon}(\mathbf{x})=B(\mathbf{x},\varepsilon)$ centered at $\mathbf{x}$ with radius $\varepsilon$. The coefficient $D_{T}(\mathbf{x},{\cal R})$ is the topological derivative of the functional at ${\mathbf{x}}$ [55]. When $D_{T}(\mathbf{x},{\cal R})$ is negative, the cost functional decreases for $\varepsilon>0$ small, that is, $J({\cal R}\setminus\overline{B_{\varepsilon}(\mathbf{x})})<J({\cal R})$. This suggests that the cost functional decreases by forming objects $\Omega_{\rm ap}$ with points where the topological derivative is negative and large [8, 20, 47]:

[TABLE]

For large enough $C>0,$ we expect $J({\cal R}\setminus\overline{\Omega}_{\rm ap})<J({\cal R}).$ In practice, to select a first approximation of the holographied objects we set ${\cal R}=\mathbb{R}^{3}$ and evaluate the topological derivative in a bounded region ${\cal R}_{\rm obs}\subset\mathbb{R}^{3}$ where objects are assumed to be located, i.e.,

[TABLE]

with $0<C_{0}<1$ arbitrarily chosen (in most of our examples we will set $C_{0}=0.15$ or $C_{0}=0.2$).

Computing the topological derivative using the expansion (18) is too costly from the computational point of view. We use explicit expressions in terms of adequate forward and adjoint fields instead. For the cost functional (15) and ${\cal R}=\mathbb{R}^{3}$:

[TABLE]

where $\beta={\mu_{e}\over\mu_{i}}$, see Appendix B.3. Notice that $\beta\sim 1$ for biological samples [40], therefore we will neglect the second term. The forward $\mathbf{E}$ and adjoint $\mathbf{P}$ fields are solutions of

[TABLE]

where $\delta_{\mathbf{x}_{j}}$ are Dirac masses concentrated at the detectors $\mathbf{x}_{j}$, $j=1,...,N,$ and $\hat{\mathbf{x}}={\mathbf{x}\over|\mathbf{x}|}.$ As shown in Appendix B.3, the equation governing the forward electric field $\mathbf{E}$ and the boundary conditions at the interface of the object determine the expression of the derivative. The specific dependence on $\mathbf{E}$ of the cost functional influences the source of the adjoint problem instead. Expressions of the form (21) were first established in [47] with a different adjoint field to account for a different cost functional.

For an incident plane wave $\mathbf{E}_{\rm inc}(\mathbf{x})=\hat{\mathbf{p}}e^{\imath k_{e}z},$ polarized in a direction $\hat{\mathbf{p}}$ orthogonal to the direction of propagation $(0,0,1)$, the solution of (24) is $\mathbf{E}=\mathbf{E}_{\rm inc}$. The conjugate of the solution of (27) is

[TABLE]

where $G_{k_{e}}(\mathbf{x})={e^{\imath k_{e}|\mathbf{x}|}\over 4\pi|\mathbf{x}|}$ is the outgoing Green function of Helmholtz equation, that is, the solution of $-\Delta G-k_{e}^{2}G=\delta$ satisfying an outgoing radiation condition. A quantifies the deviation of these 3D formulas for the topological derivatives of the holography cost functional with Maxwell constraints from the scalar Helmholtz approximations employed in [6].

Figure 2 shows results for an experimentally recorded hologram. Panel (b) represents the topological derivative constructed using the hologram in panel (a) in the vector Maxwell framework. The region where large negative values are attained would indicate the approximate location of an object. We can appreciate that this prediction is elongated and shifted towards the screen. Panel (c) displays the topological derivative obtained using the scalar formulas in A. This latter prediction improves noticeably and becomes comparable when we replace the true hologram by a synthetic hologram consistently generated using the scalar approximation for the polarized component. The increased shift is the result of neglecting the nonpolarized components, which contribute to the experimentally measured hologram. All the subsequent numerical tests in this paper use synthetically generated holograms, and devices of smaller size to reduce the computational cost of the evaluation of the involved fields in 3D regions. The distance to the hologram recording screen will be about $5\mu$m and the size of the screen $10\mu$m $\times$ $10\mu$m in the set-up sketched in Fig. 1(b).

Figure 3 illustrates the procedure for an ellipsoid when using red light of $660$ nm. The estimation of the object location improves as we place it closer to the recording screen and reduce its size. However, the elongation in the predicted shape remains and the true orientation of the ellipsoid is missed in the initial guess. The same behavior is observed for violet light of $405$ nm.

Fixing the hologram size, the distance to the hologram recording screen and the light wavelength, we have noticed that oscillations may appear in the topological derivative as the size of the object grows. Then, the objects are more clearly located tracking the peaks of a companion field [16], the topological energy $E_{T}(\mathbf{x},\mathbb{R}^{3})=|\mathbf{E}|^{2}|\overline{\mathbf{P}}|^{2},$ where $\mathbf{E}$ and $\mathbf{P}$ are the forward and adjoint fields governed by (24) and (27). This issue is further discussed in [5, 7, 26]. Here, we focus on improving the reconstructions of objects whose sizes are similar to or smaller than the light wavelength to avoid this transition.

3.2 Correction of the number of components

Once we have an initial reconstruction of the objects $\Omega_{\rm ap}$ we may improve it by a topological derivative based iteration. We construct a new approximation $\Omega_{\rm new}$ from $\Omega_{\rm ap}$ as follows [5, 8]:

[TABLE]

The positive constants $C_{1}$, $c_{1}$ in (29) are selected to ensure a decrease in the shape functional: $J({\mathbb{R}}^{3}\setminus\overline{\Omega}_{\rm new})<J({\mathbb{R}}^{3}\setminus\overline{\Omega}_{\rm ap})$. While the points where the updated topological derivative is negative and large are added to the former reconstructions, the points where it is positive and large are removed. In principle, we will set $C_{1}=c_{1}=C_{0}$ but they may be automatically reduced in case the cost functional does not decrease. This descent strategy is suggested by expansion (18) for $\mathbf{x}\in{\cal R}=\mathbb{R}^{3}\setminus\overline{\Omega}$ and its equivalent

[TABLE]

for $\mathbf{x}\in\Omega.$ Notice the change of sign when compared with expansion (18) for $\mathbb{R}^{3}\setminus\overline{\Omega}.$ This is a matter of choice. Keeping the same sign in both we should add and remove regions with large and negative topological derivative. Changing the sign we achieve a simpler and more visual interpretation, which is easier to implement. Also, this choice ensures that the global expression of the topological derivative when $\beta=1$ is continuous in $\mathbb{R}^{3}$ for scalar problems [6, 8].

In the presence of an object $\Omega=\Omega_{\rm ap}$, the topological derivative is given by

[TABLE]

when $\beta={\mu_{e}\over\mu_{i}}\sim 1$, see Appendix B.3, with forward and conjugate adjoint fields satisfying transmission Maxwell problems with object $\Omega=\Omega_{\rm ap}$

[TABLE]

where $\hat{\mathbf{n}}$ is the unit normal and $\hat{\mathbf{x}}={\mathbf{x}\over|\mathbf{x}|}.$

The procedure described in (29) allows for the creation of new components forming the objects (if missed in the previous iteration), for merging existing components (if they are close and intermediate points are indentified), for the destruction of proposed components (if they turn out to be spurious), and for the creation of holes inside them [8]. It does not rely on any object parametrization and can be used to reconstruct non convex and multiple objects [5, 6].

A delicate step in the process is the fitting of surfaces to the objects defined by (20) and (29), in such a way that the corresponding forward and adjoint fields can be computed numerically. Here we resort to star-shaped parametrizations and spectral solvers, as described in Sections 5 and 6. When we are only interested in determining the number and location of the objects, we may just fit spheres to each component and use explicit Mie series solutions (see Appendix C) as done in [6] for scalar problems. In general, we may avoid choosing specific parametrizations by fitting to (29) smooth surfaces defined by means of blobby molecules and signed distance functions, and employing coupled BEM-FEM solvers as in [5]. We could also work directly with (29) and DDA solvers [58, 59].

Figure 4 illustrates the performance of the iterative method based on (29) to detect two ellipsoids aligned in the incidence direction for the red light. The synthetic hologram is generated here solving numerically the forward problem (39) when the scatterers $\Omega$ are the true ellipsoids and adding $2\%$ random noise, fnamely, when $\eta=0.02$ in (17). The two objects are detected, with adequate approximate locations and sizes, but the process stagnates without recovering the correct orientation. The same happens if we modify the light, increase the number of detectors or change the thresholds $C_{0}$ and/or $C_{1}$ and $c_{1}$.

4 Iteratively regularized Gauss-Newton approach

Topological derivative-based methods yield predictions of the number, location and size of holographied objects. However, the objects appear to be displaced and elongated towards the hologram recording screen and their global 3D shape is not really recognizable, though some $xy$ slices may suggest it. We will show that this information can be refined by iteratively regularized Gauss-Newton methods (IRGNM) to obtain accurate reconstructions of the holographied shapes, provided we can rely on some kind of parametrization. However, IRGNM can only improve shapes, not alter the number of predicted object contours. This section adapts the IRGNM to holographic settings. A hybrid algorithm combining topological derivatives to allow for variations in the number of contours and IRGNM to correct shapes will be presented in Section 5.

4.1 General framework

We choose an algorithm studied in [31, 32] and based on the IRGN method introduced in [2], which enjoys some convergence properties [2, 31]. Let us recall the main ideas. Given two Hilbert spaces $X$, $Y$ and a Fréchet differentiable operator ${\cal F}:D({\cal F})\subset X\longrightarrow Y$, we assume that the exact data $y\in Y$ are attainable, that is, there is $x\in X$ such that ${\cal F}(x)=y,$ but only noisy data $y^{\delta}$ satisfying $\|y^{\delta}-y\|_{Y}\leq\delta$ are available. Starting from an initial guess $x_{0}$, the IRGN method [2] generates a sequence $x_{k+1}^{\delta}$ as follows. At each step, we linearize the equation at $x_{k}^{\delta}$, solve ${\cal F}(x_{k}^{\delta})+{\cal F}^{\prime}(x_{k}^{\delta})\xi=y^{\delta}$ by means of the minimization problem:

[TABLE]

and set $x_{k+1}^{\delta}=x_{k}^{\delta}+\xi_{k+1}.$ The penalty term $\alpha_{k}\|\xi\|_{X}$ would lead to the Levenberg-Marquart algorithm [39, 43]. The Tikhonov regularizing term $\alpha_{k}\|x_{k}^{\delta}+\xi-x_{0}\|_{X}^{2}$ has additional regularizing properties and facilitates convergence for specific choices of $\alpha_{k}$ and of initial guesses $x_{0}$ [31]. From the theory of linear Tikhonov regularization, the unique solution of (47) is given by

[TABLE]

where ${\cal F}^{\prime}(x_{k}^{\delta})^{*}$ denotes the adjoint of the Fréchet derivative ${\cal F}^{\prime}(x_{k}^{\delta})$, see B.1.

4.2 Variant adapted to inverse holography

We introduce next a procedure based on star-shaped parameterizations inspired in [27, 38]. We assume that the holographied object $\Omega_{\rm true}$ can be approximated by a finite number of star-shaped objects. This happens, of course, when $\Omega_{\rm true}$ is the union of star-shaped objects, but also for a wide range of 3D objects that are not star-shaped, like some “bean” or “peanut”-like objects. More complicated shapes, for instance “dolphin” or “aeroplane”-like ones cannot be well approximated by star-shaped parameterizations and the forthcoming method would fail. However, the ideas can be adapted to more general parameterizations [27] to deal with this kind of objects. This extension is out of the scope of the paper.

Let us first assume that $\Omega_{\rm true}$ consists in just one component which can be reasonably well-approximated by a star-shaped object $\Omega_{\rm star}$. The extension to multiple components is straightforward and will be discussed at the end of this section. Points on the surface of $\Omega_{\rm star}$ are located at rays emerging from a point $\mathbf{c}_{\rm star}$, at distances of this point which vary with the angles and can be referred to a spherical coordinate system by $\mathbf{r}_{\rm star}(\hat{\mathbf{x}})=r_{\rm star}(\hat{\mathbf{x}})\hat{\mathbf{x}},$ $\hat{\mathbf{x}}={\mathbf{x}\over|\mathbf{x}|}\in\mathbb{S}^{2}$. The spherical harmonics [11, 30] furnish a basis to expand functions defined in $L^{2}(\mathbb{S}^{2})$ which we use to approximate $r_{\rm star}(\mathbb{S}^{2})$, see C. Let $\mathbb{H}_{n_{max}}$ be the finite dimensional space spanned by scalar spherical harmonics with degree up to $n_{max}\in\mathbb{N}$. Given a star-shaped approximation $\Omega_{\mathbf{q}_{\rm ap}}$ of $\Omega_{\rm star}$ (namely, of $\Omega_{\rm true}$), we can describe its boundary as

[TABLE]

for $\hat{\mathbf{x}}\in\mathbb{S}^{2}.$ Our goal is to find $\mathbf{c}_{\rm ap}$ and $r_{\rm ap}$ in such a way that equation (14) is satisfied by the star-shaped object $\Omega=\Omega_{\rm\mathbf{q}_{\rm ap}}$ defined by it. Gradient methods would seek to correct the approximate parametrization $\mathbf{q}_{\rm ap}$ in the direction in which some kind of domain derivative of the functional (15) is negative. Gauss-Newton methods instead aim to correct the parametrization by solving equation (14) linearized at the available approximation $\mathbf{q}_{\rm ap}$. Both strategies may undergo stagnation in the direction of small gradient. To allow for convergence avoiding this artifact, regularizing terms are added to Gauss-Newton methods, as we have discussed in Section 4.1.

In our framework, starting from an approximate parametrization $\mathbf{q}_{0}$, the IRGNM solves the linearized equation

[TABLE]

by minimizing the regularized nonlinear least squares problem

[TABLE]

where $\|\boldsymbol{\xi}\|^{2}_{H^{s}(\mathbb{S}^{2})}=\sum_{n=0}^{n_{max}}\sum_{m=-n}^{n}(1+n^{2})^{s}|\gamma_{n,m}|^{2},$ with $\gamma_{m,n}$ as in (52). We choose $\alpha_{k}=\alpha_{0}({2\over 3})^{k}$ to ensure logarithmic convergence [31], and fix $\alpha_{0}=0.1$. We update the object shape by setting $\mathbf{q}_{k+1}=\mathbf{q}_{k}+\boldsymbol{\xi}_{k+1}$.

As indicated in Section 4.1, the minimizer of (54) is the unique solution of

[TABLE]

where ${\cal I}^{\prime}(\mathbf{q}_{k})^{*}$ represents the adjoint of the Fréchet derivative ${\cal I}^{\prime}(\mathbf{q}_{k})$.

Let us detail explicit expressions for both operators. Admissible perturbations $\mathbf{q}\in{\cal Q}$ of a given spherical parametrization $\mathbf{q}_{\rm ap}$ of an initial guess star-shaped with respect to a fixed center $\mathbf{c}_{\rm ap}$ take the form:

[TABLE]

and define deformed domains $\Omega_{\mathbf{q}(\boldsymbol{\xi})}$ such that $\partial\Omega_{\mathbf{q}(\boldsymbol{\xi})}:=\mathbf{c}_{\rm ap}+\boldsymbol{\xi}_{\rm ap}(\mathbb{S}^{2})+\boldsymbol{\xi}(\mathbb{S}^{2}).$ In this way, the function $\boldsymbol{\xi}$ is uniquely determined by $\partial\Omega_{\mathbf{q}(\boldsymbol{\xi})}$. As argued in Section B.1, the Fréchet derivative of the hologram admits the explicit formula

[TABLE]

where ${\bf x}_{1},\dots,{\bf x}_{N}$ are the screen detectors, and $\dot{\mathbf{E}}$ is the solution of

[TABLE]

with object $\Omega=\Omega_{\rm ap}$ and transmission data

[TABLE]

$\mathbf{E}=\mathbf{E}_{\Omega_{\rm ap}}$ being the solution of (12) with $\Omega=\Omega_{\rm ap}$ and $\hat{\mathbf{n}}$ the outer unit normal. Notice that this solution is continuously differentiable in a classical sense far from the obstacles [45].

Let us characterize now the adjoint. Let $\boldsymbol{h}$ be a real-valued function defined on the screen detectors $\mathbf{x}_{1},\ldots,\mathbf{x}_{N}$. The $L^{2}$-adjoint operator of ${\cal I}^{\prime}(\mathbf{q}_{\rm ap})$ is defined by

[TABLE]

and

[TABLE]

where ${\mathbf{E}^{+}_{\boldsymbol{h}}}$ is the solution of the transmission problem with object $\Omega=\Omega_{\rm ap}$ and incident field

[TABLE]

The IRGN method has never been used with intensities in electromagnetism before. Formula (69) is the same as that derived in [38] when working with measurements of the full field. However, the incident adjoint field $\mathbf{E}^{inc}_{\boldsymbol{h}}$ changes and it is given by (70). To obtain the adjoint of ${\cal I}^{\prime}(\mathbf{q}_{\rm ap})$ in $H^{s}$ it suffices to compose ${\cal I}^{\prime}(\mathbf{q}_{\rm ap})_{|L^{2}}^{*}$ with the diagonal operator

[TABLE]

where $\langle\phi_{1},\phi_{2}\rangle_{L^{2}}=\int_{\partial\Omega}\phi_{1}\cdot\overline{\phi_{2}}$. We then have

[TABLE]

see Appendix B.1 for details. Formulas (66)-(67) can also be simplified by using only $\mathbf{E}^{+}$, as done in [38], Remark 5, thanks to the transmission conditions. These explicit expressions for the Fréchet derivatives and their adjoints allow us to solve numerically equation (56) in $H^{s}(\mathbb{S}^{2})$ as explained in Section 6.

This idea can be generalized to multiple simply connected domains just by considering a set of parametrizations $\mathbf{q}=(\mathbf{q}_{\ell})_{\ell=1}^{\mathcal{L}_{\rm ap}}$ for the union of the boundaries of the components $\Omega_{\rm ap,\ell}$ of $\Omega_{\rm ap}.$ Notice that the procedure we have just described fixes the approximate center of the object $\mathbf{c}_{\rm ap}$ and seeks for an adequate spherical parametrization about it. We correct the centers by a simple procedure embedded in the algorithm described in the next section.

5 Hybrid Topological derivative/Gauss-Newton algorithm

In general, we ignore the number of components of the holographied object a priori. Therefore a scheme should allow to create, merge and destroy contours at certain stages. In this section, we propose a hybrid scheme combining topological derivative based optimization with iteratively regularized Gauss-Newton iterations for inverting holographic data which achieves that goal. The algorithm proceeds in the following steps.

*Step 1 - Observation region. * We define a meshed region $\mathcal{R}_{\rm obs}$ where we seek the objects and evaluate the different fields. In our set-up, we assume the screen is the square $[0,10]\times[0,10]\times\{10\}$. Considering $\mathcal{R}_{\rm obs}=[d,10-d]\subset[0,10]^{3}$, the 3D mesh is formed by the $(M_{d}+1)^{3}$ points $\mathbf{x}=(x_{1,m_{1}},x_{2,m_{2}},x_{3,m_{3}})$, where $x_{i,m_{i}}=d+m_{i}{{10-2d}\over M_{d}}$, $m_{i}=0,\ldots,M_{d}$, $i=1,2,3$.

Step 2 - Initial Guess.

1. (i)

Choose $0<C_{0}<1$ and define $\Omega_{\rm ap}$ as in (20) where $D_{T}(\mathbf{x},\mathbb{R}^{3})$ is given by (21) with $\beta=1$ and (24)-(28). Unless otherwise stated, in all our tests we set $C_{0}=0.15$ for the red light and $C_{0}=0.2$ for the violet one. 2. (ii)

For all simply connected components $\Omega_{{\rm ap},\ell}$ of $\Omega_{\rm ap}$, compute the centroid $\mathbf{c}_{{\rm ap},\ell}\in\mathbb{R}^{3}$ and the boundary points $\mathbf{c}_{{\rm ap},\ell}+\mathbf{r}_{{\rm ap},\ell}$. 3. (iii)

Determine star-shaped parametrizations for $\mathbf{r}_{{\rm ap},\ell}(\hat{\mathbf{x}})=r_{{\rm ap},\ell}(\hat{\mathbf{x}})\hat{\mathbf{x}}$ finding the coefficients $\gamma_{n,m}=\int_{\mathbb{S}^{2}}r_{{\rm ap},\ell}(\hat{\mathbf{x}})\overline{Y_{n,m}(\hat{\mathbf{x}})}dS_{\hat{\mathbf{x}}}$ in the expansion (52) by means of the exact quadrature rule (2.42) in [22] after interpolating $r_{{\rm ap},\ell}(\hat{\mathbf{x}})$ to obtain its values at the Gauss quadrature points, for each component.

We can store the result of checking the condition (i) on the topological derivative at each grid point as a binary vector ${\tt M}$ of $(M_{d}+1)^{3}$ entries. In practical implementations, this vector will be longer. After obtaining $\Omega_{\rm ap}$ as described in Step 2(i), we consider a small cube around each simple connected component. This cube is again meshed and the topological derivative (TD) is evaluated at these points. In the end, the domain $\Omega_{\rm ap}$ is defined by the set of all the points in this finer grid where the TD satisfies (20).

For (ii), the MATLAB command CC = bwconncomp(M) extracts the number of connected components and labels the points within them. The command boundary extracts the boundary points of each component. The property centroid of the command regionprops provides the centroids. See also [36] for simple MATLAB instructions regarding (iii).

As a result, we obtain an approximation $\Omega_{0}=\Omega_{\rm ap}$ with ${\mathcal{L}_{\rm ap}}$ components $\Omega_{0,\ell}=\Omega_{\rm ap,\ell}$ parameterized by $\mathbf{q}_{0,\ell}=\mathbf{q}_{\rm ap,\ell}=\mathbf{c}_{\rm ap,\ell}+\mathbf{r}_{\rm ap,\ell}=\mathbf{c}_{0,\ell}+\mathbf{r}_{0,\ell}$, $\ell=1,\ldots,{\mathcal{L}_{\rm ap}}$ (in principle, ${\cal L_{\rm ap}}$ does not coincide with the true number of objects ${\mathcal{L}}$). Then, for each iteration $k\geq 1$:

Step $2k+1$ - Shape Correction.

1. (i)

Update the shape of each connected component $\Omega_{k,\ell}$ of $\Omega_{k}$ by computing the solution $\boldsymbol{\xi}_{k+1,\ell}$ of (54)-(56) and setting $\mathbf{q}_{k+1,\ell}=\mathbf{q}_{k,\ell}+\boldsymbol{\xi}_{k+1,\ell}$ for $\ell=1,\ldots,{\mathcal{L}}_{\rm ap}.$ 2. (ii)

If any of the centers ${\bf c}_{k,\ell}$ is far from the gravity center of the corresponding component $\Omega_{k+1,\ell}$ (i.e. the distance is larger than $10\%$ of the diameter), then we set ${\bf c}_{k+1,\ell}$ equal to the gravity center and replace $\mathbf{q}_{k+1,\ell}$ with the parametrization that uses that center, as indicated in Step 2 (ii)-(iii). Otherwise, set ${\bf c}_{k+1,\ell}={\bf c}_{k,\ell}.$ 3. (iii)

Solve the forward problem (39) with $\Omega=\Omega_{k+1}$ and evaluate the solution $\mathbf{E}_{\Omega_{k+1}}$ and the corresponding hologram $I_{\Omega_{k+1}}$ at the detectors $\mathbf{x}_{1},\ldots,\mathbf{x}_{N}$.

Step $2k+2$ - Stop or modify contours.

1. (i)

We estimate the noise in the measured hologram $\delta=\|I_{\rm meas}^{\eta}-I_{\rm meas}\|_{2}$ and select as stopping criterion the discrepancy principle. If

[TABLE]

we end the algorithm. To ensure an accurate approximation preventing early stops, we set $\tau=1.01$. 2. (ii)

If the cost functional stagnates without fulfilling the stopping criteria, i.e.,

[TABLE]

then we automatically create and/or merge and/or destroy components using the topological derivative:

1. (a)

Solve the adjoint problem (45) with $\Omega=\Omega_{k+1}$, $\mathbf{E}=\mathbf{E}_{\Omega_{k+1}}$ and compute $D_{T}({\bf x},\mathbb{R}^{3}\setminus\overline{\Omega}_{k+1})$ using (33) for all ${\bf x}\in\mathcal{R}_{\rm obs}\setminus\overline{\Omega}_{k+1}$. 2. (b)

Replace $\Omega_{k+1}$ by

[TABLE]

If we wish to allow for the destruction of existing components too, then we use (29) to generate $\Omega_{\rm new}$ from $\Omega_{\rm ap}=\Omega_{k+1}.$ 3. (c)

Go to Step 2 (ii) and repeat the procedure with $\Omega_{\rm ap}=\Omega_{\rm new}.$ 3. (iii)

Otherwise, go to Step $2k+1$ (i) and repeat both Steps now for $\Omega_{k+1}$.

We solve the problem (54) for the parameterization correctors and the forward and adjoint problems for the topological derivatives (39) and (45) using the methods described in Section 6.

We have tested this algorithm in the set-up depicted in Fig. 1(b). After non-dimensionalizing as indicated in Section 2, the detectors are placed on the screen $[0,10]\times[0,10]\times\{10\}$. The holograms are synthetically generated on that screen on a grid with $51\times 51$ detectors by solving numerically the corresponding forward problems and adding random noise of magnitude $2\%$, i.e., such that (17) holds for $\eta=0.02$. Detectors are located at the points ${\bf x}_{k\ell}=(0.2k,0.2\ell,10)$, $k,\ell=0,\dots,50$. The direction of the incident light is $(0,0,1)$ and the polarization vector is (1,0,0). We fix the refractive indexes of the objects, so that $k_{e}=12.56$ and $k_{i}=15.12$ for a red incident light of $660$ nm and $k_{e}=20.6$, $k_{i}=24.79$ for a violet incident light of $405$ nm. As said earlier, $\beta=1$. We keep the default values for the parameters $C_{0}=0.15$ (red light), $0.2$ (violet light), $C_{1}=c_{1}=C_{0}$, $\alpha_{0}=0.1$ and $\tau=1.01$ indicated in the description of the algorithm, unless stated otherwise. We have fixed in all our tests $n_{max}=8$ for the definition of the parameterizations (52).

Figure 5 revisits the approximation obtained in Figure 4 with this technique keeping the parameters for red light. The location, size and shape of the two object components is now recovered with accuracy, setting $\Omega_{\rm ap}$ equal to the initial approximation in Fig. 4(d). The cost functional (15), which stagnated during the iterative procedure proposed in Section 3.2, decreases now after the introduction of the second component, as the location, size and orientation of both components improves. Video 1 in the Supplementary Material reproduces the whole sequence of approximations. In Figure 6 we repeat the tests for the violet light. We observe that our method is very robust with respect to the choice of the threshold $C_{0}$ in (20) and (75). It is able to find very accurate reconstructions even when the initial guess has only one component which is rather small (panels (a,b,c)), when it has only one component comparable in size with the objects to be found (panels (d,e,f)), and when the initial guess has the correct number of contours, but sizes are overestimated and orientations are incorrect (panels (g,h,i)). Comparing the decay of the cost functional for the three selected values, we find that for $C_{0}=0.3$ the cost functional decays rapidly during the first iterations because the initial guess has the correct number of components, and the number of iterations is slightly smaller because stagnation only occurs at the end of the procedure. Remarkably, when starting with a wrong number of components (panels (c,f)), the cost functional stagnates at the 5th iteration and the number of components is updated. At about the 10th iteration, the value of the cost function is almost the same for the three situations. Moreover, we have observed that a further increase in $C_{0}$ promotes initial guesses for which the cost functional increases, so that the method automatically reduces such constants. In the sequel we will set $C_{0}=0.2$ for the violet light.

We next replace the ellipsoid located at $(5,5,4)$ by a peanut. The parameterization of the boundary of a peanut centered at $(c_{x},c_{y},c_{z})$ is:

[TABLE]

with $\theta\in[0,\pi]$, $\phi\in[0,2\pi]$ and $r_{\gamma}(\theta)=\frac{1}{4\sqrt{1+\sqrt{\gamma+1}}}\sqrt{\cos(2\theta)+\sqrt{\gamma+1-\sin^{2}(2\theta)}}.$ We select $(c_{x},c_{y},c_{z})=(5,5,4)$ and $\gamma=0.5$. The parameter $\gamma>0$ determines the narrowness in the middle of the peanut shape (the closer to zero, the more constricted). Figure 7 illustrates the good performance of the hybrid algorithm, which recovers both shapes in 23 iterations for violet light (results for red light are similar in 24 iterations), see also Video 2 in the Supplementary Material.

If we replace the peanut by two spheres, we recover the three objects in $22$ iterations of the hybrid TD/IRGN method as shown in Figure 8, see also Video 3 in the Supplementary Material. Figure 9 and Video 4 illustrate convergence when the two spheres are placed asymmetrically. Notice that in this case, recovering the sphere aligned with the ellipsoid should be in principle more complicated than approximating the sphere not screened by it. However, our method overcomes the difficulty and provides an accurate description of both spheres.

Finally, in Figure 10 (see also Video 5 in the supplementary material), we illustrate the performance of the hybrid algorithm in a configuration with three objects of different sizes and shapes. One of them is a bean, a non-star shaped object described by the parameterization

[TABLE]

with $A(\theta)\!=\!\sqrt{0.64(1\!-\!0.1C(\theta))}$, $B(\theta)\!=\!\sqrt{0.64(1\!-\!0.4C(\theta))}$, $C(\theta)=\cos(\pi\cos\theta)$ for $\theta\in[0,\pi]$, $\phi\in[0,2\pi]$. The first identified object is the bean, that is closer to the screen. When the algorithm carries out a topological derivative iteration to update the number of objects, the existing approximation resembles an ellipsoid rather than a bean. After this update, only the ellipsoid is detected because is bigger than the sphere, which is not seen. Once this new component is found, the algorithm requires a new topological derivative iteration, finally determining the correct number of objects. At the 6th iteration (see panel (j)), the three objects resemble ellipsoids, but in a few iterations the true shapes are approximated with accuracy (panel (l)).

Remarkably, the quality of the previous approximations is maintained when reducing drastically the number of detectors, in agreement with the observations made in [13]. Let us consider a much coarser detector grid formed by just $121$ detectors located at the points ${\bf x}_{k\ell}=(k,\ell,10)$, $k,\ell=0,\dots,10$. Figure 11 shows the results obtained for different sets of scatterers using violet light. Comparing Figure 11(b,e,h,k) with Figure 6(e), Figure 7(l), Figure 9(c) and Figure 10(l), respectively, we observe that the quality of the approximations is similar.

6 Numerical solution of boundary value and linearized problems

The optimization techniques to address the inverse holography problem described in Sections 3-5 rely on the availability of adequate solvers for the associated forward and adjoint problems on one side, and the characterization of the Fréchet derivatives on the other. All of them can be rewritten to adopt the same mathematical structure.

6.1 Forward and adjoint fields

Expressing the total light field $\mathbf{E}$ in terms of transmitted ${\mathbf{E}}_{\rm tr}$, scattered ${\mathbf{E}}_{\rm sc}$ and incident ${\mathbf{E}}_{\rm inc}$ wave fields, and the conjugate adjoint field in a similar way where $\overline{\mathbf{P}}_{\rm inc}$ is the solution of the conjugate adjoint problem in the whole space (28) we find

[TABLE]

systems (39) and (45) governing the forward and adjoint fields become

[TABLE]

for $\mathbf{W}=\mathbf{E}$ or $\overline{\mathbf{P}}.$ In principle, any solver for transmission Maxwell problems can be used.

In our setting, for scatterers whose boundary is defined by a ${\cal C}^{2}$-parameteriza- tion (in particular, by a star-shaped one) with piecewise constant refractive indexes, we use fast boundary integral/spectral codes [22, 38]. Maxwell equations are solved by a Muller type boundary integral formulation [18, 38, 50]. The parametrization is exploited to transport the integral equations into unit spheres and the system is solved by means of Galerkin schemes using tangential vector spherical harmonics of low order [38]. This approach may be faster than implementing the discrete dipole approximation [58, 59] or combing boundary element (BEM)/ finite element (FEM) in a 3D region [5, 48, 51]. Unlike these two latter methods, it is only applicable when a specific parametrization is available. Fast multipole scattering methods would also be adequate in the presence of many particles [24].

6.2 Linearized equation to correct parameterizations

The least squares problem (54) is minimized solving numerically equation (56) in $H^{s}(\mathbb{S}^{2})$ using the conjugate gradient method with initial solution $\boldsymbol{\xi}_{k+1}=0$ and residual $10^{-2}$. In each step of the conjugate gradient method, we must apply the operators ${\cal I}^{\prime}(\mathbf{q}_{k})$ and ${\cal I}^{\prime}(\mathbf{q}_{k})^{*}$ to known vectors. Applying ${\cal I}^{\prime}(\mathbf{q}_{k})$ to a vector $\boldsymbol{\xi}$ amounts to solving system (65). System (65) has the structure (88) with transmission data (66)-(67). The techniques mentioned in Section 6.1 can be used to solve it.

To compute numerically the adjoint (71) there are two possibilities, depending on how we evaluate the Cauchy data $(\hat{\mathbf{n}}\times\,\overline{\mathbf{E}^{+}_{\boldsymbol{h}}},\hat{\mathbf{n}}\times\mathbf{curl}\,\overline{\mathbf{E}^{+}_{\boldsymbol{h}}})$:

• (i)

Either we compute the data $(\hat{\mathbf{n}}\times\,{\mathbf{E}^{+}_{\boldsymbol{h}}},\hat{\mathbf{n}}\times\mathbf{curl}\,{\mathbf{E}^{+}_{\boldsymbol{h}}})$ by solving the transmission problem with the incident field $\mathbf{E}^{inc}_{\boldsymbol{h}}$ as indicated above. Then we take the conjugate.

• (ii)

Or else, we compute the hermitian adjoint (transpose of the conjugate matrix) of the discretized approximation of the operator

[TABLE]

that provides the vector $(\hat{\mathbf{n}}\times\overline{\mathbf{E}^{+}_{\boldsymbol{h}}},\hat{\mathbf{n}}\times\mathbf{curl}\,\overline{\mathbf{E}^{+}_{\boldsymbol{h}}})$. This is our choice. Indeed, we have

[TABLE]

7 Conclusions

We have proposed a fully automatic algorithm to numerically reconstruct 3D objects from the holograms they generate. The method seeks to minimize the difference when comparing the true holograms and the holograms generated by approximate objects as predicted by a Maxwell forward model of light scattering. The topological derivative of this error functional provides an initial guess of the holographied objects in the absence of a priori information, other than the ambient refractive index and the incident wave. Working with star-shaped parameterizations, we implement an iteratively regularized Gauss-Newton method to successively correct the parametrization by solving linearized problems. Automatically combined with additional topological derivative based iterations to create or destroy objects when the decrease in the error functional stagnates, this procedure yields accurate reconstructions of a variety of objects in an experimental holography setting. This scheme is quite general and can be applied to other inverse electromagnetic problems by changing the error functional, so that it is defined in terms of the adequate data. This usually requires adjusting the sources in the adjoint problems and adjusting the Fréchet derivatives.

In general ‘nonlinear least squares’ problems, gradient methods seek to deform initial domains in the direction in which some kind of derivative of the error functional is negative, so that the functional decreases. This kind of methods includes topological derivative (TD), shape derivative (SD) and level set (LS) based optimization. Instead, Gauss-Newton type methods aim to correct the current approximation by linearizing about it and solving the resulting problem. Even when no spurious solutions are introduced by these procedures, they all may suffer from stagnation in the direction of small ‘gradient’. Here, the inclusion of Tikhonov regularizing terms in Gauss-Newton methods allows for convergence avoiding this artifact. From a technical point of view, while Gauss-Newton approaches require an explicit expression for a Fréchet derivative with respect to the domain, descent methods based on TD, SD or LS rely on explicit expressions for the shape/topological derivatives. Such expressions can usually be obtained introducing auxiliary adjoint problems without the need of characterizing the Fréchet derivative. The availability of these latter characterizations allows us to implement the IRGNM considered here.

We have focused on recovering shapes assuming the refractive indexes of the objects known. To obtain both we may combine our algorithm for shape optimization with parameter optimization, as in [6] or explore bayesian techniques [15]. Moreover, this paper considers object sizes of the same or smaller order than the employed light wavelength. The design of effective methods for larger objects and for objects whose components cannot be well approximated by star-shaped descriptions will require further developments.

In principle, the algorithm we propose could be extended to other inverse electromagnetic scattering problems by adjusting the closed-form formulas for the pertinent derivatives, and also to other inverse scattering problems, provided closed-form formulas for the derivatives as well as solvers for forward/adjoint problems are available (acoustics, elastography…).

Acknowledgments

This research has been supported by MINECO grants No. MTM2017-84446-C2-1-R (AC, MLR) and TRA2016-75075-R (MLR) as well as by sabbatical funds from Fundación Caja Madrid and the Salvador de Madariaga Program PRX18/00112 (AC). T.G. Dimiduk and A. Carpio thank V.N. Manoharan for discussions of holography and the Kavli Institute Seminars at Harvard for the interdisciplinary communication environment that initiated this work. A. Carpio thanks M.P. Brenner for hospitality while visiting Harvard University and R.E. Caflisch for hospitality while visiting the Courant Institute at NYU.

Appendix A Comparison with scalar approximations

When working with an incident wave polarized in the $x$ or $y$ directions, such as $\mathbf{E}_{\rm inc}=(E_{\rm inc},0,0)$, $E_{\rm inc}=e^{\imath k_{e}z},$ a standard approximation sets the non polarized components of the electric field equal to zero and uses a scalar Helmholtz transmission problem to approximate the polarized component $E_{1}$. In this framework, the inverse problem (14) is reformulated as finding $\Omega$ such that the solution $E_{\Omega}$ of the scalar forward problem

[TABLE]

satisfies the equations $I_{\rm meas}(\mathbf{x}_{j})=|E_{\Omega}(\mathbf{x}_{j})|^{2}$, $j=1,\dots,N$, or is a global minimum of the functional ${1\over 2}\sum_{j=1}^{N}|I_{\rm meas}(\mathbf{x}_{j})-|E_{\Omega}(\mathbf{x}_{j})|^{2}|^{2}.$ A scalar version of the topological methods explained in Section 3 was implemented in [6].

This strategy may be reasonable in the presence of isolated smooth scatterers, placed far enough from the detector screen. However, it introduces a number of errors of varying magnitude:

• •

First, $I_{\rm meas}(\mathbf{x}_{j})=E_{1}(\mathbf{x}_{j})^{2}+E_{2}(\mathbf{x}_{j})^{2}+E_{3}(\mathbf{x}_{j})^{3}$ where $\mathbf{E}=(E_{1},E_{2},E_{3})$ is the solution of the vector Maxwell system (12) for the true object $\Omega$. The boundary conditions couple the three components of the electric field at the boundary of the object. However, if we work with well separated objects and the detectors are placed far enough from them, the value of $E_{2}(\mathbf{x}_{j})^{2}+E_{3}(\mathbf{x}_{j})^{2}$ is expected to be small.

• •

The topological derivative for the scalar problem in the whole space is $D_{T}(\mathbf{x},\mathbb{R}^{3})={\rm Re}[(k_{e}^{2}-k_{i}^{2})E\overline{Q}]$ where $E=E_{\rm inc}$ and the scalar adjoint $\overline{Q}$ is an outgoing solution of

[TABLE]

Particularizing formula (21) for the vector Maxwell problem with polarized light and $\beta=1$, we find $D_{T}(\mathbf{x},\mathbb{R}^{3})={3k_{e}^{2}\over k_{i}^{2}+2k_{e}^{2}}{\rm Re}[(k_{e}^{2}-k_{i}^{2})E\overline{P}_{1}]$, where $\mathbf{P}=(P_{1},P_{2},P_{3})$ is given by (27). To quantify the difference between both formulas for topological derivatives, we need to relate $Q$ and $P_{1}$. We find a term of the form $k_{e}^{-2}\sum_{j=1}^{N}\nabla{\rm div}G_{k_{e}}*d_{j}\delta_{\mathbf{x}_{j}}$, where $G_{k_{e}}$ is the Green function of the Helmholtz equation. This error term decays with the distance to the detectors.

Indeed, conjugating system (27) we obtain:

[TABLE]

For $j=1,\ldots,N$, we take the divergence $\mathbf{curl}\,(\mathbf{curl}\,\overline{\mathbf{P}}_{j})-k_{e}^{2}\overline{\mathbf{P}}_{j}={\mathbf{d}_{j}}\delta_{\mathbf{x}_{j}}$, where $\mathbf{d}_{j}=2(I_{\rm meas}(\mathbf{x}_{j})-|\mathbf{E}(\mathbf{x}_{j})|^{2})\overline{\mathbf{E}(\mathbf{x}_{j})}$. Since ${\rm div}\,(\mathbf{curl}\,\mathbf{A})=0$ for any vector $\mathbf{A}$, we find ${\rm div}\,{\overline{\mathbf{P}}_{j}}=-{1\over k_{e}^{2}}{\rm div}\,{\mathbf{d}}_{j}\delta_{\mathbf{x}_{j}}.$ Making use of the vector identity $\mathbf{curl}\,(\mathbf{curl}\,\overline{\mathbf{P}}_{j})=\nabla({\rm div}\,\overline{\mathbf{P}}_{j})-\Delta{\overline{\mathbf{P}}_{j}}$ we have

[TABLE]

We can solve the equations by convolution with the Green function of Helmholtz equation:

[TABLE]

This expression quantifies the difference between the Green function of Maxwell equation and Green functions of Helmholtz equations. Indeed, notice that the right hand side can be rewritten as $G_{k_{e}}*{\mathbf{d}}_{j}\delta_{\mathbf{x}_{j}}+{1\over k_{e}^{2}}G_{k_{e}}*[\mathbf{curl}\,\mathbf{curl}\,{\mathbf{d}}_{j}\delta_{\mathbf{x}_{j}}+\Delta{\mathbf{d}}_{j}\delta_{\mathbf{x}_{j}}].$ Interchanging the derivatives in the convolution we find:

[TABLE]

for $\mathbf{x}\neq\mathbf{x}_{j}$. Summing over $j$, this implies that formula (28) for the adjoint holds pointwise as long as we do not reach the screen where the hologram is recorded.

• •

When we consider the topological derivatives in the presence of approximated objects $\Omega_{\rm ap}\neq\emptyset$, the differences between the vector and the scalar approach increase, since the expression for the topological derivative (33) becomes discontinuous across the boundary and the coupling of the field components at the boundary in (39) and (45) produces nonzero components for the forward and adjoint fields.

Appendix B Fréchet, shape and topological derivatives for holography

In this section we justify the formulas for the Fréchet, shape and topological derivatives for the hologram and the cost functional (15) employed throughout the paper.

Given an object $\Omega\subset\mathbb{R}^{3}$ with boundary $\Gamma=\partial\Omega$, we consider variations of $\Omega$ along vector fields $\boldsymbol{\xi}.$ For $\tau>0$, we introduce the family of deformations

[TABLE]

and consider the deformed domains $\varphi_{\tau\boldsymbol{\xi}}(\Omega)=\Omega_{\tau\boldsymbol{\xi}}.$ Let us quantify first the variations of the hologram under such deformations.

B.1 Fréchet derivative of the hologram

The Fréchet derivative extends the concept of differential to general spaces of infinite dimension. Given two Banach spaces $X$, $Y$ and a function ${\cal F}:D({\cal F})\subset X\longrightarrow Y$, its Fréchet derivative ${\cal F}^{\prime}:X\longrightarrow Y$ is a linear bounded operator satisfying ${\cal F}(x+\xi)={\cal F}(x)+{\cal F}^{\prime}(x)\xi+o(\xi)$ for $\xi\in X$ as $\|\xi\|_{X}\rightarrow 0$, for any $x\in X$. In other words, it satisfies ${\rm lim}_{\|\xi\|_{X}\rightarrow 0}{\|{\cal F}(x+\xi)-{\cal F}(x)-{\cal F}^{\prime}(x)\xi\|_{Y}\over\|\xi\|_{X}}=0.$ It is related to the directional Gateaux derivative

[TABLE]

by ${\cal F}^{\prime}(x)\xi=D_{\xi}{\cal F}(x).$

Let us introduce the operator that maps the deformation of the object boundary to the solution of the forward problem evaluated at the detectors:

[TABLE]

where $k\geq 2$, ${\mathbf{E}}_{\Omega_{\boldsymbol{\xi}}}$ is the solution of (12) with $\Omega={\Omega_{\boldsymbol{\xi}}},$ $\mathbf{x}_{j}$, $j=1,\ldots,N,$ are detectors placed on a screen. Then, we can write the corresponding hologram as ${\cal I}({\boldsymbol{\xi}})_{j}=|\mathbf{E}_{\Omega_{\boldsymbol{\xi}}}(\mathbf{x}_{j})|^{2}={\cal F}({\boldsymbol{\xi}})_{j}\cdot\overline{{\cal F}({\boldsymbol{\xi}})_{j}},$ $j=1,\ldots,N.$

**Theorem 1.**The Fréchet derivative of the operator ${\cal I}({\boldsymbol{\xi}})$ and its adjoint operator are given by (59) and (71), respectively.

Proof. The Gateaux derivative $\dot{\mathbf{E}}$ of the solution of (12) in the direction $\boldsymbol{\xi}$ is characterized as the solution of (65) in [12], Theorem 6.6. This provides the Fréchet derivative of the operator ${\cal F}$ defined by (105). Then,

[TABLE]

We deduce that

[TABLE]

where ${{\cal F}^{\prime}}^{*}(0)$ is characterized in [38], Proposition 6. The only difference with respect to the use of full measurements (phase+intensity) [38] is the multiplication of the incident adjoint field by $2\boldsymbol{h}$. $\square$

Using the same formula we get similar results inspired by [33, Eq. (1.3)] for the scalar case.

**Theorem 1bis. ***Using the scalar approximation for polarized waves, the Fréchet derivative of the operator ${\cal I}({\boldsymbol{\xi}})$ * is given by

[TABLE]

where $\dot{E}$ is the solution of

[TABLE]

with object $\Omega=\Omega_{\rm ap}$ and transmission data

[TABLE]

$E=E_{\Omega_{\rm ap}}$ being the solution of (A) with $\Omega=\Omega_{\rm ap}$ and $\hat{\mathbf{n}}$ the outer unit normal. The $L^{2}$-adjoint operator is defined by

[TABLE]

where ${E^{+}_{\boldsymbol{h}}}$ is the solution of the transmission problem (A) with object $\Omega=\Omega_{\rm ap}$ and incident field

[TABLE]

B.2 Shape derivative of the cost functional

When instead of the hologram we differentiate the cost functional (15) with respect to deformations along vector fields, we obtain the so-called shape derivative. Let us fix a vector field $\boldsymbol{\xi}$. For any region ${\cal R}$, the deformed domain $\varphi_{\tau}({\cal R})$ is the image of ${\cal R}$ by the deformation (100). Evaluating ${J}$ on the deformed regions, we obtain a scalar function $J(\tau)={J}(\varphi_{\tau}({\cal R}))$ of the deformation parameter $\tau$, which can be differentiated with respect to it. The shape derivative along the vector field $\boldsymbol{\xi}$ is precisely this derivative:

[TABLE]

Theorem 2. Let us assume that $\Omega$ is a $C^{2}$ domain and the coefficients $k_{e},k_{i},\beta$ are piecewise constant in the different components. Then the shape derivative $\langle DJ(\mathbb{R}^{3}\setminus\overline{\Omega}),\boldsymbol{\xi}\rangle$ of functional (15) is given by:

[TABLE]

*where $\mathbf{E}$ and $\overline{\mathbf{P}}$ are the solutions of the forward and conjugate adjoint problems (39) and (45). *

Proof. Given the vector field $\boldsymbol{\xi}$, we must evaluate (108). Differentiating with respect to $\tau$ the function $J(\mathbb{R}^{3}\setminus\overline{\phi_{\tau}(\Omega)})$ with $J$ given by (15) we find a vector version of the scalar formula in [6]:

[TABLE]

where $\mathbf{E}_{\tau}$ is a solution of the forward problem (39) with object $\Omega_{\tau}=\varphi_{\tau}(\Omega).$ Evaluating at $\tau=0$, we obtain

[TABLE]

since $\mathbf{E}_{0}=\mathbf{E}$ is the solution of (39) with object $\Omega$.

We set \dot{\mathbf{E}}={d\mathbf{E}_{\tau}\over d\tau}\Big{|}_{\tau=0}, which is the Gateaux derivative in the direction $\boldsymbol{\xi}$ of the solutions of the forward problems (39) with deformed domains $\Omega_{\tau}$. For the transmission Maxwell system it is characterized as the solution of (65), see [12]. Then, elliptic regularity and Sobolev’s embeddings ensure continuity of the solution and continuity of derivatives away from the interface $\partial\Omega$ [23, 25, 51]. Therefore, we can evaluate it at the detectors.

We rewrite the right hand side in identity (114) as follows. Integrating the equations for the adjoint field and using Green’s formula we obtain the identities:

[TABLE]

and

[TABLE]

Here, $B_{R}=B(\mathbf{0},R)$ is a sphere with radius $R$ large enough to contain the objects and detectors and $L_{k_{e}}$ stands for the Dirichlet to Neumann operator for Maxwell equations [49, Chapter 10]. To be able to use the transmission boundary conditions for $\dot{\mathbf{E}}$ we combine (117)-(120) to get:

[TABLE]

Following [36], we analyze the boundary terms in this identity using the boundary conditions satisfied by $\dot{\mathbf{E}}$, $\mathbf{E}$, $\overline{\mathbf{P}}$, together with classical vector identities and differential calculus relations on surfaces [51], to finally obtain the desired expression using (114), see [36] for details. $\square$

This Theorem is stated for piecewise constant coefficients and $C^{2}$ domains. However, piecewise $C^{1}$ coefficients and Lipschitz domains suffice [14, 29, 51].

Shape derivatives can be exploited to optimize the cost functional (15) since they provide vector fields $\boldsymbol{\xi}$ to deform an initial guess $\Omega_{\rm ap}$ in such a way that the cost functional decreases. The vector field $\boldsymbol{\xi}$ is selected to ensure that the shape derivative of (15) is negative. The process can be repeated, generating a sequence of reconstructions along which the cost functional diminishes.

When $\beta=1$, the shape derivative of the cost functional (15) at the scatterer $\Omega_{\rm ap}$ is given by

[TABLE]

where $\mathbf{E}$ and $\mathbf{P}$ are the solutions of the forward and adjoint problems (39) and (45) with object $\Omega_{\rm ap}$. Choosing $\boldsymbol{\xi}=\xi_{n}\hat{\mathbf{n}}$, $\hat{\mathbf{n}}$ being the unit normal vector on $\partial\Omega_{\rm ap},$ with

[TABLE]

we ensure $\langle DJ(\mathbb{R}^{3}\setminus\overline{\Omega}_{\rm ap}),\boldsymbol{\xi}\rangle<0$. A new approximation is defined as

[TABLE]

Then, $J(\mathbb{R}^{3}\setminus\overline{\Omega}_{\rm new})<J(\mathbb{R}^{3}\setminus\overline{\Omega}_{\rm ap})$ for $\tau>0$ small.

When enforcing star-shaped parameterizations we should fit a star-shaped parametrization to the resulting object. Alternatively, we can adapt an strategy to preserve the parametrization used in [9]. We define a set of directions associated to the parametrization:

[TABLE]

Then, the parameters defining the star-shaped boundary of the new approximation $\Omega_{\rm new}$ would be:

[TABLE]

with $\tau_{\ell},\tau_{n,m}>0$ small enough. We have tested both procedures, but the approximations stagnate without converging, as it happens for the topological derivative based iterations.

B.3 Topological derivative of the cost functional

Once an explicit formula for the shape derivative of a cost functional is available, we can obtain an explicit formula for the topological derivative taking an additional limit [20]. For every $\mathbf{x}\in{\cal R},$ we set $B_{\varepsilon}=B({\mathbf{x}},\varepsilon)$ and ${\cal R}_{\varepsilon}={\cal R}\setminus\overline{B_{\varepsilon}},$ we have [20]:

[TABLE]

for every $\mathbf{x}\in{\cal R}$ 555Reference [20] uses $-{4\pi\varepsilon^{2}}$. This changes the sign of $D_{T}({\mathbf{x}},{\cal R})$ at each point.. The vector field $\boldsymbol{\xi}$ is an extension to $\mathbb{R}^{3}$ of $\boldsymbol{\xi}=\hat{\bf n}({\mathbf{z}}),\,{\mathbf{z}}\in\Gamma_{\varepsilon}=\partial B_{\varepsilon}({\mathbf{x}}),$ where the normal $\hat{\bf n}({\bf z})={{\bf z}-{\bf x}\over|{\bf z}-{\bf x}|}$ points outside666Reference [20] and later work often select the inner normal vector to the ball instead of the exterior one, setting ${\bf V}=-{\bf n}({\mathbf{z}})$ at the boundary, which is equivalent. the ball, and vanishes out of a narrow neighborhood of $\partial B_{\varepsilon}$.

Theorem 3. Let us assume that $\Omega$ is a $C^{2}$ domain and the coefficients $k_{e},k_{i},\beta$ are constant. Then the topological derivative of functional (15) in $\mathbb{R}^{3}$ is given by (21) with forward and adjoint fields governed by (24) and (27).

Proof. Expressions of the form (21) for cost functionals involving the full wave field were first established in [47] using asymptotic expansions. We obtain it here exploiting the relation with shape derivatives (B.3), as in [6, 8, 20, 36], To find an expression for the topological derivative in terms of adjoint and forward fields we must first evaluate the shape derivative. We set here ${\cal R}=\mathbb{R}^{3}$. By Theorem 2 in Appendix B.2, $\langle DJ(\mathbb{R}^{3}\setminus\overline{B_{\varepsilon}}),\boldsymbol{\xi}\rangle$ is given by

[TABLE]

where $\mathbf{E}_{\varepsilon}$ and $\mathbf{P}_{\varepsilon}$ are the solutions of the forward and adjoint problems (39) and (45) with object $B_{\varepsilon}$. These solutions admit explicit forms given by the series expansions in C. From them, we obtain the following asymptotic behaviors when $|\boldsymbol{\chi}|=1$ as in [6, 36]:

[TABLE]

when $\varepsilon\rightarrow 0.$ To justify this asymptotic behavior, we use the series expansion (159) with coefficients (165) for $R=\varepsilon$ at $\mathbf{z}=\mathbf{x}+\varepsilon\boldsymbol{\chi},$ evaluating $\alpha_{n,m}$,$\beta_{n,m}$ from the series expansion of the incident wave (167). Since the spheres are centered at a point $\mathbf{x}$ instead of $\mathbf{0}$, $\mathbf{M}_{n,m}^{(1)}(k_{i},\mathbf{x})$ and $\mathbf{N}_{n,m}^{(1)}(k_{i},\mathbf{x})$ are replaced by $\mathbf{M}_{n,m}^{(1)}(k_{i},\mathbf{z}-\mathbf{x})$ and $\mathbf{N}_{n,m}^{(1)}(k_{i},\mathbf{z}-\mathbf{x})$ [10, 42, 36].

The behavior of the conjugate adjoint fields $\overline{\mathbf{P}}^{-}_{\varepsilon}$ is similar to (142), replacing $\mathbf{E}_{\rm inc}$ by $\overline{\mathbf{P}}_{\rm inc}.$ Combining the asymptotic behaviors (142) with the identities

[TABLE]

where $\boldsymbol{\chi}(\mathbf{z})={\mathbf{z}-\mathbf{x}\over|\mathbf{z}-\mathbf{x}|}$, $|\mathbf{z}-\mathbf{x}|=1$, and $\int_{\partial B_{\varepsilon}}a(\mathbf{z})dS_{\mathbf{z}}=\varepsilon^{2}\int_{\partial B_{1}}a(\mathbf{x}+\varepsilon\boldsymbol{\chi})dS_{\boldsymbol{\chi}}.$ we obtain expression (21) when taking the limit (B.3) using formula (137). $\square$

Theorem 4. Let us assume that $\Omega$ is a $C^{2}$ domain, the coefficients $k_{e},k_{i}$ are piecewise constant functions and the coefficient $\beta=1$. Then the topological derivative of functional (15) in $\mathbb{R}^{3}$ is given by (33) with forward and adjoint fields governed by (39) and (45).

Proof. We can adapt the arguments in either [6] for scalar holography or in [37] for full Maxwell measurements. Consider first points $\mathbf{x}\in{\cal R}=\mathbb{R}^{3}\setminus\Omega$. Since $\boldsymbol{\xi}$ vanishes on $\partial\Omega$, the only difference with the proof of Theorem 2 arises when justifying the limits (142). Let us first observe that (142) hold when $\mathbf{E}_{\rm inc}$ is replaced by a general field $\mathbf{W}_{\rm inc}$ which can be expanded as (167) and $\mathbf{E}_{\varepsilon}$ is replaced by the solution $\mathbf{W}_{\varepsilon}$ of (151) with $R=\varepsilon.$ We set $\mathbf{W}_{\rm inc}=\mathbf{E}$ equal to the solution of the forward problem (39) with object $\Omega$.

The proof is concluded by showing that to first order in $\varepsilon$ we have $\mathbf{W}_{\varepsilon}\sim\mathbf{E}_{\varepsilon}$ on $\partial B_{\varepsilon}(\mathbf{x}).$ To do so, let us quantify the difference between the solutions $\mathbf{E}_{\varepsilon}$ of (39) with object $\Omega\cup B_{\varepsilon}(\mathbf{x})$ and the total fields $\tilde{\mathbf{E}}_{\varepsilon}$:

[TABLE]

$\mathbf{W}_{\varepsilon}$ being the solution of (151) with $R=\varepsilon.$ The functions $\tilde{\mathbf{E}}_{\varepsilon}$ satisfy

[TABLE]

Notice that both $\mathbf{W}_{\rm inc}=\mathbf{E}$ and $\mathbf{W}_{\varepsilon,\rm sc}$ satisfy the transmission boundary conditions on $\partial\Omega$, the first one by construction and the second one by regularity. Now, $\mathbf{W}_{\varepsilon,\rm sc}(\mathbf{z})$ is given by a series of the form (158) setting in the coefficients $R=\varepsilon$ and using $\mathbf{M}^{(3)}_{n,m}(k_{e},\mathbf{z}-\mathbf{x})$, $\mathbf{N}^{(3)}_{n,m}(k_{e},\mathbf{z}-\mathbf{x})$ [10, 42]. Exploiting the asymptotic behavior of the spherical Bessel functions [11] as $\varepsilon\rightarrow 0$, $j_{n}(k\varepsilon)\sim(k\varepsilon)^{n}{2^{n}n!\over(2n+1)!}$, $y_{n}(k\varepsilon)\sim-(k\varepsilon)^{-(n+1)}{(2n-1)!\over 2^{n-1}(n-1)!}$, as well as the relations $h^{(1)}_{n}=j_{n}+\imath y_{n}$, and $z_{n}(k\varepsilon)^{\prime}=-z_{n+1}(k\varepsilon)+n(k\varepsilon)^{-1}z_{n}(k\varepsilon)$ for $z_{n}=j_{n},h_{n}^{(1)},$ we see that the coefficients $a_{n}(\varepsilon)$ and $b_{n}(\varepsilon)$ given by (165) behave like $\varepsilon^{2n+1}$. Therefore, the source term $(k_{e}^{2}-k_{i}^{2})\mathbf{W}_{\varepsilon,\rm sc}(\mathbf{z})=O(\varepsilon^{3})$ and $\tilde{\mathbf{E}}_{\varepsilon}\sim\mathbf{E}_{\varepsilon}$ at zero order in $\varepsilon$.

Let us consider now $\mathbf{x}\in\Omega$. The proof of identity (B.3) in [20] uses expansions of the form (18). If we revisit that proof using expansion (30) instead and setting ${\cal R}=\mathbb{R}^{3}\setminus\overline{\Omega}$, we find that D_{T}({\mathbf{x}},{\cal R})=-\lim_{\varepsilon\to 0}\frac{1}{4\pi\varepsilon^{2}}\,\frac{d}{d\tau}{J}(\varphi_{\tau}({\cal R}\cup B_{\varepsilon}))\big{|}_{\tau=0} for $\mathbf{x}\in\Omega$. Now, the roles of $k_{i}$ and $k_{e}$ are exchanged when computing the shape derivative \frac{d}{d\tau}{J}(\varphi_{\tau}({\cal R}\cup B_{\varepsilon}))\big{|}_{\tau=0} following the proof of Theorem 1. Therefore the formula for the shape derivative has the opposite sign and the two minus signs cancel each other, yielding (33) for $\beta=1.$ 777When $\beta\neq 1$, similar arguments show that the formula for the topological derivative inside $\Omega$ becomes $3\,{\rm Re}\left[-{k_{i}^{2}(k_{i}^{2}-k_{e}^{2}\beta^{-1})\over(k_{e}^{2}\beta^{-1}+2k_{i}^{2})}\mathbf{E}(\mathbf{x})\cdot\overline{\mathbf{P}}(\mathbf{x})+{1-\beta^{-1}\over 1+2\beta^{-1}}\mathbf{curl}\,\mathbf{E}(\mathbf{x})\cdot\mathbf{curl}\,\overline{\mathbf{P}}(\mathbf{x})\right].$ . The passage to the limit remains similar. $\square$

When $k_{i}$ and $k_{e}$ are piecewise $C^{1}$ in space, the formula in Theorem 2 persists, replacing $k_{i}$ and $k_{e}$ by $k_{i}(\mathbf{x})$ and $k_{e}(\mathbf{x})$. When taking limits in Theorems 3 and 4, we expand $k(\mathbf{x}+\varepsilon\boldsymbol{\xi})=k(\mathbf{x})+\varepsilon r(\mathbf{x})$, $r$ being a bounded function. Similar results hold replacing $k_{i}$ and $k_{e}$ by $k_{i}(\mathbf{x})$ and $k_{e}(\mathbf{x})$ again, revisiting the arguments in [8].

Appendix C Spherical harmonic expansions for spheres

In this Appendix, we obtain the series expansions for the forward and adjoint fields employed to calculate the topological derivative in Appendix B.3. We consider the transmission problems (88) when $\Omega=B_{R}=B(\mathbf{0},R)$ is a sphere centered at $(0,0,0)$ with radius $R>0$. The boundary $\partial\Omega=\partial B_{R}$ is the surface $|\mathbf{x}|=R$, with normal vector $\hat{\mathbf{n}}=\hat{\mathbf{x}}={\mathbf{x}\over|\mathbf{x}|}$. Both problems take the form:

[TABLE]

for specific choices of $\mathbf{W}_{\rm inc}.$ We set $\mathbb{S}^{2}=\partial B_{1}.$

Any solution of (151) admits a series expansion in terms of spherical harmonic functions [3, 11]. Different choices for this basis coexist in the literature. Initial calculations of Mie series for the scattering of plane waves by spheres used functions of the form $P_{n}^{m}(\cos(\theta))\cos(m\psi),$ $P_{n}^{m}(\cos(\theta))\sin(m\psi),$ $m=0,1,...$, $n=m,...$, where $P_{n}^{m}$ are associated Legendre polynomials [3]. For many computations it is useful to select a basis that is invariant by complex conjugation, such that $\overline{Y}_{n}^{m}={Y}_{n}^{-m}$. We will work with [11, 49]:

[TABLE]

where $r>0,\theta\in[0,\phi],\psi\in[0,2\phi]$ are the spherical coordinates of a point $\mathbf{x}\in\mathbb{R}^{3}$ and $P_{n}^{|m|}(\cos(\theta))$ are the associated Legendre polynomials for $m=-n,...,n$ and $n=0,1,2...$. These spherical functions constitute an orthonormal basis in $L^{2}(\mathbb{S}^{2})$. Denoting by $j_{n}$ the spherical Bessel functions of the first kind and by $h^{(1)}_{n}$ the spherical Hankel functions, the sets of functions

[TABLE]

are solutions of $\mathbf{curl}\,\mathbf{curl}\,\mathbf{W}-k^{2}\mathbf{W}=0$ in $\mathbb{R}^{3}$ for $m=-n,...,n$ and $n=0,1,2...$ and the functions

[TABLE]

are solutions of $\mathbf{curl}\,\mathbf{curl}\,\mathbf{W}-k^{2}\mathbf{W}=0$ in $\mathbb{R}^{3}\setminus\{\mathbf{0}\}$ satisfying the Silver-Müller radiation condition for $m=-n,...,n$ and $n=0,1,2...$. We can use $\mathbf{M}^{(1)}_{n,m}$ and $\mathbf{N}^{(1)}_{n,m}$ to seek series expansions of solutions of interior problems inside spheres and $\mathbf{M}^{(3)}_{n,m}$ and $\mathbf{N}^{(3)}_{n,m}$ for exterior problems outside. More precisely, the unique solution of (151) is given by [35]:

[TABLE]

with coefficients

[TABLE]

where $\alpha_{n,m}$ and $\beta_{n,m}$ are the coefficients of the series expansion of $\mathbf{W}_{\rm inc}:$

[TABLE]

These transfer coefficients are the same for the spherical harmonics used in [3], with a change of sign in $a_{n}$ and $b_{n}$ because of a sign difference in the coefficients employed in [3] for the expansion of the scattered waves. They are obtained as follows. We take the ${\bf curl}$ of the series (158)-(159) and use $\mathbf{curl\,M}^{(j)}=\imath k\mathbf{N}^{(j)}$, $\mathbf{curl\,N}^{(j)}=-\imath k\mathbf{M}^{(j)}$, $j=1,3$ to express them again in terms of the basis functions. We then impose the transmission conditions (151) on $|\mathbf{x}|=R$ computing the cross product with the normal vectors. Notice that on the surface of a sphere $|\mathbf{x}|=R$,

[TABLE]

where $\mathbf{curl}_{\mathbb{S}^{2}}$ and $\nabla_{\mathbb{S}^{2}}$ represent surface curl and gradients on the sphere $|\boldsymbol{\chi}|=1$. Similar identities hold for $\mathbf{M}^{(3)}$ and $\mathbf{N}^{(3)}$ replacing $j_{n}$ by $h_{n}^{(1)}.$ We observe that, in the transmission conditions, the terms involving ${\boldsymbol{\chi}}Y_{n}^{m}(\boldsymbol{\chi})$ vanish. The functions $\mathbf{curl}_{\mathbb{S}^{2}}Y_{n}^{m}$ and $\nabla_{\mathbb{S}^{2}}Y_{n}^{m}$ form an orthogonal basis in $L^{2}_{\mathbf{t}}(\mathbb{S}^{2})$ (tangent fields), satisfying $\hat{\mathbf{n}}\times\nabla_{\mathbb{S}^{2}}Y_{n}^{m}=-\mathbf{curl}_{\mathbb{S}^{2}}Y_{n}^{m}$ and $\hat{\mathbf{n}}\times\mathbf{curl}_{\mathbb{S}^{2}}Y_{n}^{m}=\nabla_{\mathbb{S}^{2}}Y_{n}^{m}.$ Setting the coefficients of the resulting expansions equal to zero, we obtain the expressions (165), where $a_{n},b_{n},c_{n},d_{n}$ are the solutions of the systems of equations:

[TABLE]

Thus, to obtain explicit expressions for the forward and adjoint fields in the presence of spheres we only need to determine the coefficients $\alpha_{n,m}$ and $\beta_{n,m}$ of the series expansion (167) of $\mathbf{E}_{\rm inc}$ and of $\overline{\mathbf{P}}_{\rm inc}$ given by (27).

C.1 Series expansion of forward fields for spheres

For a plane wave advancing in the direction of the axis $z$ and polarized along the axis $x$, that is, $\mathbf{E}_{\rm inc}(\mathbf{x})=(e^{\imath k_{e}z},0,0),$ the Mie series expansion for $\mathbf{E}_{\rm inc}$ is obtained from the expansion in [3] setting $\sin(m\theta)={e^{\imath m\theta}-e^{-\imath m\theta}\over 2\imath}$ and $\cos(m\theta)={e^{\imath m\theta}+e^{-\imath m\theta}\over 2}$:

[TABLE]

With these coefficients we obtain the corresponding Mie solution (158)-(167). The above formulas apply to spheres centered at $(0,0,0)$. If the sphere is centered at a point $\mathbf{y}$, the formulas must be conveniently shifted.

C.2 Series expansion of adjoint fields for spheres

The field $\overline{\mathbf{P}}_{\rm inc}$ given by (27) involves terms of the form

[TABLE]

where $\mathbf{d}_{j}=2(I_{\rm meas}(\mathbf{x}_{j})-|\mathbf{E}(\mathbf{x}_{j})|^{2})\overline{\mathbf{E}(\mathbf{x}_{j})},$ where $\mathbf{E}$ is the solution to (24). When the incident wave is a plane wave, the coefficients (167) are given in C.1. The Mie expansion of the functions (175) is given in [11], page 222 888In [11], the notation is $\mathbf{M}_{n,m}=\mathbf{M}_{n,m}^{(1)}$, $\mathbf{N}_{n,m}=\mathbf{M}_{n,m}^{(3)}$. . Summing over $j$:

[TABLE]

With these coefficients we obtain the corresponding Mie solution (158)-(167). Such formulas apply to spheres centered at $(0,0,0).$ When the sphere is centered at a point $\mathbf{y}$, we replace everywhere the functions $\mathbf{M}_{n,m}^{(1)}(\mathbf{x})$, $\mathbf{N}_{n,m}^{(1)}(\mathbf{x})$, $M_{n,m}^{(3)}(\mathbf{x}_{j})$, $\mathbf{N}_{n,m}^{(3)}(\mathbf{x}_{j})$ by $\mathbf{M}_{n,m}^{(1)}(\mathbf{x}-\mathbf{y})$, $\mathbf{N}_{n,m}^{(1)}(\mathbf{x}-\mathbf{y})$, $\mathbf{M}_{n,m}^{(3)}(\mathbf{x}_{j}-\mathbf{y})$, $\mathbf{N}_{n,m}^{(3)}(\mathbf{x}_{j}-\mathbf{y})$. The resulting formula holds when $|\mathbf{x}_{j}-\mathbf{y}|>|\mathbf{x}-\mathbf{y}|,$ for $j=1,...,N$. An alternative derivation follows [10, 42].