Towards Low-Photon Nanoscale Imaging: Holographic Phase Retrieval via Maximum Likelihood Optimization

A new algorithmic framework is presented for holographic phase retrieval via maximum likelihood optimization, which allows for practical and robust image reconstruction. This framework is especially well-suited for holographic coherent diffraction imaging in the low-photon regime, where data is highly corrupted by Poisson shot noise. Thus, this methodology provides a viable solution towards the advent of low-photon nanoscale imaging, which is a fundamental challenge facing the current state of imaging technologies. Practical optimization algorithms are derived and implemented, and extensive numerical simulations demonstrate significantly improved image reconstruction versus the leading algorithms currently in use. Further experiments compare the performance of popular holographic reference geometries to determine the optimal combined physical setup and algorithm pipeline for practical implementation. Additional features of these methods are also demonstrated, which allow for fewer experimental constraints.


Holographic CDI and Phase Retrieval
Coherent Diffraction Imaging (CDI) is a scientific imaging technique used for resolving nanoscale scientific specimens, such as macroviruses, proteins, and crystals [11]. In CDI, a coherent radiation source, often being an X-ray, is incident upon a object, whereupon diffraction occurs. The resultant diffracted wave is then incident upon a far-field detector which measures the resulting photon flux. This photon flux is (disregarding noise) approximately proportional to the squared magnitude values of the Fourier transform of the electric field within the diffraction area. Thus, in principle, the object's structure (e.g., its electron density) can be determined by solving the phase retrieval problem, the mathematical inverse problem of recovering a signal from squared magnitude values of its Fourier transform.
A popular setup for performing CDI experiments, known as holographic CDI, is when the collective object undergoing diffraction physically consists of a both a specimen (with an unknown electron density), together with a "reference", which is a portion of the object where the electric field is a priori known (see Fig. 1). The inclusion of a reference in the CDI setup is well-known both to enhance the quality of image reconstruction, and to greatly simplify the analysis and solving of the corresponding phase retrieval problem [1,7].

Figure 1:
Holographic CDI schematic. The upper portion of the diffraction area contains the imaging specimen of interest, and the lower portion consists of a known reference shape. Image courtesy of [14].
A convenient setup and mathematical notation for a specimen-reference hybrid object which we shall employ is S = [X, 0, R], (1.1) where X, 0, R ∈ R n×n and S ∈ R (3n)×n . Here, X is the "unknown" specimen, 0 is an n × n zero matrix, and R is the reference, whose entries are known. An example of this setup is shown in Fig. 2.
Thus, the photon flux data is modeled as where Y ∈ R m 1 ×m 2 and S ∈ C m 1 ×m 2 is the size m 1 × m 2 discrete Fourier transform of S, given by for (t 1 , t 2 ) ∈ {0, . . . , n − 1} × {0, . . . , 3n − 1}, and the absolute value operator |·| is understood as acting pointwise. Let the set of m 1 × m 2 measured values be denoted by M, and suppose that this set is indexed so that the low-frequency content is centered around (0, 0). Eq. (1.2) can be expressed as Given the data Y and the values of B (or, equivalently, R), obtaining the values of X amounts to solving the holographic phase retrieval problem, which can be stated as: (1.4)

Beamstop Apparatus
In CDI experiments, the central portion of diffracted radiation is often blocked from reaching the detector array by a beamstop apparatus (see Fig. 3a). This is because the low-frequency content (i.e. the Fourier transform magnitudes) of the measured data is typically much larger in magnitude than that of the higher frequencies [1] (e.g. see Fig. 3b). Thus, the low-frequency data must typically be excluded so that the range of measured values does not exceed the dynamic range of the detector sensors [9]. In this case, the data acquired is more realistically modelled as where denotes the Hadamard product (i.e. pointwise multiplication) and B ∈ C m 1 ×m 2 is of the form B ij = 0, |i| < ω 1 and |j| < ω 2 1, otherwise for some cutoff frequencies ω 1 and ω 2 .

