Motion-resolved quantitative phase imaging

: The temporal resolution of quantitative phase imaging with Diﬀerential Phase Contrast (DPC) is limited by the requirement for multiple illumination-encoded measurements. This inhibits imaging of fast-moving samples. We present a computational approach to model and correct for non-rigid sample motion during the DPC acquisition in order to improve temporal resolution to that of a single-shot method and enable imaging of motion dynamics at the framerate of the sensor. Our method relies on the addition of a simultaneously-acquired color-multiplexed reference signal to enable non-rigid registration of measurements prior to phase retrieval. We show experimental results where we reduce motion blur from fast-moving live biological samples.


Introduction
Quantitative phase imaging (QPI) [1][2][3][4][5][6][7] enables stain-free and label-free imaging of transparent biological samples in vitro [8,9]. Unlike non-quantitative phase contrast techniques (e.g. Zernike Phase Contrast [10], Differential Interference Contrast (DIC) [11]), QPI methods are able to separate out the effects of phase and absorption. However, this generally comes at a cost of lost temporal or spatial resolution due to the need for multiple measurements. Here, we implement QPI without sacrificing speed or resolution, for the specific case of coded-illumination QPI.
Quantitative Differential Phase Contrast (DPC) [3,4,12,13] recovers the complex transmittance function of a sample from several coded-illumination measurements and a phase retrieval optimization. DPC achieves spatial resolution corresponding to twice the coherent diffraction limit and is practically implemented with an LED array based coded-illumination source on a commercial microscope ( Fig. 1(a)) [3,14]. Traditional DPC measurements consist of 4 intensity images, each captured with a half-circle illumination pattern at a different rotation angle ( Fig. 1(b)). The time-multiplexed nature of the measurements requires the implicit assumption that the sample is not moving during the acquisition. Of course, live biological samples may be non-stationary (defined as moving more than one pixel during the acquisition time). When only a single measurement is required (single-shot), the exposure time of the sensor can be scaled to guarantee approximately stationary behavior; however, for multi-shot methods, the acquisition time is limited by the sensor readout time. While each individual measurement may have an appropriate exposure time to guarantee the stationary assumption, motion occurring between measurements during the multi-shot DPC acquisition will cause errors in the reconstructed complex-field.
Interferometry-based QPI techniques such as digital holographic microscopy [5] and white light diffraction phase microscopy [6] can be single-shot, but are limited in spatial resolution by the coherent diffraction limit and are sensitive to system imperfections that cause speckle. Transport of intensity equation based QPI techniques [7] can be single-shot if the chromatic aberrations of the system are great enough [15] or with additional camera hardware [16]. DIC-based QPI techniques [17] can be single-shot with the addition of specialized hardware. Other methods rely on simultaneously acquiring multiple measurements via color multiplexing [18][19][20], polarization Traditional four-image DPC acquisition with rotating half-circle sources. Because the polystyrene bead is moving, the reconstructed phase suffers from motion blur artifacts. (c) Our method, mrDPC, uses traditional DPC source patterns in the green color channel and an additional constant navigator source pattern (half-circle) in the red channel. The motion-resolved phase reconstructed corrects the effects of the sample's non-rigid motion. multiplexing [21,22], or spatial multiplexing [23]. In the case of color multiplexing, the implicit assumption is made that the sample has no chromatic dispersion and is colorless; and in the case of polarization multiplexing the implicit assumption is made that the sample is not birefringent, both of which may be difficult to guarantee when imaging biological samples. Finally, in the case of spatial multiplexing, the space-bandwidth product of the reconstructed phase will be limited by the division of the sensor into smaller non-overlapping segments.
Here, we demonstrate that non-rigid sample motion occurring between the frames of a multishot DPC acquisition can be estimated and corrected. Techniques for rigid and non-rigid motion estimation and correction have been comprehensively applied in other fields (e.g. magnetic resonance imaging [24,25], multi-frame image enhancement [26], remote sensing [27], computer vision [28,29]), but not in QPI microscopy. It is not straightforward to apply these existing methods to DPC because they make an assumption that the spatial frequency content between any two images being registered is similar [30]. This assumption is violated when estimating the motion between DPC measurements, since each coded-illumination measurement has a unique spatial frequency contrast of the sample's optical phase. Thus, estimation of motion between raw DPC measurements will fail when using traditional registration techniques.
In order to perform motion estimation for DPC images, we introduce a new method, termed motion-resolved DPC (mrDPC), that uses an additional simultaneously-acquired color-multiplexed measurement with a constant coded-illumination pattern ( Fig. 1(c)). This navigator measurement uses one color channel of the source LEDs to display a constant illumination pattern (a half-circle), thus maintaining constant spatial frequency contrast. A color camera then separates the navigator and DPC measurements (without any assumptions regarding the dispersion or color of the sample). Non-rigid motion can then be estimated from the navigator measurements and corrected in the DPC measurements prior to phase retrieval. In this way, quantitative phase images can be recovered for each time point of the captured data, resulting in temporal resolution equivalent to single-shot methods. We demonstrate proof of principle experimental results in which blurring due to live sample motion (Amoeba proteus and Caenorhabditis elegans) is reduced.

