Physics-based Learned Design: Optimized Coded-Illumination for Quantitative Phase Imaging

I Introduction

Quantitative Phase Imaging (QPI) enables stain-free and label-free microscopy of transparent biological samples in vitro . When compared with coherent methods , QPI methods that use partially coherent light achieve higher spatial resolution, more light throughput, and reduced speckle artifacts. Phase contrast may be generated using interference or defocus . More recently, coded-illumination microscopy has been demonstrated as an accurate and inexpensive QPI scheme. To realize coded-illumination, we replace a commercial microscope’s illumination unit with a light-emitting diode (LED) domed array (see Fig. 1) . This provides a flexible hardware platform for various QPI applications including super-resolution , multi-contrast , and 3D imaging .

Coded-illumination microscopy uses asymmetric source patterns and multiple measurements to retrieve 2D phase information. Quantitative Differential Phase Contrast (qDPC), for example, captures four measurements with rotated half-circle source patterns, from which the phase is computationally recovered using a partially coherent linearized model. The practical performance of qDPC is predominantly determined by how the phase information is encoded in (via coded-illumination patterns) and decoded from (via phase recovery) the intensity measurements.

The half-circle illumination designs of qDPC were derived analytically based on a Weak Object Approximation which linearizes the physics in order to make the inverse problem mathematically convenient. This linearized model enables one to derive a phase transfer function and analyze the spatial frequency coverage of any given source pattern ; however, the non-linearity of the exact model makes it impossible to predict an optimal source design without knowing the sample’s phase a priori. In addition, these types of analysis are inherently restricted to linear reconstruction algorithms and will not necessarily result in improved accuracy when the phase is retrieved via non-linear iterative methods.

Motivated by the success of deep learning for image reconstruction problems , data-driven approaches have been adopted for learning coded-illumination patterns. For instance, researchers have used machine learning to maximize the phase contrast of each coded-illumination measurement , to improve accuracy on classification tasks , and to reconstruct phase . All of these techniques learn the input-output relationship with a deep convolutional neural network (CNN) using training data. It is not straightforward to include the well-characterized system physics; hence, the CNN is required to learn both the physical measurement formation and the phase reconstruction process. This task requires training of 10s to 100s of thousands of parameters and an immense number of training examples.

Here, we introduce a new data-driven approach to optimizing the source pattern design for coded-illumination phase retrieval by directly including both the system physics and the non-linear nature of a reconstruction algorithm in the learning process. Our approach unrolls the iterations of a generic non-linear reconstruction algorithm to construct an unrolled network . Similar to CNNs, our unrolled network consists of several layers (one for each iteration); however, in our case each layer consists of well-specified operations to incorporate measurement formation and sparse regularization, instead of standard operations such as generic convolutions. The key aspects of our approach are:

incorporation of the system physics and reconstruction non-linearities in the illumination design process.

efficient parameterization of the unrolled network.

reduced number of training examples required.

We deploy our data-driven approach to learn improved coded-illumination patterns for phase reconstruction. Each layer of the unrolled network is parameterized by only a few variables (LED brightness values), enabling an efficient use of training data ($<100$ simulated training examples). We compare the QPI performance of our learned designs to previous work and demonstrate that our designs generalize well to the experimental setting with biological samples.

II Quantitative Phase Imaging

qDPC recovers a sample’s complex transmittance function from several coded-illumination measurements. The phase recovery optimization algorithm aims to minimize the Euclidean norm of the error between the measurements and the expected measurements formed with the current phase estimate. Using a gradient-based procedure, the phase estimate is iteratively updated until convergence. For a partially coherent source, the phase can be recovered with resolution up to twice the coherent diffraction limit. In this section, we describe the measurement formation process and phase recovery optimization.

A thin sample’s transmission function can be approximated as a 2D complex function, $o(\mathbf{r})=e^{j\phi(\mathbf{r})-\mu(\mathbf{r})}$, characterized by its absorption, $\mu(\mathbf{r})$, and phase, $\phi(\mathbf{r})=\frac{2\pi}{\lambda}\Delta n(\mathbf{r})d(\mathbf{r})$, where $\mathbf{r}$ are 2D spatial coordinates, $\lambda$ is the wavelength of the illumination, $d(\mathbf{r})$ is the physical thickness of the sample, and $\Delta n(\mathbf{r})$ is the change in refractive index from the background. Intensity measurements, $y(\mathbf{r})$, of the sample are a non-linear function of $o(\mathbf{r})$, mathematically described by,

