3D refractive index reconstruction from phaseless coherent optical microscopy data using multiple scattering-based inverse solvers—a study

Reconstructing 3D refractive index profile of scatterers using optical microscopy measurements presents several challenges over the conventional microwave and RF domain measurement scenario. These include phaseless and polarization-insensitive measurements, small numerical aperture, as well as a Green’s function where spatial frequencies are integrated in a weighted manner such that far-field angular spectrum cannot be probed and high spatial frequencies that permit better resolution are weighed down. As a result of these factors, the non-linearity and the ill-posedness of the inverse problem are quite severe. These limitations have imposed that inverse scattering problems in the microscopy domain largely consider scalar wave approximations and neglect multiple scattering. Here, we present first inverse scattering results for optical microscopy setup where full-wave vectorial formulation and multiple scattering is incorporated. We present (a) how three popular inverse scattering solvers from microwave domain can be adapted for the present inverse problem, (b) the opportunities and challenges presented by each of these solvers, (c) a comparative insight into these solvers and contrast with the simpler Born approximation approach, and (d) potential routes to improve the performance of these solvers for the hard inverse problem of optical microscopy.

Reconstructing 3D refractive index profile of scatterers using optical microscopy measurements presents several challenges over the conventional microwave and RF domain measurement scenario.These include phaseless and polarization-insensitive measurements, small numerical aperture, as well as a Green's function where spatial frequencies are integrated in a weighted manner such that far-field angular spectrum cannot be probed and high spatial frequencies that permit better resolution are weighed down.As a result of these factors, the non-linearity and the ill-posedness of the inverse problem are quite severe.These limitations have imposed that inverse scattering problems in the microscopy domain largely consider scalar wave approximations and neglect multiple scattering.Here, we present first inverse scattering results for optical microscopy setup where full-wave vectorial formulation and multiple scattering is incorporated.We present (a) how three popular inverse scattering solvers from microwave domain can be adapted for the present inverse problem, (b) the opportunities and challenges presented by each of these solvers, (c) a comparative insight into these solvers and contrast with the simpler Born

Introduction
Reconstruction of 3D refractive index (RI) or relative permittivity profile of scatterers given far-field measurements is a standard inverse scattering problem [1].It has been of interest for radio frequency wave, microwave, and millimeter wave range applications involving urban [2,3], geological [4,5] and biological [6][7][8] applications.However, interest in this inverse problem in the optical domain is relatively less explored [9][10][11].The primary reasons are (a) limited numerical aperture (NA) available for illumination of the sample and the collection of scattered light, and (b) optical range detectors being intensity-only measurement devices.Both these problems make optical domain 3D RI reconstruction highly ill-posed and non-linear.
The two popular optical systems where 3D RI reconstruction has been performed are optical holography and optical tomography systems [12,13].The inversion process in tomography generally involves two steps: 2D phase retrieval for each angle and 3D RI reconstruction by tomography, such as projection algorithm [14,15] which models the light as a ray and diffraction tomography [16,17] where weak scattering is assumed.Such methods require interferometric systems [18][19][20] to obtain phase information.
Tomographic reconstruction methods based on intensity-only data have also been proposed [21,22] which allow direct reconstruction of RI without phase retrieval.Beyond Born [23] and Rytov [24] approximation, forward models considering multiple scattering along the propagation direction for thick samples have been proposed and used for inversion.In such methods, the sample is decomposed of thin slices along optical axis and the field is calculated layer by layer with some linear approximation, such as multi-slice beam propagation methods [25][26][27] and multi-layer Born methods [28].These forward models have been used with tomography for 3D sample reconstruction.
An important aspect of all of these methods is that they use scalar approximations of far field in order to simplify the inverse problem and make it tractable.Full-wave vectorial problem and one that incorporates high orders of multiple scattering have not been addressed in optical tomography and holography, nor in high NA coherent microscopy systems with multiple angles of illumination.
Coherent wide-field optical microscopes image the intensity of electric far-field pattern on the camera which maps the light scattered by sample after it passes through a two-lens system.For such systems, the positions of the detectors or camera pixels have a one-to-one correspondence with a conjugate plane in the sample region, i.e. such a system is an imaging system and the measurements represent a qualitative image of the conjugate sample plane.In addition to the conventional challenges for optical regime discussed before, a big impediment is that the dyadic Green's function of such optical system is significantly more complex than usually used in non-optical inverse scattering problems and scalar tomography problems.Now, appropriate formulation of Green's functions have been derived [29,30] incorporating the imaging system as well as reflections from substrate.An accurate forward model based on Lippmann-Schwinger equation has also been developed and experimentally verified [31].Therefore, it is opportune to consider contemporary full-wave multiple scattering based inverse solvers for 3D RI reconstruction.

The problem statement
The microscope setup under investigation is sketched in figure 1(a).The biological samples are placed on a substrate, characterized by RI n sub and immersed in a homogeneous background with n obj .The values of n sub and n obj are known a priori.The sample has an inhomogeneous RI distribution n(r).
The sample is illuminated by transverse electric (TE) polarized plane waves from different angles which propagate parallel to the immersion-substrate interface successively, as shown in figure 1(b).The entire sample is illuminated, and the light scattered by the sample is collected by the objective lens with a specified NA and further focused by the tube lens onto the camera placed at the focal plane of the tube lens.Since the illumination direction is perpendicular to the optics axis and the NA of the objective lens does not overlap the illumination, the illumination components do not reach the camera.In the microscopy parlance, this implies a dark-field configuration and in the inverse scattering parlance, the measurements are the intensity of the scattered far field mapped to the camera pixels through the optical systems.Notably, neither the sample nor the camera is moved, which means that only one lateral image is taken per illumination.The problem is to reconstruct n(r) given these measurements and corresponding illumination conditions.

