Magazine Cover
Download article (PDF)

This article is published open access.

MAG125pp. 4–12

Seeing the invisible: Digital holography

  • Ana Carpio

    Universidad Complutense de Madrid, Spain
For the past years there has been an increasing interest in developing mathematical and computational methods for digital holography. Holographic techniques furnish noninvasive tools for high-speed 3D live cell imaging. Holograms can be recorded in the millisecond or microsecond range without damaging samples. A hologram encodes the wave field scattered by an object as an interference pattern. Digital holography aims to create numerical images from digitally recorded holograms. We show here that partial differential equation constrained optimization, topological derivatives of shape functionals, iteratively regularized Gauss–Newton methods, Bayesian inference, and Markov chain Monte Carlo techniques provide effective mathematical tools to invert holographic data with quantified uncertainty. Holography set-ups are particularly challenging because a single incident wave is employed. Similar tools could be useful in inverse scattering problems involving other types of waves and different emitter/receiver configurations, such as microwave imaging or elastography, for instance.

1 Introduction

Experimental sciences have traditionally been a source of challenging mathematical problems with a double edge: while mathematical theories are created, technology moves fast and industry develops. Imaging sciences provide a remarkable example. Typical imaging systems, such as radar [28 S. W. McCandless and C. R. Jackson, Principles of synthetic aperture radar. In SAR Marine Users Manual, NOAA (2004) ], magnetic resonance tomography, ultrasound, echography [25 A. Maier, S. Steidl, V. Christlein, J. Hornegger, Medical imaging systems: An introductory guide. Springer (2018) ], and seismic imaging [34 J. Tromp, Seismic wavefield imaging of Earth’s interior across scales. Nature Reviews Earth & Environment 1, 40–53 (2020) ], pose inverse scattering problems with a similar mathematical structure. In all of them, waves generated by a set of emitters interact with a medium under study and the wave field resulting from the interaction is recorded at a set of receivers [10 D. Colton and R. Kress, Inverse acoustic and electromagnetic scattering theory. Appl. Math. Sci. 93, Springer, Berlin (1992) ]. Different imaging systems resort to different types of waves and arrange emitters and receivers according to varied geometries. The nature of the employed waves depends on factors such as the size of the specimens under study, the contrast between components, and the damage caused to the sample during the imaging procedure. Knowing the emitted and recorded waves, we aim to infer the structure of the medium.

Approximating the solutions of inverse scattering problems is a challenging task because such problems are severely ill posed [10 D. Colton and R. Kress, Inverse acoustic and electromagnetic scattering theory. Appl. Math. Sci. 93, Springer, Berlin (1992) ]. Given arbitrary data, the problem under study may not admit a solution, the solution may not be unique, or it may not depend continuously on the given data. This means that small errors may lead to a solution different from the searched one. In view of the relevant technological applications in a host of fields, such as medicine, security, geophysics, or materials testing, to mention a few, there is a need of even better mathematical techniques for classical imaging problems, as well as a need of new ideas to tackle new imaging set-ups.

We focus here on recent developments in digital holography, summarizing work done during the past 10 years in collaboration with experimentalists designing holographic microscopes. This collaboration started in 2012 thanks to the interdisciplinary communication environment created at the Harvard University’s Kavli Institute seminars. Since then, we have developed analytical and computational tools to handle inverse problems arising in digital holography, in collaboration with researchers from Harvard University and Tesla, Universidad Complutense de Madrid, Universidad Politécnica de Madrid, Universidad de Oviedo, Université de Technologie de Compiègne, and New York University.

Digital in-line holography is a noninvasive tool for accelerated three-dimensional imaging of soft matter and live cells [23 S. H. Lee, Y. Roichman, G. R. Yi, S. H. Kim, S. M. Yang, A. van Blaaderen, P. van Oostrum and D. G. Grier, Characterizing and tracking single colloidal particles with video holographic microscopy. Optics Express 15, 18275–18282 (2007) , 16 J. Fung, R. P. Perry, T. G. Dimiduk and V. N. Manoharan, Imaging multiple colloidal particles by fitting electromagnetic scattering solutions to digital holograms. J. Quant. Spectroscopy Radiative Transfer 113, 212–219 (2012) , 26 P. Marquet, B. Rappaz, P. J. Magistretti, E. Cuche, Y. Emery, T. Colomb and C. Depeursinge, Digital holographic microscopy: a noninvasive contrast imaging technique allowing quantitative visualization of living cells with subwavelength axial accuracy. Optics Letters 30, 468–478 (2005) , 37 A. Yevick, M. Hannel and D. G. Grier, Machine-learning approach to holographic particle characterization. Optics Express 22, 26884–26890 (2014) ] that achieves high spatial (nanometers) and temporal (microseconds) resolution without the need of toxic fluorescent markers or stains. In this context, a hologram is a two-dimensional light interference pattern encoding information about the optical and geometrical properties of a set of objects [35 T. Vincent, Introduction to holography. CRC Press (2012) ]. Shining a properly chosen light beam back through the hologram we can recreate the original three-dimensional image. Instead, digital holography is designed to produce numerical reconstructions of the objects in an automatic way, which amounts to solving computationally an inverse scattering problem. We will show next that optimization schemes with partial differential equation constraints, analysis of the topological derivative of objective functions, regularized Gauss–Newton iterations, and Bayesian inference are effective tools to invert holographic data in the presence of noise while quantifying uncertainty.

Figure 1. Formation of an in-line hologram. A laser beam hits an object. The scattered and undiffracted beams form an interference pattern on a screen, which is recorded at a mesh of detectors. Laser lights have wavelengths varying from about 405 nm (violet light) to about 660 nm (red light). Typical object sizes are in the micron range (1 µm = 10610^{-6} m, 1 nm = 10910^{-9} m). 🅭🅯 CC BY 4.0

2 The forward problem

The forward problem is a mathematical model of how a hologram is generated. Figure 1 illustrates how an in-line hologram is formed, though more complicated set-ups are possible. First, a laser light beam interacts with a sample. Then, interference of the scattered light field with the undiffracted beam generates the hologram on a detector screen past the object [23 S. H. Lee, Y. Roichman, G. R. Yi, S. H. Kim, S. M. Yang, A. van Blaaderen, P. van Oostrum and D. G. Grier, Characterizing and tracking single colloidal particles with video holographic microscopy. Optics Express 15, 18275–18282 (2007) ]. The light wave field obeys the Maxwell equations. Typically, the emitted laser beams are time harmonic, that is, Einc(x,t)=Re[eıωtEinc(x)]\mathcal{E}_{\mathrm{inc}}(\mathbf{x},t)=\operatorname{Re}[e^{-\imath\omega t}{\mathbf{E}}_{\mathrm{inc}}(\mathbf{x})]. The resulting wave field is also time harmonic, namely, EΩ,κ(x,t)=Re[eıωtEΩ,κ(x)]\mathcal{E}_{\Omega,\kappa}(\mathbf{x},t)=\operatorname{Re}[e^{-\imath\omega t}{\mathbf{E}}_{\Omega,\kappa}(\mathbf{x})], with complex amplitude EΩ,κ(x)\mathbf{E}_{\Omega,\kappa}(\mathbf{x}) governed by the stationary Maxwell equations. The resulting forward problem is

curl(1μecurlE)κe2μeE=0in R3Ω,curl(1μicurlE)κi2μiE=0in Ω,n^×E=n^×E+on Ω,1μin^×curlE=1μen^×curlE+on Ω,limxxcurl(EEinc)×xxıκe(EEinc)=0,\begin{gathered}\begin{aligned} \operatorname{\mathbf{curl}}\biggl(\frac{1}{\mu_{e}}\operatorname{\mathbf{curl}}\mathbf{E}\biggr)-\frac{\kappa_{e}^{2}}{\mu_{e}}\mathbf{E}&=0&&\textrm{in}\ \mathbb{R}^{3}\setminus\overline{\Omega},\\ \operatorname{\mathbf{curl}}\biggl(\frac{1}{\mu_{i}}\operatorname{\mathbf{curl}}\mathbf{E}\biggr)-\frac{\kappa_{i}^{2}}{\mu_{i}}\mathbf{E}&=0&&\textrm{in}\ \Omega,\\ \hat{\mathbf{n}}\times\mathbf{E}^{-}&=\hat{\mathbf{n}}\times\mathbf{E}^{+}&&\textrm{on}\ \partial\Omega,\\ \frac{1}{\mu_{i}}\hat{\mathbf{n}}\times\operatorname{\mathbf{curl}}\mathbf{E}^{-}&=\frac{1}{\mu_{e}}\hat{\mathbf{n}}\times\operatorname{\mathbf{curl}}\mathbf{E}^{+}&&\textrm{on}\ \partial\Omega,\end{aligned}\\ \lim_{\lvert\mathbf{x}\rvert\to\infty}\lvert\mathbf{x}\rvert\biggl|\operatorname{\mathbf{curl}}(\mathbf{E}-\mathbf{E}_{\mathrm{inc}})\times\frac{\mathbf{x}}{\lvert\mathbf{x}\rvert}-\imath\kappa_{e}(\mathbf{E}-\mathbf{E}_{\mathrm{inc}})\biggr|=0,\end{gathered}