The Poisson Shot Noise Model for CDI
The well-known Poisson distribution is a discrete probability distribution, which depends on a parameter λ > 0. A discrete random variable with this distribution -denoted by Z ∼ Pois(λ)has as probability mass function given by As a consequence of the quantum dynamics of photon emission in a radiation source, the photon flux measured at a CDI detector follows a Poisson shot noise model. Specifically, let Y denote the average value of Y (averaged over the number of detector pixels), and N p be the average photon flux per pixel. The data measured at the detector, at each pixel location (i, j) ∈ M, is modeled as [1]: (a) Holographic CDI setup with beamstop apparatus. (Image courtesy of [9]).
(b) Typical acquired CDI data (which is the squared Fourier magnitude data corresponding to the setup shown in Fig. 2). Low-frequency content is the largest by several orders of magnitude.

Figure 3
where A random variable Z ∼ Pois(λ) is well-known to have both its mean and variance equal to λ, from which it follows that Thus, as the average photon flux N p decreases, the corrupting effects of the Poisson shot noise increase.

Prior Art
Holographic phase retrieval dates back to the seminal work by Gabor [5], and has recently gained much attention due to the emerging technology of CDI. The most popular reconstruction algorithm for holographic phase retrieval to date is known as inverse filtering, which is a deconvolution procedure by which exact reconstruction of an image may be achieved given noiseless data [1] . While straightforward, this method does not incorporate any denoising. A popular variant on this algorithm, known as Wiener filtering, incorporates denoising into the inverse filtering algorithm [8]. However, this denoising is based on an additive noise model -an assumption which does not hold true for the Poisson shot noise encountered in CDI. Thus, both these methods are not well-suited and perform poorly given highly noise-corrupted (e.g. low-photon) CDI data.
The advent of a maximum likelihood approach for phase retrieval has previously been suggested for phase retrieval reconstruction in CDI in the non-holographic setting [15,2]. However, this method alone without the usage of a holographic reference is insufficient for high-quality imaging in the low-photon regime (see Fig. 13 and Fig. 14), and has not found practical application in CDI.

Our Contributions
This work proposes and studies the advent of maximum likelihood optimization together with the usage of a holographic reference as a new framework for low-photon CDI imaging. Algorithms are presented by which this optimization setup can be robustly implemented in practical CDI applications. Through extensive numerical experiments, we demonstrate the superior image reconstruction achieved by this method given CDI data in the low-photon regime. This surpasses the performance of the leading popular algorithms to date, while as well not requiring several of these algorithms' constraints. We further experiment with various reference objects to determine the optimal experimental setup. These results thus provide a viable pathway towards the advent of low-photon nanoscale imaging via CDI.

The HoloML Objective Function
Following the general strategy of maximum likelihood estimation, we seek to determine the image X which maximizes the probability of obtaining the measured data Y , given the probability distribution corresponding to the Poisson shot noise model. From Eq. (1.7), Using the standard assumption that measured pixel values are independent (e.g. see [17]) and substituting in Eq. (1.6), it then follows that the probability of obtaining the set of measured data Y as a function of X is given by Since the function log(·) is monotonically increasing, the global maximizer of g(X) is equal to the global minimizer of the corresponding negative log-likelihood function, i.e. of the function

The value of
Np Y does not significantly vary with X (see [1]), and thus this term can be viewed as being approximately constant and not affecting the global minimizer. We can thus formulate the holographic phase retrieval problem via the solution of the following maximum-likelihood optimization problem: We shall term this the HoloML objective function.

Optimization Methods
The optimization problem of Eq. (2.1) is nonconvex [3]. Likewise, the leading algorithms for classical (i.e. non-holographic) phase retrieval can be formulated as nonconvex optimization problems [18]. In the classical setting, these optimization problems are effectively optimized by alternating projection type schemes exclusively, and are known to not be amenable to first-order methods [12].
In contrast, Eq. (2.1) is observed to be effectively minimized via the usage of both first-order and alternating projection type algorithms. For the advent of first-order methods, the gradient function is given by where F † denotes the adjoint of F. An alternating projection scheme is also subsequently presented, being the popular ADMM optimization algorithm. Remarkably, as shown in Section 3, both these algorithm types produce similar solutions that are of high quality. This underscores that the effectiveness of the HoloML method is in the objective function, and is not tied to a particular optimization algorithm. In turn, this allows for flexible implementations that can be tailored to best suit particular CDI applications.