where $|\cdot|^{2}$ denotes squared absolute value, $*$ denotes convolution, $\odot$ denotes elementwise multiplication, $s(\mathbf{r})$ is the illumination’s complex-field at the sample plane and $p(\mathbf{r})$ is the point spread function (PSF) of the microscope. The illumination from each LED is approximated as a tilted plane wave, $s(\mathbf{r})=e^{\frac{j}{\lambda}\mathbf{u}_{pos}^{T}\mathbf{r}}$, with tilt angle, $\mathbf{u}_{pos}$, determined by the physical position of the LED relative the microscope .

Because the measured image in Eq. 1 is non-linear with respect to the sample’s transmission function, recovering phase generally requires non-convex optimization. However, biological samples in closely index-matched fluid have a small scatter-scatter term. This means that a weak object approximation can be made; linearizing the measurement formation model such that phase recovery requires only a linear deconvolution of the measurements with their respective weak object transfer functions (WOTFs) . Further, unstained biological samples are predominantly phase objects since they are only weakly absorbing (i.e. $\mu(\mathbf{r})$ is small). With these approximations, we can express each intensity measurement as a linear system with contributions from the background and phase contrast. In Fourier space,

where $\widehat{\cdot}$ denotes Fourier transform, $\mathbf{u}$ are 2D spatial-frequency coordinates, $B$ is the measurement’s background energy concentrated at the DC and $h(\mathbf{u})$ is the phase WOTF. The phase WOTFs are a function of the illumination source and the pupil distribution of the microscope . For a single LED the WOTF is:

In , multiple LEDs are turned on simultaneously to increase signal-to-noise (SNR) and improve phase contrast. Because the fields generated by each LED’s illumination are spatially incoherent with each other, the measurement from multiple LEDs will simply be the weighted sum of each LED’s individual measurement, where the weights correspond to the LEDs’ brightness values. The phase WOTF for illumination by multiple LEDs will also be the weighted sum of the single-LED phase WOTFs. Mathematically,

where $\mathcal{W}$ is the set of LEDs turned on and $c_{w}\geq 0$ are the LEDs’ brightness values.

Following common practice , we discretize the 2D spatial distributions and format them as vectors (bold lower case) (e.g. $\widehat{\mathbf{h}}$ represents the transfer function’s 2D spatial-frequency distribution and $\boldsymbol{\phi}$ represents the 2D spatial phase distribution). The measurementsIn practice, $\mathbf{y}$ typically refers to the so-called flattened image, where the background energy in (2) is removed via background subtraction. are described in Fourier space as $\widehat{\mathbf{y}}=\mathbf{A}\widehat{\boldsymbol{\phi}}$ with system function $\mathbf{A}=diag(\widehat{\mathbf{h}})$.

Based on this model, we define $\mathbf{Y}\in\mathds{R}^{M\times S}$ as the Fourier transform of $S$ single LED measurements, $\widehat{\mathbf{y}}$, along the columns. Then, $\mathbf{C}\in\mathds{R}^{S\times K}$ is defined as the $S$ single-LED weights for each of $K$ measurements, and $\mathbf{c}_{k}\in\mathds{R}^{S}$ is the $k^{th}$ column of $\mathbf{C}$. The product $\widehat{\mathbf{y}}_{k}=\mathbf{Y}\mathbf{c}_{k}$ simulates the $k^{th}$ multiple-LED measurement. Similarly, we define $\mathbf{H}\in\mathds{R}^{N\times S}$ as $S$ single LED phase WOTFs, $\widehat{\mathbf{h}}$ along the columns, such that the product $\mathbf{A}_{k}=diag(\mathbf{H}\mathbf{c}_{k})$ gives the corresponding multiple-LED phase WOTF for the $k^{th}$ measurement.

II-B Phase Recovery

