[figure]font=footnotesize

Sensitivity fields and parameter estimation
from dielectric objects

Abstract

The quantitative phase image formation process is posed as a problem of parameter estimation from intensity measurements. This approach is inclusive of traditional pixel-oriented imaging, where the sought parameters are the pixel values. The resulting optimization process to find the parameters is then seen to depend on the Sensitivity Field: this is the gradient of the scattered field with respect to the parameters, and it turns out to obey a scattering relationship that is analogous to that of the original scattered field. Examples are given from several regimes of scattering strength.

1     Introduction

“Quantitative Phase Imaging” has become, not least thanks to the tireless work of our late colleague and friend Gabi Popescu [1], whose memory this special issue honors, a standard and popular sub-discipline of Imaging Science. It is the topic of dedicated technical meetings and special sessions within bigger conferences, and of course numerous theses and government grants. Taken at face value, the topical term suggests as its ultimate goal the numerically accurate reconstruction of the spatial details of electromagnetic field’s phase.

The question of phase, and of its meaning and utility for imaging, dates back to at least eighty years ago. This was when the origins of present optical, x-ray and electron microscopy methods were being set by Lawrence Bragg and others. In the words of Denis Gabor [2]:

“[…] Bragg’s method, in which a lattice is reconstructed by diffraction from an X-ray diffraction pattern, can be applied only to a rather exceptional class of periodic structures. It is customary to explain this by saying that diffraction diagrams contain information on the intensities only, but not on the phases. The formulation is somewhat unlucky, as it suggests at once that since phases are unobservables, this state of affairs must be accepted. In fact, not only that part of the phase which is unobservable drops out of conventional diffraction patterns, but also the part which corresponds to geometrical and optical properties of the object, and which in principle could be determined by comparison with a standard reference wave.”

Already in this introductory paragraph, and many times subsequently in this remarkable paper, Gabor makes the point that the phase is not a goal in its own right, but rather as the means to recover more information about the object of interest.

Of course, one important difference between Gabor’s time and ours is that then image interpretation was done entirely manually; therefore, it was reassuring for the human processors, scientists and their assistants alike, to be operating on visually detailed images whose fidelity to the true scene would be somehow guaranteed. Nowadays, image exploitation is largely automated, and yet we still treat the processing algorithms “anthropomorphically”; that is, in equal need of spatial richness and accuracy. Even more importantly, the intervening development of Information Theory—with considerable contributions by Gabor himself—has taught us that the information content of an object is limited by its statistical variation. Spatial variability that is irrelevant to a specific information retrieval task is, at best, unnecessarily resource demanding and, in some cases, it may be outright harmful.

As an example, consider the task of localizing a blob-like object whose blob shape is already known. The anthropomorphic way would be to image the blob as accurately as possible, and then find its center of mass. In the 1960’s this would have involved calculators and rulers. Nowadays, with a computer one might first try to “improve” the image as much as possible with regularizing filters and the like; and then apply a centroid or similar algorithm. This seems like a lot of work just to get out a single coordinate—the center of mass. Moreover, even though my students and I followed this exact same approach before and showed that at least our specific instance was accurate [3, 4], in general there is no guarantee that the image “improvement algorithm” is bias-free.

On the other hand, it is possible to imagine methods of determining the center of mass that are impervious to shape specifics. In fact, one of the earliest “optical computing” methods, the van der Lugt correlator [5, sec. 8.4], does just that. Thus, it is evident that a gap exists between the traditional imaging mentality that prioritizes exquisite spatial detail even if only a few parameters are to be extracted from it; and the information theoretic-driven view that acknowledges the sparsity of the sought-after description and orients the retrieval method accordingly.

My purpose in the present paper is to bridge this gap by recasting spatial imaging as a problem of parameter estimation. This includes fully spatially resolved images as special cases. When the parameters of interest are significantly fewer than the number of pixels, certain efficiencies become possible for imaging system design: it is no accident that the parametric view is popular in modalities deprived of the benefit of high pixel-count spatial reconstructions, e.g. the geosciences.

The core development of the parametric approach is in Section 2. It consists of formulating an optimization functional on the detected intensity—which is always the case, even if the estimate of a phase-like quantity is to be sought after later as an intermediate product. Taking gradients on the functional reveals the need for a Sensitivity Field, which is simply defined as the gradient of the detected field to the sought-after parameters. The Sensitivity Field is an artificial, mathematical construct; however, it turns out to obey propagation rules similar to the physical fields whose gradients it represents.

Section 3 derives the Sensitivity Fields for three regimes that are commonly assumed in light-matter interaction: the thin-film approximation, the first Born (or Wolf) kernel, and strong scattering approaches such as the Lippmann-Schwinger equation and the Beam Propagation Method. Interestingly, the Sensitivity Fields themselves scatter similar to and driven by the electromagnetic fields they originated from. These results are analogous, albeit not entirely identical to sensitivity and adjoint analysis of dynamical systems. [6, 7]

Before diving into the technical details, let me note what the paper is not attempting to do. First, it is not meant to be an exhaustive analysis of electromagnetic scattering models, and it could not credibly be in such a short space. Several excellent textbooks cover this topic in detail [8]. The purpose of using a few such scattering models is to show how their Sensitivity Functions are derived.

Second, my goal here is not to compare the validity or fidelity of different parametric descriptions or scattering models. The former requires a careful formulation of detection noise statistics, which is beyond the scope of the present paper. The basics of the latter are in the aforementioned textbooks, and still subject to active research. [9] Lastly, in Section 2 I only use a simple quadratic—that is, optimal for additive Gaussian statistics only—loss function to derive the sensitivity. Examining the behavior of the Sensitivity Field with alternative loss functions, e.g. adversarial, and more generally adoption of modern optimization methods would be interesting topics for follow-up work.