Methods
Our proposed method, mrDPC, captures each DPC measurement sequentially, while simultaneously capturing a color-multiplexed navigator, using an LED array microscope. We program the green color channel of the LED array with rotating half-circle patterns (for DPC) and the red color channel with a constant half-circle pattern (for the navigator), illustrated in Fig. 2(a). The signal is measured on a color camera and separated via demosaiking with a precalibrated spectral sensitivity matrix (see App. color multiplexing) into DPC measurements and navigator measurements ( Fig. 2(a)).
To achieve motion correction between the four DPC measurements, we need to register each image to the others. Because the different DPC illumination patterns result in different contrast, they cannot be directly registered to each other. The navigator measurements circumvent this problem; they can be registered to each other and the resulting motion estimates can then be applied to the DPC measurements. Specifically, we estimate the motion between three pairs of measurements (between t 0 and t 2 , t 1 and t 2 , and t 3 and t 2 ) as outlined in Sec. 2.1. The reference time point can be any of the four measurement time points; we chose T = t 2 . The motion estimates are plotted as vector fields in Fig. 2(b), where the arrows' magnitude corresponds to the amount of the sample's local displacement and the arrows' orientation corresponds to the direction of the sample's local displacement. The three motion estimates are applied to the raw DPC measurements at times t 0 , t 1 , t 3 , respectively, to register them to the reference measurement. The registration is performed by resampling with linear interpolation. Using the physics model in Sec. 2.2, we then linearly deconvolve the motion-corrected DPC measurements (Eq. 7) to recover the sample's absorption and optical phase (Fig. 2(c)).