where μi\mu_{i}, εi\varepsilon_{i} and κi=ωi2εiμi\kappa_{i}=\omega_{i}^{2}\varepsilon_{i}\mu_{i} are the permeabilities, permittivities and wavenumbers of the imaged objects Ω\Omega, while μe\mu_{e}, εe\varepsilon_{e} and κe\kappa_{e} correspond to the ambient medium [3 C. F. Borhen and D. R. Huffman, Absorption and scattering of light by small particles. Wiley Sciences, John Wiley & Sons, Berlin (1998) ] and are known. In biomedical applications, μiμeμ0\mu_{i}\sim\mu_{e}\sim\mu_{0}, μ0\mu_{0} being the vacuum permeability. The upper signs - and ++ represent limit values from inside and outside Ω\Omega, respectively, and n^\hat{\mathbf{n}} denotes the outer unit normal vector. Incident waves are polarized in a direction p^\hat{\mathbf{p}} orthogonal to the direction of propagation d^\hat{\mathbf{d}}, that is, Einc(x)=E0p^eıκed^x\mathbf{E}_{\mathrm{inc}}(\mathbf{x})=E_{0}\hat{\mathbf{p}}\,e^{\imath\kappa_{e}\hat{\mathbf{d}}\cdot\mathbf{x}}, where E0E_{0} stands for the magnitude of the incident field.

For any smooth region ΩR3Ω\Omega^{\prime}\subset\mathbb{R}^{3}\setminus\overline{\Omega} and any real κe>0\kappa_{e}>0, system (1) has a unique solution [31 J.-C. Nédélec, Acoustic and electromagnetic equations. Appl. Math. Sci. 144, Springer, New York (2001) ] in the Sobolev space H2,0(Ω)={EH2(Ω),divE=0}H^{2,0}(\Omega^{\prime})=\{\mathbf{E}\in H^{2}(\Omega^{\prime}),\operatorname{div}\mathbf{E}=0\} that is continuous in Ω\Omega^{\prime} (see [19 P. Grisvard, Elliptic problems in nonsmooth domains. Classics Appl. Math. 69, SIAM, Philadelphia (2011) ]). For collections of spheres and piecewise-constant κi\kappa_{i}, one can calculate Mie series solutions [3 C. F. Borhen and D. R. Huffman, Absorption and scattering of light by small particles. Wiley Sciences, John Wiley & Sons, Berlin (1998) ]. Starshaped object parametrizations with piecewise-constant μi\mu_{i} allow for fast spectral solvers [20 H. Harbrecht and T. Hohage, Fast methods for three-dimensional inverse obstacle scattering problems. J. Integral Equations Appl. 19, 237–260 (2007) , 24 F. Le Louër, A spectrally accurate method for the direct and inverse scattering problems by multiple 3D dielectric obstacles. ANZIAM J. 59, E1–E49 (2018) ]. Coupled BEM/FEM formulations [29 S. Meddahi, F.-J. Sayas and V. Selgás, Nonsymmetric coupling of BEM and mixed FEM on polyhedral interfaces. Math. Comp. 80, 43–68 (2011) , 31 J.-C. Nédélec, Acoustic and electromagnetic equations. Appl. Math. Sci. 144, Springer, New York (2001) ] are convenient for more general parametrizations, while discrete dipole approximations [36 A. Wang, T. G. Dimiduk, J. Fung, S. Razavi, I. Kretzschmar, K. Chaudhary and V. N. Manoharan, Using the discrete dipole approximation and holographic microscopy to measure rotational dynamics of non-spherical colloidal particles. J. Quant. Spectroscopy Radiative Transfer 146, 499–509 (2014) , 38 M. A. Yurkin and A. G. Hoekstra, The discrete-dipole-approximation code ADDA: Capabilities and known limitations. J. Quant. Spectroscopy Radiative Transfer 112, 2234–2247 (2011) ] solve the problem avoiding the use of parametrizations.

In principle, the hologram is obtained evaluating the solution of the forward problem (1) at detectors placed on the screen: IΩ,κi=Einc+Esc,Ω,κi2=EΩ,κi2I_{\Omega,\kappa_{i}}=\lvert\mathbf{E}_{\mathrm{inc}}+{\mathbf{E}}_{\mathrm{sc},\Omega,\kappa_{i}}\rvert^{2}=\lvert\mathbf{E}_{\Omega,\kappa_{i}}\rvert^{2}. In practice, the measured holograms Imeas\mathbf{I}_{\mathrm{meas}} are corrupted by noise.

3 Deterministic inverse problem

Given a hologram Imeas\mathbf{I}_{\mathrm{meas}} measured at screen points xj\mathbf{x}_{j}, j=1,,Nj=1,\ldots,N, the inverse holography problem seeks objects Ω==1LΩ\Omega=\bigcup_{\ell=1}^{L}\Omega_{\ell} and functions κi ⁣:ΩR+\kappa_{i}\colon\Omega\to\mathbb{R}^{+} such that

Imeas(xj)=EΩ,κi(xj)2,j=1,,N,I_{\mathrm{meas}}(\mathbf{x}_{j})=\lvert\mathbf{E}_{\Omega,\kappa_{i}}(\mathbf{x}_{j})\rvert^{2},\quad j=1,\ldots,N,

where EΩ,κi=Einc+Esc,Ω,κi\mathbf{E}_{\Omega,\kappa_{i}}=\mathbf{E}_{\mathrm{inc}}+\mathbf{E}_{\mathrm{sc},\Omega,\kappa_{i}} is the solution of the forward problem (1) with an object Ω\Omega and the wavenumber κi\kappa_{i} (see [5 A. Carpio, T. G. Dimiduk, F. Le Louër and M. L. Rapún, When topological derivatives met regularized Gauss–Newton iterations in holographic 3D imaging. J. Comput. Phys. 388, 224–251 (2019) ]). Since the measured data are not exact, in practice one seeks shapes Ω\Omega and functions κi\kappa_{i} for which the error between the recorded hologram and the synthetic hologram that would be generated solving (1) for the proposed objects and wavenumbers is as small as possible.

3.1 Constrained optimization

We recast the inverse problem as an optimization problem with a partial differential equation constraint: find Ω\Omega and κi\kappa_{i} minimizing the cost functional

J(Ω,κi)=12j=1NIΩ,κi(xj)Imeas(xj)2,J(\Omega,\kappa_{i})=\frac{1}{2}\sum_{j=1}^{N}\lvert I_{\Omega,\kappa_{i}}(\mathbf{x}_{j})-I_{\mathrm{meas}}(\mathbf{x}_{j})\rvert^{2},

where IΩ,κi=EΩ,κi2I_{\Omega,\kappa_{i}}=\lvert\mathbf{E}_{\Omega,\kappa_{i}}\rvert^{2} and EΩ,κi\mathbf{E}_{\Omega,\kappa_{i}} is the solution of (1). Here Ω\Omega and κi\kappa_{i} are the design variables and the stationary Maxwell system (1) is the constraint. For exact data, the true objects would be a global minimum at which the functional (2) vanishes. In general, spurious local minima may arise.

3.2 Topological derivative based approximations

A topological study of the shape functional (2) for κi\kappa_{i} fixed provides first guesses of the imaged objects without a priori information on them. The topological derivative of a shape functional [32 J. Sokołowski and A. Żochowski, On the topological derivative in shape optimization. SIAM J. Control Optim. 37, 1251–1272 (1999) ] quantifies its sensitivity to removing and including points in an object. Given a point x\mathbf{x} in a region R\mathcal{R}, we have the expansion

J(RBε(x))=J(R)+43πε3DT(x,R)+o(ε3),ε0,J(\mathcal{R}\setminus\overline{B_{\varepsilon}(\mathbf{x})})=J(\mathcal{R})+\frac{4}{3}\pi\varepsilon^{3}D_{\mathrm{T}}(\mathbf{x},\mathcal{R})+o(\varepsilon^{3}),\quad\varepsilon\to 0,