2     Parameter estimation from spatial images

2.1 Intensity detection and the sensitivity field

Figure 1: General scattering geometry with a single intensity detection on a photoelectric transducer.

The geometry to be used in the subsequent analysis is shown in Figure [$\rightarrow$ Make.] Let ${\psi}_{\text{i}}(\boldsymbol{r})$ denote the incident electromagnetic field as function of the Cartesian spatial coordinate $\boldsymbol{r}$; and let $n(\boldsymbol{r};\bm{\theta})$ denote the refractive index of the object whose parameters $\bm{\theta}$ are of interest. The refractive index is assumed to be real, i.e. signifying a pure dielectric object. The parameters $\bm{\theta}$ may be thought of as “controlling” the spatial distribution of the refractive index; some examples are given in Section 2.5 below.

The resulting scattered field is denoted as ${\psi}_{\text{s}}(\boldsymbol{r};\bm{\theta})$ and transduced at $M$ pixel locations $\boldsymbol{r}_{m},$ $m=1,\ldots,M,$ where the intensities

$$g_{m}(\bm{\theta})=\left|{\psi}_{\text{s}}(\boldsymbol{r}_{m};\bm{\theta})\right|^{2}$$ (1)

are recorded. The vector

$$\bm{\mathrm{g}}(\bm{\theta})=\left(\rule[-2.15277pt]{0.0pt}{12.91663pt}\begin{array}[]{ccccc}g_{1}(\bm{\theta})&\cdots&g_{m}(\bm{\theta})&\cdots&g_{M}(\bm{\theta})\end{array}\right)^{\dagger}$$ (2)

is “the ideal measurement,” where $\dagger$ denotes the transpose. The right-hand side in (1) is the the physical model for light propagation through the object for any choice of parameters $\bm{\theta}$. It is also known as the “forward operator.” In optical systems usually the $\boldsymbol{r}_{m}$’s are arranged on a planar 2D lattice of pixels, but in other modalities, e.g. in geo-imaging, they may be randomly placed.

In a noiseless measurement and with perfect knowledge of $\bm{\theta}$, eq. 1 holds with strict equality. In practice, a noisy measurement $\bm{\mathrm{\tilde{g}}}$ will be instead recorded; and the parameter vector $\bm{\theta}$ is to be determined. To carry out this latter task, define a loss function such as

$$E(\bm{\theta})=\frac{1}{4M}\>\sum_{m=1}^{M}\left[\rule[-2.15277pt]{0.0pt}{10.76385pt}\tilde{g}_{m}-\left|{\psi}_{\text{s}}(\boldsymbol{r}_{m};\bm{\theta})\right|^{2}\right]^{2}.$$ (3)

It is useful to reiterate that inside the square bracket the first term is the noisy measurement, whereas the second term is the forward model, parameterized according to $\bm{\theta}$. The quadratic loss function is optimal in the maximum likelihood sense for additive Gaussian detection statistics. Other classical choices, e.g. for Poisson (shot-noise) statistics, or more modern ones such as adversarial loss are possible, but I will not work them out in detail in this paper.

The optimal parameter vector is found from solving the minimization problem

$$\bm{\hat{\theta}}=\>\stackrel{{\scriptstyle[}}{{\bm{}}}{\theta}]{}{\text{argmin}}E(\bm{\theta}).$$ (4)

If additional prior information on the parameters can be captured in the form or a regularizer $\Phi(\bm{\theta})$, then the minimization problem is modified as the unconstrained

$$\bm{\hat{\theta}}=\>\stackrel{{\scriptstyle[}}{{\bm{}}}{\theta}]{}{\text{argmin}}\left\{\rule[-2.15277pt]{0.0pt}{10.76385pt}E(\bm{\theta})+\kappa\>\Phi(\bm{\theta})\right\};$$ (5)

or the inequality-constrained

$$\bm{\hat{\theta}}=\>\stackrel{{\scriptstyle[}}{{\bm{}}}{\theta}]{}{\text{argmin}}\Phi(\bm{\theta})\quad\text{subject to}\quad E(\bm{\theta})<\kappa.$$ (6)

In either option, $\kappa$ is acting as a regularization parameter. In the interest of brevity, the regularized forms (5-6) are largely left unexploited in this paper.

Existence and uniqueness of minima need to be analyzed on a case-by-case basis, depending on the forward problem and its parameterization. If a minimum does exist, it may be found in a number of ways: e.g. genetic algorithms, gradient-based methods, etc. The simplest among the latter, the gradient descent, proceeds iteratively as

$$\bm{\hat{\theta}}[t+1]=\bm{\hat{\theta}}[t]-\eta[t]\,\frac{\partial E\left(\bm{\hat{\theta}}[t]\right)}{\partial\bm{\theta}}.$$ (7)

Here, $t$ denotes the iteration index and $\eta$ is the learning rate. The gradient term is obtained from (3) as

$$\frac{\partial E\left(\bm{\hat{\theta}}[t]\right)}{\partial\bm{\theta}}=\frac{1}{M}\sum_{m=1}^{M}\text{Re}\left\{\,{\psi}_{\text{s}}^{*}(\boldsymbol{r}_{m})\>\frac{\partial{\psi}_{\text{s}}(\boldsymbol{r}_{m})}{\partial\bm{\theta}}\,\right\}\>\left[\rule[-2.15277pt]{0.0pt}{10.76385pt}\tilde{g}_{m}-\left|{\psi}_{\text{s}}(\boldsymbol{r}_{m};\bm{\theta})\right|^{2}\right].$$ (8)

