pyoptics.phase_retrieval

Phase retrieval: recover a phase from intensity measurements.

Two problems are covered:

1. Support-based phase retrieval (coherent diffractive imaging): an object g of limited extent (its support) is to be found from the magnitude |FT(g)| of its far field alone. ER() implements Fienup’s error reduction, HIO() the hybrid input-output algorithm, both with support and optionally positivity constraints and a shrink-wrap support update (ShrinkWrap, Marchesini et al. 2003). With a distance z, the detector sits in the near field instead (Fresnel coherent diffractive imaging): the Fourier transforms are replaced by a pyoptics propagator to z and back.

2. Pupil phase retrieval from through-focus PSFs (phase diversity): the pupil phase of an imaging system is found from intensity PSFs measured at known defocus values. phase_diversity() offers a multi-plane Gerchberg-Saxton (Misell) iteration and a least squares fit of Fringe Zernike coefficients using the PSF model of imaging.get_psf().

Conventions

Fields vary as exp(i(k.r - omega*t)). Far field magnitudes are those of fft.FT() (centered, unnormalized): magnitude = abs(FT(g)). Pupils live on the centered frequency grid of the image plane grid x, y, as in imaging, with normalized pupil coordinates (P_x, P_y) = (k_x, k_y)/(k*NA); the PSF is inv_FT() of the pupil function. Pupil phases are in radians, Zernike coefficients (Fringe numbering, R_norm = 1 on the unit pupil) in waves, as in litho.

Ambiguities

abs(FT(g)) does not change if g is shifted cyclically by whole pixels, multiplied by a constant phase factor or replaced by its twin conj(g(-x, -y)). Support-based retrieval therefore can recover g only up to these operations; align() removes them for a comparison with a reference. In the near field (z given), only the constant phase factor remains: the curvature of the Fresnel kernel distinguishes the object from its twin, and a shift of the object shifts the measured magnitude. Likewise, a single in-focus PSF cannot distinguish a pupil phase phi(P) from -phi(-P), so the sign of the even Zernikes (defocus, astigmatism, spherical aberration, …) is undetermined. A known defocus applied to some of the measurements (phase diversity) removes this ambiguity.

References

R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures”, Optik 35, 237-246 (1972).

J. R. Fienup, “Reconstruction of an object from the modulus of its Fourier transform”, Opt. Lett. 3, 27-29 (1978).

J. R. Fienup, “Phase retrieval algorithms: a comparison”, Appl. Opt. 21, 2758-2769 (1982).

S. Marchesini, H. He, H. N. Chapman, S. P. Hau-Riege, A. Noy, M. R. Howells, U. Weierstall and J. C. H. Spence, “X-ray image reconstruction from a diffraction pattern alone”, Phys. Rev. B 68, 140101 (2003).

G. J. Williams, H. M. Quiney, B. B. Dhal, C. Q. Tran, K. A. Nugent, A. G. Peele, D. Paterson and M. D. de Jonge, “Fresnel coherent diffractive imaging”, Phys. Rev. Lett. 97, 025506 (2006).

D. L. Misell, “An examination of an iterative method for the solution of the phase problem in optics and electron optics: I. Test calculations”, J. Phys. D 6, 2200-2216 (1973).

R. A. Gonsalves, “Phase retrieval and diversity in adaptive optics”, Opt. Eng. 21, 829-832 (1982).

Classes

RetrievalResult	Result of an iterative phase retrieval.
ShrinkWrap	Shrink-wrap support update after Marchesini et al., Phys.

Functions

ER	Fienup's error reduction algorithm.
HIO	Fienup's hybrid input-output algorithm.
align	Remove the trivial ambiguities of support-based phase retrieval.
autocorrelation_support	Initial support estimate from the autocorrelation of the object.
defocus_to_waves	Fringe Zernike Z4 coefficient (in waves) of a paraxial defocus dz.
diversity_psfs	Intensity PSFs at a series of known defoci.
fourier_error	Normalized rms deviation of abs(FT(obj)) from the measured magnitude.
phase_diversity	Retrieve the pupil phase from intensity PSFs at known defoci.
twin	Twin image conj(g(-x, -y)) on a centered grid; it has the same far field magnitude as g.

class pyoptics.phase_retrieval.RetrievalResult(obj: ndarray, support: ndarray, errors: ndarray, errors_object: ndarray = None, phase: ndarray = None, zernikes: dict = <factory>)