for any ball Bε(x)=B(x,ε)B_{\varepsilon}(\mathbf{x})=B(\mathbf{x},\varepsilon) centered at x\mathbf{x} with radius ε\varepsilon. The factor DT(x,R)D_{\mathrm{T}}(\mathbf{x},\mathcal{R}) is the topological derivative of the functional at x{\mathbf{x}} (see [32 J. Sokołowski and A. Żochowski, On the topological derivative in shape optimization. SIAM J. Control Optim. 37, 1251–1272 (1999) ]). If DT(x,R)D_{\mathrm{T}}(\mathbf{x},\mathcal{R}) is negative, J(RBε(x))<J(R)J(\mathcal{R}\setminus\overline{B_{\varepsilon}(\mathbf{x})})<J(\mathcal{R}) for ε>0\varepsilon>0 small. We expect the cost functional to decrease by forming objects Ωap\Omega_{\mathrm{ap}} with points below a large enough negative threshold [9 A. Carpio and M.-L. Rapún, Solving inhomogeneous inverse problems by topological derivative methods. Inverse Problems 24, Article ID 045014 (2008) , 14 G. R. Feijoo, A new method in inverse scattering based on the topological derivative. Inverse Problems 20, 1819–1840 (2004) , 27 M. Masmoudi, J. Pommier and B. Samet, The topological asymptotic expansion for the Maxwell equations and some applications. Inverse Problems 21, 547–564 (2005) ]:

Ωap{xRDT(x,R)<C0},C0>0.\Omega_{\mathrm{ap}}≔\{\mathbf{x}\in\mathcal{R}\mid D_{\mathrm{T}}(\mathbf{x},\mathcal{R})<-C_{0}\},\quad C_{0}>0.

When μe=μi\mu_{e}=\mu_{i}, R=R3\mathcal{R}=\mathbb{R}^{3} and Einc(x)=p^eıκez\mathbf{E}_{\mathrm{inc}}(\mathbf{x})=\hat{\mathbf{p}}\,e^{\imath\kappa_{e}z}, asymptotic expansions yield the formula [27 M. Masmoudi, J. Pommier and B. Samet, The topological asymptotic expansion for the Maxwell equations and some applications. Inverse Problems 21, 547–564 (2005) ]

DT(x,R3)=3Re[κe2(κe2κi2)(κi2+2κe2)E(x)P(x)],xR3,D_{\mathrm{T}}(\mathbf{x},\mathbb{R}^{3})=3\operatorname{Re}\biggl[\frac{\kappa_{e}^{2}(\kappa_{e}^{2}-\kappa_{i}^{2})}{(\kappa_{i}^{2}+2\kappa_{e}^{2})}\mathbf{E}(\mathbf{x})\cdot\overline{\mathbf{P}}(\mathbf{x})\biggr],\quad\mathbf{x}\in\mathbb{R}^{3},

where E=Einc\mathbf{E}=\mathbf{E}_{\mathrm{inc}} and

P(x)=j=1Ncurlcurl(2κe2Gκe(xxj)(Imeas(xj)Einc(xj)2)Einc(xj))\overline{\mathbf{P}}(\mathbf{x})=\sum_{j=1}^{N}\operatorname{\mathbf{curl}}\operatorname{\mathbf{curl}}\biggl(\frac{2}{\kappa_{e}^{2}}G_{\kappa_{e}}(\mathbf{x}-\mathbf{x}_{j})\bigl(I_{\mathrm{meas}}(\mathbf{x}_{j})-\lvert\mathbf{E}_{\mathrm{inc}}(\mathbf{x}_{j})\rvert^{2}\bigr)\overline{\mathbf{E}_{\mathrm{inc}}(\mathbf{x}_{j})}\biggr)

with Gκe(x)=14πxeıκexG_{\kappa_{e}}(\mathbf{x})=\frac{1}{4\pi\lvert\mathbf{x}\rvert}e^{\imath\kappa_{e}\lvert\mathbf{x}\rvert} denoting the outgoing Green function of the Helmholtz equation [31 J.-C. Nédélec, Acoustic and electromagnetic equations. Appl. Math. Sci. 144, Springer, New York (2001) ]. Once Ωap\Omega_{\mathrm{ap}} is constructed, we fit aparametrized contour qap\mathbf{q}_{\mathrm{ap}} to its boundary. Starshaped parametrizations are typical choices. Figure 2 exemplifies the procedure. The method is robust to noise, in the sense that perturbations of the data with random 10 % or 20 % noise, for instance, produce similar results. Notice that the value of κi\kappa_{i} enters through a factor that we may scale out in (5) and it is not really needed to localize the object. Similar results are obtained using the topological energy [6 A. Carpio, T. G. Dimiduk, M. L. Rapún and V. Selgas, Noninvasive imaging of three-dimensional micro and nanostructures by topological methods. SIAM J. Imaging Sci. 9, 1324–1354 (2016) ]

ET(x,R3)=E(x)2P(x)2,E_{\mathrm{T}}(\mathbf{x},\mathbb{R}^{3})=\lvert\mathbf{E}(\mathbf{x})\rvert^{2}\lvert\mathbf{P}(\mathbf{x})\rvert^{2},

which does not involve κi\kappa_{i} at all. No knowledge of κi\kappa_{i} is needed to construct a first guess of the objects.

Figure 2. Slice y=0y=0 of the topological derivative computed using expression (5) for holographic data Imeas\mathbf{I}_{\mathrm{meas}} corresponding to a sphere of radius 0.45 µm illuminated by polarized light of wavelength 520 nm and placed at a distance 28 µm of a CMOS screen. Axis units are microns. The red contour marks the location of the true object, while the cyan contour represents the approximation. Redrawn from [5]. 🅭🅯 CC BY 4.0

3.3 Regularized Gauss–Newton iterations

Fast methods to improve our knowledge of the objects starting from an initial guess are based on the following result. Let us consider two Hilbert spaces XX, YY and a Fréchet differentiable operator F ⁣:D(F)XY\mathcal{F}\colon D(\mathcal{F})\subset X\to Y. Assuming that the exact data yYy\in Y are attainable (that is, there is xXx\in X such that F(x)=y\mathcal{F}(x)=y), but only noisy data yδy^{\delta} verifying yδyYδ\lVert y^{\delta}-y\rVert_{Y}\leq\delta are accessible, the iteratively regularized Gauss–Newton (IRGN) method [1 A. B. Bakushinskii, On a convergence problem of the iterative-regularized Gauss–Newton method. Zh. Vychisl. Mat. i Mat. Fiz. 32, 1503–1509 (1992) ] constructs a sequence xk+1δx_{k+1}^{\delta} as follows. We linearize the equation at xkδx_{k}^{\delta} at each step, approximate the solution of F(xkδ)+F(xkδ)ξ=yδ\mathcal{F}(x_{k}^{\delta})+\mathcal{F}^{\prime}(x_{k}^{\delta})\xi=y^{\delta} through the minimization problem

ξk+1=ArgminξXF(xkδ)+F(xkδ)ξyδY2+αkxkδ+ξx0X2\begin{split}\xi_{k+1}=\operatorname*{Argmin}_{\xi\in X}&\lVert\mathcal{F}(x_{k}^{\delta})+\mathcal{F}^{\prime}(x_{k}^{\delta})\xi-y^{\delta}\rVert_{Y}^{2}\\[-8.6pt] &\quad+\alpha_{k}\lVert x_{k}^{\delta}+\xi-x_{0}\rVert_{X}^{2}\end{split}

and set xk+1δ=xkδ+ξk+1x_{k+1}^{\delta}=x_{k}^{\delta}+\xi_{k+1}. The Tikhonov term αkxkδ+ξx0X2\alpha_{k}\lVert x_{k}^{\delta}+\xi-x_{0}\rVert_{X}^{2} has regularizing properties and promotes convergence for specific choices x0x_{0} and αk\alpha_{k} (see [21 T. Hohage, Logarithmic convergence rates of the iteratively regularized Gauss–Newton method for an inverse potential and an inverse scattering problem. Inverse Problems 13, 1279–1299 (1997) ]). The theory of linear Tikhonov regularization guarantees that

ξk+1=(F(xkδ)F(xkδ)+αkI)1[F(xkδ)(F(xkδ)yδ)+αk(xkδx0)],\begin{split}\xi_{k+1}=-\bigl(\mathcal{F}^{\prime}(x_{k}^{\delta})^{*}\mathcal{F}^{\prime}(x_{k}^{\delta})+\alpha_{k}I\bigr)^{-1}[&\mathcal{F}^{\prime}(x_{k}^{\delta})^{*}(\mathcal{F}(x_{k}^{\delta})-y^{\delta})\\[-3.0pt] &\quad+\alpha_{k}(x_{k}^{\delta}-x_{0})],\end{split}

