Model-independent noise-robust extension of ptychography

A. P. Konijnenberg; W. M. J. Coene; H. P. Urbach

doi:10.1364/OE.26.005857

1. Introduction

Ptychography [1] has recently been proven to be a very promising technique to image objects ranging from biological samples [2] to semiconductor structures [3] in both the visible light regime [4,5], as well as in the extreme ultraviolet (EUV) [6], X-ray [3,7], and electron regimes [8]. In this method, multiple intensity patterns are recorded, from which the object is reconstructed computationally using an iterative algorithm. Since the object is recorded from a series of intensity patterns, it is important that these patterns are recorded with a sufficiently high signal-to-noise ratio (SNR).

For both regular [1] and Fourier ptychography [5], the effects of noise in the intensity measurements on the object reconstruction have been studied, and to make the reconstruction algorithm more noise robust several methods have been proposed, most of which rely on choosing the right cost function to minimize [9], sometimes combined with a regularization scheme [10]. One choice for the cost-function is the log-likelihood cost function which can be derived from a maximum-likelihood principle [10,11]. In this scheme, one needs to assume a certain noise model and maximize the likelihood of obtaining the measured intensity patterns by minimizing the negative logarithm of the likelihood function. A very good initial guess is required, or some regularization to prevent divergence is necessary [10,11].

An alternative method uses variance stabilizing transforms [12–14]. Also in this scheme a noise model needs to be assumed. In the noise model, the variance of the measured intensity depends on its expectation value, but by applying a variance stabilizing transform, the variance becomes approximately independent of the noise-free intensity value. It can be shown that the cost function using the variance stabilizing transform closely approximates the log-likelihood cost function [14]. However, the update function derived from the log-likelihood cost function can diverge, which would therefore require some additional regularization, whereas the cost function using a variance-stabilizing transformation for Poisson noise does not suffer from such instabilities. Another way to reliably increase noise-robustness is to simply reduce the step size of the update function [15].

We propose an extension to the ptychographic algorithm that can be applied to any of the aforementioned algorithms [10–14], does not require an assumed noise model, and works by adapting the intensity constraints. We test the algorithm with a proof-of-principle experiment in the visible light regime where we use a phase-only Spatial Light Modulator (SLM) to create and shift the object. The proposed algorithm is shown to compare favorably to regular ePIE [16] with reduced step size combined with a position correction scheme [17]. Note that the idea behind this algorithm is not necessarily restricted to ptychography, but could also be applied to other reconstruction schemes where a field is reconstructed from multiple intensity measurements, such as for example a through-focus scan or similar scans [18].

2. Noise robust ptychographic algorithms

Previous proposals for noise robust ptychographic algorithms mainly revolved around assuming a known noise model, and using a maximum-likelihood optimization scheme [10,11] or a variance-stabilizing scheme [12–14] to construct a cost function which, when minimized, yields a superior reconstruction, because prior knowledge of the noise model has been incorporated. In this section, we discuss the methods in more detail to see where there is room for improvement in these reconstruction schemes, which is what motivated the development of the algorithm that is proposed in this article.

2.1. Maximum Likelihood scheme

Let us start by describing the problem posed in ptychography when noise is present. We have an unknown complex-valued object O(x) that we want to estimate with O_est(x). Here, x is a two-dimensional position vector in object space. The object is illuminated with a known probe P(x) that is shifted to different positions X. We denote the noise-free diffraction pattern intensities by

m_{X} (u) = {| ℱ {O (x) P (x - X)} (u) |}^{2},

and the estimated diffraction patterns by

z_{X} (u) = {| ℱ {O_{est} (x) P (x - X)} (u) |}^{2} .

Here, ℱ{·} denotes the Fourier transform, and u denotes a two-dimensional position vector in diffraction space. We assume that u is discrete due to the pixelation of the detector. We denote the measured noisy diffraction patterns by y_X(u). The problem that needs to be solved is: how do we find an estimated object O_est(x) that matches the noise-free diffraction patterns m_X(u), while only noisy diffraction patterns y_X(u) are given?

The reasoning behind the maximum-likelihood optimization scheme is the following. We assume that y_X(u) obeys a certain noise model. For example, if we assume Poissonian noise, then the probability that we measure a certain intensity y given a noise-free intensity m is

P (y | m) = \frac{m^{y} e^{- m}}{y!} .

The total probability that we measure the data set y_X(u) is given by

P_{tot} [m_{X} (u)] = \prod_{X, u} P (y_{X} (u) | m_{X} (u)),

assuming that the noise at each pixel is uncorrelated. We then try to find O_est such that P_tot[z_X(u)] is maximized, or equivalently, such that − log P_tot[z_X(u)] is minimized. Thus, we define the cost function

\begin{array}{l} L [O_{est} (x)] & = - log P_{tot} [z_{X} (u)] \\ = \sum_{X, u} z_{X} (u) - y_{X} (u) log z_{X} (u) + log y_{X} (u)! . \end{array}