It is evidently convenient at this point to introduce the Sensitivity Field

$${\chi}_{\text{s}}(\boldsymbol{r};\bm{\theta})\equiv\frac{\partial{\psi}_{\text{s}}(\boldsymbol{r};\bm{\theta})}{\partial\bm{\theta}},$$ (9)

which is simply the gradient of the scattered field with respect to the chosen parameterization; and the local error term

$$\Delta_{m}(\bm{\theta})\equiv\tilde{g}_{m}-\left|{\psi}_{\text{s}}(\boldsymbol{r}_{m};\bm{\theta})\right|^{2}.$$ (10)

In terms of these two newly defined quantities, the loss function gradient is written more compactly as

$$\frac{\partial E}{\partial\bm{\theta}}=\frac{1}{M}\sum_{m=1}^{M}\text{Re}\left\{\,\rule[-2.15277pt]{0.0pt}{10.76385pt}{\psi}_{\text{s}}^{*}(\boldsymbol{r}_{m};\bm{\theta})\>{\chi}_{\text{s}}(\boldsymbol{r}_{m};\bm{\theta})\,\right\}\>\Delta_{m}(\bm{\theta}).$$ (11)

2.2 Multiple exposures

Figure 2: General scattering geometry with several intensity detection steps, each time with the illumination and photoelectric transducer at a different orientation or rotation.

For the sake of completeness, I develop this notation a little bit further for the numerous methods that employ multiple exposures onto the same object. Popular examples are confocal microscopy, transport of intensity, ptychography, computed (Radon) tomography, optical diffraction tomography, optical coherence tomography, etc. A drawing suggestive of rotations between exposures, as in tomography, is shown in Figure 2. Let $l$ denote the counter of exposures, entered as superscript in the other notations. Thus, $\boldsymbol{r}_{m}^{(l)}$ is the coordinate of the $m$-th detector pixel at the $l$-th exposure, ${\psi}_{\text{i}}^{(l)}(\boldsymbol{r})$ is the illuminatioun at the $l$-th exposure, etc. The loss function is then expressed as

$$E(\bm{\theta})=\frac{1}{4LM}\sum_{l=1}^{L}w^{2}_{l}\sum_{m=1}^{M}\left[\rule[-2.15277pt]{0.0pt}{10.76385pt}\tilde{g}^{(l)}_{m}-\left|{\psi}_{\text{s}}^{(l)}(\boldsymbol{r}^{(l)}_{m};\bm{\theta})\right|^{2}\right]^{2}.$$ (12)

Here, for full generality we have weighed different exposures with the positive numbers $w^{2}_{l}$. These express relative confidence, for example noisier exposures might receive lower weights.

All the derivations from the previous section follow similarly, finally obtaining

$$\frac{\partial E}{\partial\bm{\theta}}=\frac{1}{LM}\sum_{l=1}^{L}w^{2}_{l}\sum_{m=1}^{M}\text{Re}\left\{\,\rule[-2.15277pt]{0.0pt}{10.76385pt}{\psi}_{\text{s}}^{(l)\,*}(\boldsymbol{r}^{(l)}_{m};\bm{\theta})\>{\chi}_{\text{s}}^{(l)}(\boldsymbol{r}^{(l)}_{m};\bm{\theta})\,\right\}\>\Delta^{(l)}_{m}(\bm{\theta}).$$ (13)

This formulation is appealing because it is compatible with tensorial formulations in computing platforms such as Python and Julia.

2.3 Holographic detection and the sensitivity field

Figure 3: General scattering geometry with holographic detection.

In the holographic configuration, rather than detecting the intensity scattered field itself, one interferes it with a reference field ${\psi}_{\text{r}}(\boldsymbol{r})$. Thus, the intensity at the $m$-th pixel is the interferogram

$$g_{m}=\left|{\psi}_{\text{r}}(\boldsymbol{r})+{\psi}_{\text{s}}(\boldsymbol{r};\bm{\theta})\right|^{2}$$ (14)

It is also common in experiments for the reference power to be much stronger than the scattered field; that is, $|{\psi}_{\text{r}}|\gg|{\psi}_{\text{s}}|$ everywhere. We then obtain

$$\frac{\partial E}{\partial\bm{\theta}}\approx\frac{1}{M}\sum_{m=1}^{M}\text{Re}\left\{\,\rule[-2.15277pt]{0.0pt}{10.76385pt}{\psi}_{\text{r}}^{*}(\boldsymbol{r}_{m};\bm{\theta})\>{\chi}_{\text{s}}(\boldsymbol{r}_{m};\bm{\theta})\,\right\}\>\Delta_{m}(\bm{\theta}).$$ (15)

Whether (15) is an improvement over the reference-free version (11) depends on the nature of the scattered field. For slowly-varying fields, the reference term acts as an amplification factor on the sensitivity, which should be beneficial. On the other hand, the high frequencies in the scattered field would make the interference vary too rapidly, and that might lead to local minima in the loss function. A detailed comparison of direct vs. holographic intensity detection is beyond the scope of this paper, but certainly very interesting for future work.

2.4 Partially coherent illumination

So far the illumination has been assumed to be temporally and spatially coherent. To account for spatial partial coherence, one should replace the field in (1) with

$$g_{m}(\bm{\theta})={\Gamma}_{\text{I}}(\boldsymbol{r}_{m},\boldsymbol{r}_{m};\bm{\theta}).$$ (16)