Phase recovery using the forward model in Sec. II-A can be formulated as a regularized linear inverse problem,

where $\boldsymbol{\phi}^{\star}$ is the recovered phase, $K$ is the number of measurements acquired, $\mathbf{\widehat{y}}_{k}$ is the Fourier transform of the $k^{th}$ measurement and $\mathcal{P}(\cdot)$ is a user-chosen regularizer. We solve this optimization problem efficiently using the accelerated proximal gradient descent (APGD) algorithm by iteratively applying an acceleration update, a gradient update and a proximal update . The algorithm is detailed in Alg. 1, where $\alpha$ is the gradient step size, $N$ is the number of iterations, $\mathbf{s}$ and $\mathbf{z}$ are intermediate variables, $\mu^{(n)}$ is the acceleration parameter derived by the recursion, $\mu^{(n)}=\frac{1+\sqrt{1+4\mu^{(n-1),2}}}{2}$ , and $\text{prox}_{\mathcal{P}}(\cdot)$ is the proximal operator corresponding to the user-chosen regularizer $\mathcal{P}(\cdot)$ .

III Physics-Based Learned Design

Given the phase recovery algorithm in Sec. II-B, we now describe our main contribution of learning the coded-illumination designs for a given reconstruction algorithm and training set.

Traditionally, DNNs contain many layers of weighted linear mixtures and non-linear activation functions . Here, we consider specific linear functions which capture the system physics of measurement formation and specific non-linear activation functions which promote sparsity. Starting from Alg. 1, we treat each iteration as a layer such that when unrolled they form a network of $N$ layers, denoted $\mathcal{R}$ (Fig. 2). Each layer of $\mathcal{R}$ contains a module for each of the iterative algorithm’s updates (i.e. an acceleration module, a gradient module (incorporates system physics), and a proximal module (incorporates sparsity)). The regularization and step size parameters specified for Alg. 1 are fixed. The network’s inputs comprise $(\widehat{\mathbf{y}}_{k})_{k=1}^{K}$ and the network’s output is $\widehat{\boldsymbol{\phi}}^{(N)}$. The design parameters of the network, which will be learned, govern the relative brightness of the LEDs and are incorporated in the measurement formation and the system WOTFs.

III-B Learning Objective

Our learning objective is to minimize the phase reconstruction error of the training data over the space of possible LED configurations, subject to constraints that enforce physical feasibility and eliminate degenerate and trivial solutions:

Here, $(\mathbf{Y}_{l},\boldsymbol{\phi}^{\prime}_{l})_{l=1}^{L}$ are $L$ training pairs for which $\mathbf{Y}_{l}$ is a matrix of the Fourier transform of single-LED measurements for the $l^{th}$ sample with optical phase, $\boldsymbol{\phi}^{\prime}_{l}$. $\odot$ is the elementwise product operator, $\mathbf{m}_{k}$ is a geometric constraint mask for the $k^{th}$ measurement, and $\mathbf{0}$ is the null vector.

The non-negativity constraint (Eq. 9) prevents non-physical solutions by enforcing the brightness of each LED to be greater than or equal to zero. This is enforced by projecting the parameters onto the set of non-negative real numbers. The scale constraint (Eq. 10) enforces that each coded-illumination design must have weights with sum equal to 1, in order to eliminate arbitrary scalings of the same design. This is enforced by scaling the parameters for each measurement such that their sum is one. The geometric constraint (Eq. 11) enforces that the coded-illumination designs do not use conjugate-symmetric LED pairs to illuminate the sample within the same measurement, since these would also result in degenerate solutions (e.g. two symmetric LEDs produce opposite phase contrast measurements that would cancel each other out). To prevent this, we force the source patterns for each measurement to reside within only one of the major semi-circle sets (e.g. top, bottom, left, right). This constraint is enforced by setting the LED brightnesses outside the allowed semi-circle to zero.

We solve Eq. 8 iteratively via accelerated projected gradient descent (Alg. 2). At each iteration, the coded-illumination design for each measurement is updated with the analytical gradient, projected onto the constraints (denoted by $\mathcal{B}(\cdot)$) and updated again with a contribution from the previous iteration (weighted by $\beta^{(t)}$). $\mathcal{B}(\cdot)$ enforces the constraints in the following order: non-negativity, geometric, and scale.

