Multispectral total-variation reconstruction applied to lens-free microscopy

: Lens-free microscopy multispectral acquisitions are processed with an inverse problem approach: a multispectral total variation criterion is deﬁned and minimized with the conjugate gradients method. Reconstruction results show that the method is eﬃcient to recover the phase image of densely packed cells.


Introduction
Lens-free microscopy [1] is a minimalist microscopy technique, which consists in acquiring diffraction patterns of a thin sample and computing (reconstructing) the optical properties of an object which diffraction pattern would match the measurements. By its simple design, the system is cheap, robust and incubator-proof. It can easily provide videos with a large field of view. Moreover, biological samples can be imaged without staining [2].
Lens-free acquisition is however problematic. While a complete description of light requires complex value fields, only their absolute values (intensities) are measured with a lens-free setup featuring a basic sensor. The consequence of this lack of information is to create artifacts in reconstructed images known as twin images [3].
A first strategy to reduce such artifacts consists in using priors on the sample, such as non negativity [4], or priors related to sparsity [3,5]. Another strategy consists in enriching the measurements set by acquiring diffraction images at multiple object-to-detector distances [6,7], or at multiple wavelengths [2,8].
This article relies on both strategies, i.e. the use of an enriched measurements set (multispectral acquisitions) associated with physically valid priors of the unknown sample. Hence, we developed an algorithm that minimizes a multispectral total variation (MS-TV) criterion. The algorithm uses the MS-TV criterion to reconstruct phase images that show shape similarities along the different acquisition wavelengths.
The paper is organized in three parts. The method part describes the models used (i.e. the optical field, its propagation and its measurement), the multispectral total variation criterion, the computation of its gradient and an overview of the reconstruction process. The material part describes the 3-chromatic RGB acquisition system, which sequentially acquires three diffraction images of biological samples and the results part shows an example of a reconstruction obtained from a culture of A549 mammalian cells. The reconstructed time-lapse movie of this cell culture is provided (field of view of 29.4 mm 2 , 3.7 days) in supplementary material.

Method
Light is described as a complex scalar field and its propagation from one plane to another is obtained using a convolution kernel (h). Therefore, if A 0 is the complex field after the sample, the field A Z in the detector plane at the distance Z from the object is: where '*' is the convolution operator. Here, we used a kernel h Z derived from Fresnel theory, but other propagation models such as the angular spectrum [9] approach could be employed with no difficulties. In the following, Z is assumed to be known. In practice, the focusing distance has to be estimated with a precision of approximately 10 µm and it can be evaluated by inspecting Z-stack of back-propagated measurements or with an automated method as described in [10].
The propagation kernel is h λ Z (ì r) = (1/iλZ) exp (iπì r 2 /λZ), i being the imaginary unit. Note that the quantities A Z , A 0 , h Z , as well as the measured intensity I Z , the phase in the detector plane ϕ Z and the objective function gradient ∇ε (that are introduced in the following), are functions of the wavelength λ and the spatial position ì r in a plane. Hence, for example, the convolution of Eq. (1) can be explicitly written as: where ì r is the coordinate in the sample plane and ì r is the one in the detector plane. Comments about spatial discretization of Eq. (1) and Eq. (2) and computation using Fourier transforms can be found in Appendix A.
The quantity measured by a standard detector is the intensity of the light field: I Z = | A Z | 2 . At this stage, we see that the modulus of A Z is known but its phase is undetermined. Whatever the phase map ϕ Z , A Z = √ I Z . exp (i.ϕ Z ) and the corresponding field A 0 at the sample plane, would perfectly match the data. The approach presented here consists in finding the phase field ϕ Z so that A 0 (ϕ Z ) minimizes the MS-TV criterion (ϕ Z ): where ε is a small valued coefficient (10 −4 ) used to make Eq. (3) differentiable and prevent division by 0 in the gradient derivation presented afterwards. MS-TV is a L 1 -norm applied to the 2N λ fields (∂ A λ 0 /∂ x and ∂ A λ 0 /∂ y) where N λ is the number of wavelengths used (here N λ = 3). A comment about discretization of operators ∂/∂ x and ∂/∂ y can be found in Appendix A. As illustrated in Fig. 1, MS-TV criterion is lower when high values of these 2N λ fields are colocalized. With this consideration plus the fact that L 1 -norm tends to promote sparsity, we can see that minimization of the MS-TV criterion leads to results with similar sharp edges at all wavelengths. Contrary to [2,8] which present another way to data process multi-spectral acquisitions, our MS-TV criterion is not designed to produce sample fields A 0 with the same values for different wavelengths (see Appendix D) but favors colocalization of edges across wavelengths.
The optimization is performed using the non-linear conjugate gradients iterative method as described by the following pseudocode: .D (k) 5/ test of convergence, if not reached go back to 2/ For step 3/, the gradient ∇ (k) is analytically computed as detailed in Appendix C. The advancement direction D (k) is obtained by conjugating the gradient ∇ (k) with the previous direction of advancement D (k−1) , using Hestenes-Stiefel formula [11]. The advancement step σ (k) (a scalar) is obtained by developing the MS-TV criterion at order 2 in the advancement direction (Hessian computation) and by minimizing a second-order polynomial of variable σ (k) . For the test of step 5/, we simply stop the algorithm when the number of iterations reaches 100.

