PaperlayerView exact revision on arXiv

Source-linked AI summary

Experimental robustness of Fourier Ptychography phase retrieval algorithms

Li-Hao Yeh, Jonathan Dong, Jingshan Zhong, Lei Tian, Michael Chen, Gongguo Tang, Mahdi Soltanolkotabi, Laura Waller

arXiv:1511.02986v2physics.opticscs.CV

TL;DR

Fourier ptychography must recover high-resolution complex fields from many angle-dependent measurements despite noise and model mismatch. This paper compares cost functions and algorithms using theory, simulations, and experiments, finding that amplitude-based objectives are more robust and that Newton-type methods offer strong reconstruction performance with computational trade-offs.

  • Problem

    Existing Fourier ptychography algorithms require comparison under experimental imperfections, especially systematic model mismatch from aberrations and LED misalignment.

  • Method

    The paper develops a maximum-likelihood framework, compares cost functions and algorithms using simulations and experiments, and introduces a global Newton’s method.

  • Results

    Amplitude-based algorithms are more experimentally robust than intensity-based algorithms, while global Newton’s method gives the best reconstruction and sequential Gauss-Newton offers the best trade-off for large-scale applications.

  • Takeaways & Limitations

    Cost functions should account for intensity-dependent noise, and algorithm choice must balance robustness, reconstruction quality, and computational cost.

  • Takeaways & Limitations

    The FPM implementation of Wirtinger flow deviates from the original theory and is not guaranteed to converge to the global minimum.

Abstract

from arXiv · show

Fourier ptychography is a new computational microscopy technique that provides gigapixel-scale intensity and phase images with both wide field-of-view and high resolution. By capturing a stack of low-resolution images under different illumination angles, a nonlinear inverse algorithm can be used to computationally reconstruct the high-resolution complex field. Here, we compare and classify multiple proposed inverse algorithms in terms of experimental robustness. We find that the main sources of error are noise, aberrations and mis-calibration (i.e. model mis-match). Using simulations and experiments, we demonstrate that the choice of cost function plays a critical role, with amplitude-based cost functions performing better than intensity-based ones. The reason for this is that Fourier ptychography datasets consist of images from both brightfield and darkfield illumination, representing a large range of measured intensities. Both noise (e.g. Poisson noise) and model mis-match errors are shown to scale with intensity. Hence, algorithms that use an appropriate cost function will be more tolerant to both noise and model mis-match. Given these insights, we propose a global Newton's method algorithm which is robust and computationally efficient. Finally, we discuss the impact of procedures for algorithmic correction of aberrations and mis-calibration.

1 Introduction

Fourier ptychography reconstructs high-resolution complex fields from multiple low-resolution images, but experimental imperfections make algorithm robustness crucial. The paper compares algorithms and identifies cost function, noise, and model mismatch as key determinants of performance.

  • Fourier ptychography combines multiple illumination angles to achieve gigapixel-scale quantitative phase images with wide field-of-view and high resolution.
  • Experimental robustness is limited mainly by measurement noise, aberrations, and LED misalignment, which produce model mismatch.
  • Amplitude-based, intensity-based, and Poisson-likelihood algorithms show similar performance within each cost-function class, whereas Gerchberg-Saxton produces artifacts.
  • Second-order methods reduce low-frequency artifacts and often need fewer iterations, while sequential and global implementations show little performance difference for a given algorithm.
  • The cost function is a critical predictor of experimental performance because intensity-based methods suffer artifacts from noise and model mismatch, while amplitude-based methods are more robust.
  • Brightfield images can have several orders of magnitude higher intensity than darkfield images, producing substantially different noise levels.

2 Algorithm Formulation

The paper formulates Fourier ptychography through a forward model and maximum-likelihood noise models, then analyzes controlled experimental imperfections. Amplitude, intensity, and Poisson-likelihood costs encode different assumptions about how measurement errors depend on intensity.

  • 2.1 Forward problem for Fourier ptychography: Tilted plane-wave illumination shifts the object’s Fourier spectrum, and the objective lens passes the shifted spectrum through its pupil.
  • 2.2 Possible noise and simulated dataset: Controlled simulations isolate ideal data, Poisson noise, aberrations, and LED misalignment as separate sources of deviation from the forward model.
  • 2.3 Optimization problem based on different noise models: The intensity-based cost function treats pixel noise as identical and independent of measured intensity.
  • 2.3 Optimization problem based on different noise models: The amplitude-based cost function minimizes differences between estimated and measured amplitudes and underlies projection-based and Gerchberg-Saxton algorithms.
  • 2.3 Optimization problem based on different noise models: The white-Gaussian-noise formulation assumes equal variance across all pixels and images, an assumption that does not hold for Fourier ptychography datasets.
  • 2.3 Optimization problem based on different noise models: Poisson noise grows with intensity, so the Poisson-likelihood cost function assigns smaller weights to brighter, noisier pixels.
  • 2.4 Vectorization Notation: Vectorization raster-scans measured images into vectors and represents the higher-space-bandwidth object and cost functions in linear-algebra form.