III-C Gradient Update

The gradient of the loss function (Eq. 8) with respect to the design parameters has contributions at every layer of the unrolled network through both the measurement terms, $\widehat{\mathbf{y}}_{k}$, and the phase WOTF terms, $\mathbf{A}_{k}$, for each measurement $k\in\{1...K\}$. Here, we outline our algorithm for updating the coded-illumination design weights via a two-step procedure: backpropagating the error from layer-to-layer and computing each layer’s gradient contribution. For simplicity, we outline the gradient update for only a single training example, $l$, as the gradient for all the training examples is the sum of their individual gradients.

Unlike pure gradient descent, where each iteration’s estimate only depends on the previous’, accelerated methods like Alg. 1 linearly combine the previous two iteration’s estimates to improve convergence. As a consequence, backpropagating error from layer-to-layer requires contributions from two successive layers. Specifically, we compute the error at all $N$ layers with the recursive relation,

where each partial gradient constitutes a single step in Alg. 1 (fully derived in the supplement).

With the backpropagated error at each layer, we compute the gradient of the loss function with respect to $\mathbf{C}$ as,

Here, $\left(\partial\widehat{\boldsymbol{\phi}}^{(n)}\middle/\partial\mathbf{z}^{(n)}\right)$ backpropagates the error through the proximal operator and other partials with respect to $\mathbf{C}$ relate the backpropagated error at each layer to the changes in $\mathbf{C}$. Derivations of these partial gradients are included in the supplementary material. In Alg. 3, we unite these two steps to form a recursive algorithm which efficiently computes the analytic gradient for a single training example. Alternatively, general purpose auto-differentiation included in learning libraries (e.g. PyTorch, TensorFlow) can be used to perform the gradient updates.

IV Results

where $\tau=1\text{e}^{-3}$ is set to trade off the TV cost with the data consistency cost and $D_{i}$ is the first-order difference operator along the $i^{th}$ image dimension. We efficiently implement the proximal operator of Eq. 17 in closed form via parallel proximal method (details in supplement).

Traditional qDPC uses 4 measurements to adequately cover frequency space. Our learned designs are more efficient and may require fewer measurements; hence, we show learned designs for the cases of 4, 3 and 2 measurements. The designs and their corresponding phase WOTFs are shown in Fig. 3.

Comparing our learned designs with previous work, Fig. 4 shows the phase reconstruction for a single simulated test example using 4, 3 and 2 measurements. The ground truth phase is compared with the phase reconstructed using traditional qDPC designs , annular illumination designs , condition number optimized designs , A-optimal designs , and our physics-based learned designs. Table I reports the peak SNR (PSNR) statistics (mean and standard deviation) for the phase reconstructions from $\mathcal{R}$ evaluated on our set of testing examples. Our learned designs give significant improvement, recovering both the high and low frequencies more accurately.

IV-B Experimental Validation

To demonstrate that our learned designs generalize well in the experimental setting, we implement our method on an LED array microscope. A commercial Nikon TE300 microscope is equipped with a custom quasi-Dome illumination system (581 programmable RGB LEDs: $\lambda_{R}=625$ nm, $\lambda_{G}=532$ nm, $\lambda_{B}=450$ nm) and a PCO.edge 5.5 monochrome camera ($2560\times 2160$, $6.5\mu m$ pixel pitch, 16 bit). We image two samples: a USAF phase target (Benchmark Technologies) and fixed 3T3 mouse fibroblast cells (prepared as detailed in the supplement). In order to validate our method, we compare results against phase experimentally estimated via pupil-corrected Fourier Ptychography (FP) with equivalent resolution. FP is expected to have good accuracy, since it uses significantly more measurements (69 single-LED measurements) and a non-linear reconstruction process.

Using the USAF target, we compare phase reconstructions from FP with traditional qDPC and our learned design measurements (Fig. 5). Traditional qDPC reconstructions consistently under-estimate the phase values. However, phase reconstructions using our learned design measurements are similar to phase estimated with FP. As the number of measurements is reduced, the performance quality of the reconstruction using traditional qDPC degrades, while the reconstruction using the learned design remains accurate.