The differentiating features of this phaseless inverse scattering problem
Unlike the measurement system considered in traditional microwave imaging, in a microscope setup, the scattered wave propagates through the objective lens and tube lens before collected by the camera.Thus, this problem presents some important and interesting features that differentiates it from the conventional phaseless inverse scattering setup.
An obvious difference is that not only the phase information but also the polarization information is lost in such a set-up.This means that the entire electric far field of the light scattered by the sample, containing complex valued azimuthal and elevation components, is reduced to a single intensity measurement, i.e. the loss of information is four-fold.
Due to the existence of the objective lens with a limited size, only part of the propagating plane-wave components of the scattered wave can reach camera, i.e. there is an upper limit to the transverse spatial frequency that gets measured by the system that is determined by the NA of the microscope.Therefore, the far-field Green's function works as low-pass filter.Moreover, it is easy to see in the derivations of Green's functions for such two-lens systems [29,31,32] that the lowest transverse spatial frequencies are assigned the highest weights in the integrand, implying suppression of high spatial frequency content vital for resolution of 3D RI distribution.
In a conventional lens-less setup, the electric field due to a point source (such as used for Green's function formalism) at an arbitrary measurement point in the far-field comprises of a single wave-vector.This is not the case in a two-lens system such as considered here.The first lens, i.e. the objective lens, preserves the wave vectors of the wave coming from a point source in the sample region, effectively functioning as band-limited Fourier transform of the field incident on it.However, the second lens, i.e. the tube lens, integrates over all the wavevectors originating from a point source in the sample region and passing the objective lens, effectively functioning as a band-limited inverse Fourier transform.This implies that the spectral separation of the information in the scattered light is lost.
Lastly, even though it appears in the first glance that the presence of several pixels in a camera presents an opportunity, there is a subtle but important detail.Since this system is an imaging system, the light scattered from one point is mapped to the camera as a diffractionlimited blob whose main lobe spans only a small region in the camera plane.The side lobes are usually small enough to be negligible in practical scenario.With sufficiently large magnification, the size of this blob implies that only a few pixels span the main lobe, restricting the number of measurements per point scatterer to be in the range 5-50.

Inverse scattering solver candidates
Inspired by the success of full-wave inversion method in microwave imaging, here we formulate the problem 3D RI reconstruction using coherent optical microscope as a full-wave inverse scattering problem and test the performance of the existing inverse algorithms that have been widely used.
Various full-wave nonlinear inverse scattering algorithms that incorporate multiple scattering have been proposed.These include distorted Born iterative method [33], contrast source inversion (CSI) method [34,35], subspace-based optimization method (SOM) [36] and its variants [37,38], contraction integral equation (CIE) based inversion [39], which have been widely applied to both 2D [40] and 3D [41,42] imaging problems.Their adaptions to the phaseless measurement have also been presented [32,43,44].However, all of these adaptations assume that polarization-separated measurements are available.On the other hand, our problem pertains absence of both phase and polarization-separated measurements.
We choose three algorithms that have been widely used or proposed for highly-nonlinear problem, namely, CSI, CIE based inversion and SOM to solve this specific problem.These have also been adapted for phaseless inversion, therefore they provide a good case study for investigation of their suitability for phaseless, polarization-suppressed measurements using a microscopy system.

Outline
The rest of the paper is organized as follows.A brief description of the imaging system and forward formulation is given in section 2. The chosen inverse algorithms and their adaptations to the concerned problem are shown in section 3 with their performance on the numerical examples is given in section 4. We discussed on the present algorithms in section 5 and conclusion is drawn in section 6.

Forward model
The microscope setup has already been presented in figure 1(a) and discussed briefly in section 1.1.We continue our discussion of forward model from thereon.We define D as the domain of interest (DoI) with all samples enclosed which has an inhomogeneous RI distribution n(r).They are illuminated by TE-polarized plane waves.Time dependence exp(−iωt) is discarded for simplicity in the following formulation.The scattered wave due to the presence of the samples will propagate through the objective lens and tube lens, and is collected by the camera (S).With N inc illumination, N inc camera images are recorded.
Two local Cartesian coordinate systems are built with the origins being the focal points of the objective lens and tube lens.z-axis is chosen as the optical axis.For pth illumination, it generates vectorial incident field E inc p (r) which is decomposed into x, y and z components (z component E inc p,z (r) = 0 with TE wave for all p).The relation between the total field and the sample contrast can be described by Lippmann-Schwinger integral equation where contrast is defined as χ(r) = (n(r)/n obj ) 2 − 1 for nonmagnetic objects.Here we assume that the substrate is thick enough and G d (r, r ′ ) is dyadic half-space Green's function [45].
The scattered wave reaching camera pixel r cam can be calculated as where G s (r cam , r ′ ) is the far-field Green's function giving the field on the camera generated by a dipole source under certain microscope configuration.The expression of far-field Green's function for our configuration has been derived in [29] and used in [31] for demonstrating the forward solver.Since the incident wave propagates perpendicular to the optical axis, only scattered field is captured by the camera and the camera measures the total intensity of the field as To solve the problem numerically, the method of moments [46] with pulse basis function and point matching is used to discretize the equations.The DoI is discretized into N = N x × N y × N z cells with the center of each cell located at r n , n = 1, . . ., N. The image plane consists of M = M x × M y square pixels, centered at r cam m , m = 1, . . ., M. The above equations can be written in matrix vector form where

Inversion with intensity-only images
With the forward model, we can retrieve the 3D RI profile of the sample with inverse algorithms given multiple 2D camera images, each corresponding to a different illumination direction.The image recorded by the microscope camera contains no phase and polarization information of the scattered field, which increases the ill-posedness and nonlinearity of the inverse problem.Due to the fact that the z-axis is the optical axis, the z-component of the scattered field arriving at the camera is much smaller than transverse counterparts.Besides, since TE illumination is used, the z-components of all the illuminations are zero.As a consequence, the z-component in the contrast source distribution is induced only by the multiple-scattering effect, so that it is also smaller than J x and J y .Therefore, to ease the computational burden, we discard E sca z and J z in the following inversion process, except for one situation which we discuss later.Thus, in the description below G In the inverse scattering solver design, the design of cost function, optimization variables, and optimization strategy are the core differentiating features among the different methods.The cost function typically contains at least a term called the data equation error, which corresponds to the misfit between the actual measurements and the equivalent quantities computed using an estimate of the contrast.In addition, methods that incorporate multiple scattering typically include an additional state equation error term in their cost function.The state equation error term usually computes a misfit between two estimates of a physical quantity in the DoI in which one estimate directly uses the current estimate of the contrast.The state equation error is largely responsible for incorporating the multiple scattering information in near field, which may not propagate to the far field where measurements are conducted.The total electric field in the DoI or the contrast source distribution in the DoI are examples of physical quantities which may be used for defining the state equation error.The specific expressions of the data and state equation errors and their relative contribution to the cost function may vary across methods and their adaptations.
Next, we briefly present our adaptions of the inverse scattering algorithms which we use to solve this specific problem and conduct our comparative study.These include Born approximation (BA) as it is the most accepted and used approach for 3D RI reconstruction in the optical domain, and three popular multiple scattering based inverse scattering solvers namely CSI, CIE, and SOM.For all of them, we use the simplification discussed in the previous paragraph relating to neglecting the z-components of the electric field and the contrast source.However, for CIE which presented relative computational simplicity, we considered also the case where the z-components of the electric field and the contrast source are included in the inversion.

BA
In BA, the size and the contrast of the object is assumed to be small enough such that multiple scattering can be ignored and the total field E tot (r) can be replaced by the incident field E inc (r) in the DoI.BA uses only the data equation error as the cost function.We define the data equation error for BA in a very conventional manner as follows: with η S p being the inverse of the normalization factor and expressed as η S p = 1/∥I p ∥ 2 .Even though BA is used, the problem is still nonlinear due to the lack of phase and polarization information.We use gradient-based optimization to find a solution.The gradient of the cost function with respect to the contrast is with scattered field Here superscript ' †' means conjugate transpose and diag{a} generates a diagonal matrix with the diagonal vector being a.

CSI
In CSI, contrast source J(r) = χ(r)E tot (r) is regarded as an intermediate optimization variable for the data equation error and as the physical quantity on which the state equation error is defined.The data and state equation errors being used in our adaptation of CSI are as follows: (8a) with η D p being the inverse of the normalization factor for state equation and we choose p is the same as used in BA.Note that in comparison to the data equation error for BA, here the contrast source is used which incorporates multiple scattering through E tot (r).Further, the term X 2 G d J p incorporates multiple scattering in the state equation error.The cost function is (F S + F D ).
Alternate update of J and χ is performed to minimize the cost function.The gradient of the cost function with respect to the contrast source is with After the contrast source J is updated, the contrast at nth voxel element in DoI can be obtained explicitly written as (10) with superscript * represents conjugate and

CIE
CIE based inversion was proposed to solve highly nonlinear inverse scattering problem.Modified contrast R(r) is defined as R(r) = βχ(r)/(βχ(r) + 1) with regularization parameter β to control the ratio between local and global contribution to multiple scattering.Regularization technique proposed in [38] is employed to constrain the contrast source in a stable subspace, where a parameter M F is used to determine the size of Fourier basis utilized to represent contrast source J p .Accordingly, J p is represented using α p which are the bandlimited Fourier coefficients of J p where the band-limit is determined by M F .Therefore, the size of α p is determined by M F .With increasing M F , the bandwidth for representing J p increases, higher spatial frequency components are included and the number of entries in Fourier coefficients α p increases.As the optimization proceeds, β becomes smaller and M F becomes larger.
Combining work related to CIE for 3D problem assuming both phase and intensity measurements [47] and for phaseless 2D problem [44], we define the data and state equation errors for given values of β and M F as (11a) is the same as used in BA.F −1 {•} represents the inverse Fourier transform.Note that here we operate such transform to x, y and z components respectively and then concatenate the results as a vector.The cost function is (F S + 100F D ), where the multiplicative coefficient is empirically chosen.
Similar with CSI, alternate minimization is used to find a solution.However, in CIE the Fourier coefficient vector α p and the modified contrast R are used as the optimization variables instead of J p and χ, respectively.The gradient with respect to α p is expressed as The modified contrast R at the nth voxel in the DoI is updated as with

SOM
Before discussing SOM, we bring to the notice of the readers that while we adapted each method above for our problem through custom definitions of cost functions, our adaptation of SOM is significantly different from the original SOM, not just in the cost function definition, but also in other aspects that will be appropriately emphasized in the text below.SOM utilizes the spectral features of the mapping from the contrast source to the scattered field given by singular value decomposition (SVD).With , where U and V represent the basis set of left and right singular vectors, the contrast source can be represented by right singular vectors as J = Vα.Please note that α here is different from the vector α in CIE.In CIE α represents the coefficients of Fourier basis vectors and in SOM α represents the coefficients of right singular basis set V. SOM splits the contrast source J into a deterministic components J + and an ambiguous component J − , i.e.J = J + + J − .This is accomplished by splitting V into two sets V + and V − using a user-specified parameter L which defines the number of basis vectors used to represent the deterministic part.Accordingly, and the deterministic and ambiguous parts of the contrast source are given as , respectively.Further, as obvious, α + represents the coefficients of deterministic part and α − represent the coefficients of ambiguous part of the contrast source, such that α is a concatenation of α + and α − .
For phaseless data, similar with [43], the deterministic part of the contrast source J + in signal space is obtained with optimization method by considering data misfit, unlike original method where it is determined directly from known complex scattered field.We follow the same strategy as [43] for reconstructing the deterministic part.The reconstruction of the ambiguous part of the contrast source is extremely computationintensive for 3D RI reconstruction problem due to the large dimensionality of V − .Therefore, instead of following the conventional practice in SOM of using the entire V − for reconstruction from the beginning, we gradually increase the number of vectors in V − used in reconstruction over the progression of optimization.
The whole optimization therefore consists of two steps.In the first step, the coefficients α + are solved by optimizing The gradient with respect to the coefficients α + is with In the second step, we use the first L δ vectors of V − , referred to as , to approximate the ambiguous part of the current.For this, we define ], and the corresponding coefficient vector α δ which also includes coefficients α + .We optimize α δ directly by considering both data and state equation errors expressed as (16a) η S p is the same as used in BA.The cost function in this step (F S + 0.01F D .The factor 0.01 is chosen empirically, but the choice is motivated by the fact that the data error F S has already been reduced due to step 1, and a small coefficient for F D helps in balancing the contribution of data and state error equations in the optimization.The gradient with respect to the coefficients α δ is The expression for the contrast is the same as (10) with J p = V δ α δ p .We use the same optimization strategy as CIE where multiple rounds are run and increasingly more L δ singular vectors are used to form the basis in each round.

Summary of the methods and their adaptations
BA uses only data equation and reconstructs the contrast directly.Our problem being phaseless, our BA formulation is still non-linear and we use gradient based optimization for the contrast.
CSI uses both data and state error equations, and it employs the contrast source as an intermediate optimization variable.A two step update scheme is used in which first the contrast source is optimized, and the contrast is updated using it.Our adaptation of CSI features one main difference from the CSI formulation for phaseless problem, which is the absence of polarization-separated measurements.This difference also applies to CIE and SOM.
CIE uses a modified contrast function whose parameter β can be varied through the progression of optimization.Further, CIE uses Fourier coefficients of the contrast source as an intermediate optimization variable, and allows the number of Fourier coefficients of the contrast source to change with the progression of optimization.Beside these important differences, our adaptation of CIE generally follows a CSI-like scheme for the cost function and the optimization process.
SOM uses SVD of G s to obtain its right singular vectors.These right singular vectors are used as a basis set for the contrast source.The coefficients of these right singular vectors are considered as intermediate optimization variable.Furthermore, instead of optimizing all the coefficients at the same time, SOM optimizes a subset of coefficients representing the deterministic part of the contrast source separately using the data equation error alone, and the complementary subset of coefficients representing the ambiguous part of the contrast source is optimized separately using both data and state equation errors.
Beside these characteristics, some adaptations of SOM including ours may follow a scheme similar to CIE for the cost function and optimization process.In our adaptation of SOM, there are two key differences.The first difference is that in the first step we optimize the coefficients of the deterministic part alone and using the data equation error only, and in the second step we optimize the coefficients of both the deterministic and ambiguous parts of the contrast source and using both the data and state error equations.The contrast is optimized using the contrast source in each iteration of the second step.The second difference is that we optimize only a subset of the ambiguous part of contrast source in the second step.The size of this subset is allowed to vary as the optimization progress.
Further, for all of the above, we have neglected the z-components of electric fields and contrast source, which presents computational simplification without any loss of generality.Lastly, we have included experiments on CIE where we consider also the z-components.We refer to these results using CIE(z).

Results
We test the above inverse algorithms with several numerical examples under the same microscope setup and simulation configuration.Here we will mainly focus on the reconstruction of the contrast χ.Once the contrast is obtained, the 3D RI profile of the object can be easily calculated with known RI of the immersion medium.
Microscope setup.The glass substrate with n sub = 1.515 is located at z sub = 0.The samples are immersed in air n obj = 1.0 and illuminated by plane waves with wavelength λ = 638 nm from N inc = 8 angles which are evenly distributed in [0, 2π).The image of the sample is obtained by a 40 × 0.65 NA microscope.The size of a single camera pixel is 6.5 µm and each image consists of 18 × 18 pixels.Simulation Configuration.Different grid sizes are used for forward and inverse problems to avoid inverse crime.The DoI has a size of 1.5 µm × 1.5 µm × 0.6 µm and is divided into cubic cells with side length of 20 nm to compute synthetic images on the camera.The camera images are obtained by solving the forward problem with Bi-CGSTAB method [48].Here we test the algorithms with clean data such that no noise is added to the generated data.For inversion, the cell size is set as 30 nm, leading to 50 × 50 × 20 cells.With such cell side length (∼λ/21), the reconstructed object resolutions are in the deep subwavelength regime.Initial guess.Since the problem is highly nonlinear, it is important to find a good initial guess.Fortunately microscope can provide images which give blurred structures of the sample.A simple practice is to combine the normalized camera images as I C = N inc p=1 I p / max(I p ) with I p being the 2D camera image of the pth illumination and max(I p ) being the maximum intensity in the image I p .Projection of this image into sample region cylindrically based on the magnification of the microscope can be used as the initial guess of the contrast.The contrast at each cell is proportional to the projected value and a maximum contrast is set manually according to Table 1.Implementation details of the inverse algorithms with total intensity data and difference with their origins with full-wave data.

Methods Our implementation
Origin with full-wave data BA data equation ( 6) the prior information about the sample.In this way we can get the initial contrast distribution χ (0) in DoI.Our choice of initial guess for the contrast source or its coefficients is zero in CSI and CIE.For the coefficients used in SOM, the initial guess of the coefficients for the deterministic part of the contrast source cannot be zero, and we set them as α For the new coefficients added in the following rounds, their values start from 0.
Regularization and optimization.We impose weighted L 2 total variation (TV) regularization F TV defined in [49] on all the algorithms.It is incorporated as an additive term for BA and as a multiplicative term for the other algorithms.Thus, we have the cost function F BA = F S + γF TV and F CSI/CIE/SOM = (F S + F D )F TV .Update of the contrast source is not affected except for BA where the gradient should consider the contribution from the regularization term.Conjugate gradient method [50] with Polak-Ribière search direction v is used to update the contrast in BA, contrast source in CSI, and the coefficients of the contrast source in CIE and SOM.Step size d is obtained by solving d = min F(a + dv) where F is the cost function and a is the current variable to be updated.With multiplicative regularization, we employ the update method described in chapter 4.6 of [51].The contrast is calculated analytically first as given above and then updated with two Jacobi iterations to impose the regularization effect.The implementation details are summarized in table 1 and the difference between the adaptions of these algorithms from their origins are also displayed.
Hyperparameters, termination criterion, and other notes.For BA, the regularization coefficient γ = 0.001 and δ 2 in the expression of TV is chosen as δ 2 = F S .This parameter is set δ 2 = F D in CSI, CIE and SOM.In CIE, three rounds are run with β = [6, 3, 1] and M F = [5,8,13].In SOM, for the first step, L = 36 right singular vectors are used to form the basis for optimization of the deterministic part of the contrast source.In the second step, number of vectors in V − involved in optimization starts from 100 with an increment of 100 in each round.The number of rounds for the second step is not set explicitly.Instead it is determined by either the termination condition being achieved or L δ reaching a value equal or more than 2000.
The optimization for all algorithms stops at nth iteration if the relative change in the cost function F (n) − F (n−5)∼(n−1) /F (n) is smaller than a threshold.Here, F (n−5)∼(n−1) is the mean of the cost function errors for (n − 5)th to (n − 1)th iterations.The threshold is chosen as 0.0005 for BA and CSI.In CIE, such thresholds for its three rounds are chosen as 0.002, 0.001 and 0.0005.For SOM, the threshold is 0.01 in the first step to get an estimate of J + and 0.001 in the second step for each round.While performing reconstruction using SOM, we found that there is artefact with much higher contrast on the four corners of the top surface compared with low estimated contrast of the object.To suppress such artefacts, we constrained that the contrast on all the boundary surfaces of the DoI to be zero except the bottom surface.The bottom surface corresponds to the interface between the substrate and the sample, and since the sample is not free floating and rests on the bottom surface, not all the voxels on this surface can have zero contrast.We incorporated this condition for all the solvers, including SOM.We do note that such constraint led to a degradation in the reconstruction of the sample for all the solvers and especially for SOM.In the future, it may be desirable to design more suitable strategy to deal with this artefact.
Assessment criteria.To assess the result quantitatively, the relative reconstruction error is defined as where χ gt is the groundtruth contrast and χ is the reconstruction result.We also test the errors for contrast source if it is reconstructed in the algorithms as Such error is defined for J x , J y and J z respectively.Lastly, where relevant, we have included the structural similarity index (SSIM) [52] and multiscale structural similarity index (MS-SSIM) metrics [53] to measure the structure similarity between the 3D reconstruction result of the contrast and the ground truth.Here both the result and ground truth are normalized by their maximum value such that both the value ranges are [0, 1].The dynamic range in the parameter of SSIM is set 1 and default values are used for other parameters.

Test example 1
In the first example, two beads with radius of a = 0.2 µm are placed along the y-axis and touching each other, as shown in figure 2(a).The RI of the beads is n = 1.05, corresponding to a contrast of χ = 0.1.The simulated camera images with illumination k 1 -k 8 are displayed in figure 2(c).We can distinguish these two beads from almost all the camera images except the ones obtained with illumination k 1 and k 5 where the scattering of the two beads are in phase and that constructive interference occurs.The normalized intensity image I C (figure 2(b)) also shows clear separation of these two beads when all the normalized camera images are summed up, and it is used to generate the initial guess of the contrast.
As mentioned before, the maximum initial contrast is set manually based on prior information about the samples.Here we test two values, 0.2 and 0.05, to see the influence of the initial guess on the reconstruction results in view that prior knowledge may be not accurate.The reconstruction results of all the algorithms are shown in figures 3 and 4 respectively.The SSIM, MS-SSIM and the reconstruction errors are summarized in table 2. Effect of the value of initial contrast on the solvers.In the qualitative sense as observed in figures 3 and 4, the initial value of the contrast does not affect the results of any solver significantly.Visually, one can distinguish the two beads in the results of each solver irrespective of the initial contrast value.The invariance along z-axis in the initial guess has also been well amended by each solver.Each presents a decent estimate of the shape of the two beads, although some subtle differences are interesting to note.CSI reconstructs shapes and sizes of the beads closer to the ground truth when higher initial value of contrast is employed.CIE(z) and CIE present better visual separability and better indication of shape of the beads when lower initial value is used, but better estimate of size of the beads is obtained when higher initial value is used.BA and SOM present reconstructions which seem robust and unaffected by the initial value.Quantitatively as noted in table 2, BA shows quite similar results for both initial contrast values, in the sense of both relative error and structural similarity.SOM also presents quantitative metrics not significantly affected by the initial contrast values.This is consistent with the qualitative observations of their robustness against the initial contrast value.In contrast, there is an obvious improvement in the performance of CSI and CIE when the optimization starts with a smaller contrast.This effect is more prominent for CIE, irrespective of the inclusion or exclusion of the z-components of the scattered field and contrast source.
Comparison across algorithms.Most algorithms give a reconstruction result with a higher contrast than the actual contrast of the beads.An exception is SOM which gives a significantly low contrast, about one tenth of the ground truth.In terms of size, all algorithms present results that indicate a size smaller than the actual size of the beads.
Qualitatively as noted from figures 3 and 4, the shapes and sizes of the beads are best appreciated in the result of CSI, closely followed by the results of CIE and CIE(z).Interestingly, quantitatively as noted in table 2, when using low initial contrast value, CIE and CIE(z) present the best MS-SSIM values and the lowest reconstruction errors.They also present shapes and sizes relatively comparable to CSI.Incorporating z-components of the scattered field and contrast source in CIE does not show evident improvement in this test example for both cases, which may be due to the increased computational complexity.Our adaptation of SOM gives the worst result.The relative error is high, which can be inferred from the reconstructed contrast of the objects being significantly smaller than the actual value (see in both figures 3 and 4).When we remove the influence of the contrast value by normalization in SSIM, the result is still unsatisfactory, which may be also observed from the figure that the beads reconstructed by SOM are most squeezed along z-axis.
Lastly, we note that BA does not perform too poorly in comparison to the other solvers.Its worst performance is quantitatively better than the best results of SOM and CSI, and only marginally poorer than the best results of CIE.
Convergence.The evolution of the reconstruction error e(χ) is shown in figure 5. Besides, the sum of data equation error and state equation error, denoted as F ′ , for each algorithm is also displayed.Readers are referred to table 1 for the expressions of the cost function of these algorithms in terms of the data and state equation errors.BA has the best convergence speed among all the solvers being compared although some oscillations are observed in the cost function during the optimization process.This advantage of BA is attributed to less nonlinearity of the inverse problem due to neglecting the multiple scattering effect and the least number of unknowns.CSI converges very slowly.It is also interesting that although CSI gives a smaller total cost than BA, the reconstruction error of CSI is higher, which may imply that CSI has converged to a wrong local minimum.CIE converges faster than CSI, which is potentially because the regularization parameters β and M F in CIE control the local effect of the contrast source in a subspace and therefore converge to more accurate reconstruction by changing the size of the region of local effect.Also, it is notable that when lower initial contrast value is used in CIE(z), the convergence clearly shows the transition between the first and second rounds at iteration number 242.This indicates the change of the size of local effect affecting the convergence directly.For the presented SOM, the reconstruction error almost does not change after the first update even though the cost function decreases slowly.We attribute this to our adaptation being sub-optimal.This is discussed in more detail in section 5.

Test example 2
This example consists of two ellipsoids with different sizes and contrasts.The center of the outer one is located at (0, 0, 0.2) µm with semiaxis lengths a x = 0.6 µm, a y = 0.6 µm and a z = 0.2 µm.The inner ellipsoid is located at (0.1, 0.05, 0.2) µm and the semiaxis lengths are 0.2 µm, 0.2 µm and 0.06 µm.The RI of the outer ellipsoid is n = 1.05 and the inner one has From figure 6, we can see that image I C simply obtained by summing up all the normalized camera images cannot distinguish the inner ellipsoid.Here we test the reconstruction results with initial contrast calculated from I C to see if the inner ellipsoid can be well distinguished from its surroundings and the contrast of both the objects can be estimated with sufficient accuracy.The maximum initial contrast value is chosen as 0.1.Figure 7 shows the reconstruction results and the assessment metrics are listed in table 3.
We see that all the algorithms can distinguish the inner ellipsoid from the outer one laterally even though it is not clear in initial guess.Most algorithms can also identify well that the inner ellipsoid is enclosed by the outer one axially.An exception to this observation is CSI, which gives the result such that both ellipsoids are elongated along z-axis.Compared with other algorithms, SOM gives a quite low contrast (as also noted in test example 1) and an irregular shape in the xy-plane.We recall that we introduced a constraint that the contrast on the boundary surfaces is zero so that SOM did not present high contrast artifacts on the boundary.We noted in results not displayed here that if this constraint is removed, SOM reconstructs the shape of this example much better than that seen in figure 7, however the boundary artifacts as significantly so pronounced that the reconstructed sample is almost invisible.Quantitatively, we see in table 3 that CIE and CIE(z) present the best similarity with the structure as well as the lowest reconstruction error, performing marginally better than BA.A practical challenge-inaccurate knowledge about sample-substrate interface.Experimentally, despite using the best optomechanical components and good calibration protocols, a situation is often encountered where we cannot precisely estimate the relative position of the substrate surface to the focal plane of the objective lens.In simple words, we cannot ensure that z sub is known accurately.This affects the expression of the far-field Green's function G s , and therefore the entire reconstruction.Using the test example 2, we test the performance of these inverse algorithms when wrong estimate of the substrate surface location is used, and as a consequence the wrong far-field Green's function G s is used.Apart from the correction location z sub = 0, another two hypothetical incorrect estimates of locations are also tested, namely z sub = 0.5 µm and z sub = −0.5 µm.New far-field Green's functions are calculated accordingly.Figure 8 shows the reconstructed results by the algorithms with the same initial guess but using the two incorrect locations.The defined errors are summarized in table 3.
When the location of the substrate is not correctly estimated with an error of 500 nm and wrong far-field Green's function is used for inversion, we can see that the object can be still reconstructed in general with the main structure identified.However, deformation appears especially significant for BA and CSI with case estimated z sub = 500 nm.CIE is the one that has the least change when the far-field Green's function changes and the reconstruction results look similar.We also plot the reconstruction error e(χ) and the cost F ′ as a function of iterations in figure 9 to study convergence of the different methods for the challenge scenarios.We see that the error decreases smoothly for BA, CSI and CIE with correct G s , however, when incorrect G s is used, the monotonic decrease may be interrupted and the error may start to rise even though the total cost is decreasing.We conjecture that this is because the data equation has low fidelity in the case of incorrect G s and so minimizing the data equation error may lead to a wrong solution.CIE performs well perhaps due to our choice of the cost function where the state equation error is amplified by 100 times.Consequently, the optimization is mainly guided by the state error equation instead of the data equation.
From this example, we can see that the inverse algorithms are robust to the chosen substrate location to some extent.With an error of about 0.8λ, there is a degradation in the performance but the main structures of the object are reconstructed.
An investigation into the consequence of polarization unavailability.Compared with phaseless inverse scattering problems that have been widely investigated, the present problem is more complicated due to the absence of the polarization information in the measurement data.Here we use test example 2 to check the consequence of the polarization unavailability.When the polarization information is available, the data equation error can be rewritten as Here, z component is ignored (as also elsewhere) and the state equation remains unchanged.Figure 10 shows the reconstruction result with and without polarization sensitive measurements while using the same initial contrast.The reconstruction error and data equation error are summarized in table 4.
For the ease of comparison, the results reconstructed from data with no polarization (figure 7) are reshown here in figure 10(b).There is no obvious qualitative difference between the results when considering whether the intensity data contains polarization information, except for SOM where the shape of the reconstructed object is evidently improved.Most algorithms show a smaller reconstruction error if the provided data contains polarization information, except CIE.We also noted that the convergence rate is not obviously improved, except for SOM where the number of iterations is decreased from 576 to 109 when polarized data is employed for reconstruction.For the cases where the scatterer is isotropic and not strong scattering, the relation between the polarization components of the contrast source almost follows the one between the components of the incident field which is already known.Therefore, in this case, the lack of polarization information will not affect the inversion results much.

Test example 3
Here we show an example where all the inverse algorithms fail.This example contains two spheres with radius a = 0.1 µm, located on y-axis with a center-to-center distance of 0.3 µm, as shown in figure 11, also with the RI n = 1.05.The size of the object is smaller than diffraction spot of the imaging system (spot size of a dipole source).From most camera images, we cannot   identify that there are two objects except the ones obtained with k 3 k 7 where destructive interference occurs.
The maximum initial contrast is set 0.05 and the reconstruction results are shown in figure 12.We see that all the inverse algorithms fail to identify these two beads and underestimate the contrast a lot.Almost all the algorithms show that the objects are on the bottom of the DoI, except CIE where object with a large height is reconstructed.The performance of CIE relies on choice of proper regularization parameters and there may exist a better set of parameters that can give satisfactory result.Reconstruction results of plane z = 105 nm by various algorithms with the center-to-center distance between the two beads equal to 300 nm (top), 350 nm (middle) and 400 nm (bottom).When the distance is 300 nm, no approach can resolve the two beads.When the distance increases to 350 nm, all the algorithms can resolve them but some only barely resolve.In contrast, all the algorithms are able to clearly identify the two beads when the distance is 400 nm.
We see that the inverse algorithms are not able to identify closely-located small objects, which may be a problem if we want to investigate small structures.
An investigation into resolution.In the example 3, we have shown that all the algorithms are unable to resolve two small beads with a center-to-center distance of 300 nm.Here, based on this example, we made a further investigation into the resolution performance of the algorithms.Figure 13 gives the reconstruction results of the two beads placed on y-axis with an increasing distance from 300 nm to 400 nm with an increment of 50 nm.All the algorithms are unable to resolve the two beads when the distance is 350 nm, while with more 50 nm, the two beads are resolvable with all the algorithms, however, CSI and CIE can only barely resolve the objects.When the center-to-center distance becomes 400 nm, the two beads are identified by all the algorithms clearly.In all, with the given example, BA and SOM show superior resolution to CSI and CIE, although the estimated contrast is smaller.

Summary and discussion of results
The good performance of BA and the questions it opens.BA gives comparable results with other solvers which take multiple scattering effect into consideration.It shows best computational efficiency and convergence rate due to the reduced nonlinearity and the least number of unknowns to reconstruct.The reconstruction results are not sensitive to the maximum initial contrast value with the given initial contrast distribution.Due to the fact that the reconstruction result of BA totally relies on the data equation, a degradation in performance is observed when the location of the substrate is wrongly estimated, which introduces error into far-field Green's function G s and further results in additional error in the data equation.Therefore, it is sensitive to the estimated location of the substrate, as can be seen from the case where the immersion-substrate interface is estimated 0.5 µm above the real location.
For the situations where the Green's function is not erroneous, the fact that BA competes well with the other solvers is quite and begs more investigation.One may argue that the contrast considered in these test examples is quite small, indicating weak scattering case.However, even though locally the contrast is small at any location, the fact that the sample objects are comparable to the wavelength implies that the net multiple scattering effect is not negligible.In this event, BA provides competitive performance.It essentially indicates that for such a microscopy setup, the reconstruction is largely governed by the far-field measurements, and it is difficult to use the near-field multiple scattering phenomenon for improving the reconstruction.This is validated further by the fact that CIE is the only solver that outperforms BA, but only after we amplified the state equation error by a multiple of 100.While this may suggest that we could blindly endorse using an amplified contribution of the state error equation in the cost function, this is a risky endeavor.In SOM, for example, we compared the result with three different cost functions, i.e. (F S + 0.01F D ) (reported here), (F S + 1F D ) and (F S + 100F D ), and we found that the latter two cost functions did not perform as well as the first one.In fact, even with the first cost function, the value of 0.01F D always exceeded F S , indicating a severe mismatch in the near-field despite a good match in the far-field.This also indicates that the blind amplification of the state error equation is not necessarily useful.
Based on these observations, it is worthwhile to consider then the following questions in the future investigations: (a) does the nature of G s somehow restrict the reconstruction quality by limiting the induced current profile that can satisfy the data equation error?(b) can we modify the measurement system to better endow G s with qualities that support better reconstruction?(c) could we design either the solvers or the relative contribution of the state error equation to be adaptable in a manner such that state equation error can play a powerful role in the reconstruction process without introducing instability or misconvergence?CSI-qualitatively best results and other potential representations.CSI gives a good shape reconstruction qualitatively in our numerical tests.It is not affected by the initial contrast much in the qualitative sense, however quantitatively, we see improved performance when lower initial contrast value is used.Further, the reconstruction performance of CSI is influenced by the estimated location of the substrate.With wrong far-field Green's functions arising from incorrect estimate of the location of the substrate, we observed an increase in the reconstruction error as the optimization proceeds even while the total cost decreases continuously.In general, CSI seems to provide a result that has overfitted the problem, potentially due to many unknowns present in the cost function.And this effect is further more severe if the G s is incorrect.In this sense, is quite motivating to use other representations of the contrast source, such as Fourier representation in CIE and orthonormal basis set such as right singular vectors of G s in SOM.It may be interesting to consider other representations which may in some sense be better.Potentially considering separate basis sets for the contrast sources along the different axes may be a good solution to explore in the future.Or perhaps, instead of looking for representations of the contrast source, we may consider alternative representations of the shapes, differing from the cubic sampling of the DOI in such a manner that the number of unknowns is smaller and the physical coordinates of these unknowns are quite informative.To this end, it may be interesting to consider differential topology for this problem.
CIEs-the most robust-physically grounded for the application domain.CIE shows the best potential among the multiple-scattering based inverse solvers which we tested.It is the only solver that surpassed BA.However, this is possible only when we use a cost function with amplified contribution of the state equation error.
CIE shows the most robust performance to the various location estimation of the substrate for the test result, but this may also be due to that the state equation error has been magnified.When a wrong substrate location is used, additional error is introduced into data equation but the state equation is still correct.
We also observed that when z-components of scattered field and contrast source are considered, there is no obvious advantage over the one where the z-components are ignored in the test examples, which means the balance between the accuracy of the model and computational complexity needs to be considered.In our test, the hyperparameters in CIE are same for all the examples.Properly chosen parameters may give a better result but such work is out of the scope of this paper.
One of the most exciting features of CIE is that it is more intuitive and physically grounded for the microscopy community.This is because of two reasons.The first reason is the use of Fourier representation of the contrast source, since Fourier domain analysis is quite conventional and popular in this community.The second reason is that the value of β in the modified contrast indicates the local region of influence of multiple scattering, which is more intuitive to understand and incorporate.
Subspace based optimization-our ineffective adaptation-failure analysis and future considerations.We noted in section 3 that our adaptation of SOM is the most drastic divergence from the original version.Our results indicate that our adaptation has been quite ineffective, and here we present a breakdown of our investigation into the reasons for this.
We highlighted in section 3 that as opposed to the conventional practice in SOM of using the entire ambiguous current for reconstruction, we only use a subset of ambiguous current for reconstruction.We were driven by the fact that singular values corresponding to the singular vectors indicate the strength of contribution of the associated features in the measurements.However, our results inspired us to look more closely into this hypothesis.Here we take test 2 as an example. 2 , which is the mean energy in the part of the contrast source used for reconstruction relative to the energy in the ground truth contrast source.If the value of Γ is close to 1, the chosen part of the contrast source is sufficient to represent the contrast source.The plot of Γ as a function of L δ is given in figure 14.It is seen that even when up to 2000 coefficients are used to represent the contrast source, it represents only ∼33% of the energy of the actual contrast source.It is also seen that we need to use almost the full basis set for representing the contrast source.
Despite the unsatisfactory performance of our adaptation of SOM, we still see a potential in SOM itself.For instance, with α + obtained from the first step, the first update of contrast actually gives good structural estimate of the object in spite of low contrast, as shown in figure 15.Interestingly, we note that the initial guess of the contrast was set to be invariant along the z-axis.SOM gets rid of this invariance along the axial direction in the first step itself, and retains the lateral structure with good structural details.On the other hand, the optimization of the contrast along that direction is very slow in all the other solvers.We have tested hyperparameter L in the range (12, 60) and the results are similar, but in some cases the results are blurred laterally if L > 60 is used.Yet we did not investigate further but such result can be at least used as initial guess for other inverse algorithms.
These factors indicate that a suitable adaptation of SOM in the future can deliver the performance of SOM for microscopy as well.In this regard, we would like to investigate if retaining the concept of deterministic and ambiguous contrast sources, but using another suitable basis set can be useful.This suggestion is based on the fact that using the entire basis set is computationally expensive and that Γ does have flat zones where the relative increase in the energy is small.This indicates that there is a possibility to represent the contrast source sparsely, however not so in the current basis set.We also think that potentially separating the basis set for contrast source components along the different directions may be a useful step since inherently the basis set for each of these components will be smaller than the basis set for all of them put together.

Conclusion
To solve the problem of 3D RI reconstruction with label-free microscopy by inverse scattering theory, we have adapted several inverse algorithms that incorporate multiple scattering to work with measured data without phase and polarization information for the configuration of a darkfield microscope.
From the numerical tests, we see that it is possible to utilize the algorithms for 3D sample reconstruction with only a few measurements.The algorithms show promising reconstruction results for extended objects, however, yet the performance is not good for objects smaller than the diffraction limit.Our results also show that the inverse algorithms are kind of robust to the choice of maximum initial guess and far-field Green's function of the microscope to some extent, in view that the prior knowledge of the sample and the substrate location may be not accurate.Still, degradation in the reconstruction performance exists when wrong G s is used.
From the numerical results, we see that BA usually gives a good reconstruction result for the test samples despite their sizes being comparable to the wavelength.It is seen that the results of BA are often comparable or better than the inverse scattering solvers that incorporate multiple scattering.So far the numerical examples under investigation here have a low contrast and a comparable size with the wavelength (the largest object has a dimension of ∼ 2λ × 2λ × 0.6λ).In real application, biological samples are usually low-contrast but may have a larger size.In this case, multiple-scattering effect may be further stronger, thus increases the nonlinearity of the problem.Also the increased dimension will lead to more unknowns and pose challenge to the computation process.Thus, more investigation is needed for problem formulation which can decrease the nonlinearity and the improvement of optimization process.We consider this as an open challenge for the inverse problems community to propose better practical microscope designs and/or inverse solvers that can facilitate decoding of the multiple scattering effect for better quantitative accuracy and qualitative reliability of 3D reconstruction of RI profiles of samples.Apart from the investigation in the algorithms, we are also in the process of developing an experimental system such that the algorithms can be validated with experimental data in the future.

Figure 1 .
Figure 1.Configuration (a) and illumination for case N inc = 8 (b) of the microscope under investigation.DoI is on the substrate which encloses all the objects.The illumination propagates parallel with the substrate from various angles.For each illumination, the intensity of the scattered field at different camera pixels is measured which forms a 2D image.
way.Contrast source is expressed as J = [J x ; J y ] and similarly for all field quantities.

Figure 2 .
Figure 2. Configuration (a) and camera images (c) of the first test example.I C is obtained by combining all 8 images and is used to get the initial guess of the contrast.(a) Axis units: µm.(b), (c) Scale bar: 1 µm.

Figure 3 .
Figure 3. Reconstruction results of the test example 1 with maximum contrast of the initial guess set as 0.2.Three cross sections are shown from top to bottom: z = 195 nm, y = −15 nm and x = −15 nm.Axis units: µm.

Figure 4 .
Figure 4. Reconstruction results of the test example 1 with maximum contrast of the initial guess set as 0.05.Same layout as above.Axis units: µm.

Figure 5 .
Figure 5. Reconstruction Error e(χ) and cost function F ′ versus iterations for test 1.The effect of the contrast value of initial guess is shown in (a), (b).Maximum 1000 iterations shown.

Figure 6 .
Figure 6.Configuration (a) and camera images (c) of the second test example.(b) shows the intensity maps obtained by combining all 8 camera images described in paragraph initial guess.(a) Axis units: µm.(b), (c) Scale bar: 1 µm.

Figure 7 .
Figure 7. Reconstruction results of the test example 2 with correct estimation of the substrate location.Three cross sections are shown from top to bottom: z = 195 nm, y = 15 nm and x = −15 nm.Axis units: µm.

Figure 8 .
Figure 8. Illustration of DoIs and reconstruction results of example 2 using G s calculated with (top panel) z sub = 0.5 µm (0.5 µm offset along positive z direction) and (bottom panel) with z sub = −0.5 µm (0.5 µm offset along negative z direction).Axis units: µm.

Figure 9 .
Figure 9.Reconstruction error e(χ) and cost function F ′ versus iterations for test example 2, where the position of the interface may not be known correctly.Maximum 700 iterations shown.

Figure 10 .
Figure 10.Reconstruction results of example 2 with polarized measurements (a) and without polarized measurements (b).Three cross sections are shown from top to bottom: z = 255 nm, y = 15 nm and x = 15 nm.Axis units: µm.

Figure 11 .
Figure 11.Configuration (a) and camera images (c) of the test example 3. The final intensity image is obtained by combining all 8 images and is used to get the initial guess of the contrast.(a) Axis units: µm.(b), (c) Scale bar: 1 µm.

Figure 12 .
Figure 12.Reconstruction results of two beads with a center-to-center distance of 0.3 µm.Three cross sections are shown from top to bottom: z = 105 nm, y = 15 nm and x = 15 nm.Axis units: µm.

Figure 13 .
Figure13.Reconstruction results of plane z = 105 nm by various algorithms with the center-to-center distance between the two beads equal to 300 nm (top), 350 nm (middle) and 400 nm (bottom).When the distance is 300 nm, no approach can resolve the two beads.When the distance increases to 350 nm, all the algorithms can resolve them but some only barely resolve.In contrast, all the algorithms are able to clearly identify the two beads when the distance is 400 nm.

Figure 14 .
Figure 14.Energy in the portion of contrast source considering up to L δ coefficients relative to the actual energy in the contrast source, Γ, plotted as a function of L δ .

Figure 15 .
Figure 15.Contrast reconstruction after the first step of SOM (using only the deterministic part of contrast source) for the test examples 1, 2 and 3 (from left to right).Axis units: µm.

Table 2 .
Comparison of errors between inversion algorithms for example 1.

Table 3 .
Comparison of errors between inverse algorithms for example 2. Here, z sub = 0 corresponds to the correct location.

Table 4 .
Comparison of errors between cases with and without polarization (abbreviated polar. in the columns) sensitive measurement system.