Here, ${\Gamma}_{\text{I}}(\boldsymbol{r}_{1},\boldsymbol{r}_{2})$ is the spatial correlation function of the scattered field, more usually referred to as mutual intensity, at a general pair of spatial coordinates $\boldsymbol{r}_{1}$, $\boldsymbol{r}_{2}$. This may also be expressed in terms of the coherent modes [$\rightarrow$ Ref] as

$${\Gamma}_{\text{I}}(\boldsymbol{r}_{1},\boldsymbol{r}_{2};\bm{\theta})=\sum_{q=1}^{\infty}\beta_{q}(\bm{\theta})\>\zeta_{q}(\boldsymbol{r}_{1};\bm{\theta})\>\zeta^{*}_{q}(\boldsymbol{r}_{2};\bm{\theta}).$$ (17)

The ideal measurement then is

$$g_{m}(\bm{\theta})=\sum_{q=1}^{\infty}\beta_{q}(\bm{\theta})\left|\zeta_{q}(\boldsymbol{r}_{m};\bm{\theta})\right|^{2}.$$ (18)

Expression (16) or (18) is then to be substituted in the loss function (3) and subsequent formulae. In the interest of keeping the present paper confined, a full analysis of parameter estimation under partially coherent illumination is left for future work.

2.5 On the choice of parameterization

The following three examples elucidate how the parameterization choice influences the formulation of the estimation problem.

• Traditional pixel-oriented spatial imaging. Conceding that at best one may only hope to recover a bandlimited version of the refractive index, let

  $$\left\{\boldsymbol{r}_{j}\right\}_{j=1,\ldots,J}$$

  denote a sampling lattice at the object space, and

  $$n(\boldsymbol{r})=\sum_{j=1}^{J}f_{j}\>\delta(\boldsymbol{r}-\bm{\mathrm{r}}_{j}).$$ (19)

  The coefficients $f_{j}$ are what traditionally one might call the spatial pixel values of the (quantitative) phase in 2D, or of the refractive index reconstruction in 3D. By assembling these values into the parameter vector

  $$\bm{\theta}=\left(\rule[-2.15277pt]{0.0pt}{12.91663pt}\begin{array}[]{ccccc}f_{1}&\cdots&f_{j}&\cdots&f_{J}\end{array}\right)^{\dagger},$$

  one can see that the minimum of (3), if it exists, will correspond to the traditional spatial “quantitative phase image.”

  

  Figure 4: Geometry for estimating the displacement $\boldsymbol{r}_{0}$ of a known blob shape $f(\boldsymbol{r})$.

  Note that if $LM>J$, then the problem is overdetermined, i.e. one has more measurements than unknowns; whereas, if $LM<J$ then it is ill-posed, or underdetermined. Alternatively, (3-4) may be thought of as solving a set of $M$ nonlinear equations with $J$ unknowns.

• Localization of a known shape. Returning to the localization example from the Introduction section, suppose that

  $$n(\boldsymbol{r};\boldsymbol{r}_{0})=f(\bm{\mathrm{r}}-\boldsymbol{r}_{0}),$$ (20)

  where $f(\boldsymbol{r})$ is a known shape, and it is desired to localize it; i.e., the only unknown is the displacement $\boldsymbol{r}_{0}$. This geometry is shown in Figure 4. The parameterization with respect to displacement is

  $$\bm{\theta}=\boldsymbol{r}_{0}.$$

  If, moreover, the imaging system is spatially shift-invariant, then (11) reduces to

  $$\frac{\partial E}{\partial\bm{\theta}}=\frac{1}{2M}\sum_{m=1}^{M}\frac{\partial\left|{\psi}_{\text{s}}(\boldsymbol{r}-\bm{\theta})\right|^{2}}{\partial\boldsymbol{r}}\>\Delta_{m}(\bm{\theta}).$$ (21)

  The problem has become one of finding the best fit displacement coordinate vector to the forward model for the $LM$ intensity measurements. As long as $LM>3$ (or $2$ for planar blobs) this estimation problem is overdetermined.

• Spatial representation in terms of basis functions. Compressive Sensing relies on the existence of a set of basis functions $c_{j}(\boldsymbol{r})$ where the object can be represented efficiently. For the refractive index as object, this is expressed as

  $$n(\boldsymbol{r})=\sum_{j=1}^{J}f_{j}\,c_{j}(\boldsymbol{r})$$ (22)

  where most of the coefficients $f_{j}$ should turn out to be zero. Such sparse bases are discovered by unsupervised learning, e.g. dictionaries [$\rightarrow$ Ref] The pixel-sampled parameterization (19) is seen to be a special albeit not necessarily sparse version of (22). The parameter vector is similarly chosen as

  $$\bm{\theta}=\left(\rule[-2.15277pt]{0.0pt}{12.91663pt}\begin{array}[]{ccccc}f_{1}&\cdots&f_{j}&\cdots&f_{J}\end{array}\right)^{\dagger}.$$

  Now a regularized approach such as (5) or (6) suggests itself, with a suitable choice of sparsity-promoting regularizer, e.g. $\Phi(\bm{\theta})=\left|\!\left|\bm{\theta}\right|\!\right|_{1}.$

3     Sensitivity Fields for some scattering models

The Sensitivity Field (9) in general may be obtained numerically by automatic differentiation. Instead, it is worthwhile to derive explicit expressions for a few common cases of scattering strengths: the thin-film object, the single-event scatterer or First Born approximation, and the multiple scattering regime. The derivations are simple and yield results with a certain intuitive appeal. Moreover, they may provide insights about convergence and fidelity in future real-life applications.