To demonstrate our method with live biological samples, we repeated the experiments with 3T3 mouse fibroblast cells. Figure 6 shows that phase reconstructions from traditional qDPC again consistently under-estimate phase values, while phase reconstructions using learned design measurements match the phase estimated with FP well.

V Discussion

Our proposed experimental design method efficiently learns the coded-illumination designs by incorporating both the system physics and the non-linear nature of iterative phase recovery. Learned designs with only 2 measurements can efficiently reconstruct phase with quality similar to Fourier Ptychography ($69$ measurements) and better than qDPC ($4$ measurements), giving an improvement in temporal resolution by a factor of 2$\times$ over traditional qDPC and far fewer than FP. Additionally, we demonstrate (Table I) that the performance of our designs on a set of testing examples is superior to previously-proposed coded-illumination designs. Visually, our learned design reconstructions closely resemble the ground truth phase, with both low-frequency and high-frequency information accurately recovered.

By parameterizing our learning problem with only a few weights per measurement, our method can efficiently learn an experimental design with a small simulated dataset. This enables fast training and reduces computing requirements significantly. Obtaining large experimental datasets for training may be difficult in microscopy, so it is important that our method can be trained on simulated data only. Experimental results in Sec. IV-B show similar quality to simulated results, with both using the designs learned from simulated data only.

Finally, phase recovery with the learned designs’ measurements are trained with a given number of reconstruction iterations (e.g. determined by a CPU budget). This makes our method particularly well-suited for real-time processing. qDPC can also be implemented in real-time, but limiting the compute time for the inverse problem (by restricting the number of iterations) limits convergence and causes low-frequency artifacts. Our learned designs incorporate the number of iterations (and hence processing time) into the design process, producing high-quality phase reconstructions within a reasonable compute time.

VI Outlook

Our method is general to the problem of experimental design. Similar to QPI, many fields (e.g. Magnetic resonance imaging (MRI), fluorescence microscopy) use physics-based non-linear iterative reconstruction techniques to achieve state-of-the-art performance. With the correct model parameterization and physically-relevant constraints, our method could be applied to learn optimal design for these applications (e.g. undersampling patterns for compressed sensing MRI , PSFs for fluorescence microscopy ).

Requirements for applying our method are simple: the reconstruction algorithm’s updates must be differentiable (e.g. gradient update and proximal update) so that analytic gradients of the learning loss can be computed with respect to the design parameters. Of practical importance, the proximal operator of the regularizer should be chosen so that it has a closed form. While this is not a strict requirement, if the operator itself requires an additional iterative optimization, error will have to be backpropagated through an excessive number of iterations. Here, we choose to penalize anisotropic TV, whose proximal operator can be approximated in closed form . Further, including an acceleration update improves the convergence of gradient-based reconstructions. As a result, the unrolled network can be constructed using fewer layers than its unaccelerated counterpart. This will reduce both computation time and training requirements.

VII Conclusion

We have presented a general framework for incorporating the non-linearities of regularized reconstruction and known system physics to learn optimal experimental design. Here, we have applied this method to learn coded-illumination source designs for quantitative phase recovery. Our coded-illumination designs can improve the temporal resolution of the acquisition and enable real-time processing, while maintaining high accuracy. We demonstrated here that our learned designs achieve high-quality reconstructions experimentally without the need for retraining.

Funding Information

This work was supported by STROBE: A National Science Foundation Science & Technology Center under Grant No. DMR 1548924 and by the Gordon and Betty Moore Foundation’s Data-Driven Discovery Initiative through Grant GBMF4562 to Laura Waller (UC Berkeley). Laura Waller is a Chan Zuckerberg Biohub investigator. Michael R. Kellman is additionally supported by the National Science Foundation’s Graduate Research Fellowship under Grant No. DGE 1106400. Emrah Bostan’s research is supported by the Swiss National Science Foundation (SNSF) under grant P2ELP2 172278.

Acknowledgment

The authors would like to thank Professor Michael Lustig for his guidance and advice.

References