2.1

The inverse scattering problem

The inverse scattering problem consists in retrieving the characteristics of an unknown object by probing it with known incident waves and measuring the scattered waves. In this paper, we consider a transparent object of unknown refractive index n(r) (r = (x, y, z)). The object is described by its the scattering potential defined as F(r) = k²(n(r)² – n₀²)/4π, with n₀ the background refractive index and k the wave vector of the incident light. F is the function that we aim at retrieving. The object is illuminated with a sequence of plane waves at different angles. For every incidence angle, the total field (incident and scattered) is recorded using digital holography, so that we gain knowledge of the amplitude and phase of the field. The field u_s scattered by the object can be calculated using the Lippmann-Schwinger equation:

where u is the total optical field, u_i the incident field.

2.2

Linear scattering

Equation (1) can be linearized if the scattered field u_s is much smaller than the incident field u_i. In this case, u can be replaced by u_i in the integrand, which is referred to as the first Born approximation, and we get:

The assumption that the scattered field is much smaller that the incident field implies that the field resulting from the re-scattering of u_s will be totally negligible. The first Born approximation is thus equivalent to the assumption of single scattering. With the linearized equation (2), it is possible to solve the inverse scattering problem (i.e. retrieve F from the collection of measurements of u) using the direct inversion formula proposed by E. Wolf [8], which we shall hence call the Wolf transform:

where f is the scattering amplitude, k the wave vector in the direction of observation and k_i the wave vector of the incident field. F̃_3D denotes the three-dimensional Fourier transform of the scattering potential. The scattering amplitude is related to the scattered field by u_s^B(r) = f(k) exp(jk|r|)/|r|, for |r| sufficiently large. Note that equation (3) is valid for plane wave illumination. For more complex incident field, we have to linearly add the scattering amplitude corresponding to each Fourier mode of the incident field.

As we will see in the results section, the Wolf transform produce reconstructions that are strongly distorted, typically close to the edges for spheroidal objects like beads or biological cells. The above-mentioned distortions are suppressed efficiently by using the Rytov approximation, which is an alternate way of linearizing equation (1) that is widely used in optical tomography and was initially proposed by Devaney in that context [13]. The underlying idea is to express the field u in term of a complex phase ψ: u = e^jkψ, which is nothing more than a change of variable. We further define u_i = e^jkψ_i and δψ = ψ –ψ_i so that we have u = u_i e^jkδψ. When inserted into Helmholtz equation, this lead to a Riccati equation in term of δψ. The Ryotv approximation consists in retaining only the first term of the Neumann expansion of the Riccati equation [13], which results in a linear equation that happens to have the same form as Helmholtz equation. For that reason, the Wolf transform, which normally applies to the field u, can now be applied to δψ directly. We can write:

For the Born approximation, we take the digital two-dimensional Fourier transform of the field u measured in the plane of the detector, which is related to f by the following equation:

where the tilde denotes the two-dimensional Fourier transform and f^B the scattering amplitude under Born approximation. Note that, this formula is valid for a plane wave incident field: u_i = exp(jk_i·r). For Rytov approximation, we make use of equation (4):

where f^R is the scattering amplitude under Rytov approximation and 𝓕 the Fourier transform. At this point, the procedure is the same for both approximations, f^B or f^R is used to populate the three-dimensional Fourier space of the object according to equation (3).

For each incident plane wave (indexed by p up to the total number of views N), the scattering amplitude f_p is used to populate the Fourier domain of the object function F̃ according to the Wolf Transform (3):

The scattering potential is obtained by taking the inverse three-dimensional Fourier transform of the previous equation:

where Δz is the discretization step in direction z. Finally, the refractive index distribution reconstructions is obtained from F:

according to Born and Rytov approximations respectively.

2.3

Learning tomography

We briefly review here the basic concepts used in LT. A digital map n(x, y, z) of the refractive index (the estimate) of the unknown object is generated in the computer and a forward model H simulating the physics of light propagation through the sample is used to predict the set of measurements m̂_i that would be obtained from the current estimate and a set of incident illumination waves u_{inc, i}. This prediction is compared to the actual measurement m_i performed experimentally on the real sample. The goal of LT is to ‘learn’ the digital map of the refractive index that minimizes the residual r between the actual measurement and the prediction, i.e. find the minimum of the following cost function C:

where D is the data fidelity term and R contains additional constraints such as a regularizing operator. The parameter τ determines the relative weight given to the additional constraints. The estimate of the refractive index will then be given by:

The operator H that maps the refractive index distribution to the measured field is nonlinear. This can be understood from Eq. (1) in which the scattered field is a nonlinear function of the scattering potential since doubling the strength of the scattering potential does not double the strength of the scattered field. Due to the nonlinearity, the cost function is generally not convex. Additionally, the inverse problem of recovering the refractive index from the light measurements is ill posed, which is addressed by regularizing the solution with the isotropic total-variation (TV) of the refractive index contrast and positivity constraints:

where D_j represent the discrete variation in direction j. The TV regularization favors piecewise constant solutions. Indeed, we expect biological tissues to consist mainly of constant refractive index parts separated by sharp boundaries, such as organelles in a cell.

In order to solve the problem posed in Eq. (11), we use the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [14]. The FISTA does not require computing the gradient of the total cost function C, but only the gradient of the data fidelity term D. The regularization is taken care of by solving the following sub-problem at each iteration k:

We use the Spilt-Step Fourier (SSF) Beam Propagation Method (BPM) [11, 12] as a forward model because the gradient of the data term can then be computed efficiently. At each iteration, we randomly select a subset of N_eff (typically 8) views from the set of all the views. We compute the gradient of the data fidelity term D for the m_j in that subset of views with respect to n and we update the refractive index map accordingly:

where n^k is the refractive index estimate at iteration k and s the step size. In addition to the TV regularization, we impose a positivity constraint on the refractive index contrast, i.e. we set n ← max(n, n₀) at each iteration, where n₀ is the background refractive index.

The initial guess, i.e. n^k=0, is taken to be constant and equal to n₀ for weakly scattering objects. For high index contrast, the nonlinearity of H is stronger and the cost function generally has many local minima. In order to converge. Therefore is desirable that we start with an initial guess that is close enough to the solution to avoid the LT algorithm to fall in such minima. The generation of the initial guess is discussed in the next section.

2.4

Linear iterative model

A fair comparison between linear and nonlinear models can only stand if both are subjects to the same type of regularization and positivity constraints. In order to produce such an iterative linear algorithm, we simply take the LT algorithm in which we replace the nonlinear forward model, i.e. H in Eq. (10), by a linear model A. Specifically, we replace the BPM, which is nonlinear, by the Wolf transform, which is linear under the Born or Rytov approximations. These steps are detailed below. Let’s define the indicator function M of the set of tomographic measurements in the Fourier space:

where k_d is the detection direction vector, k_i the incident wave vector, and K the coordinate in the (three-dimensional) Fourier domain of the sample. For convenience, the sample is described here by its scattering potential function F. According to the Wolf transform, all the tomographic measurements lie on the locus of point where M = 1, i.e. the union of all the subset of the Ewald spheres that correspond to the measurements m_i. In that context, we refer to the measurement m as being the collection of the projections of m_i(k_d, k_i) onto the Ewald spheres, i.e. m_i(K). Note that some measurements will fall on the same location K (e.g. if k_d1-k_i1= k_d2-k_i2), in this case, we simply average all the contributions at the same location. The linear forward operator A is a three-dimensional Fourier transform of the scattering potential F followed by a masking operation using the function M:

The cost function is now defined as:

The gradient of the data term D is then:

where r is the residual now only evaluated where M = 1. From this point, the procedure is exactly the same as for the LT algorithm. We then use equations (13) and (14) as in LT in order to enforce the TV and positivity constraints.

2.5

Experimental method

Experimental apparatus

The experimental apparatus is shown in Fig. 1. It consists in a Mach-Zehnder interferometer. The sample is mounted between two 150 micron thick cover slips and sits between two Olympus oil immersion PlanApo 100x NA1.4 objectives in the signal arm of the interferometer. The light source is a continuous wave diode laser at 406nm, with a coherence length of approximately 100 microns. The beam is spatially filtered, expanded and focused in the back focal plane of the illumination objective (IO) thus projecting a plane wave on the sample. A couple of galvo-mirrors are placed before the illumination objective, one steering the beam in the x dimension and one in the y dimension. The surfaces of the galvo-mirror lie in conjugate image planes and are imaged in turn onto the sample. By changing the angles of the galvo- mirrors, we can change the incidence angle of the plane wave illuminating the sample. The sample is imaged with a lateral magnification of 111 through the collection objective (CO) and a relay lens on a detector array (Andor Neo EMCCD, pixel size 6.5 microns), where a hologram is formed through interference with the reference beam. For each view, we measure the total field u with the object in place. We then move the sample together with the cover slips to the side so that the field of view becomes clear of any object. We record then a second field, which is in fact a measurement of the incident field u_i itself. This technique allows us to compensate aberration present in the optical system.

Fig. 1:

Experimental apparatus. The light source is a continuous-wave laser diode at 406nm. The beam is spatially filtered, expanded and split into a reference wave and a signal wave. The signal wave incidence angle can be steered with two galvo-mirrors. The sample is placed between two 100x microscope objectives and imaged on a CCD camera, where a digital hologram is captured by recombining the reference beam.