The starting point is the spatially modulated scalar Helmholtz wave equation

$$\nabla^{2}\psi(\boldsymbol{r};\bm{\theta})+k^{2}\,n^{2}(\boldsymbol{r};\bm{\theta})\,\psi(\boldsymbol{r};\bm{\theta})=0,$$ (23)

where $k$ is the free-space wavenumber and $\psi(\boldsymbol{r})$ the appropriate element of the electric field vector. Eq. 23 already neglects polarization effects due to propagation at large numerical apertures and optical anisotropy.

3.1 Point-spread functions

Before developing the Sensitivity Fields, one important detail needs to be taken care of. In the following, ${\psi}_{\text{G}}(\boldsymbol{r})$ denotes the optical system’s Green function, also known as impulse response and point-spread function (PSF). Some typical Green’s functions are listed below and may be replaced in any of the scattering formulae of subsequent Sections, as appropriate.

• The free-space Green’s function

  $${\psi}_{\text{G}}(\boldsymbol{r})=\frac{\exp\left\{ik\left|\boldsymbol{r}\right|\right\}}{ik\left|\boldsymbol{r}\right|},$$ (24)

  which is an ideal spherical wavefront. Since an ideal point source does not exist, this expression does not represent a physically realizable solution to the Helmholtz equation (23) for electromagnetic fields; closest to a realistic solution is the dipole field [$\rightarrow$ Ref.] Nevertheless, for numerical apertures less than $60^{\circ}$ the ideal spherical wave approximation works quite well.

• The paraxial free-space Green’s function, also referred to as the Fresnel kernel

  $${\psi}_{\text{G}}(\boldsymbol{r})=\exp\left\{ik\,z+ik\,\frac{x^{2}+y^{2}}{2z}\right\}.$$ (25)

• The ideal in-focus diffraction-limited system PSF

  $${\psi}_{\text{G}}(x,y)=\text{jinc}\left(\frac{\sqrt{x^{2}+y^{2}}}{{R}_{\text{DL}}}\right),$$ (26)

  where ${R}_{\text{DL}}$ is the radius of the diffraction-limited spot and $\text{jinc}(\cdot)$ is the same function as in [$\rightarrow$ Ref specific equation in Goodman].

3.2 “Phase objects” and the thin-film approximation

Loosely speaking, a phase object is generally understood as a spatial material distribution that interacts with light in such a way as to impart phase delay only on the wavefield; no attenuation is to be imparted anywhere. A good approximation for such an object is an infinitessimally thin film yet with a thickness along the optical axis that is modulated as $d(x,y;\bm{\theta})$ as function of the lateral spatial coordinates $x,y$; and whose index of refraction $n(x,y;\bm{\theta})$ is purely real, and also modulated. This configuration is depicted in Figure 5.

It is customary to define the optical path difference (OPD)

$$\phi(x,y;\bm{\theta})=k\,\left[\rule[-2.15277pt]{0.0pt}{10.76385pt}n(x,y;\bm{\theta})-1\right]\,d(x,y;\bm{\theta})$$ (27)

and say that the wavefront itself is modulated along the lateral dimensions $(x,y)$ according to the OPD. That is, the field emanated directly after the rightmost edge of the thin transparency, at the vertical dashed line in Figure 5, is

$${\psi}_{\text{s}}(\boldsymbol{r};\bm{\theta})={\psi}_{\text{i}}(\boldsymbol{r})\,\exp\left\{\rule[-2.15277pt]{0.0pt}{10.76385pt}i\,\phi(x,y;\bm{\theta})\right\}.$$ (28)

A brief proof of (28) including the list of conditions that have to be met for it to be acceptable is in the Appendix.

Figure 5: The thin transparency is an object that imparts pure phase modulation onto an incident field.

Past the sample, the field propagates either in free space or through an optical system till it reaches the detector. In both cases, we denote the point-spread function as ${\psi}_{\text{G}}(\boldsymbol{r})$, as discussed in Section 3.1. The final field on the detector is

$${\psi}_{\text{d}}(\boldsymbol{r};\bm{\theta})=\iint{\psi}_{\text{i}}(x^{\prime},y^{\prime},z)\,\exp\left\{\rule[-2.15277pt]{0.0pt}{10.76385pt}i\,\phi(x,y;\bm{\theta})\right\}{\psi}_{\text{G}}(x-x^{\prime},y-y^{\prime})\>\text{d}x^{\prime}\text{d}y^{\prime}.$$ (29)

Note that for applications utilizing the thin film approach, such as imaging of biological cells, it may be necessary to parameterize both refractive index and thickness. This is even though they cannot be decoupled from a single intensity measurement and in the absence of other priors, as is evident from (27).

Leaving this issue aside, the Sensitivity Field is readily obtained as

$${\chi}_{\text{d}}(\boldsymbol{r};\bm{\theta})=\iint{\chi}_{\text{i}}(x^{\prime},y^{\prime},z;\bm{\theta})\,\exp\left\{\rule[-2.15277pt]{0.0pt}{10.76385pt}i\,\phi(x,y;\bm{\theta})\right\}{\psi}_{\text{G}}(x-x^{\prime},y-y^{\prime})\>\text{d}x^{\prime}\text{d}y^{\prime},$$ (30)

where the “Sensivity illumination” field is

$${\chi}_{\text{i}}(\boldsymbol{r};\bm{\theta}) ={\psi}_{\text{i}}(\boldsymbol{r})\>\frac{\partial\phi(x,y;\bm{\theta})}{\partial\bm{\theta}}$$ (31)