where F(xkδ)\mathcal{F}^{\prime}(x_{k}^{\delta})^{*} denotes the adjoint of the Fréchet derivative F(xkδ)\mathcal{F}^{\prime}(x_{k}^{\delta}). The noise level δ\delta affects the stopping criterion, the so-called discrepancy principle.

In a holography set-up, the map F\mathcal{F} is the operator that to each parametrization of objects q\mathbf{q} assigns the synthetic hologram I(q)\mathbf{I}(\mathbf{q}) generated by solving the forward problem for those objects. Starshaped parametrizations are a standard choice for simple objects. They describe each object by a few parameters: its center and a radius function represented by a finite combination of spherical harmonics [5 A. Carpio, T. G. Dimiduk, F. Le Louër and M. L. Rapún, When topological derivatives met regularized Gauss–Newton iterations in holographic 3D imaging. J. Comput. Phys. 388, 224–251 (2019) , 20 H. Harbrecht and T. Hohage, Fast methods for three-dimensional inverse obstacle scattering problems. J. Integral Equations Appl. 19, 237–260 (2007) ]. Given a starshaped parametrization qk\mathbf{q}_{k} and a recorded hologram Imeas\mathbf{I}_{\mathrm{meas}} with a level of noise δ\delta, the IRGN method first solves the linearized equation

I(qk)+I(qk)ξ=Imeas\mathbf{I}(\mathbf{q}_{k})+\mathbf{I}^{\prime}(\mathbf{q}_{k})\mathbb{\xi}=\mathbf{I}_{\mathrm{meas}}

by addressing the nonlinear least squares problem

ξk+1=Argminξ{ImeasI(qk)I(qk)ξ22+αkqk+ξqapHs(S2)2},\begin{split}\mathbb{\xi}_{k+1}=\operatorname{Argmin}_{\mathbb{\xi}}\bigl\{&\lVert\mathbf{I}_{\mathrm{meas}}-\mathbf{I}(\mathbf{q}_{k})-\mathbf{I}^{\prime}(\mathbf{q}_{k})\mathbb{\xi}\rVert^{2}_{2}\\ &\quad+\alpha_{k}\lVert\mathbf{q}_{k}+\mathbb{\xi}-\mathbf{q}_{\mathrm{ap}}\rVert^{2}_{H^{s}(\mathbb{S}^{2})}\bigr\},\end{split}

where Hs(S2)H^{s}(\mathbb{S}^{2}), s>0s>0, is an adequate Sobolev space [5 A. Carpio, T. G. Dimiduk, F. Le Louër and M. L. Rapún, When topological derivatives met regularized Gauss–Newton iterations in holographic 3D imaging. J. Comput. Phys. 388, 224–251 (2019) ], and then sets qk+1=qk+ξk+1\mathbf{q}_{k+1}=\mathbf{q}_{k}+\mathbb{\xi}_{k+1}. The initial parametrization q0=qap\mathbf{q}_{0}=\mathbf{q}_{\mathrm{ap}} represents the first guess of the objects constructed by topological methods. The updated objects Ωk\Omega_{k} correspond to the parametrizations qk\mathbf{q}_{k}. The stopping criterion for the noise level δ\delta is as follows. If the synthetic hologram calculated numerically for the current approximation of the objects I(qk)\mathbf{I}(\mathbf{q}_{k}) satisfies

I(qk)Imeas2τδ,\lVert\mathbf{I}(\mathbf{q}_{k})-\mathbf{I}_{\mathrm{meas}}\rVert_{2}\leq\tau\delta,

we stop the algorithm, τ>0\tau>0 being a parameter adjusted to guarantee a reasonable approximation while preventing early stops.

Figure 3. For the hologram in Figure 4 (redrawn from [5]): (a) True geometry. (b) Slice x=5x=5 of the topological derivative (5). Red contours are true objects. (c) Initial guess defined by (4). (d) Approximation after 4 steps of the IRGN method. (e) Slice of the topological derivative (6) at step 4. Cyan regions are approximate objects. (f) Approximation after creating an object at step 4 and applying once the IRGNM. (g) Slice of the topological derivative (6) at step 5. (h) Approximation after creating an object at step 5 and applying once the IRGNM. (i) Final approximation. Axis units are µm. (j) Decrease in the cost. 🅭🅯 CC BY 4.0

Figure 4. Hologram generated by the three objects represented in Figure 3(a), obtained with violet light having a wavelength of 405 nm emitted at z=0z=0 and recorded at z=10z=10. Axis units are microns. Redrawn from [5]. 🅭🅯 CC BY 4.0

Figures 3 and 4 illustrate the process. Figure 4 depicts the hologram generated by the configuration with three objects shown in Figure 3 (a). We use the topological derivative (5) to spot a first dominant object at the top and locate an object there, see panel (b). Then we apply the IRGN method, see panels (c) and (d). At step 4 the cost functional, depicted in panel (j), stagnates without fulfilling the stopping criteria. This suggests that more objects should be created. This can be done by hybrid methods, as we explain next.

3.4 Topologically informed IRGN methods

Approaches that use initial object parametrization as reference have a drawback: the initial guess of the number of objects may be wrong. To overcome it, we have developed hybrid algorithms [5 A. Carpio, T. G. Dimiduk, F. Le Louër and M. L. Rapún, When topological derivatives met regularized Gauss–Newton iterations in holographic 3D imaging. J. Comput. Phys. 388, 224–251 (2019) ] combining topological derivatives and regularized Gauss–Newton iterations [5 A. Carpio, T. G. Dimiduk, F. Le Louër and M. L. Rapún, When topological derivatives met regularized Gauss–Newton iterations in holographic 3D imaging. J. Comput. Phys. 388, 224–251 (2019) ]. We fit an initial parametrization qap\mathbf{q}_{\mathrm{ap}} to the first guess of the objects constructed by topological methods. Then, we apply the IRGN method and check that the cost (2) decreases. When the cost stagnates without fulfilling the stopping criteria, we reset Ωap\Omega_{\mathrm{ap}} equal to the current guess of the objects Ωk\Omega_{k} for the last parametrization obtained qk\mathbf{q}_{k} and calculate the topological derivative of the cost for R=R3Ωap\mathcal{R}=\mathbb{R}^{3}\setminus\overline{\Omega_{\mathrm{ap}}}. This is given by (3) if xR=R3Ω\mathbf{x}\in\mathcal{R}=\mathbb{R}^{3}\setminus\overline{\Omega} and its equivalent

J(R3(ΩBε(x)))=J((R3Bε(x))Ω)=J(R3Ω)43πε3DT(x,R3Ω)+o(ε3)\begin{split}&J\bigl(\mathbb{R}^{3}\setminus(\overline{\Omega\setminus\overline{B_{\varepsilon}(\mathbf{x})}})\bigr)\\ &\qquad=J\bigl((\mathbb{R}^{3}\cup B_{\varepsilon}(\mathbf{x}))\setminus\overline{\Omega}\bigr)\\[-3.0pt] &\qquad=J(\mathbb{R}^{3}\setminus\overline{\Omega})-\frac{4}{3}\pi\varepsilon^{3}D_{\mathrm{T}}(\mathbf{x},\mathbb{R}^{3}\setminus\overline{\Omega})+o(\varepsilon^{3})\end{split}

if xΩ\mathbf{x}\in\Omega. Asymptotic calculations yield the formula [9 A. Carpio and M.-L. Rapún, Solving inhomogeneous inverse problems by topological derivative methods. Inverse Problems 24, Article ID 045014 (2008) , 5 A. Carpio, T. G. Dimiduk, F. Le Louër and M. L. Rapún, When topological derivatives met regularized Gauss–Newton iterations in holographic 3D imaging. J. Comput. Phys. 388, 224–251 (2019) ]

DT(x,R3Ω)={3Re[κe2(κe2κi2)(κi2+2κe2)E(x)P(x)],xR3Ω,3Re[κi2(κe2κi2)(κe2+2κi2)E(x)P(x)],xΩ,\begin{split}&D_{\mathrm{T}}(\mathbf{x},\mathbb{R}^{3}\setminus\overline{\Omega})\\ &\quad=\begin{cases}3\operatorname{Re}\biggl[\dfrac{\kappa_{e}^{2}(\kappa_{e}^{2}-\kappa_{i}^{2})}{(\kappa_{i}^{2}+2\kappa_{e}^{2})}\mathbf{E}(\mathbf{x})\cdot\overline{\mathbf{P}}(\mathbf{x})\biggr],&\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Omega},\\[6.45pt] 3\operatorname{Re}\biggl[\dfrac{\kappa_{i}^{2}(\kappa_{e}^{2}-\kappa_{i}^{2})}{(\kappa_{e}^{2}+2\kappa_{i}^{2})}\mathbf{E}(\mathbf{x})\cdot\overline{\mathbf{P}}(\mathbf{x})\biggr],&\mathbf{x}\in\Omega,\end{cases}\end{split}