Bases: object

Result of an iterative phase retrieval.

Attributes:

objarray

The retrieved object (ER(), HIO()) or pupil function (phase_diversity()).

supportarray of bool

Final support (ER(), HIO()) or pupil (phase_diversity()).

errorsarray

Error metric per iteration (see the function that returned the result).

errors_objectarray

Object domain error per iteration (ER(), HIO() only): the fraction of the energy of the Fourier projected object that violates the object constraints.

phasearray

Pupil phase in radians on the frequency grid, zero outside the pupil (phase_diversity() only).

zernikesdict

Fringe Zernike coefficients of the pupil phase in waves, {j: c_j} for j = 2, …, n_zernike (phase_diversity() only).

pyoptics.phase_retrieval.fourier_error(obj, magnitude)

Normalized rms deviation of abs(FT(obj)) from the measured magnitude.

E_F = sqrt(sum((abs(FT(obj)) - magnitude)**2)/sum(magnitude**2)).

Parameters:

obj2d array

Object estimate.

magnitude2d array

Measured Fourier magnitude (centered, as FT() returns it).

Returns:

errornumber

pyoptics.phase_retrieval.autocorrelation_support(magnitude, threshold=0.04, sigma=1.0)

Initial support estimate from the autocorrelation of the object.

The inverse transform of magnitude**2 is the autocorrelation of the object, whose support contains the support of the object (it is twice as wide). It is a safe, loose start for ShrinkWrap.

Parameters:

magnitudearray

Far field magnitude abs(FT(g)).

thresholdnumber, optional

Fraction of the maximum of the (smoothed) autocorrelation above which a pixel belongs to the support.

sigmanumber, optional

Width in pixels of a Gaussian smoothing applied before thresholding.

Returns:

supportarray of bool

class pyoptics.phase_retrieval.ShrinkWrap(every=20, sigma=3.0, sigma_min=1.0, decay=0.95, threshold=0.2)

Bases: object

Shrink-wrap support update after Marchesini et al., Phys. Rev. B 68, 140101 (2003).

Every every iterations, the modulus of the current object estimate is smoothed by a Gaussian of width sigma (in pixels) and thresholded at threshold times its maximum; this is the new support. sigma is multiplied by decay after each update until it reaches sigma_min, so the support tightens gradually around the object.

Parameters:

everyint, optional

Number of iterations between two support updates.

sigmanumber, optional

Initial width of the Gaussian smoothing in pixels.

sigma_minnumber, optional

Final width of the Gaussian smoothing.

decaynumber, optional

Factor applied to sigma after every update.

thresholdnumber, optional

Fraction of the maximum of the smoothed modulus.

pyoptics.phase_retrieval.ER(magnitude, support, N_iter=200, start=None, positivity=False, real=False, shrinkwrap=None, rng=None, callback=None, z=None, x=None, y=None, wavelength=None, n=1.0, propagator=None)

Fienup’s error reduction algorithm.

Alternates the projection onto the measured Fourier magnitude, g’ = inv_FT(magnitude*exp(i*angle(FT(g)))), with the object constraints: g_new = g’ where g’ meets the constraints (inside the support and, with positivity, Re(g’) >= 0) and 0 elsewhere. The Fourier error never increases, but ER stagnates easily; it is mostly used to polish the result of HIO().

Parameters:

magnitudearray

Measured far field magnitude abs(FT(g)) (the square root of the intensity), centered as returned by fft.FT().

supportarray of bool or None

Pixels in which the object may be nonzero. None: start from autocorrelation_support() (use with shrinkwrap).

N_iterint, optional

Number of iterations.

startarray, optional

Start object, e.g. the result of a previous run. Defaults to random values inside the support.

positivitybool, optional

Require a real, nonnegative object.

realbool, optional

Require a real object (implied by positivity).

shrinkwrapShrinkWrap, True or None, optional

Update the support during the iteration; True uses ShrinkWrap() with default parameters.

rngint or numpy.random.Generator, optional

Seed for the random start object.

callbackcallable, optional

Called after every iteration as callback(iteration, estimate, support), with the iteration counted from 1 and the current estimate (the Fourier projection with the object constraints applied), e.g. to record snapshots.

znumber, optional