3 Derivation of algorithms

The algorithm-derivation section compares existing Fourier ptychography solvers and introduces a second-order global Newton’s method under multiple cost functions, using a common initialization.

  • Existing algorithms and a new second-order global Newton’s method are compared by how they update the object estimate while optimizing different cost functions.
  • All algorithms use the amplitude of the on-axis LED image as their initialization.

3.1 First-order methods

The first-order methods include sequential amplitude-based updates and Wirtinger flow, with implementation choices determining their optimization behavior. Gerchberg-Saxton applies alternating projections, while the studied Wirtinger-flow implementation uses an intensity-based cost and non-guaranteed initialization.

  • Sequential amplitude-based method: The sequential amplitude-based algorithm optimizes one measurement at a time and is equivalent to the Gerchberg-Saxton alternating-projection approach.
  • Sequential amplitude-based method: The algorithm cycles through all measurements sequentially, updating the object after each measurement.
  • Sequential amplitude-based method: Each sub-iteration performs an amplitude replacement in the real domain followed by a Fourier-domain projection onto the updated region.
  • Wirtinger flow algorithm: Wirtinger flow uses gradient-descent-like updates with a specially chosen initialization and step size for theoretical convergence.
  • Wirtinger flow algorithm: The FPM implementation of Wirtinger flow uses an intensity-based cost function and differs from the original theory through additive-noise handling and low-resolution-image initialization.
  • Wirtinger flow algorithm: Because of these deviations, the studied FPM implementation is essentially gradient descent with a special step size and lacks a guarantee of convergence to the global minimum.
  • Wirtinger flow algorithm: Wirtinger flow can be implemented with amplitude-based, intensity-based, or Poisson-likelihood-based gradients for comparison.

3.2 Second-order methods

Second-order methods use curvature information to improve Fourier ptychography optimization, while the proposed global Newton method derives updates for multiple cost functions.

  • 3.2 Second-order methods: Second-order methods use both first and second derivatives to improve convergence speed and stability for nonlinear, non-convex problems.The Hessian supplies curvature information for more direct updates, but its computational cost can scale quadratically with image size.
  • 3.2.1 Sequential Gauss-Newton method [19]: The Gauss-Newton method approximates the Hessian using the Jacobian and regularizes its inversion by dropping off-diagonal terms and adding a constant regularizer.The regularized approximation is used because direct inversion is computationally expensive and may be ill-posed.
  • 3.2.1 Sequential Gauss-Newton method [19]: Applying the approximated Hessian inversion yields a Newton update for the object estimate, with α^(i) controlling the descent step.When P is fixed to binary values, the method reduces to sequential gradient descent with tunable regularization; simultaneous P updates make second-order optimization more important.
  • 3.2.2 New algorithm implementing a global Newton’s method: The proposed algorithm is a global second-order Newton method designed to compare amplitude-, intensity-, and Poisson-likelihood-based cost functions.The paper derives gradients and Hessians directly for the FPM problem rather than relying on approximations in the derivation.
  • 3.2.2 New algorithm implementing a global Newton’s method: The proposed Newton implementation solves Hessian inversions with conjugate gradients and selects step sizes using backtracking line search.The exact cost function and Hessian depend on whether the algorithm is amplitude- or intensity-based.

3.3 Convex-based methods

Convex-based methods lift phase retrieval into a higher-dimensional matrix problem, then use factorization and relaxations to make computation feasible, at the cost of losing the original guarantees.

  • 3.3 Convex-based methods: PhaseLift reformulates the non-convex object problem using the rank-1 matrix X = OO† and estimates X instead of O.The lifted formulation is intended to make the optimization convex and guarantee convergence to the global solution.
  • 3.3.1 PhaseLift algorithm [23–27]: The convex relaxation replaces rank minimization subject to I = A(X) with trace minimization, while the implementation uses intensity-based or related cost functions.Intensity- and Poisson-likelihood-based costs are convex on estimated intensity and can therefore be transformed into convex functions of X.
  • 3.3.1 PhaseLift algorithm [23–27]: The lifted matrix X has dimensions n^2 × n^2, making direct storage and eigenvalue decomposition impractical on a normal computer.The implementation therefore factorizes X as RR†, with R an n^2 × k matrix.
  • 3.3.1 PhaseLift algorithm [23–27]: The factorized problem is solved with an augmented Lagrangian multiplier, using tunable multiplier parameters to optimize R iteratively.The factorization avoids directly handling the oversized lifted matrix.
  • 3.3.1 PhaseLift algorithm [23–27]: After factorization, minimizing with respect to R makes the problem non-convex and removes the original provable guarantee.The paper identifies this loss of convexity as a consequence of the computationally practical reformulation.
  • 3.3.1 PhaseLift algorithm [23–27]: PhaseLift can be represented using stacked image vectors and a linear operator A that maps X to estimated intensities |g|^2.The operator combines the transformation from the lifted variable to predicted measurements.
  • 3.3.1 PhaseLift algorithm [23–27]: When σ is large relative to components of yℓ and O, factorized rank-1 PhaseLift becomes equivalent to intensity-based optimization and uses L-BFGS.L-BFGS approximates Hessian inversion from previous gradients.
  • 3.3.1 PhaseLift algorithm [23–27]: With a Poisson-likelihood cost, rank-1 PhaseLift is equivalent to the paper’s global Newton method for the same reason.This equivalence is stated for the rank-1 case.