Motion estimation
The task of removing motion artifacts can be formulated as a blind deconvolution problem [31,32] where the unblurred image and the blur kernel are jointly estimated; however, this does not account for non-rigid motion. In this work, we use our navigator measurements to correct for the sample's non-rigid motion via image registration, enabling a wider array of biological applications.
To model non-rigid sample motion, our proposed method estimates a deformable mapping between pairs of images, for which there exist many algorithms [28,29,[33][34][35]. We chose the Symmetric Normalization (SyN) method [33,34] for its state-of-the-art performance [36] and open-source availability [37]. The method is called symmetric because it is commutative with respect to the ordering of the two input images and therefore does not over-fit the deformable mapping estimate to either image. This is particularly important to us, so that we can arbitrarily choose the reference time point without biasing our results.
The SyN algorithm solves an optimization problem (Eq. 1) to estimate the deformable mapping between two images, I 0 (r) and I 1 (r), such that a similarity metric, S(I 0 , I 1 ), is maximized and the deformable mapping is spatially smooth. The deformable mapping, g(r, t), is a function of space and time, where r denotes 2D spatial coordinates and t ∈ [0, 1] denotes a dimensionless time coordinate. At time t = 0, g(r, 0) maps I 0 (r) to itself, while at time t = 1, g(r, 1) maps I 0 (r) to I 1 (r). The SyN method achieves symmetry by jointly estimating a forward deformable mapping g 0 (r, t) between I 0 and I 1 and a backwards deformable mapping, g 1 (r, t), between I 1 and I 0 . Mathematically, the algorithm can be written as: max g 0 (r,t),g 1 (r,t) S(I 0 (g 0 (r, 0.5)), subject to, where our similarity metric is the normalized cross-correlation, defined as S(I a , I b ) = I a ,I b I a I b . The mappings' spatial smoothness is achieved by penalizing the term R(v(r, t)) = ∫ 0.5 t=0 Lv(r, t) 2 dt for each map. Here, L = ∇ 2 + I is the linear differential operator, where ∇ is the first-order difference, I is the identity operator and v(r, t) is the velocity field corresponding to g(r, t). This correspondence is enforced with the Lagrangian-Euler constraint [38] in Eq. 2. The SyN optimization is solved via gradient descent [33,34] and implemented in Dipy [37].

Phase retrieval
After using the navigator measurements to estimate the motion, we correct for motion in the DPC measurements, which can then be used as input to a phase retrieval algorithm. Generally, the relationship between an object's 2D complex transmittance function and measured intensity is non-linear, so recovery of phase requires non-linear optimization and an iterative solver. For in vitro biological samples, the "scatter-scatter" term is small and so we can make a weak object approximation, thus enabling phase recovery by simple linear deconvolution with weak object transfer functions (WOTFs) [3,4,7,39]. This linearization decouples the contributions of absorption and phase and allows us to express intensity in terms of linear contributions from: background, absorption and phase contrast. In the Fourier domain, where y is the intensity measurement and B is the DC term. Here,· denotes Fourier transform and u are 2D spatial frequency coordinates. H µ (u) is the WOTF for the sample's absorption and H φ (u) is the WOTF for the sample's phase. These terms are derived in [3], where P(u) is the complex pupil function, S(u) is the illumination source distribution, and denotes cross-correlation. In traditional DPC, the illumination sources are four rotating half-circles with radius N A ob j oriented right, bottom, left, and top (see Fig. 1(b)).
Here, λ µ and λ φ are regularization parameters, which are set to trade off data consistency and penalties for the low-frequencies in φ and the high-frequencies in µ. The necessity of the regularization comes from the phase WOTF's low-sensitivity to low-frequencies and absorption WOTF's low-sensitivity to high-frequencies. The optimization in Eq. 7 can be reformulated as a least-squares problem, to yield a closed form solution for the motion-resolved absorption,μ * , and motion-resolved quantitative phase images,φ * (· denotes the conjugate operator):
Validating our method's motion correction ability is challenging because of the difficulty in obtaining ground truth phase for comparison. To address this, we start by capturing full framerate videos (50×, NA=0.55) of a slow-moving sample Amoeba proteus (Carolina Biological Supply) and assume the traditional DPC reconstruction to be ground truth, since we expect negligible motion between frames. We then decimate the dataset in time by a factor of 8× to emulate faster motion, reconstruct with our method and traditional DPC, and compare results to the ground truth DPC result. As can be seen in Fig. 3(a), the bright water vacuoles and well-defined wall edges in the amoeba's nucleus and contractile vacuole (gold arrows in Fig. 3) appear blurred in the time-decimated traditional DPC reconstruction, but not in our method's reconstruction.
We next demonstrate our method with an even faster sample, live C. elegans (12.5×, NA=0.25), which generates significant non-rigid motion between the raw measurements ( Fig. 4(a)). From these measurements, we reconstruct and compare traditional DPC with our mrDPC method ( Fig. 4(b)). Insets highlight mrDPC's correction of distortion around the head region (pink insets in Fig. 4(c)) and blurring of internal features (blue insets in Fig. 4(c)).
By capturing a continuous video at the full framerate of the sensor, we can reveal biological motion dynamics of the C. Elegans. Reconstructions are performed on a sliding window of measurements such that each window has a full set of four DPC measurements and each window is offset by one measurement. The motion-resolved absorption and quantitative phase video reconstructions are compared with traditional DPC in Supplementary Material Visualization 1.
Finally, we compare our method to color-multiplexed DPC [18] (colorDPC), which is a single   shot QPI method (Fig. 5). ColorDPC is able to encode the information required for reconstruction into a single measurement using color multiplexing of the RGB LEDs and a color camera, under the assumption that the sample is non-dispersive (an assumption which is not required for mrDPC). Since colorDPC is a single-shot method, it will not suffer inter-frame motion blur, but uses fewer measurements than mrDPC and traditional DPC for each reconstruction and thus will have lower reconstruction SNR. In addition, the three color-encoded measurements for colorDPC will have different bandwidths, each defined by their respective encoding wavelength. As a result, the final reconstruction has orientation-varying high-frequency contrast, while traditional DPC and mrDPC do not (highlighted in Fig. 5).

Discussion
Our proposed method, mrDPC, can correct artifacts due to sample motion that is fast enough to cause motion blur across the four captured DPC images, but does not cause motion blur within each measurement. In the case of Amoeba proteus, the sample motion is slow enough over the duration of the multi-shot DPC acquisition that it can be assumed stationary and no motion correction is necessary. In the case of C. elegans, the stationarity assumption is violated, but each individual captured image is unblurred, so mrDPC helps to resolve the motion between frames.
Most generally, the sample is non-stationary in each single measurement and motion-induced blur is present. To address this, strobed illumination could be used to effectively shorten the capture time of each measurement to ensure stationarity. This strategy is analogous to our time-decimated validation in Sec. 3 (measurements acquired with delays between them). Design of the navigator pattern also affects performance. We make the assumption that the structural motion between the navigator measurements is the same motion that occurs between DPC measurements. This should hold true, since the navigator measurement is acquired simultaneously with its corresponding DPC measurement and its illumination is similar to the DPC illumination pattern in terms of bandwidth. This ensures that the highest resolution features' motion in the DPC measurements will be captured in the navigator measurements. Further, the navigator measurement must have sufficient SNR and gradient information [40] to perform motion estimation. While SNR isn't rigorously measured here, the power and number of LEDs in the navigator pattern is equal to that of the DPC pattern so that when spectrally unmixed neither contribute much additional noise to each others' measurements.
The work of Phillips et al. [18] discusses that only three half-circle coded-illumination measurements are required to perform the quantitative DPC reconstruction. We could incorporate this by only performing our method on three rather than four measurements; however, performance might degrade similar to in Fig. 5c. Our method is not limited to a fixed number of measurements; if more measurements give improved reconstruction results (e.g. increased SNR), we can incorporate them simply by registering additional measurements to the reference.
One limitation of the present method is that of using the weak object approximation, which only applies to samples with relatively weak phase and absorption. Since the motion correction is independent of the phase retrieval method, nonlinear methods can be used when the weak object approximation is violated. In that case, the linear reconstruction can serve as a good initialization to a non-linear phase retrieval optimization [4].

Conclusion
We present a computational method, motion-resolved DPC, that achieves similar reconstruction quality to that of traditional DPC's quantitative phase images, but corrects for the blurring caused by sample motion during the four-image acquisition. Validation of our method's navigator-based non-rigid motion estimation and correction of live Amoeba proteus sample motion is performed. Furthermore, we motion resolve even faster live C. elegans and reveal motion dynamics at the frame rate of the camera with video reconstruction.

Appendix: color multiplexing
Color sensors with a traditional Bayer filter [41] spatially multiplex color filters to capture spatial-spectral information. However, these color filters are not perfectly selective to a single spectra, but rather are sensitive to overlapping distributions of spectra. This cross-talk makes it necessary to calibrate the pixels' sensitivity relative to the spectrum of the illumination source, so that the desired spectral response can be demixed from the acquired measurements.
For our method, we estimate the spectral sensitivity of our sensor's red and green pixels to our illumination source's red (625nm) and green (532nm) LEDs. This is accomplished by spatially averaging each of the color channels' intensity response to red-only and green-only illumination [18]. The results form the entries in a matrix, C, that can be used by applying its pseudo inverse, C † , to unmix future color multiplexed measurements.
Here, we use the green pixels (I g1 , I g2 ) and green LEDs to encode the DPC signal, I dpc , and we use the red pixels (I r ) and red LEDs to encode the navigator signal, I nav .