Low-Photon Counts Coherent Modulation Imaging via Generalized Alternating Projection Algorithm

Abstract

Phase contrast imaging is advantageous for mitigating radiation damage to samples, such as biological specimens. For imaging at nanometer or atomic resolution, the required flux on samples increases dramatically and can easily exceed the sample damage threshold. Coherent modulation imaging (CMI) can provide quantitative absorption and phase images of samples at diffraction-limited resolution with fast convergence. When used for radiation-sensitive samples, CMI experiments need to be conducted under low illumination flux for high resolution. Here, an algorithmic framework is proposed for CMI involving generalized alternating projection and total variation constraint. A five-to-ten-fold lower photon requirement can be achieved for near-field or far-field experiment dataset. The work would make CMI more applicable to the dynamics study of radiation-sensitive samples.

1. Introduction

The past 20 years saw the rapid advance of coherent diffraction imaging (CDI) techniques [1,2,3]. CDI disposes of the imaging lens and its theorical resolution is only determined by the wavelength of illumination radiation and the maximum angles that the diffracted light can be reliably detect. CDI is suitable for any coherent radiations, including synchrotron X-rays [4,5,6] or pulsed free electron lasers X-rays [7] and cold emission electron beams [8]. Various variants of algorithms, mostly derived from the Gerchberg–Saxton (GS) and the Fienup algorithm [9,10,11,12,13], were proposed. The oversampling smoothness algorithm in [9] properly introduces a spatial frequency filter based on the hybrid input–output (HIO) [3]. This algorithm makes full use of the smooth constraint of the region outside the support in the iterative process, and is proven to be robust when there is Poisson noise in the diffraction pattern. Single-shot CDI, however, requires isolated samples, and the sample phase variation is also restricted [14,15].

Ptychography ameliorates the drawbacks of CDI by recording multiple diffraction patterns with illumination overlapping. Ptychography became a main imaging modality at synchrotron facilities worldwide [15,16,17,18,19,20,21,22,23,24], and also gained popularity in the electron microscopy community [24,25,26,27]. In the presence of the noise of light and detectors, some optimization frameworks are proposed for ptychography and an order of magnitude of reconstruction quality is achieved by combining the conjugate gradients with preconditioning strategies, as well as regularization accounting for the statistical noise model of real experiments [28,29].

For imaging of biological samples, the ultimate limit of resolution is set by the radiation damage, and thus the performance of the imaging method under low illumination flux is essential [30,31,32,33]. Previous researchers also demonstrated that electron ptychography is a promising phase-contrast imaging technique with the advantage of high dose efficiency [30,34,35,36]. According to the dose factorization theorem [37], there is a strong hope that ptychography can keep its high performance under low light conditions. Using two-dimensional MoS2 crystal samples, it was demonstrated that highly quantitative phase can be obtained under low electron dose conditions [34].

The requirement of multiple measurements in ptychography makes it unsuitable for fast sample dynamics investigation. CMI is an alternative solution in addressing the weakness of CDI, but in a single-shot manner [38,39]. In CMI, a modulator is placed between the object and the detector to strengthen the coding of the phase of diffracted waves in a diffraction pattern. CMI is shown to work well under illumination flux greater than ${10}^7$ photons. Nevertheless, collected diffraction patterns can be severely interfered by Poisson noise as the illumination weakens, which affects the stability and accuracy of the algorithm. Therefore, it is worthy of retrieving phase information from contaminative diffraction patterns in biological specimens with efficient algorithm and ideal illumination photons. For lower flux, the current algorithm starts to encounter convergence difficulty and image artifacts in reconstruction. Recently, the PhENN neural network was proposed to improve the low-flux performance of CDI [40], as well as CMI [41]. Their training requirement can be computationally expensive and still there is concern about their suitability on general samples. Untrained learning methods also exist [42], but they are yet to be demonstrated for the low photon scenario.

In this manuscript, we report a computationally efficient and robust algorithm, based on generalized alternating projections (GAP), that preserves the advantages of CMI to recover the phase information from low signal-to-noise ratios data. GAP consists of two sub-problems, including measurement formation and statistical prior regularization. The paper is organized as follows: Section 2 details the theory and metric of the proposed method. In Section 3, the results and quantitative analysis of the simulation are presented. Experimental results of visible light and X-rays are described in Section 4 before a summary of our method is drawn in Section 5.

2. Methods

Generally, phase retrieval is a non-convex optimization problem, due to the intensity measurement of wavefields. Unlike convex situations, non-convex optimization problems are at risk of stagnation at local minima and saddle points [43]. Current phase retrieval algorithms usually introduce additional regularizing priors to the optimization function to reduce the scope of the solution space.