Material
Measurements were carried out with a Cytonote (Iprasense) lens-free video microscope ( Fig.  2(b)) which geometry is illustrated in Fig. 2(a). Workable diffraction measurements must exhibit fringes. Therefore, the optical system is designed to provide partial coherence. The source is a 3-chromatic RGB-LED spatially filtered by a pinhole illuminating a thin sample on a microscope slide at 50 mm of the source. A CMOS detector of 6.4 × 4.6 mm 2 with 3840 × 2748 pixels (pixel pitch of 1.67 µm) is used to measure diffraction patterns at a distance of ≈ 1 mm from the sample. In the configuration used, optical magnification is close to 1. For the aim of observing living cells during a long period of time without significant disturbance, 1 frame per 5 minutes is the maximum sampling rate achievable on this setup. Indeed, due to its short cell-to-sensor distance, the heat of the detector generated during a frame acquisition is partially transferred to the biological sample. A faster acquisition rate would lead to the overheating of the sample and may cause cell death. The system spatial resolution limitations are discussed in Appendix B. We show that the limiting physical factor in this system is the temporal coherence of the source.

Results
To illustrate the capabilities of the system, a video of the diffraction of a biological sample (A549 cells observed during 3.7 days, 1065 frames with a sampling of 1 frame per 5 minutes) was acquired. A typical diffraction measurement is presented on Fig.3(d) and Fig.3(g). Since such samples are transparent (non-absorbing), most of the contrast is expected to be due to phase variations.
An optimal distance Z (1340 µm) was first determined by visual inspection of retropropagation (A 0 = √ I Z * h −Z ) of measurements for various Z. This distance is used for the entire set of 1065 frames. We then performed the reconstruction with the proposed method (100 iterations) using either tri-chromatic acquisition, or only one acquisition channel (red). An example of reconstruction obtained with the proposed method from a tri-chromatic acquisition is presented on Fig.3(a), Fig. 3(b) and Fig.3(e). The reconstruction time of one frame is about 40 seconds with a script running on Matlab R2017b and a desktop PC equipped with a Nvidia Titan XP graphic card. The reconstruction results for all the 1065 frames of the acquisition set are presented as a movie in supplementary material (see Visualization 1 ). For comparison, the reconstruction obtained with the proposed method but using only the red channel acquisition, is presented on Fig. 3(c) and Fig. 3(f) (18 seconds runtime).
We note that the algorithm is successful in retrieving sharp results with the lens-free setup even when fringes of various cells are intermixed ( Fig.3(g)). Fig. 3(c) and Fig.3(f) show that the proposed algorithm used with only one wavelength is efficient in recovering a valid phase image of the sample even in the difficult case of densely packed biological cells. Furthermore, the detailed image of reconstructed phase (Fig.3(e)) shows the benefit of using wavelengths diversity since results from tri-chromatic acquisitions are less impaired by twin image residues (see arrows in Fig.3

(e) and (f), as well as the flatter area void of cells).
In Appendix D, reconstructions performed from the same set of acquisition data but using the standard multispectral method as described in [8] are presented for comparison purpose. We show that the two methods provide results with the same spatial resolution. However the MS-TV algorithm allows to reduce twin-image residuals and artifacts on this example.

Conclusion
In this work, we developed an algorithm for processing multispectral lens-free microscopy acquisitions. This technique is well adapted for visualizing transparent samples such as unstained cell cultures since it relies on acquisitions of diffraction images. The algorithm is based on the definition of a novel multispectral total variation criterion. Its minimization is performed through a standard non-linear conjugated gradient method, which uses analytical computation of the criterion gradient also presented in the paper.
Acquisitions from a densely packed culture of A549 mammalian cells are presented. In such a case, algorithms based on the determination of the support of the object would have failed since the later occupies most of the field of view. On the contrary, the proposed algorithm is able to reconstruct valid phase images of the sample over a wide field of view. Detailed inspection of reconstructions shows that the data-processing of multispectral acquisitions is advantageous compared to one channel acquisitions since it provides results with less artifacts.