$$={\psi}_{\text{i}}(\boldsymbol{r})\>\left\{\rule[-4.30554pt]{0.0pt}{12.91663pt}\frac{\partial n(x,y;\bm{\theta})}{\partial\bm{\theta}}\,d(x,y;\bm{\theta})+\left[\rule[-2.15277pt]{0.0pt}{10.76385pt}n(x,y;\bm{\theta})-1\right]\frac{\partial d(x,y;\bm{\theta})}{\partial\bm{\theta}}\right\}.$$ (32)

Evidently, the Sensitivity Field is the same as the field produced by the same object when the illumination is modulated by the OPD sensitivity $\partial\phi/\partial\bm{\theta}$. This interesting symmetry will recur in Section 3.3 with the first Rytov approximation (43) and in Section 3.4 with the Lippmann-Schwinger integral equation (45).

Alternatively, the terms within the integral may be regrouped as

$${\chi}_{\text{d}}(\boldsymbol{r};\bm{\theta})=\iint{\psi}_{\text{i}}(x^{\prime},y^{\prime},z;\bm{\theta})\,\exp\left\{\rule[-2.15277pt]{0.0pt}{10.76385pt}i\,{\chi}_{\text{t}}(x,y;\bm{\theta})\right\}{\psi}_{\text{G}}(x-x^{\prime},y-y^{\prime})\>\text{d}x^{\prime}\text{d}y^{\prime},$$ (33)

where now the phase transmissivity is

$${\chi}_{\text{t}}(x,y;\bm{\theta})=\phi(x,y;\bm{\theta})-i\log\left(\frac{\partial\phi(x,y;\bm{\theta})}{\partial\bm{\theta}}\right).$$ (34)

Due to the imaginary part in this expression, this is no longer a pure “phase object;” it includes an amplitude modulation. An analogous result occurs in Section 3.3 with the first Born approximation (38).

3.3 First Born and First Rytov approximations

The Born series is a rigorous, if not necessarily stable numerically method to obtain a solution for the scattered field in (23). Keeping the first term only implies that multiple scattering is negligible. This results from the substitution

$$\psi(\boldsymbol{r};\bm{\theta})={\psi}_{\text{i}}(\boldsymbol{r})+{\psi}_{\text{s}}(\boldsymbol{r};\bm{\theta})$$ (35)

for the total field, and the approximation

$${\psi}_{\text{s}}(\boldsymbol{r};\bm{\theta})\approx\int{\psi}_{\text{i}}(\boldsymbol{r})\,v(\boldsymbol{r^{\prime}};\bm{\theta})\,{\psi}_{\text{G}}(\boldsymbol{r}-\boldsymbol{r^{\prime}})\,\text{d}\boldsymbol{r^{\prime}}.$$ (36)

Here,

$$v(\boldsymbol{r};\bm{\theta})\equiv\frac{1}{2}\left[n^{2}(\boldsymbol{r};\bm{\theta})-1\right]$$ (37)

is the scattering potential. Result (36) is known as the First Born approximation or Born integral and more recently it has started being referred to as the Wolf transform. [$\rightarrow$ Ref Wolf 1967]

It is straightforward to obtain the Sensitivity Field for the First Born approximation as

$${\chi}_{\text{s}}(\boldsymbol{r};\bm{\theta})=\int{\psi}_{\text{i}}(\boldsymbol{r^{\prime}})\>{\chi}_{\text{v}}(\boldsymbol{r^{\prime}};\bm{\theta})\>{\psi}_{\text{G}}(\boldsymbol{r}-\boldsymbol{r^{\prime}})\,\text{d}\boldsymbol{r^{\prime}},$$ (38)

where the modified potential

$${\chi}_{\text{v}}(\boldsymbol{r};\bm{\theta})=\frac{\partial n(\boldsymbol{r};\bm{\theta})}{\partial\bm{\theta}}\>n(\boldsymbol{r};\bm{\theta}).$$ (39)

Here, as in (33), the Sensitivity Field is produced by the same illumination as the original scattered field but with a modified scattering potential.

One alternative to the Born series is the Rytov perturbation expansion, which starts by the transformation ${\psi}_{\text{s}}=\text{e}^{{\phi}_{\text{s}}}$. To a first-order approximation then the scattered field is obtained as

$$\psi(\boldsymbol{r}) \approx{\psi}_{\text{i}}(\boldsymbol{r})\>\exp\left\{\rule[-2.15277pt]{0.0pt}{10.76385pt}\phi_{1}(\boldsymbol{r})\right\},\qquad\text{where}\qquad$$ (40)

$$\phi_{1}(\boldsymbol{r};\bm{\theta}) =\int\frac{{\psi}_{\text{i}}(\boldsymbol{r^{\prime}})}{{\psi}_{\text{i}}(\boldsymbol{r})}\>v(\boldsymbol{r^{\prime}};\bm{\theta})\>{\psi}_{\text{G}}(\boldsymbol{r}-\boldsymbol{r^{\prime}})\>\text{d}\boldsymbol{r^{\prime}}.$$ (41)

This result is known as the Rytov approximation or Rytov integral. Consequently, its Sensitivity Field is

$${\chi}_{\text{s}}(\boldsymbol{r};\bm{\theta}) ={\chi}_{\text{i}}(\boldsymbol{r};\bm{\theta})\>\exp\left\{\rule[-2.15277pt]{0.0pt}{10.76385pt}\phi_{1}(\boldsymbol{r})\right\},\qquad\text{where}\qquad$$ (42)