In the following, we briefly sketch the basics of the proposed CMI-GAP algorithm. As is common in machine learning [44], dictionary learning [45], compressive sensing [46], as well as sparse coding [47], the phase retrieval problem in CDI can be represented as an objective function

$$\underset{x}{\min }\frac{1}{2}{‖\left|Ax\right|-y‖}_2^2+R(x).$$

(1)

Here, $x$ is the complex object wavefield to be recovered. $y$ is the measured diffraction data. $A$ stands for the free space propagation operator. In the near-field, forward propagation is defined as

$$P_{near}=As\left\{M\right\{As(x\cdot s)\left\}\right\}.$$

(2)

$M$ is the transmission function of the modulator and $s$ is the support function. $As$ denotes the angular spectrum propagation method. Additionally, in far-field, forward propagation is defined as

$$P_{far}=F\left\{M\right\{As(x\cdot s)\left\}\right\},$$

(3)

where $F$ denotes the Fourier transform.

In Equation (1), ${‖\left|Ax\right|-y‖}_2^2$ is a smooth fidelity term, which can be minimized under a standard gradient descent method, while $R\left(x\right)$ is a non-differentiable regularization term. The proximal gradient algorithm can be adopted for such optimization problems; the so-called proximal operator is defined as

$$pro{x}_g(x,\lambda )=\underset{x}{\arg \min }\frac{1}{2}{‖\left|Ax\right|-y‖}_2^2+\lambda g(x).$$

(4)

The phase solution x is calculated by minimizing the combination of the fidelity term ${‖\left|Ax\right|-y‖}_2^2$ and regularization term $g\left(x\right)$. Constant $\lambda$ adjusts the weighting of the regularization term relative to the fidelity term.

Iterative phase retrieval algorithms for CDI involve projections between two spaces and imposition of constraints, which can generate satisfying results under certain limitations [48]. For CMI, the propagation operator $A$ also includes wavefront modulation. We observed that a forward operator $A$ with a complex matrix structure could enhance the encoding effectiveness of the sample phase and accelerate the retrieval process greatly. CMI shows good performance in X-ray and optical experiments when the measured data are not subject to heavy noise [38,39].

The GAP framework is one popular method that can retrieve high fidelity images from contaminated data by incorporating sample priors. Compared with other algorithms, such as ADMM [49] and TwIST [50], the GAP algorithm uses only Euclidean projection in updating the unknown $x$. PnP-GAP exhibits global convergence for data recorded in the scheme of snapshot compressive imaging [51,52].

In this work, we adopt the GAP strategy to suppress the noise and enhance the performance of CMI in a low lighting scenario. According to the GAP framework, Equation (4) can be reformulated as

$$(x,z)=\arg \min \frac{1}{2}{‖y-\left|Ax\right|‖}_2^2+\lambda g(z).s.t.x=z.$$

(5)

Here, a new variable $z$ is introduced into the objective function. The conditional constraint term $x=z$ is then replaced with a penalty term, given by

$$(x,z)=\arg \min \frac{1}{2}{‖y-\left|Ax\right|‖}_2^2+{‖x-z‖}_2^2+\lambda g(z)$$

(6)

where $x$ and $z$ are the variables that need to be optimized. The optimum may be effectively reached by updating $x$ and $z$ in an alternating fashion. In updating $x$, we are inclined to get more information from the diffraction pattern $y$ to form the wavefront; while in updating z, the resultant $x$ is conditioned to fit more accurately with the real samples. The specific details of implementation are as follows:

Updating $x$: assume $z$ is fixed, Equation (6) can be simplified to

$$x=\arg \min \frac{1}{2}{‖y-\left|Ax\right|‖}_2^2+{‖x-z‖}_2^2$$

(7)

and be solved as

$$x^{k+1}=z^k+A^T{\left(A{A}^T\right)}^{-1}\left(y-A{z}^k\right).$$

(8)

Updating $z$: assume $x$ is fixed and $g\left(x\right)$ the total variation (TV) regulation, Equation (6) can be rewritten as

$$z=\arg \min {‖x-z‖}_2^2+\lambda {‖z‖}_{TV}$$

(9)

$$z^{k+1}=\mathrm{TV}(x^{k+1}).$$

(10)

The TV regulation term features a good edge-preserving property, and Equation (9) forms a denoiser. The anisotropic TV defined as

$${\mathrm{TV}}_{anisotropic}(z)={‖D_{x}z‖}_1+{‖D_{y}z‖}_1={\displaystyle \sum _{i=1}^N\sqrt{{\left(D_{x}z\right)}_i^2}}+\sqrt{{\left(D_{y}z\right)}_i^2}$$