when μe=μi\mu_{e}=\mu_{i}, with forward and conjugate adjoint fields satisfying transmission Maxwell problems with object Ω=Ωap\Omega=\Omega_{\mathrm{ap}}:

curl(curlE)κe2E=0in R3Ω,curl(curlE)κi2E=0in Ω,n^×E=n^×E+on Ω,n^×curlE=n^×curlE+on Ω,\displaystyle\begin{aligned} \operatorname{\mathbf{curl}}(\operatorname{\mathbf{curl}}\mathbf{E})-\kappa_{e}^{2}\mathbf{E}&=0&&\textrm{in}\ \mathbb{R}^{3}\setminus\overline{\Omega},\\ \operatorname{\mathbf{curl}}(\operatorname{\mathbf{curl}}\mathbf{E})-\kappa_{i}^{2}\mathbf{E}&=0&&\textrm{in}\ \Omega,\\ \hat{\mathbf{n}}\times\mathbf{E}^{-}&=\hat{\mathbf{n}}\times\mathbf{E}^{+}&&\textrm{on}\ \partial\Omega,\\ \hat{\mathbf{n}}\times\operatorname{\mathbf{curl}}\mathbf{E}^{-}&=\hat{\mathbf{n}}\times\operatorname{\mathbf{curl}}\mathbf{E}^{+}&&\textrm{on}\ \partial\Omega,\end{aligned}
limxxcurl(EEinc)×x^ıκe(EEinc)=0,\displaystyle\lim_{\lvert\mathbf{x}\rvert\to\infty}\lvert\mathbf{x}\rvert\lvert\operatorname{\mathbf{curl}}(\mathbf{E}-\mathbf{E}_{\mathrm{inc}})\times\hat{\mathbf{x}}-\imath\kappa_{e}(\mathbf{E}-\mathbf{E}_{\mathrm{inc}})\rvert=0,
curl(curlP)κe2P=2j=1N(ImeasE2)Eδxjin R3Ω,curl(curlP)κi2P=0in Ω,n^×P=n^×P+on Ω,n^×curlP=n^×curlP+on Ω,\displaystyle\begin{aligned} \operatorname{\mathbf{curl}}(\operatorname{\mathbf{curl}}\overline{\mathbf{P}})-\kappa_{e}^{2}\overline{\mathbf{P}}&=2\sum_{j=1}^{N}(\mathbf{I}_{\mathrm{meas}}-\lvert\mathbf{E}\rvert^{2})\overline{\mathbf{E}}\delta_{\mathbf{x}_{j}}&&\textrm{in}\ \mathbb{R}^{3}\setminus\overline{\Omega},\\ \operatorname{\mathbf{curl}}(\operatorname{\mathbf{curl}}\overline{\mathbf{P}})-\kappa_{i}^{2}\overline{\mathbf{P}}&=0&&\textrm{in}\ \Omega,\\ \hat{\mathbf{n}}\times\overline{\mathbf{P}}^{-}&=\hat{\mathbf{n}}\times\overline{\mathbf{P}}^{+}&&\textrm{on}\ \partial\Omega,\\ \hat{\mathbf{n}}\times\operatorname{\mathbf{curl}}\overline{\mathbf{P}}^{-}&=\hat{\mathbf{n}}\times\operatorname{\mathbf{curl}}\overline{\mathbf{P}}^{+}&&\textrm{on}\ \partial\Omega,\end{aligned}
limxxcurlP×x^ıκeP=0,\displaystyle\lim_{\lvert\mathbf{x}\rvert\to\infty}\lvert\mathbf{x}\rvert\lvert\operatorname{\mathbf{curl}}\overline{\mathbf{P}}\times\hat{\mathbf{x}}-\imath\kappa_{e}\overline{\mathbf{P}}\rvert=0,

where n^\hat{\mathbf{n}} is the unit outer normal, x^=x/x\hat{\mathbf{x}}={\mathbf{x}}/\lvert\mathbf{x}\rvert and δxj\delta_{\mathbf{x}_{j}} are Dirac masses concentrated at the detectors xj\mathbf{x}_{j}, j=1,...,Nj=1,...,N.

We create a new approximation Ωnew\Omega_{\mathrm{new}} from Ωap\Omega_{\mathrm{ap}} by removing the points in Ωap\Omega_{\mathrm{ap}} at which the topological derivate surpasses a positive threshold cnewc_{\mathrm{new}} and adding the points outside Ωap\Omega_{\mathrm{ap}} at which the topological derivate falls below a negative threshold Cnew-C_{\mathrm{new}}, see [9 A. Carpio and M.-L. Rapún, Solving inhomogeneous inverse problems by topological derivative methods. Inverse Problems 24, Article ID 045014 (2008) , 6 A. Carpio, T. G. Dimiduk, M. L. Rapún and V. Selgas, Noninvasive imaging of three-dimensional micro and nanostructures by topological methods. SIAM J. Imaging Sci. 9, 1324–1354 (2016) ]:

Ωnew{xΩapDT(x,R3Ωap)<cnew}{xR3ΩapDT(x,R3Ωap)<Cnew}.\begin{split}\Omega_{\mathrm{new}}&≔\{\mathbf{x}\in\Omega_{\mathrm{ap}}\mid D_{\mathrm{T}}(\mathbf{x},\mathbb{R}^{3}\setminus\overline{\Omega}_{\mathrm{ap}})<c_{\mathrm{new}}\}\\ &\qquad\cup\{\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Omega}_{\mathrm{ap}}\mid D_{\mathrm{T}}(\mathbf{x},\mathbb{R}^{3}\setminus\overline{\Omega}_{\mathrm{ap}})<-C_{\mathrm{new}}\}.\end{split}

The constants CnewC_{\mathrm{new}}, cnewc_{\mathrm{new}} are selected to ensure a decrease in the cost functional (2) keeping κi\kappa_{i} fixed. Once Ωnew\Omega_{\mathrm{new}} is constructed, we fit a parametrization qnew\mathbf{q}_{\mathrm{new}} to its contour and restart the IRGN procedure for qap=qnew\mathbf{q}_{\mathrm{ap}}=\mathbf{q}_{\mathrm{new}}. The procedure stops when the changes in the cost and the parametrizations fall below selected thresholds.

Let us revisit the example studied in Figures 3 and 4. At step 4 of the IRGN method the cost stagnates without fulfilling the stopping criteria. We calculate the topological derivative (6) of the cost for the current approximation of the objects, illustrated in Figure 3 (e). A new region where the topological derivative attains large negative values appears. We create a new object there and update the parametrization, see panel (f). Then we apply the IRGN method again. Since the cost functional still stagnates without fulfilling the stopping criteria, we recalculate the topological derivative (6) for the available object approximation. Panel (g) suggests the creation of a third object. We update the IRGN method using this new configuration, and evolve the resulting object configuration, represented in panel (h), until the stopping criterion is met at panel (i) after 24 steps. Panel (j) illustrates stagnation and decrease of the cost as new objects are added to the parametrization using topological information and the updated IRGN method evolves, in a logarithmic scale. These simulations assume κi\kappa_{i} known and fixed. Once first guesses for κi\kappa_{i} are available, we can implement this procedure considering constant values for κi\kappa_{i} at each component of the parametrization. Obtaining first guesses for κi\kappa_{i} that are reliable enough is a hard task [7 A. Carpio, T. G. Dimiduk, V. Selgas and P. Vidal, Optimization methods for in-line holography. SIAM J. Imaging Sci. 11, 923–956 (2018) ] and the optimization procedure can encounter difficulties. Bayesian approaches provide alternative procedures that can handle these difficulties while quantifying uncertainty associated to noise and missing information.

4 Bayesian inverse problem

Bayesian formulations consider all unknowns in the inverse problem as random variables. Given a recorded hologram Imeas\mathbf{I}_{\mathrm{meas}}, we seek a finite-dimensional vector of parameters ν\mathbb{\nu} characterizing the imaged objects. When we assume the presence of LL objects, ν\mathbb{\nu} is formed by LL blocks, one per object. Using Bayes’ formula [22 J. Kaipio and E. Somersalo, Statistical and computational inverse problems. Appl. Math. Sci. 160, Springer, New York (2005) , 33 A. Tarantola, Inverse problem theory and methods for model parameter estimation. Society for Industrial and Applied Mathematics (SIAM), Philadelphia (2005) ]