$${\chi}_{\text{i}}(\boldsymbol{r};\bm{\theta}) =\int\frac{{\psi}_{\text{i}}(\boldsymbol{r^{\prime}})}{{\psi}_{\text{i}}(\boldsymbol{r})}\>{\chi}_{\text{v}}(\boldsymbol{r^{\prime}};\bm{\theta})\>{\psi}_{\text{G}}(\boldsymbol{r}-\boldsymbol{r^{\prime}})\>\text{d}\boldsymbol{r^{\prime}}$$ (43)

is the sensitivity illumination field and ${\chi}_{\text{v}}$ is again given by (39). In this case, as in (30), the Sensitivity Field results from a modified illumination acting on the same scattering potential as in the original scattered field.

3.4 Multiple scattering via the Lippmann-Schwinger equation

The Lippmann-Schwinger integral equation captures multiple scattering orders and is formally equivalent with, albeit usually better behaved numerically than the Born series. The scattered field is obtained from

$${\psi}_{\text{s}}(\boldsymbol{r})={\psi}_{\text{i}}(\boldsymbol{r})+\int{\psi}_{\text{s}}(\boldsymbol{r^{\prime}})\,v(\boldsymbol{r^{\prime}})\,{\psi}_{\text{G}}(\boldsymbol{r}-\boldsymbol{r^{\prime}})\,\text{d}\boldsymbol{r^{\prime}}.$$ (44)

The scattering potential $v(\boldsymbol{r})$ is the same as in (37).

Taking gradients on both sides, the Sensitivity Field for the Lippmann-Schwinger integral equation is found to satisfy

$${\chi}_{\text{s}}(\boldsymbol{r};\bm{\theta})={\chi}_{\text{i}}(\boldsymbol{r})+\int{\chi}_{\text{s}}(\boldsymbol{r^{\prime}};\bm{\theta})\,v(\boldsymbol{r^{\prime}};\bm{\theta})\,{\psi}_{\text{G}}(\boldsymbol{r}-\boldsymbol{r^{\prime}})\,\text{d}\boldsymbol{r^{\prime}},$$ (45)

where the Sensitivity illumination field is

$${\chi}_{\text{i}}(\boldsymbol{r};\bm{\theta})=\int{\psi}_{\text{s}}(\boldsymbol{r^{\prime}};\bm{\theta})\>{\chi}_{\text{v}}(\boldsymbol{r};\bm{\theta})\>{\psi}_{\text{G}}(\boldsymbol{r}-\boldsymbol{r^{\prime}})\,\text{d}\boldsymbol{r^{\prime}}.$$ (46)

Once again one arrives at a situation like (32) and (43), where the Sensitivity Field is scattered by the same potential as the original scattered field with an illumination that results from the original scattered field modulated by the gradient of the scattering potential.

3.5 Multiple scattering via the Beam Propagation Method

If back scattering may be neglected, then a popular alternative is the Beam Propagation Method. In [9] it was shown that the Beam Propagation Method approximates the Lippmann-Schwinger equation solution well for weakly scattering objects, and the error grows as the fourth power of the maximum scattering angle, i.e. the spatial bandwidth of the scatterer.

In BPM, the object is divided into multiple slices oriented perpendicular to the directional axis $z$. Each slice is treated like a thin transparency, and subsequently the field is forward-propagated as if through free space to the next slice. That is,

$$\psi(x,y,z_{l+1})=\int\psi(x^{\prime},y^{\prime},z_{l})\,\exp\left\{\rule[-2.15277pt]{0.0pt}{10.76385pt}ik\,\varepsilon(x^{\prime},y^{\prime};z_{l})\right\}\\ {\psi}_{\text{G}}(x^{\prime}-x,y^{\prime}-y,z^{\prime}-z)\,\text{d}x^{\prime}\text{d}y^{\prime}.$$ (47)

The recursion is initialized with

$$\psi(x,y,z_{1})={\psi}_{\text{i}}(x,y,z_{1}),$$ (48)

the incident illumination. Since the field from a given slice can reach the next slice at a considerably large angle, usually the ideal spherical wave kernel (24) is used in BPM.

The sensitivity analysis on (47) was carried out in [$\rightarrow$ Ref Kamilov] leading to an elegant analogy between the sensitivity equation (11) for BPM and backpropagation in neural networks. Even though these papers considered only the pixel-oriented spatial imaging parameterization, #1 in Section 2.5, it is straightforward to generalize them to other parameterizations.

4     Conclusions

There are two main messages from this paper. The first is that treating the quantitative phase imaging problem as one of parameter estimation allows some interesting generealizations that are less “anthropomorphic.” Thus, perhaps accuracy improvements and efficiency savings may result. The parametric approach is inclusive of the traditional spatial imaging way of thinking, as the reader has seen in case #1 of Section 2.5.

The second main point is the Sensitivity Field. It results in straightforward fashion from taking gradients on the loss function, eqs. 3-9, with respect to the sought parameters. Superficially, thanks to automatic differentiation it may appear that there is no reason to carry out all the derivations of Sensitivity Fields in Section 3. I beg to disagree. The results reveal an interesting intuition: the Sensitivity Fields are either produced by the same given scattering potential but with a modified Sensitivity Illumination field that is computed from the original scattered field, as in eqs. 30, 42, 45; or by the same illumination field but instead with a scattering potential that is the gradient of the original potential, as in eqs. 33 and 38. Such analogies are well known in the community of feedback control [6, 7] and have led to interesting innovations. [10]

Acknowledgments