(11)

was used in this work, and the accelerated implementation method in [53] was adopted. The symbols $D_x$ and $D_y$ denote the first derivatives along the horizontal and vertical directions, respectively.

Figure 1 shows the overall flowchart of the CMI-GAP algorithm.

Both numerical simulations and optical experiments with similar parameters were carried out to evaluate the performance of the proposed algorithm. For both cases, the image quality was quantitively evaluated via root-mean-square (RMS) [54],

$$RMS=\sqrt{1-\frac{{\left|{\displaystyle \sum _{m,n}F\left. [\varphi \left(m,n\right)\right]\cdot F{\left. [\phi \left(m,n\right)\right]}^{*}}\right|}^2}{\left\{{\displaystyle \sum _{m,n}{\left|F\left. [\varphi \left(m,n\right)\right]\right|}^2\cdot {\displaystyle \sum _{m,n}{\left|F\left. [\phi \left(m,n\right)\right]\right|}^2}}\right\}}}$$

(12)

where $\varphi$ denotes the reconstructed wavefront, $\phi$ denotes the ground truth, $F\left. [\cdot \right]$ denotes the Fourier transform operator. The peak signal-to-noise ratio (PSNR) was also used.

3. Numerical Simulations

In this section, simulations with a different number of illumination photons and experiment settings are presented. The sample illumination was formed by a 1.5 mm pinhole irradiated by a collimated 632.8$\mathrm{nm}$ laser beam. Images of printed circuit boards were used for sample amplitudes and phases. The code was written with MATLAB and run on a computer equipped with 16 GB RAM and a GTX1060 6 GB GPU.

3.1. Near-Field Simulation

The convergency requirement of the CMI method determined other distance parameters; a large sample (pinhole) to modulator distance increases the effectiveness of the support constraint but at the cost of degraded image resolution in the reconstruction. A value between 20 and 40 mm was found to give the best balance. Therefore, the separations between the sample, modulator, and detector were set to $22.1$ mm and $21.5$ mm according to our experimental settings. Wave propagation between the three planes was calculated with the angular spectrum algorithm. The sampling intervals were 5.04 mm at all planes.

The performance of CMI-GAP was evaluated with the illumination photon number of $N$ being ${10}^8,{10}^7,{10}^6,{10}^5$. Reconstructed object amplitudes and phases, as well as the PSNR and RMS, are shown in Figure 2 and Figure 3, respectively. The iteration number of the traditional algorithm is about 500 and 50 with CMI-GAP. Overall, one can see the CMI-GAP algorithm shows better performance than CMI. When $N$ is above ${10}^7$, their reconstruction quality is comparable; as $N$ decreases, the RMS curve of CMI increases rapidly. When $N$ drops below ${10}^6$, the CMI-GAP algorithm can still retrieve details in most of the regions.

The convergence rate dependence of the CMI-GAP algorithm on the illumination photons was then studied. Figure 4 shows the comparison of the convergence rate for the photon counts of ${10}^6$, ${10}^7$, and ${10}^8$. Notably, the proposed algorithm converges rapidly in the first ten iterations when the number of photons is ${10}^7$. The RMS value reaches 0.08 before it stagnates, while the RMS using the CMI algorithm converges relatively slowly and reaches an RMS of 0.15 after 30 iterations. As $N$ decreases further, the performance difference between the two algorithms becomes more apparent.

3.2. Far-Field Simulation

We also validated the effectiveness of the CMI-GAP algorithm for the far-field geometry. In the simulations, a $632.8 \mathrm{nm}$ laser source was assumed. The sample to modulator distance was $30.8 \mathrm{mm}$. The modulator was placed on the front focal plane of a Fourier transform lens with a focal length of $30 \mathrm{mm}$. At the back focal plane of the lens, the detector was located.

Reconstructed amplitudes and phases are shown in Figure 5. Figure 6 shows the RMS and PSNR against the illumination photon numbers of $N$. The iteration number of the traditional algorithm is about 500 and 50 with CMI-GAP. A clear deviation in the image quality can be seen when $N$ is below ${10}^6,$ similar to the near-field case. In contrast, CMI-GAP produces significant enhancement in reconstruction quality, with as much as 20 dB and three times improvement in PSNR and RMS, respectively. The RMS bears a similar behavior as the near-field geometry. Figure 7 demonstrates the convergence curve when $N$ are ${10}^6$ and ${10}^7$. Generally speaking, the CMI-GAP algorithm has a relatively low RMS value and a faster convergence. The vibration of the convergence curve of the CMI algorithm was due to the algorithm switching, i.e., the HIO and DM algorithm are switched at specified iterations to promote convergence.