ADMM for HoloML Optimization
We derive the the Alternating Direction Method of Multipliers (ADMM) algorithm for optimizing the HoloML objective function of Eq. (2.1). This proceeds based on the key observation that this optimization problem can be equivalently formulated as 3) The augmented Lagrangian corresponding to Eq. (2.3) is where ·, · and · denote the Frobenius inner product and norm, respectively, and ρ > 0 is a fixed constant. The variable updates given by the ADMM algorithm are thus given by The explicit formulas for these minimizations are derived as follows.
Updating U : By completing the square, the subproblem for updating U can be reformulated as This problem can be decoupled into minimizing each of the n 1 × n 2 entries in the matrix U separately, i.e. for each (i, j) we seek to minimize the subproblem The first two terms depend only on the magnitude of U ij , and the last term implies that an optimizer must takes the phase which is minimized by Thus, altogether the update for U is given by

Updating X:
By similarly completing the square, it follows that the update of X is given by the minimizer of the subproblem This is given by

Performance
Numerical experiments were conducted comparing the results of optimization algorithms applied to minimizing the HoloML objective function versus the current leading holographic phase retrieval algorithms. These experiments were conducted using size 256 × 256 pixel test images of biophysical specimens -the mimivirus [6], embryo [4], oocytes [4], S. pistillata [16], salmonella [13], and sifA protein [13]. For the experiments in the following two subsections, the setup used is that given by Eq. (1.1) and shown in Fig. 2, where the test image, zero region, and reference object (being the URA reference) are each of size 256 × 256 pixels, altogether forming a hybrid object of size 256 × 768. Two-times oversampling was implemented when generating the corresponding Fourier transform magnitudes, which thus produced a data array of size 512 × 1536. This data was then subject to Poisson shot noise, whose average photon flux value is subsequently denoted as N p . For each setting, the resulting data was considered both without any beamstop occlusion and with a beamstop of size 25 × 25 pixels zeroing out the lowest-frequency data. Two algorithms were applied towards optimizing the HoloML objective function. The first such method is the conjugate gradient algorithm (a standard first-order optimization algorithm) which is implemented given the HoloML objective function of Eq. (2.1) and its corresponding gradient given by Eq. (2.2), and was implemented in Matlab using the Manopt optimization package. We term this algorithm HoloML-CG. The second method is an implementation of the ADMM algorithm derived in Section 2.2.1, and is termed HoloML-ADMM. We compare these new methods to the leading holographic phase retrieval algorithms in use to date, namely the inverse filtering [1] and Wiener filtering methods [8].
The relative error used for these experiments is given by where Y 0 denotes the experimental data, and Y = |F(X) + B| 2 is the data corresponding to the reconstructed image X. 1 This error metric is standard and naturally applicable for CDI data, for which only the measured data Y 0 is known.

Low-Photon Imaging Experiments
We consider the comparative performance of the HoloML algorithms on various biophysical specimen test images that are corrupted with Poisson shot noise with an average photon flux value of N p = 1. The images reconstructed from these simulations, with no beamstop occluding data, are shown in Fig. 4, and the corresponding relative error values are shown in Fig. 6. The analogous images and data for experiments given data occluded by a beamstop are shown in Fig. 5 and Fig. 7, respectively. In all these experiments, the HoloML algorithms clearly produce the best image reconstruction, as well as the smallest relative error.

Varying the Photon Flux
The comparative performance of these algorithms at varying average photon flux values, specifically for N p = 1000, 100, 10, 1, 0.1, is studied. The resulting reconstructed images and corresponding  Figure 4: Reconstruction of biophysical test images from low-photon data using various holographic phase retrieval algorithms, for data with no beamstop occlusion. The HoloML algorithms provide superior results.
relative errors are shown in Fig. 8 and Fig. 10, respectively, for reconstruction of the mimivirus image given data with no beamstop, and given data with beamstop occlusion in Fig. 9 and Fig. 11, respectively. It is observed that as the value of N p decreases, the image reconstruction quality is dramatically improved using the HoloML algorithms.

Reference Object Performance
Experiments were performed in which the reference object was varied (while using the otherwise same simulation parameters) to study the performance of the most popular CDI reference choices towards low-photon image reconstruction. Specifically, the reference choices implemented are the pinhole reference, block reference [1], and uniformly redundant array (URA) reference [10] as shown in Fig. 12, as well as no reference (i.e. a region consisting of all zero values). In these simulations, given various biophysical test images and for each of these reference choices, experiments were conducted for which data was subject to Poisson shot at N p = 1 and the HoloML-CG algorithm was applied. Fig. 13 and Fig. 14 show the results of these simulations, for data with and without beamstop occlusion, respectively. It is observed that the URA reference consistently produces the best image reconstruction, while the block reference performs the best from amongst references with  Figure 5: Reconstruction of biophysical test images from low-photon data using various holographic phase retrieval algorithms, with low-frequency data occluded by a beamstop. The HoloML algorithms provide superior results.  simpler geometries. This is consistent with prior theoretical and experimental results comparing reference performance [1], [10].