This research was funded by the U. S. Air Force Office of Scientific Research through the Multi-University Research Initiative (MURI) “Searching for what’s new: the systematic development of dynamic x-ray microscopy,” grant No. FA9550-23-1-0284; and by Singapore’s National Research Foundation through the Intra-Create thematic grant “Retinal Analytics through Machine learning aiding Physics” (RAMP), grant No. NRF2019-THE002-0006. The author is grateful to Steven G. Johnson of the MIT Mathematics Department; Chris Jacobsen of the Northwester University Physics Department; Shen Zuowei and Ji Hui of the National University of Singapore’s Department of Mathematical Sciences; Fabian Bräu, Bingyao Tan, Michael Girard, Leopold Schmetterer and Aung Tin of the Singapore Eye Research Institute; Sandip Mondal of SMART; and Qihang Zhang of Tsinghua University for helpful discussions.

Appendix: Derivation of the thin film approximation

Starting with the space-variant Helmholtz equation (23), we elect a directional optical beam along the $z$ axis as

$$\psi(\boldsymbol{r})=\tilde{\psi}(\boldsymbol{r})\>\text{e}^{ikz}.$$ (49)

Here, $\tilde{\psi}(\boldsymbol{r})$ is a slow-varying modulation , i.e.

$$\left|\frac{\partial\tilde{\psi}}{\partial w}\right|\ll k\quad\text{for\ }w=x,y,z.$$ (50)

Substituting (49) into (23) we obtain

$$\frac{\partial^{2}\tilde{\psi}}{\partial{x}^{2}}+\frac{\partial^{2}\tilde{\psi}}{\partial{y}^{2}}+\frac{\partial^{2}\tilde{\psi}}{\partial{z}^{2}}+2ik\,\frac{\partial\tilde{\psi}}{\partial z}+k^{2}\,\left[\rule[-2.15277pt]{0.0pt}{10.76385pt}n^{2}(\boldsymbol{r})-1\right]\,\tilde{\psi}(\boldsymbol{r})=0.$$ (51)

We now make the strong assumption

$$\left|\frac{\partial^{2}\tilde{\psi}}{\partial{w}^{2}}\right| \ll k\>\frac{\partial\tilde{\psi}}{\partial z},\quad\text{for\ }w=x,y,z;\qquad\text{and}\qquad$$ (52)

$$\frac{\partial n}{\partial z} =0.$$ (53)

In Geometrical and Hamiltonian optics, condition (53) is also known as a “guide medium.” The remaining terms in (51) then have the solution

$$\tilde{\psi}(x,y,z)=\tilde{\psi}_{0}(x,y)\,\exp\left\{ik\,\frac{n^{2}(x,y)-1}{2}\,z\right\},$$ (54)

where $\tilde{\psi}_{0}(x,y)\equiv\tilde{\psi}(x,y,0)$ is the incident field. To bring this into an even more recognizable form, we further assume a weak refractive index modulation

$$n(x,y)=1+\varepsilon(x,y),\qquad\text{where}\qquad\left|\varepsilon(x,y)\right|\ll 1.$$ (55)

Then (54) becomes

$$\tilde{\psi}(x,y,z)=\tilde{\psi}_{0}(x,y)\,\exp\left\{\rule[-2.15277pt]{0.0pt}{10.76385pt}ik\,\varepsilon(x,y)\,z\right\}.$$ (56)

In the exponent, we recognize the optical path difference

$$\text{OPD}(x,y)=\varepsilon(x,y)\,z.$$ (57)

Clearly, the thin dielectric object under these approximations imparts a pure phase modulation on the incident field. This is a key consequence and justification for the term “phase object.”

To summarize, the phase object or thin film approximation (56) relies on all of the following assumptions to be true:

• The refractive index is independent of the axial variable $z$, according to (53).

• The refractive index modulation is small relative to the surrounding medium, according to (55).

• The following three conditions that are necessary for (56) to be consistent with (52):

  $$\left|\frac{\partial^{2}\tilde{\psi}_{0}(x,y)}{\partial{w}^{2}}\right| \ll 2k^{2}\varepsilon(x,y),$$ (58)

  $$\left|\frac{\partial\tilde{\psi}_{0}}{\partial x}(x,y)\frac{\partial\varepsilon(x,y)}{\partial w}z\right| \ll k\,\varepsilon(x,y),$$ (59)

  $$\left|\tilde{\psi}_{0}(x,y)\frac{\partial^{2}\varepsilon(x,y)}{\partial{w}^{2}}z\right| \ll 2k\,\varepsilon(x,y),$$ (60)

where the spatial derivatives are with respect to $w=x,y$. Without going into a full analysis of these expressions, one might just point out that as the depth $z$ increases, so does the potential of (58-60) becoming violated. The same holds for the lateral bandwidth of the incident illumination $\tilde{\psi}_{0}(x,y)$ and the modulation $\varepsilon(x,y)$.

Remarkably, biological cells seem to satisfy these smoothness conditions most of the time. This is probably why QPI has been fairly succesful in visualizing them and identifying correlates to biological functions. On the other hand, it is also cutomary to apply the thin film approximation to certain artificial objects that seem to break it, e.g. if they include sharp edges. This is the case for almost all binary optics and 3D-printed dielectric objects, for example. Even more remarkably, the thin-film approximation even then seems to predict the output of the overall optical system in good agreement with experiments.

It is possible that such objects in their physical realization are not actually as sharp as their models assume. Moreover, usually the light scattered from the highly sloped edges is directed outside the aperture of the optical system. Its absence from the detected intensity would appear to cause no damage other than smoothening; which may just be ignored or corrected by use of prior knowledge, i.e. an edge-enhancing regularizer such as Total Variation (TV).