ppt(ν)p(νImeas)=p(Imeasν)p(Imeas)ppr(ν),p_{\mathrm{pt}}(\mathbb{\nu})≔p(\mathbb{\nu}|\mathbf{I}_{\mathrm{meas}})=\frac{p(\mathbf{I}_{\mathrm{meas}}|\mathbb{\nu})}{p(\mathbf{I}_{\mathrm{meas}})}p_{\mathrm{pr}}(\mathbb{\nu}),

where ppr(ν)p_{\mathrm{pr}}(\mathbb{\nu}) represents the prior probability of the variables, which incorporates our previous knowledge on them, while p(Imeasν)p(\mathbf{I}_{\mathrm{meas}}|\mathbb{\nu}) is the conditional probability or likelihood of observing Imeas\mathbf{I}_{\mathrm{meas}} given ν\mathbb{\nu}. The solution of the Bayesian inverse problem is the posterior probability ppt(νImeas)p_{\mathrm{pt}}(\mathbb{\nu}|\mathbf{I}_{\mathrm{meas}}) of the parameters given the data. Sampling the posterior distribution, we obtain statistical information on the most likely values of the object parameters with quantified uncertainty.

4.1 Likelihood choice

Assuming additive Gaussian measurement noise, the measured hologram and the synthetic hologram obtained for the true object parameters are related by Imeas=I(νtrue)+ε\mathbf{I}_{\mathrm{meas}}=\mathbf{I}(\mathbb{\nu}_{\mathrm{true}})+\mathbb{\varepsilon}, where the measurement noise ε\mathbb{\varepsilon} is distributed as a multivariate Gaussian N(0,Γn)\mathcal{N}(0,\mathbb{\Gamma}_{\mathrm{n}}) with zero mean and covariance matrix Γn\mathbb{\Gamma}_{\mathrm{n}}. A possible choice for the likelihood p(Imeasν)p(\mathbf{I}_{\mathrm{meas}}|\mathbb{\nu}) is [8 A. Carpio, S. Iakunin and G. Stadler, Bayesian approach to inverse scattering with topological priors. Inverse Problems 36, Article ID 105001 (2020) ]

p(Imeasν)=1(2π)N/2Γnexp(12I(ν)ImeasΓn12)p(\mathbf{I}_{\mathrm{meas}}|\mathbb{\nu})=\frac{1}{(2\pi)^{N/2}\sqrt{\lvert\mathbb{\Gamma}_{\mathrm{n}}\rvert}}\exp\Bigl(-\frac{1}{2}\lVert\mathbf{I}(\mathbb{\nu})-\mathbf{I}_{\mathrm{meas}}\rVert^{2}_{\mathbb{\Gamma}_{\mathrm{n}}^{\smash{-1}}}\Bigr)

with vΓn12=vtΓn1v\lVert\mathbf{v}\rVert_{\mathbb{\Gamma}_{\mathrm{n}}^{-1}}^{2}=\overline{\mathbf{v}}^{\mathrm{t}}\mathbb{\Gamma}_{\mathrm{n}}^{-1}\mathbf{v}. Here, I(ν)\mathbf{I}(\mathbb{\nu}) represents the synthetic hologram obtained solving the forward problem (1) for objects characterized by parameters ν\mathbb{\nu}, see Section 2.

4.2 Topological priors

A typical choice for the prior distribution is a multivariate Gaussian

ppr(ν)=1(2π)n/21Γprexp(12(νν0)tΓpr1(νν0))p_{\mathrm{pr}}(\mathbb{\nu})=\frac{1}{(2\pi)^{n/2}}\frac{1}{\sqrt{\lvert\mathbb{\Gamma}_{\smash[b]{\mathrm{pr}}}\rvert}}\exp\Bigl(-\frac{1}{2}(\mathbb{\nu}-\mathbb{\nu}_{0})^{t}\mathbb{\Gamma}_{\mathrm{pr}}^{-1}(\mathbb{\nu}-\mathbb{\nu}_{0})\Bigr)

if ν\mathbb{\nu} is “admissible”, and ppr(ν)=0p_{\mathrm{pr}}(\mathbb{\nu})=0 when ν\mathbb{\nu} is “not admissible”, that is, it does not satisfy known constraints on the parameter set, see [8 A. Carpio, S. Iakunin and G. Stadler, Bayesian approach to inverse scattering with topological priors. Inverse Problems 36, Article ID 105001 (2020) ] for details. Here, Γpr\mathbb{\Gamma}_{\mathrm{pr}} is the covariance matrix and nn is the total number of parameters characterizing the objects. The mean ν0\mathbb{\nu}_{0} is typically a set of parameter values characterizing an initial guess of the objects. Sharp priors are obtained fitting parametrizations to first guesses of the objects obtained from the study of topological fields associated to deterministic shape costs, as explained in Section 3.2.

4.3 Markov chain Monte Carlo sampling

Combining (7), (8) and (9), the posterior probability becomes (neglecting normalization constants)

ppt(ν)exp(12I(ν)ImeasΓn1212νν0Γpr12)p_{\mathrm{pt}}(\mathbb{\nu})\propto\exp\Bigl(-\frac{1}{2}\lVert\mathbf{I}(\mathbb{\nu})-\mathbf{I}_{\mathrm{meas}}\rVert^{2}_{\mathbb{\Gamma}_{\mathrm{n}}^{\smash{-1}}}-\frac{1}{2}\lVert\mathbb{\nu}-\mathbb{\nu}_{0}\rVert_{\mathbb{\Gamma}_{\mathrm{pr}}^{\smash{-1}}}^{2}\Bigr)

when ν\mathbb{\nu} is admissible, and ppt(ν)=0p_{\mathrm{pt}}(\mathbb{\nu})=0 otherwise. Markov chain Monte Carlo (MCMC) methods provide tools to sample unnormalized posteriors. Classical MCMC methods, such as Hamiltonian Monte Carlo or Metropolis–Hastings [30 R. M. Neal, MCMC using Hamiltonian dynamics. In Handbook of Markov Chain Monte Carlo, edited by S. Brooks, A. Gelman, G. L. Jones and X. L. Meng, Chapman & Hall/CRC Handb. Mod. Stat. Methods, CRC Press, Boca Raton, 113–162 (2011) ] construct a chain of nn-dimensional states ν(0)ν(1)ν(k)\mathbb{\nu}^{(0)}\to\mathbb{\nu}^{(1)}\to\cdots\to\mathbb{\nu}^{(k)}\to\nobreak\cdots which evolve to be distributed in accordance with the target distribution ppt(ν)p_{\mathrm{pt}}(\mathbb{\nu}). After sampling an initial state ν(0)\mathbb{\nu}^{(0)} from the prior distribution (9), the chain advances from one state ν(k)\mathbb{\nu}^{(k)} to the next ν(k+1)\mathbb{\nu}^{(k+1)} by means of a transition operator that varies with the method employed [30 R. M. Neal, MCMC using Hamiltonian dynamics. In Handbook of Markov Chain Monte Carlo, edited by S. Brooks, A. Gelman, G. L. Jones and X. L. Meng, Chapman & Hall/CRC Handb. Mod. Stat. Methods, CRC Press, Boca Raton, 113–162 (2011) ]. More recent ensemble MCMC samplers [13 M. M. Dunlop and G. Stadler, A gradient-free subspace-adjusting ensemble sampler for infinite-dimensional Bayesian inverse problems, preprint, arXiv:2202.11088v1 (2022) , 18 J. Goodman and J. Weare, Ensemble samplers with affine invariance. Commun. Appl. Math. Comput. Sci. 5, 65–80 (2010) ] draw WW initial states from the prior distribution (the “walkers” or “particles”) and transition to new states while mixing the previous ones to generate several chains. This approach allows for parallelization and can handle multimodal posteriors [8 A. Carpio, S. Iakunin and G. Stadler, Bayesian approach to inverse scattering with topological priors. Inverse Problems 36, Article ID 105001 (2020) ].

Figure 5. For the two-dimensional object depicted in red, with violet light having a wavelength of 405 nm emitted at y=5y=-5 and recorded at y=5y=5, we present statistical information obtained from the samples generated by MCMC sampling. A contour projection of a two-dimensional histogram represents the probability of belonging to the object, compared to the contours of the true object, the prior mean and the MAP (maximum a posteriori) approximation. The probability is multimodal, as evidenced by additional two-dimensional histograms representing the probability of being the center of mass of the object, the most likely values for the largest/smallest object radii and their orientation, the most likely areas, and the deviation from a spherical shape. Two main peaks are identified. The main mode corresponds to a majority of samples wrapping around the object, while the second mode represents the contribution of additional large samples elongated in the direction of incidence of the incoming wave and represents an aberration of this imaging system, which uses only one incident wave. This effect may not be observed for other shapes, sizes, or light wavelengths – it depends on the geometry. Axis units are microns. Redrawn from [8]. 🅭🅯 CC BY 4.0