Here, the term

\sum_{X, u} log y_{X} (u)!

is independent of z_X(u), so it can be omitted. The ptychographic algorithm is then constructed by minimizing L[O_est(x)] using an optimization scheme, which is commonly the steepest-descent scheme. In the steepest-descent scheme, O_est(x) is updated using the Wirtinger derivative of the cost function

O_{est} (x) : = O_{est} (x) - μ \frac{d L}{d O_{est} {(x)}^{*}} .

The derivative of Eq. (5) may diverge. One way to remedy that is to assume that z_X(u) ≈ y_X(u) when we have a good object estimate (i.e. O_est(x) ≈ O(x)), so that we can Taylor expand Eq. (5) in terms of

\sqrt{z_{X} (u)}

around the point

\sqrt{z_{X} (u)} = \sqrt{y_{X} (u)}

, which gives [11]

z - y log z \approx y - y log y + 2 {(\sqrt{z} - \sqrt{y})}^{2} .

Thus, ignoring irrelevant additive and multiplicative constants, we can define the cost function

L [O_{est} (x)] = \sum_{X, u} {(\sqrt{z_{X} (u)} - \sqrt{y_{X} (u)})}^{2} .

This is the amplitude-based cost function that is used to derive regular PIE and ePIE, which supports the claim that regular ePIE is optimized to deal with Poisson noise.

2.2. Variance stabilization scheme

Another way to look at the problem is from the perspective of variance stabilization. For example, consider the intensity-based cost function

L [O_{est} (x)] = \sum_{X, u} {(z_{X} (u) - y_{X} (u))}^{2} .

If we again assume that y_X(u) is drawn from a Poisson distribution, then for larger y_X(u) we expect a larger error (m_X(u) − y_X(u))², and therefore the cost function weighs the pixels with high intensity more strongly than pixels with low intensity. If we want each pixel to have equal weight regardless of its measured intensity, we need to apply a transformation so that the expected error is independent of the measured intensity, which is a variance-stabilizing transformation. For Poisson noise, this transformation is taking the square-root, which leads once more to the amplitude-based cost functional of Eq. (8). It has been proposed to use other variance-stabilizing transforms, such as the Anscombe transform [14].

We have seen that the amplitude-based cost function in Eq. (8) can be derived either by using a Maximum-Likelihood approach that aims to find an object for which the probability of obtaining the measured intensity patterns is maximized, or by using a variance-stabilizing approach that aims to make the variance the same for each pixel regardless of its measured intensity. Even though these two approaches appear to be very different, it can be shown that the Maximum Likelihood approach is also a way to make the probability distribution the same for each pixel regardless of its measured intensity value. If P(y|m) denotes the probability we measure a noisy value y given a noise-free value m, then we can use Bayes’ rule to calculate the probability P(m|y) that m is the noise-free value if we measure a noisy value y [19]

P (m | y) = \frac{P (y | m) P (m)}{P (y)} .

We now want to find a transformation T_y such that T_y(m) is normally distributed with mean 0 and standard deviation 1, independently of the measured intensity value y

P (m | y) = \frac{e^{- T_{y} (m^{2}) / 2}}{\sqrt{2 π}} .

The cost function we then want to minimize is

L [O_{est} (x)] = \sum_{X, u} T_{y} {(z_{X} (u))}^{2} .

We can solve for T_y(m)²