Distance of the detector from the object. None (default): far field, magnitude = abs(fft.FT(g)). Otherwise magnitude is the modulus of the field propagated by z on the same grid, magnitude = abs(propagator(g, x, y, z, wavelength, n)[0]), and the projection propagates back by -z (Fresnel coherent diffractive imaging). The support is then required.

x, yarrays, optional

Linear grid coordinates (needed with z).

wavelengthnumber, optional

Vacuum wavelength (needed with z).

nnumber, optional

Refractive index between object and detector.

propagatorcallable, optional

Grid preserving pyoptics propagator with the signature propagator(field, x, y, z, wavelength, n) -> (field, x, y), used for +z and -z, e.g. propagators.fresnel_TF or functools.partial(propagators.propagate_padded, propagators.rayleigh_sommerfeld_I_TF). Defaults to the angular spectrum method propagators.rayleigh_sommerfeld_I_TF, which needs a grid without evanescent frequencies, dx >= wavelength/(2*n), so that the back propagation is the inverse of the forward one.

Returns:

resultRetrievalResult

obj: retrieved object (real if positivity or real), support: final support, errors: normalized rms deviation of the modulus in the detector plane from magnitude, before each projection (in the far field this is fourier_error()), errors_object: object domain error.

Notes

In the far field, the object is found only up to a cyclic shift by whole pixels, a constant phase factor and the twin image conj(g(-x, -y)); see align(). In the near field, only the constant phase factor remains.

pyoptics.phase_retrieval.HIO(magnitude, support, beta=0.9, N_iter=500, start=None, positivity=False, real=False, shrinkwrap=None, rng=None, callback=None, z=None, x=None, y=None, wavelength=None, n=1.0, propagator=None)

Fienup’s hybrid input-output algorithm.

Like ER(), but outside the valid region the iterate is not set to zero: g_new = g’ where g’ meets the object constraints and g - beta*g’ elsewhere (Fienup, Appl. Opt. 21, 2758 (1982)). This negative feedback pushes the iterate out of the local minima in which ER stagnates. The Fourier error of HIO does not decrease monotonically; a few ER() iterations with start=result.obj usually finish the reconstruction.

Parameters:

magnitude, support, N_iter, start, positivity, real, shrinkwrap, rng,

callback, z, x, y, wavelength, n, propagator

See ER().

betanumber, optional

Feedback parameter, typically 0.5 to 1.

Returns:

resultRetrievalResult

See ER(). obj is the last Fourier projection with the object constraints applied (the iterate g itself is not an estimate of the object outside the support).

pyoptics.phase_retrieval.twin(obj)

Twin image conj(g(-x, -y)) on a centered grid; it has the same far field magnitude as g.

pyoptics.phase_retrieval.align(obj, reference, allow_twin=True)

Remove the trivial ambiguities of support-based phase retrieval.

Finds the cyclic shift by whole pixels, the constant phase factor and, if allow_twin, the choice between obj and its twin() that bring obj closest to reference.

Parameters:

objarray

Retrieved object.

referencearray

True object.

allow_twinbool, optional

Also try the twin image.

Returns:

alignedarray

Shifted, phase corrected (and possibly twinned) obj.

errornumber

Normalized rms difference norm(aligned - reference)/norm(reference).

pyoptics.phase_retrieval.defocus_to_waves(dz, NA, k)

Fringe Zernike Z4 coefficient (in waves) of a paraxial defocus dz.

Observing the PSF a distance dz behind the focal plane (dz > 0 further downstream, in the direction of propagation) multiplies the pupil by exp(i*k_z*dz); paraxially k_z = k - k*NA**2*P**2/2 with the normalized pupil radius P. Up to piston this is the Fringe Zernike Z4 = 2*P**2 - 1 with the coefficient

c_4 = -k*dz*NA**2/(8*pi)  (waves).

Parameters:

dznumber or array

Defocus distance.

NAnumber

Numerical aperture as used in imaging (the pupil edge is at |k_perp| = k*NA).

knumber

Wavenumber in the image space.

Returns:

c_4number or array

Z4 coefficient in waves.

pyoptics.phase_retrieval.diversity_psfs(x, y, NA, k, defocus_waves=None, pupil_mask=None, phase=None, defocus_z=None)

Intensity PSFs at a series of known defoci.

Parameters:

x, yarrays

Linear image plane coordinates.

NAnumber