4 Performance analysis of various algorithms

Amplitude-based and Poisson-likelihood-based algorithms are more robust to noise and model mismatch than intensity-based methods. Second-order optimization improves reconstruction and convergence, while global updates generally offer little benefit over sequential ones and can cost more computationally.

  • Second-order optimization generally performs better than first-order optimization, with faster convergence and greater stability under imperfect data.The second-order method without pupil estimation can outperform first-order reconstruction with pupil estimation in the aberrated-data example.
  • Intensity-based algorithms produce high-frequency amplitude artifacts and phase-wrapping artifacts when data contain noise or model mismatch.These artifacts resemble those observed in experimental reconstructions.
  • Sequential methods generally match global methods while requiring less computation, and sequential Gauss-Newton is preferred in practice for robustness and runtime.Global Newton’s method requires substantial time to compute the full Hessian, whereas sequential Gauss-Newton combines robustness with fast runtime and reasonable complexity.
  • Amplitude-based and Poisson-likelihood-based algorithms tolerate Poisson noise, aberrations, and LED misalignment better than intensity-based algorithms.Their assumed error model matches intensity-dependent errors in FPM data.
  • Intensity-based cost functions overweight low spatial frequencies, causing the high-frequency artifacts seen in reconstruction.Amplitude-based and Poisson-likelihood-based costs have similar gradients and Hessians and yield similar reconstruction quality.
  • Pupil estimation improves complex-field reconstruction, with second-order Gauss-Newton plus pupil estimation producing the best result in the aberrated-data experiment.The pupil-estimation step approximately doubles runtime because it takes about as long as object reconstruction.

5 Conclusion

The paper formulates FPM phase retrieval through maximum-likelihood cost functions, evaluates correction methods, and identifies reconstruction and efficiency trade-offs among algorithms.

  • Maximum-likelihood analysis classifies FPM algorithms into amplitude-based methods associated with Poisson noise and intensity-based methods associated with white Gaussian noise.
  • A new Poisson-likelihood algorithm is derived to better address measurement imperfections.
  • Second-order pupil-estimation methods give the best results among the compared pupil-estimation algorithms.
  • Simulated annealing can correct misaligned LEDs and estimate their actual positions.
  • Global Newton’s method gives the best reconstruction but may have high computational cost.
  • Sequential Gauss–Newton provides the best trade-off between robustness and computational efficiency for large-scale applications.

Gradient:

The gradient formulations derive update directions for amplitude-, intensity-, and Poisson-likelihood-based objectives, with residuals carrying measurement information into corresponding Fourier regions.

  • The Poisson-likelihood gradient is equivalent to the intensity-based gradient with an added 1/|gℓ|2 weight for each image component.
  • The gradient-descent update uses a line-search step size and either an intensity-based or amplitude-based cost function.
  • Residual terms compare estimates with measurements and provide the information used to update the previous estimate.
  • The Q†ℓdiag(¯P) term maps measurement updates back to their corresponding spatial-frequency regions.
  • For the amplitude-based objective, the residual includes real-domain amplitude replacement, making gradient descent resemble projection-based phase retrieval.

Hessian:

The Hessian captures curvature for complex-valued FPM objectives, and Newton’s method uses its inverse to obtain faster and more stable updates than first-order methods.

  • Hessian entries encode curvature information across dimensions of the multivariate cost function.
  • For convex problems, a positive-semidefinite Hessian reflects nonnegative cost-function curvatures and supports optimization toward a global minimum.
  • The Hessian for a real function of complex variables is represented as a 2n^2 × 2n^2 matrix.
  • The amplitude-, intensity-, and Poisson-likelihood-based Hessians are obtained by differentiating their respective gradients again.
  • Newton’s method inverts the Hessian and generally offers faster convergence and greater stability than gradient descent.
Loading 1511.02986v2…