Figure 5 illustrates the results in a two-dimensional geometry, to reduce the computational cost in the tests. A few million samples were generated, which requires solving an identical number of forward problems. In two-dimensional set-ups we replace the stationary transmission problem for the Maxwell equations by a transmission problem posed for the Helmholtz equation [8 A. Carpio, S. Iakunin and G. Stadler, Bayesian approach to inverse scattering with topological priors. Inverse Problems 36, Article ID 105001 (2020) ]. Assuming κi\kappa_{i} is piecewise constant, we resort to fast boundary elements to solve the Hemholtz transmission problems in two dimensions [12 V. Domínguez, S. Lu and F.-J. Sayas, A fully discrete Calderón calculus for two dimensional time harmonic waves. Int. J. Numer. Anal. Model. 11, 332–345 (2014) ]. Once a large enough collection of samples is generated [17 A. Gelman and D. B. Rubin, Inference from iterative simulation using multiple sequences. Statist. Sci. 7, 457–472 (1992) ], we extract statistical information describing the imaged object: the most likely shapes, sizes, locations, as well as uncertainty in the predictions. While starshaped two-dimensional objects can be reasonably characterized with 10–20 parameters, three-dimensional objects require 80–90. Full characterization of the posterior probability by MCMC sampling becomes more expensive as the number of parameters and the time required to solve forward problems increase.

4.4 Laplace approximation

The full characterization of the posterior probability is a challenging and costly probability problem for moderate- and high-dimensional parameters ν\mathbb{\nu}. Low-cost approximations of the posterior distribution often rely on finding the maximum a posteriori (MAP) point, that is, the set of parameters that maximize the posterior probability. Upon taking logarithms, maximizing the posterior probability of the parameter set ν\mathbb{\nu} given the data Imeas\mathbf{I}_{\mathrm{meas}} is equivalent to minimizing the regularized cost functional [2 C. M. Bishop, Pattern recognition and machine learning. Springer, New York (2006) ]

J(ν)12I(ν)ImeasΓn12+12νν0Γpr12.J(\mathbb{\nu})≔\frac{1}{2}\lVert\mathbf{I}(\mathbb{\nu})-\mathbf{I}_{\mathrm{meas}}\rVert^{2}_{\mathbb{\Gamma}_{\mathrm{n}}^{\smash{-1}}}+\frac{1}{2}\lVert\mathbb{\nu}-\mathbb{\nu}_{0}\rVert^{2}_{\mathbb{\Gamma}_{\mathrm{pr}}^{\smash{-1}}}.

This is a nonlinear least-squares problem of the form previously considered in deterministic inversion, including regularization terms provided by the prior knowledge. We can solve it efficiently by using an adapted Levenberg–Marquardt–Fletcher iterative scheme [15 R. Fletcher, Modified Marquardt subroutine for non-linear least squares. Technical report 197213 (1971) ]. Starting from ν0=ν0\mathbb{\nu}^{0}=\mathbb{\nu}_{0}, we set νk+1=νk+ξk+1\mathbb{\nu}^{k+1}=\mathbb{\nu}^{k}+\mathbb{\xi}^{k+1}, where ξk+1\mathbb{\xi}^{k+1} is the solution of

(HλkGN(νk)+ωkdiag(HλkGN(νk)))ξk+1=gλk(νk).\bigl(\mathbf{H}^{\mathrm{GN}}_{\lambda_{k}}(\mathbb{\nu}^{k})+\omega_{k}\operatorname{diag}(\mathbf{H}^{\mathrm{GN}}_{\lambda_{k}}(\mathbb{\nu}^{k}))\bigr)\mathbb{\xi}^{k+1}=-\mathbf{g}_{\lambda_{k}}(\mathbb{\nu}^{k}).

Here, HGN\mathbf{H}^{\mathrm{GN}} is the Gauss–Newton approximation to the Hessian of the functional (10) and g\mathbf{g} is its gradient, while λk\lambda_{k} is a scaling factor for Γpr1\mathbb{\Gamma}_{\mathrm{pr}}^{-1} that balances the different orders of magnitude of the two terms defining the cost in the first iterations, and becomes equal to 1 at a certain point. At each step, the adjustable parameter ωk>0\omega_{k}>\nobreak 0 increases until the cost J(νk)J(\mathbb{\nu}^{k}) decreases, and decreases otherwise, making the iteration closer to Gauss–Newton or gradient schemes as required.

Linearization about the resulting MAP point νMAP\mathbb{\nu}_{\mathrm{MAP}} (the so-called Laplace approximation) provides an approximation of the posterior distribution by a Gaussian with mean νMAP\mathbb{\nu}_{\mathrm{MAP}} and posterior covariance Γpo=HGN(νMAP)1\mathbb{\Gamma}_{\mathrm{po}}=\mathbf{H}^{\mathrm{GN}}(\mathbb{\nu}_{\mathrm{MAP}})^{-1}. Sampling this Gaussian, we extract statistical information representing the dominant mode at a much lower computational cost, see Figure 6. Reaching νMAP\mathbb{\nu}_{\mathrm{MAP}} takes about 20 steps of scheme (11). The whole process, sampling included, is finished in a few minutes, instead of a few days.

Figure 6. Counterpart of Figure 5 using the Laplace approximation of the posterior density. The main mode corresponding to the true object is captured. Axis units are microns. Redrawn from [8]. 🅭🅯 CC BY 4.0

We have considered κi\kappa_{i} fixed and known in these tests. In case it is constant and unknown, it becomes an additional parameter included in ν\mathbb{\nu}. In the end, we obtain additional histograms reflecting uncertainty about the value with highest probability [8 A. Carpio, S. Iakunin and G. Stadler, Bayesian approach to inverse scattering with topological priors. Inverse Problems 36, Article ID 105001 (2020) ].

5 Perspectives

Digital holography poses challenging inverse problems which provide an opportunity to develop and test a variety of analytical and computational tools. First guesses of imaged objects are obtained by calculating the topological derivative of misfit functionals comparing the true hologram and the synthetic holograms that would be generated for different object configurations according to the selected forward model. Such guesses are robust to noise in the data. To reduce dimensionality, one can characterize the imaged objects by means of starshaped parametrizations. In a deterministic framework, we have shown that hybrid schemes combining iteratively regularized Gauss–Newton methods with topological derivative initializations and updates lead to good reconstructions of simple object configurations in a few steps, using stopping criteria that take into account the expected level of noise in the data. We are able to quantify uncertainty in such predictions by resorting to Bayesian formulations with topological priors. In two dimensions, Markov chain Monte Carlo methods provide a complete characterization of the posterior probability of the observed hologram being that is generated by a few starshaped objects. Three-dimensional tests are affordable for very simple shapes, such as a sphere or a cylinder [11 T. G. Dimiduk and V. N. Manoharan, Bayesian approach to analyzing holograms of colloidal particles. Optics Express 24, 24045–24060 (2016) ]. Handling high-dimensional parametrizations, in three dimensions or just irregular shapes, requires the introduction of strategies to reduce the computational cost. Laplace approximations based on optimizing to find the highest probability parameter set and then linearizing the posterior probability about it to obtain a multivariate Gaussian distribution are useful tools for uncertainty quantification when there is a single dominant mode. Developing fast sampling methods which are robust as dimension grows would be an important step forward to handle more general situations.

Holography set-ups are particularly challenging due to the fact that a single incident wave is used. We have focused here on light imaging, though acoustic waves can also be used to resolve at different scales. We expect similar techniques to be useful in inverse scattering problems involving other types of waves and different emitter/receiver configurations, such as microwave imaging or elastography, for instance.

Acknowledgements. Thanks to Krzysztof Burnecki for the kind invitation to write this article. The author is indebted to the Kavli Institute Seminars at Harvard for the interdisciplinary communication environment that initiated this work, M. P. Brenner for hospitality while visiting Harvard University, R. E. Caflisch for hospitality while visiting the Courant Institute at New York University, V. N. Manoharan for discussions of light holography, and D. G. Grier for an introduction to acoustic holography. The work summarized here has been done in collaboration with Thomas G. Dimiduk (Harvard University and Tesla), María Luisa Rapún (Universidad Politécnica de Madrid), Virginia Selgas (Universidad de Oviedo), Frédérique Le Louër (Université de Technologie de Compiègne), Georg Stadler (New York University), and Sergei Iakunin (Universidad Complutense de Madrid, now at the Basque Center for Applied Mathematics). This research was partially supported by the FEDER/MICINN–AEI grants PID2020-112796RB-C21, MTM2017-84446-C2-1-R, Fundación Caja Madrid Mobility grants, and MICINN “Salvador de Madariaga” Mobility grant PRX18/00112.