Breaking Classical Algorithm Barriers
Both the inverse filtering and Wiener filtering algorithms are hampered by two strict constraints. They require both a minimum oversampling ratio for the Fourier magnitude data, as well as a minimum separation distance between the specimen and reference objects [8,1]. Given a ground-truth image of size n 1 × n 2 (which includes both the specimen and reference portions) and corresponding Fourier magnitude data of size m 1 × m 2 , the oversampling ratio in the xand y-axis directions are defined as Often in CDI imaging setups these two ratios will coincide and both be referred to collectively as the oversampling ratio. For both the inverse filtering and Wiener filtering methods, there is a strict lower bound on the number of data measurements required for reconstruction. Specifically, it is necessary that m 1 ≥ 2n 1 − 1 and m 2 ≥ 2n 2 − 1 [8]. Thus, these methods essentially require an oversampling ratio of at least two.
In contrast, there is no minimal oversampling ratio necessary for implementing the HoloML methods. Numerical simulations demonstrate that these methods are capable of image recovery given data at a far smaller oversampling ratio, without loss of reconstruction quality, as shown in Fig. 15 and Fig. 16. Since a smaller oversampling ratio allows for fewer data measurements, these results imply that the HoloML algorithms allow for practical CDI imaging using far fewer CCD detectors (as each detector corresponds to a data measurement [9]). Inverse Filtering. Wiener Filtering HoloML-CG HoloML-ADMM For reconstruction via the inverse filtering and Wiener filtering algorithms, the specimen and reference setup must as well satisfy the holographic separation condition, which is precisely stated as follows. Given a specimen and reference both of size n × n, there must be a zero-separation of at least size n × n between them [8].
In contrast, the HoloML methods do not require any minimum separation between the specimen and reference objects. Numerical simulations demonstrate that for the HoloML methods, the quality of image reconstruction is not signficiantly impacted by the distance between the specimen and reference objects, including when this distance is zero. This is illustrated in Fig. 17 and Fig. 18, which shows the reconstruction of the mimivirus image using the HoloML-CG algorithms given a specimen and reference input which are separated by a distance d, and given data at various photon flux levels N p . (The case where d = n coincides with the previous experimental setups, and is shown in Fig. 2.) By allowing for a smaller separation distance between the specimen and reference objects, or none at all, the HoloML methods thus allow for more robust and flexible design of CDI experiments.

Discussion
A new algorithmic framework for holographic phase retrieval given data subject to Poisson shot noise is introduced via maximum likelihood optimization, and is termed HoloML. This method is implemented using practical optimization algorithms, namely conjugate gradient and the alternating direction method of multipliers (ADMM), on simulated holographic CDI data that is highly corrupted by shot noise. The results obtained demonstrate the far improved image reconstruction quality achieved by HoloML algorithms given noisy (i.e. low-photon) data. As well, it is shown that the implementation of a reference object in the experimental setup is necessary for obtaining superior image quality, and amongst popular reference object choices the URA reference achieves optimal results. Additional advantages of the HoloML methods -namely that they require less oversampling and no specimen-reference separation -further distinguish them from popular existing methods.
Based on these algorithmic and numerical results, an experimental realization of the HoloML method with a suitable reference object setup (e.g. a URA or block reference) thus has the potential  to produce superior imaging results, and with fewer experimental constraints. Future research in this direction would be of great interest towards impactful advances in the field of low-photon nanoscale imaging.   Figure 13: Reconstructed images from low-photon data, with no beamstop, acquired from setups with various reference objects. The best recovered image quality is provided by the URA reference, followed by the block reference.
spurred the exploration of statistical methods for holographic phase retrieval. D.B. is also very grateful to several colleagues at the Flatiron Institute for much helpful discussion throughout the  Figure 14: Reconstructed images from low-photon data, with low-frequency data occluded by a beamstop, acquired from setups with various reference objects. The best recovered image quality is provided by the URA reference, followed by the block reference.