4. Verification with Experimental Data

Optical experiments were performed to verify the CMI-GAP algorithm in a low light flux scenario. A stabilized 632.8 $\mathrm{nm}$ He-Ne laser was used to form a quasi-monochromatic and collimated illumination beam. The sample was a positive USAF resolution target stuck on a $1.5 \mathrm{mm}$ hole. Scattered waves from the sample were modulated by a downstream thin modulator. An sCMOS camera (Thorlabs CS2100) with a pixel size of 5.04$\mathsf{\mu }\mathrm{m}$ was used to collect the diffraction data. A ptychography dataset was recorded with an illumination probe flux of ${10}^8$ photons. The retrieved modulator will be used in our CMI and CMI-GAP reconstruction.

In the near-field experiment, the distance between the sample and the modulator, and the distance between the modulator and the detector, were $22.1 \mathrm{mm}$ and $21.5 \mathrm{mm},$ respectively. CMI data with four levels of photon flux were acquired. The reconstruction results are shown in Figure 8. The iteration number of the traditional algorithm is about 500 and 50 with CMI-GAP. The overall degradation of quality is believed to be caused by residual reflections and detector noise. From Figure 8, one can see that the amplitude and phase images can be accurately reconstructed when $N$ equals to ${10}^8$ or ${10}^7$. In contrast, the reconstruction from the CMI-GAP algorithm has fewer artifacts and a higher resolution. The quantitative comparison, in terms of RMS, is given in

Table 1. RMS corresponding to the near-field CMI.

Number of Photons	${10}^8$	${10}^7$	$5\times {10}^6$	$3\times {10}^6$
RMS of CMI-GAP	0.17	0.31	0.45	0.62
RMS of CMI	0.32	0.50	0.76	0.85

. The ground truth images were obtained by ptychography.

When $N$ drops to $5\times {10}^6,$ the traditional CMI algorithm fails to produce a recognizable image, while the proposed algorithm can still reconstruct most of the sample information. Both algorithms failed when the average photon count per pixel dropped to 13.

For far-field geometry, we performed optical experiments using setup parameters the same as the simulations. A similar improvement as the near-field scenario was obtained. Here we show the results when applied to the far-field X-ray experiment data [39]. In these experiments, a $6.2 \mathrm{keV}$ X-ray was used to form a coherent beam. The illumination probe, carrying about $3.6\times {10}^6$ photons, was formed by a $2.2 \mathsf{\mu }\mathrm{m}$ pinhole. The modulator was placed at $2.5 \mathrm{mm}$ downstream of the object, while the distance between the modulator and the detector with a $172 \mathrm{mm}$ pixel size was about $7.2 \mathrm{m}$.

The photon flux was about $3.6\times {10}^6$ for this X-ray dataset. The reconstruction results obtained by CMI and CMI-GAP are shown in Figure 9. With the original CMI, the amplitude recovered is drowned by noise. In contrast, reconstructions with fewer artifacts and higher contrast can be seen clearly using the proposed CMI-GAP algorithm. For the X-ray experiment, the detector had a 172 mm pixel, leading to the insufficient sampling of the diffraction patterns. A detector pixel subdivision was then used to better model the detector sampling process. With a 2 × 2 subdivision, the quality is significantly improved for both CMI and CMI-GAP reconstructions, but the CMI-GAP images have a much clean and smooth appearance.

5. Conclusions

We proposed a reconstruction algorithm for low-photon count coherent modulation imaging. The algorithm uses the total variation constraint under the GAP framework. We compared the newly proposed CMI-GAP algorithm with the current algorithm for CMI using extensive simulations and experiments of various levels of illumination photons. The CMI-GAP algorithm still preserves the fast convergence speed and the high robustness against noise of CMI. Our results show the method outperforms the current CMI algorithm in convergence and image quality when the flux of illumination is lower than ${10}^7$. When the number of photons was reduced below ${10}^6$, CMI failed completely, while some sample features could still be retrieved with the CMI-GAP algorithm. CMI-GAP is robust to noise, with enhancement as much as 20 dB on PSNR and three times on RMS.

The CMI-GAP method could be further extended. First, the algorithmic parameters can be adaptively adjusted to the photon flux level, for which the Adam optimizer could be adopted. Second, the dataset statistics could be incorporated into the algorithm. Third, the continuity of diffraction patterns may be used as prior. It is necessary to mention that the convergence rate will decrease substantially when the computation array gets larger. The possibility of using the ADMM algorithm to accelerate the convergence is worth investigating.