Ana Carpio graduated in numerical analysis from Universidad del País Vasco in Spain. She holds a PhD in mathematics from Laboratoire Jacques Louis Lions (Université Paris VI, now Paris Sorbonne), and has been a postdoctoral fellow at the Oxford Centre for Industrial and Applied Mathematics. She is a recipient of the SEMA (Spanish Society of Applied Mathematics) Prize to Young Researchers. Since 2006, she is a professor of applied mathematics at Universidad Complutense de Madrid and a member of the Gregorio Millán Barbany Institute for Modelling and Simulation in Fluid Dynamics, Nanoscience and Industrial Mathematics at Universidad Carlos III de Madrid. Currently, she serves as a Spanish representative in the ECMI (European Consortium for Mathematics in Industry) Council. Her main topics of research nowadays are inverse problems and data driven computational models in biomedicine. ana_carpio@mat.ucm.es

    References

    1. A. B. Bakushinskii, On a convergence problem of the iterative-regularized Gauss–Newton method. Zh. Vychisl. Mat. i Mat. Fiz. 32, 1503–1509 (1992)
    2. C. M. Bishop, Pattern recognition and machine learning. Springer, New York (2006)
    3. C. F. Borhen and D. R. Huffman, Absorption and scattering of light by small particles. Wiley Sciences, John Wiley & Sons, Berlin (1998)
    4. T. Bui-Thanh, O. Ghattas, J. Martin and G. Stadler, A computational framework for infinite-dimensional Bayesian inverse problems Part I: The linearized case, with application to global seismic inversion. SIAM J. Sci. Comput. 35, A2494–A2523 (2013)
    5. A. Carpio, T. G. Dimiduk, F. Le Louër and M. L. Rapún, When topological derivatives met regularized Gauss–Newton iterations in holographic 3D imaging. J. Comput. Phys. 388, 224–251 (2019)
    6. A. Carpio, T. G. Dimiduk, M. L. Rapún and V. Selgas, Noninvasive imaging of three-dimensional micro and nanostructures by topological methods. SIAM J. Imaging Sci. 9, 1324–1354 (2016)
    7. A. Carpio, T. G. Dimiduk, V. Selgas and P. Vidal, Optimization methods for in-line holography. SIAM J. Imaging Sci. 11, 923–956 (2018)
    8. A. Carpio, S. Iakunin and G. Stadler, Bayesian approach to inverse scattering with topological priors. Inverse Problems 36, Article ID 105001 (2020)
    9. A. Carpio and M.-L. Rapún, Solving inhomogeneous inverse problems by topological derivative methods. Inverse Problems 24, Article ID 045014 (2008)
    10. D. Colton and R. Kress, Inverse acoustic and electromagnetic scattering theory. Appl. Math. Sci. 93, Springer, Berlin (1992)
    11. T. G. Dimiduk and V. N. Manoharan, Bayesian approach to analyzing holograms of colloidal particles. Optics Express 24, 24045–24060 (2016)
    12. V. Domínguez, S. Lu and F.-J. Sayas, A fully discrete Calderón calculus for two dimensional time harmonic waves. Int. J. Numer. Anal. Model. 11, 332–345 (2014)
    13. M. M. Dunlop and G. Stadler, A gradient-free subspace-adjusting ensemble sampler for infinite-dimensional Bayesian inverse problems, preprint, arXiv:2202.11088v1 (2022)
    14. G. R. Feijoo, A new method in inverse scattering based on the topological derivative. Inverse Problems 20, 1819–1840 (2004)
    15. R. Fletcher, Modified Marquardt subroutine for non-linear least squares. Technical report 197213 (1971)
    16. J. Fung, R. P. Perry, T. G. Dimiduk and V. N. Manoharan, Imaging multiple colloidal particles by fitting electromagnetic scattering solutions to digital holograms. J. Quant. Spectroscopy Radiative Transfer 113, 212–219 (2012)
    17. A. Gelman and D. B. Rubin, Inference from iterative simulation using multiple sequences. Statist. Sci. 7, 457–472 (1992)
    18. J. Goodman and J. Weare, Ensemble samplers with affine invariance. Commun. Appl. Math. Comput. Sci. 5, 65–80 (2010)
    19. P. Grisvard, Elliptic problems in nonsmooth domains. Classics Appl. Math. 69, SIAM, Philadelphia (2011)
    20. H. Harbrecht and T. Hohage, Fast methods for three-dimensional inverse obstacle scattering problems. J. Integral Equations Appl. 19, 237–260 (2007)
    21. T. Hohage, Logarithmic convergence rates of the iteratively regularized Gauss–Newton method for an inverse potential and an inverse scattering problem. Inverse Problems 13, 1279–1299 (1997)
    22. J. Kaipio and E. Somersalo, Statistical and computational inverse problems. Appl. Math. Sci. 160, Springer, New York (2005)
    23. S. H. Lee, Y. Roichman, G. R. Yi, S. H. Kim, S. M. Yang, A. van Blaaderen, P. van Oostrum and D. G. Grier, Characterizing and tracking single colloidal particles with video holographic microscopy. Optics Express 15, 18275–18282 (2007)
    24. F. Le Louër, A spectrally accurate method for the direct and inverse scattering problems by multiple 3D dielectric obstacles. ANZIAM J. 59, E1–E49 (2018)
    25. A. Maier, S. Steidl, V. Christlein, J. Hornegger, Medical imaging systems: An introductory guide. Springer (2018)
    26. P. Marquet, B. Rappaz, P. J. Magistretti, E. Cuche, Y. Emery, T. Colomb and C. Depeursinge, Digital holographic microscopy: a noninvasive contrast imaging technique allowing quantitative visualization of living cells with subwavelength axial accuracy. Optics Letters 30, 468–478 (2005)
    27. M. Masmoudi, J. Pommier and B. Samet, The topological asymptotic expansion for the Maxwell equations and some applications. Inverse Problems 21, 547–564 (2005)
    28. S. W. McCandless and C. R. Jackson, Principles of synthetic aperture radar. In SAR Marine Users Manual, NOAA (2004)
    29. S. Meddahi, F.-J. Sayas and V. Selgás, Nonsymmetric coupling of BEM and mixed FEM on polyhedral interfaces. Math. Comp. 80, 43–68 (2011)
    30. R. M. Neal, MCMC using Hamiltonian dynamics. In Handbook of Markov Chain Monte Carlo, edited by S. Brooks, A. Gelman, G. L. Jones and X. L. Meng, Chapman & Hall/CRC Handb. Mod. Stat. Methods, CRC Press, Boca Raton, 113–162 (2011)
    31. J.-C. Nédélec, Acoustic and electromagnetic equations. Appl. Math. Sci. 144, Springer, New York (2001)
    32. J. Sokołowski and A. Żochowski, On the topological derivative in shape optimization. SIAM J. Control Optim. 37, 1251–1272 (1999)
    33. A. Tarantola, Inverse problem theory and methods for model parameter estimation. Society for Industrial and Applied Mathematics (SIAM), Philadelphia (2005)
    34. J. Tromp, Seismic wavefield imaging of Earth’s interior across scales. Nature Reviews Earth & Environment 1, 40–53 (2020)
    35. T. Vincent, Introduction to holography. CRC Press (2012)
    36. A. Wang, T. G. Dimiduk, J. Fung, S. Razavi, I. Kretzschmar, K. Chaudhary and V. N. Manoharan, Using the discrete dipole approximation and holographic microscopy to measure rotational dynamics of non-spherical colloidal particles. J. Quant. Spectroscopy Radiative Transfer 146, 499–509 (2014)
    37. A. Yevick, M. Hannel and D. G. Grier, Machine-learning approach to holographic particle characterization. Optics Express 22, 26884–26890 (2014)
    38. M. A. Yurkin and A. G. Hoekstra, The discrete-dipole-approximation code ADDA: Capabilities and known limitations. J. Quant. Spectroscopy Radiative Transfer 112, 2234–2247 (2011)

    Cite this article

    Ana Carpio, Seeing the invisible: Digital holography. Eur. Math. Soc. Mag. 125 (2022), pp. 4–12

    DOI 10.4171/MAG/99
    🅭🅯
    This open access article is published by EMS Press under a CC BY 4.0 license, with the exception of logos and branding of the European Mathematical Society and EMS Press, and where otherwise noted.