\begin{array}{l} T_{y} {(m)}^{2} & = - 2 (log \sqrt{2 π} + log P (m | y)) \\ = - 2 (log \sqrt{2 π} + log P (y | m) + log (P (m) - log P (y)) . \end{array}

This expression for T_y can be plugged into Eq. (12). If we assume we have no prior information about m (so P(m) is independent of m), and we ignore all irrelevant additive constants that are independent of z, and we ignore global multiplicative constants, we find the cost function

L [O_{est} (x)] = - \sum_{X, u} log P (y_{X} (u) | z_{X} (u)),

which is the negative log-likelihood function as used in the maximum-likelihood method.

2.3. Discussion

We have seen that the maximum-likelihood schemes and variance-stabilizing schemes aim to construct a cost function that gives equal weight to pixels with low intensity and pixels with high intensity. It has been demonstrated that this is important in Fourier ptychography because the dark-field images which carry information about the high spatial frequencies of the object have a very low intensity compared to the bright-field images [9]. However, in regular ptychography such a distinction is absent (although there may be a large dynamic range in each single measurement, depending on what illuminating probe function one uses), and it has even been suggested to give pixels with a higher signal-to-noise ratio (which for Poisson noise are the pixels with higher intensity) more weight [20].

There is still some room for improvement in these schemes because all these cost functions push the algorithm towards the solution for which z_X(u) = y_X(u) which we know almost certainly to be wrong. After all, if y_X(u) is randomly distributed, the probability that y_X(u) = m_X(u) for all X and u is extremely small, so we expect that for the true solution O_est(x) = O(x) that L ≠ 0, even though the aforementioned algorithms try to achieve L = 0.

Another way to think about it is that in Eq. (13) we assumed that we have no prior information P(m_X(u)) about the probability that a certain noise-free intensity value is m_X(u). However, with all the overlap constraints in ptychography, there is in fact some information about this probability. For example, consider a probe position X and all its adjacent probe positions. If we measure certain intensity patterns for the adjacent probe positions, then the intensity pattern we measure for probe position X is not completely arbitrary, which means P(m_X(u)) is not independent of m_X(u). More specifically, one could create an object reconstruction using the intensity patterns for all probe positions except X, and then use that object estimate to give a rather accurate prediction of m_X(u) that can be used to update P(m_X(u)). By including this term in the cost function as prescribed by Eq. (13), the algorithm is not necessarily pushed to the solution z_X(u) = y_X(u) anymore. Note that P(m_X(u)) is independent of the assumed noise model.

In the extension that we propose in this article, we aim to find a cost function for which L = 0 can be achieved, making use of the overlap constraints that are present in ptychography.

3. The proposed extension

In the noise-free case, one would obtain O_est(x) by minimizing for each probe position X a cost functional such as

L_{X} [O_{est} (x)] = \sum_{u} {(\sqrt{z_{X} (u)} - \sqrt{m_{X} (u)})}^{2} .

In this ideal case, ePIE would be able to obtain a good reconstruction of the object. However, since we only have noisy measurements y_X(u), we cannot define this cost functional. What we can do, however, is to estimate m_X(u) from y_X(u) and all the additional information we have, such as the probe and the probe positions. So by running the ePIE algorithm we would be using all our available information to estimate m_X(u) from y_X(u) using z_X(u). Inspired by the observation that ePIE becomes more noise robust when decreasing the step size [15], one can come up with the following scheme to update the estimate for m_X(u):

Choose the initial estimate m_X,est,0(u) of m_X(u) to be y_X(u) or some denoised version of y_X(u).
Run the PIE algorithm with a small step size for a certain number of iterations, using m_X,est,k(u) as the intensity constraints. We denote the resulting estimated diffraction patterns as z_X,k(u)
Set $m_{X, est, k + 1} (u) = μ z_{X, k} (u) + (1 - μ) m_{X, est, k} (u) .$ Here, μ is a step size that should be chosen much smaller than 1 but larger than 0.
Repeat steps 2 and 3 until z_X,k(u) ≈ m_X,est,k(u).

Adapting the intensity constraints has also been considered when combining the Hybrid Input-Output (HIO) algorithm [21] with PIE [22]. Note that if we choose μ = 0 this scheme is identical to the regular ePIE algorithm, and in case the intensity measurements are corrupted by noise, the value of the cost functional L_X[O_est(x)] can in general not be minimized to be 0. By contrast, if we choose μ > 0 (but still small), then the algorithm stagnates when z_X,k(u) ≈ m_X,est,k (u) which means that L_X[O_est(x)] ≈ 0. However, it is important to realize that if we obtain L = 0, the error of the reconstructed object need not be 0 as well.

4. Method

To create and shift a phase object on which to test our proposed algorithm experimentally, we use a phase-only Spatial Light Modulator (SLM) to which we assign an image that we want to reconstruct. The SLM in which the object is created is a reflective liquid crystal phase-only PLUTO SLM by Holoeye, with a resolution of 1920 × 1080 pixels, and a pixel pitch of 8.0μm. With a lens with a focal length of 15cm we create the Fourier transform of the field that is reflected by the SLM. To reduce the dynamic range of the diffraction patterns [4,23], a fixed rapidly varying phase pattern as shown in Fig. 1(c) is added with the SLM on top of the shifted object, which more or less defines the probe P(x). The illuminated area of the SLM (which corresponds to the size of the probe) is a circle with a radius of 250 pixels. The object is shifted along a 7× 7 square grid with a period of 50 pixels with some random offsets to reduce the raster grid pathology in the reconstruction [24]. The images are recorded with an 8-bit SVS-VISTEK eco204MVGE CCD camera with a resolution of 1024 × 768 pixels and a pixel size of 4.65μm × 4.65μm. For each intensity measurement, we take the average of 50 pictures. The object is then reconstructed using ePIE [16] and a probe position correction scheme [17]. Another reconstruction is then performed which uses, in addition to ePIE and the probe position correction scheme, the proposed noise-robust scheme where the intensity constraints are updated as in (16). The experimental setup is shown in Fig. 2.

Fig. 1 The phase objects and the illuminating probe used for the ptychography experiment.

Abstract

1. Introduction

2. Noise robust ptychographic algorithms

2.1. Maximum Likelihood scheme

2.2. Variance stabilization scheme

2.3. Discussion

3. The proposed extension

4. Method

5. Results

6. Simulations

7. Phase retrieval using focus variation

8. Conclusion

References and links

Cited By

Figures (15)

Equations (19)

Optics Express