pyoptics.TIE

Routines for TIE

Convention: fields vary as exp(i(kz - omega*t)) and are written as sqrt(I)*exp(i*phase). The transport of intensity equation then reads

k dI/dz = -div(I grad(phase)),

so a converging wave (phase decreasing away from the axis, e.g. a thin lens exp(-i k r**2/(2f))) has dI/dz > 0 on axis.

The module contains no propagator of its own: to compute defocused images, use pyoptics.propagators, e.g. rayleigh_sommerfeld_I_TF(), which follows the same convention.

References

M. R. Teague, “Deterministic phase retrieval: a Green’s function solution”, J. Opt. Soc. Am. 73, 1434-1441 (1983).

N. Streibl, “Phase imaging by the transport equation of intensity”, Opt. Commun. 49, 6-10 (1984).

D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light”, Phys. Rev. Lett. 80, 2586-2589 (1998).

C. Zuo et al., “Transport of intensity equation: a tutorial”, Opt. Lasers Eng. 135, 106187 (2020).

Functions

RMS
desymmetrize_image	Upper left quarter of an image, the inverse of symmetrize_image().
desymmetrize_stack	Inverse of symmetrize_stack() (not implemented yet, returns None).
energy_normalize_stack	Scale all images of a stack to the energy of one of them.
enhanced_tie_solve	Solve the TIE iteratively, correcting the residual flux.
get_image_energy	Energy of an intensity image, i.e. the sum of its pixels.
get_k_grid	Angular frequency grid for centered FFTs.
inv_laplacian	Inverse Laplacian with Fourier methods and Tikhonov regularization.
laplacian	Laplacian of f with Fourier derivatives.
low_pass	Remove spatial frequencies that the imaging system cannot transmit.
low_pass_stack	Apply low_pass() to each image of a stack.
scaled_laplace_like	Compute div(c*grad(f)) with Fourier derivatives.
smoothed_deriv	Derivative along the stack axis from finite difference weights.
solve_tie	Solve the transport of intensity equation for the phase.
symmetrize_image	Antisymmetric extension of an image to twice its size.
symmetrize_stack	Apply symmetrize_image() to each image of a stack.
tie	Forward transport of intensity equation.

pyoptics.TIE.low_pass(img, x, y, sigma_cutoff, NA, lam)

Remove spatial frequencies that the imaging system cannot transmit.

All frequencies beyond sigma_cutoff times the cut-off k*NA are set to zero (hard filter); the result is rescaled to the energy (sum) of the input.

Parameters:

img2d array

Intensity image.

x, y1d arrays

Coordinates of the image (equidistant).

sigma_cutoffnumber

Cut-off in units of the coherent cut-off frequency k*NA.

NAnumber

Numerical aperture.

lamnumber

Wavelength.

Returns:

filtered2d array

Filtered (real) image with the energy of img.

pyoptics.TIE.low_pass_stack(stack, x, y, sigma_cutoff, NA, lam)

Apply low_pass() to each image of a stack.

Parameters:

stack3d array

Image stack of shape (Ny, Nx, N_images).

x, y, sigma_cutoff, NA, lam

See low_pass().

Returns:

filtered3d array

Filtered stack, same shape as stack.

pyoptics.TIE.smoothed_deriv(stack, x, x0, order)

Derivative along the stack axis from finite difference weights.

The weights for the derivative of the given order at x0 from the samples at x are computed with sympy (optional dependency) and applied to all pixels.

Parameters:

stack3d array

Image stack of shape (Ny, Nx, len(x)).

xsequence of numbers

Positions (e.g. defocus values z) of the images in the stack.

x0number

Position at which to evaluate the derivative.

orderint

Order of the derivative.

Returns:

deriv2d array

Derivative at x0, shape (Ny, Nx).

pyoptics.TIE.get_k_grid(x, y)

Angular frequency grid for centered FFTs.

Parameters:

x, y1d arrays

Equidistant coordinates.

Returns:

KX, KY2d arrays

k_x and k_y (2*pi times the spatial frequencies) on a meshgrid, centered as returned by fftshift().

pyoptics.TIE.scaled_laplace_like(f, c, x, y)

Compute div(c*grad(f)) with Fourier derivatives.

For constant c this is c times the Laplacian of f.

Parameters:

f2d array

Scalar function.

c2d array or number

Scalar weight applied to the gradient.

x, y1d arrays