Appendix A: Comments on discretization and the Fourier transform
For the sake of simplifying the text, fields are written as functions of the continuous spatial variable ì r, but, since data acquisitions are performed on discrete grids corresponding to the detector sampling (on the example, grid of size 3840 × 2748 with a pixel pitch of 1.67 m), actual computations are discrete operations.
Integrals as the one of Eq. (2) are to be read as a discrete sum; for example the integral of a field f is: i and j being indexes of the matrix f . For the same reason, partial derivatives, as appearing in Eq. (3), are discretized. For example the partial derivative with respect to x of a field f is: with S x 0 = 1 and S x 1 = −1.
For computational performance reasons, operations involving convolution with Fresnel kernel, as appearing in Eq. (1), are performed in the Fourier space ( using the Fast Fourier Transform) as described in Eq.(5) of reference [12]. We utilize the fact that the Fourier transformh Z of h Z is: where the Fourier transform definition used isf ( ì µ) = ∫ dì r f (ì r) exp(−2iπ ì µ.ì r). This way to compute Eq.1 is valid as long as the Nyquist theorem is not violated (see Eq. 51 of [13]); it is the case of our experimental conditions since Z < N∆x 2 /λ, where Z is the sample to detector distance (1340 µm), ∆x the detector pitch (1.67 µm), N the number of pixels in the smallest dimension of the detector (N = 2748) and λ is around 0.5 µm.

Appendix B: Comments on system resolution due to partial coherence
Coherent wave formulation is used throughout this article, although the system source of light are RGB LEDs, which are only partially coherent. This lack of coherence leads to a decrease of the spatial resolution of the system and may result in specific artifacts due to modeling errors. The three main factors that limit the spatial resolution are the detector pixel pitch and the source temporal and spatial coherences as discussed in the article [14]. Using Eq. (4), Eq.(5), Eq.(6) and Eq.(7) of this reference, we compute the half pitch resolution ∆x of these three limiting factors; we get ∆x = 1.67 µm due to the detector pixel pitch, ∆x = 1.12 µm due to source spatial coherence and ∆x = 2.40 µm due to source temporal coherence. (Results obtained with the values: λ = 0.5 µm, source-sample distance=50 mm, sample-detector distance=1340 µm, source size=50 µm, source spectral width=20 nm). Therefore, the limiting factor of our system is the temporal coherence of the source and we expect a complete loss of contrast for lines with a width around 2.40 µm on phase test targets.
Such a loss of contrast is observed on a standard USAF1951 phase test target (SiO 2 material, 320 nm-height pillars) as shown in Fig.4. The algorithm described in the present article was applied to a tri-chromatic acquisition of this test target; group 6 is perfectly well contrasted whereas lines of group 7 are gradually less visible with a complete loss of contrast for element 4 of group 7. As experimentally evaluated on this phase test target, the system half pitch resolution is around 2.8 µm, which is good agreement with the theoretical value.
The main artifact of the reconstruction is the non homogeneity of the recovered phase values visible on large structures (such as the biggest square in Fig.4), which is caused by the fact that intensity measurements are not sensitive to small spatial frequencies of phase samples.

Appendix C: Derivation of the gradient of the problem
The goal of this annex is to give mathematical details on how the gradient of the objective function is calculated. As explained in the part 2, for any detector phase ϕ Z , a multispectral TV criterion (ϕ Z ) of the reconstructed object field A 0 can be computed. The operations chain linking ϕ Z and are the following Eq.7-10: A Z = I Z .exp(i.ϕ Z ) (detector light field) (8) To obtain the gradient of , we use the differential notion. If we have a differentiable function f : A → B, for every a and h of A, there exists a linear function D a f called differential of f around a, such that: From now on, we adopt usual notations in physics to describe the differential; we write the above formula in the following syntax, where a is implied and where the "progress" h (noted δa) of the variable and the differential are noted with the symbol δ: Note that this syntax can be used for a scalar a or a field such as ϕ Z . Since is a scalar which depends on ϕ Z , a field of indexes λ and spatial coordinates ì r, the differential of is a linear form that can be written as: ∇ is called the gradient and has the same dimension as ϕ Z , the variable of . To derive its analytical formula, let us differentiate Eq.7-10 to obtain an equation of the shape of Eq.13.