Numerical aperture.

knumber

Wavenumber.

defocus_wavessequence of numbers, optional

Fringe Z4 coefficient (waves) added to the pupil phase for each PSF, see defocus_to_waves() (paraxial defocus).

pupil_mask, phasearrays or callables, optional

See imaging.pupil_function(); phase is the unknown pupil phase in radians.

defocus_zsequence of numbers, optional

Instead of defocus_waves: distances of the observation planes from the focal plane (positive downstream). The pupil is multiplied by the exact angular spectrum propagator exp(i*(k_z - k)*dz), so the PSFs equal propagators.rayleigh_sommerfeld_I_TF() of the in-focus PSF.

Returns:

psfsarray

Intensity PSFs |inv_FT(P*D_m)|**2 with the diversity factors D_m, shape (M, Ny, Nx), each normalized to a total of 1.

pyoptics.phase_retrieval.phase_diversity(psfs, defocus_waves, x, y, NA, k, n_zernike=15, method='fit', pupil_mask=None, N_iter=None, start=None, tol=1e-12, solver='trf', defocus_z=None)

Retrieve the pupil phase from intensity PSFs at known defoci.

Parameters:

psfsarray

Measured intensity PSFs, shape (M, Ny, Nx), sampled on x, y with the PSF at the grid center as in imaging.get_psf(). Their scale does not matter; each is normalized to the energy of the model PSF.

defocus_wavessequence of M numbers or None

Known defocus of each PSF as Fringe Z4 coefficient in waves (see defocus_to_waves()), i.e. paraxial defocus. None if defocus_z is given. At least one value should be nonzero: with a single in-focus PSF, the sign of the even aberrations is ambiguous. Values of about +-0.5 to +-1 wave work well.

x, yarrays

Linear image plane coordinates.

NAnumber

Numerical aperture.

knumber

Wavenumber.

n_zernikeint, optional

Fringe Zernikes 2, …, n_zernike describe the phase (piston is irrelevant).

method{“fit”, “gs”}, optional

“fit” (default): least squares fit of the Zernike coefficients to the measured intensities (scipy.optimize.least_squares with the analytic Jacobian). Fast and accurate if the start is within the basin of attraction; noise is handled in the least squares sense. “gs”: multi-plane Gerchberg-Saxton (Misell) iteration, which cycles through the planes and replaces the modulus of the PSF by the measured one. It needs no basis, so it also finds phases that are not described by n_zernike Zernikes. Its phase is unwrapped on the pupil (unwrap.unwrap_quality_guided()) and returned as is; the Zernike coefficients are fitted to it. Convergence is slower than “fit”, and a “gs” result can serve as start of a “fit”.

pupil_maskarray or callable, optional

Known pupil amplitude, see imaging.pupil_function().

N_iterint, optional

“gs”: number of iterations (default 200). “fit”: maximum number of function evaluations (default 200).

startdict or array, optional

Start phase: Fringe Zernike coefficients {j: c_j} in waves or a phase array in radians on the frequency grid. Defaults to zero.

tolnumber, optional

“fit” only: tolerances ftol, xtol and gtol of least_squares.

solver{“trf”, “lm”, “dogbox”}, optional

“fit” only: algorithm of scipy.optimize.least_squares. The trust region default “trf” is robust; Levenberg-Marquardt (“lm”) is somewhat faster on well posed problems but can run away on ill posed ones (e.g. a single in-focus PSF).

defocus_zsequence of M numbers, optional

Instead of defocus_waves: the distances of the PSF planes from the focal plane (positive downstream, in the units of x and 2*pi/k), e.g. the positions of a focus stage. The defocus is applied exactly (non-paraxially) as the angular spectrum propagator, see diversity_psfs().

Returns:

resultRetrievalResult

obj: pupil function, support: pupil (bool), phase: pupil phase in radians (zero outside the pupil, piston removed), zernikes: Fringe Zernike coefficients in waves {j: c_j}, errors: normalized rms deviation of the model amplitudes sqrt(I) from the measured ones, per iteration (“gs”) or per function evaluation (“fit”).

Notes

The PSF grid has to sample the intensity: dx <= wavelength/(4*NA). The pupil is sampled with about L*NA/wavelength pixels across its radius (L: width of the window); a PSF window much larger than the PSF itself thus resolves the pupil better, as in imaging.