Equidistant coordinates.

Returns:

U2d array

div(c*grad(f)) (real part).

pyoptics.TIE.laplacian(f, x, y)

Laplacian of f with Fourier derivatives.

Parameters:

f2d array

Scalar function.

x, y1d arrays

Equidistant coordinates.

Returns:

L2d array

Laplacian of f (real part).

pyoptics.TIE.inv_laplacian(f, x, y, alpha=0.0)

Inverse Laplacian with Fourier methods and Tikhonov regularization.

Divides the spectrum by -(k_x**2 + k_y**2 + alpha). Without regularization, the zero frequency is undefined and set to zero, i.e. the mean free solution is returned.

Parameters:

f2d array

Scalar function.

x, y1d arrays

Equidistant coordinates.

alphanumber, optional

Regularization parameter (in units of k**2).

Returns:

L_inv2d array

Inverse Laplacian of f (real part).

pyoptics.TIE.solve_tie(I, dI_dz, x, y, lam, alpha)

Solve the transport of intensity equation for the phase.

Inverts k dI/dz = -div(I grad(phase)) with two inverse Laplacians (see the module docstring for the sign convention).

Parameters:

I2d array

Intensity in focus.

dI_dz2d array

Axial derivative of the intensity.

x, y1d arrays

Equidistant coordinates.

lamnumber

Wavelength.

alphanumber

Regularization parameter, see inv_laplacian().

Returns:

phase2d array

Recovered phase in radians (mean free).

pyoptics.TIE.tie(phase, I, x, y, lam)

Forward transport of intensity equation.

Compute dI/dz = -div(I grad(phase))/k from I and a phase (see module docstring for the sign convention).

Parameters:

phase2d array

Phase in radians.

I2d array

Intensity.

x, y1d arrays

Equidistant coordinates.

lamnumber

Wavelength.

Returns:

dI_dz2d array

Axial derivative of the intensity.

pyoptics.TIE.enhanced_tie_solve(N, I, dI_dz, x, y, lam, alpha, true_phase=None)

Solve the TIE iteratively, correcting the residual flux.

Starts with solve_tie(), then repeats N times: compute the flux of the current phase with tie(), solve for the phase of the residual flux and add it to the estimate.

Parameters:

Nint

Number of refinement iterations.

I, dI_dz, x, y, lam, alpha

See solve_tie().

true_phase2d array, optional

Known phase; if given, the RMS error is printed after each iteration.

Returns:

phase2d array

Refined phase in radians.

pyoptics.TIE.get_image_energy(image)

Energy of an intensity image, i.e. the sum of its pixels.

Parameters:

imagearray

Intensity image.

Returns:

energynumber

pyoptics.TIE.energy_normalize_stack(stack, ind)

Scale all images of a stack to the energy of one of them.

Parameters:

stack3d array

Image stack of shape (Ny, Nx, N_images).

indint

Index of the reference image.

Returns:

normalized3d array

Copy of the stack with equal energies.

pyoptics.TIE.symmetrize_image(image)

Antisymmetric extension of an image to twice its size.

The image is mirrored with a sign change across its right and lower edges (upper left: image, upper right and lower left: -mirror, lower right: double mirror), which removes the jumps at the borders for FFT based solvers.

Parameters:

image2d array

Image of shape (Ny, Nx).

Returns:

symmetrized2d array

Image of shape (2*Ny, 2*Nx).

pyoptics.TIE.desymmetrize_image(image)

Upper left quarter of an image, the inverse of symmetrize_image().

Parameters:

image2d array

Symmetrized image of shape (2*Ny, 2*Nx).

Returns:

image2d array

Image of shape (Ny, Nx).

pyoptics.TIE.symmetrize_stack(stack, x, y)

Apply symmetrize_image() to each image of a stack.

Parameters:

stack3d array

Image stack of shape (Ny, Nx, N_images).

x, y1d arrays

Equidistant coordinates of the images.

Returns:

symmetrized3d array

Stack of shape (2*Ny, 2*Nx, N_images).

x_new, y_new1d arrays

Coordinates of the symmetrized images, starting at x[0] and y[0].

pyoptics.TIE.desymmetrize_stack(stack, x, y)

Inverse of symmetrize_stack() (not implemented yet, returns None).

Parameters:

stack3d array

Symmetrized stack.

x, y1d arrays

Coordinates of the symmetrized stack.