pyoptics.holography

Digital holography: recording, reconstruction, refocusing and autofocus.

A hologram is the intensity of an object wave O and a reference wave R on a camera,

H = |O + R|**2 = |O|**2 + |R|**2 + O*conj(R) + conj(O)*R.

The first two terms are the zero order, O*conj(R) is the object (+1 order) and conj(O)*R its twin (-1 order). Recording and reconstruction:

• hologram(), phase_shifted_holograms(), reference_wave(): simulate off-axis and in-line holograms of a given object field (add camera noise with processing.CameraNoiser).

• off_axis_angle_range(), max_object_bandwidth(): the sampling condition that separates the orders on the camera.

• reconstruct_off_axis(): single shot reconstruction by sideband filtering.

• reconstruct_phase_shifting(): in-line or off-axis phase-shifting holography from N frames with known reference phase shifts (least squares, via interferometry.psi_n_step()).

• reconstruct_gabor(): in-line (Gabor) holography with plane wave illumination, optionally with iterative twin image suppression.

• suppress_zero_order(): remove the zero order from a hologram.

• find_carrier(): locate the +1 order in the hologram spectrum.

Working with the reconstructed field:

• refocus(): numerical refocusing with the angular spectrum or Fresnel propagator (padded).

• focus_metric(), autofocus(): sharpness metrics (amplitude variance, Tamura coefficient, gradient, Laplacian) and a search for the best focus.

• quantitative_phase(), remove_tilt(), optical_path_difference(): unwrapped phase maps (with pyoptics.unwrap), background compensation and OPD.

Conventions: fields vary as exp(i(k.r - omega*t)). The reference wave tilted by the angles (theta_x, theta_y) is R = A*exp(i*(k_x*x + k_y*y + delta)) with (k_x, k_y) = k*(sin(theta_x), sin(theta_y)), its “carrier”. The object term O*conj(R) thus sits at -(k_x, k_y) in the spectrum FT(H) and the twin at +(k_x, k_y). All reconstructions return the object field O in the hologram plane (scaled to the reference amplitude). An object at the distance d in front of the camera is brought back into focus with refocus(O, x, y, -d, …): propagation by z > 0 multiplies each plane wave by exp(+i*k_z*z), so the camera plane is reached from the object with +d and the object from the camera with -d. Wavelengths are vacuum wavelengths; n is the refractive index of the medium between object and camera.

References

4. Gabor, “A new microscopic principle”, Nature 161, 777-778 (1948).

E. N. Leith and J. Upatnieks, “Reconstructed wavefronts and communication theory”, J. Opt. Soc. Am. 52, 1123-1130 (1962).

U. Schnars and W. Jüptner, “Direct recording of holograms by a CCD target and numerical reconstruction”, Appl. Opt. 33, 179-181 (1994).

E. Cuche, P. Marquet and C. Depeursinge, “Simultaneous amplitude-contrast and quantitative phase-contrast microscopy by numerical reconstruction of Fresnel off-axis holograms”, Appl. Opt. 38, 6994-7001 (1999).

E. Cuche, P. Marquet and C. Depeursinge, “Spatial filtering for zero-order and twin-image elimination in digital off-axis holography”, Appl. Opt. 39, 4070-4075 (2000).

I. Yamaguchi and T. Zhang, “Phase-shifting digital holography”, Opt. Lett. 22, 1268-1270 (1997).

T. Latychevskaia and H.-W. Fink, “Solution to the twin image problem in holography”, Phys. Rev. Lett. 98, 233901 (2007).

P. Memmolo, C. Distante, M. Paturzo, A. Finizio, P. Ferraro and B. Javidi, “Automatic focusing in digital holography and its application to stretched holograms”, Opt. Lett. 36, 1945-1947 (2011).

21. Schnars and W. Jüptner, “Digital Holography”, Springer (2005).

M. K. Kim, “Digital Holographic Microscopy: Principles, Techniques, and Applications”, Springer (2011).

Functions

autofocus	Find the best focus of a reconstructed field.
find_carrier	Locate the off-axis +1 order in the spectrum of a hologram.
focus_metric	Sharpness metric of the amplitude of a field.
hologram	Record the hologram |O + R|**2 of an object field O.
max_object_bandwidth	Largest object bandwidth for which an off-axis hologram on this grid can separate the orders.
off_axis_angle_range	Range of reference tilt angles that separate the orders of an off-axis hologram.
optical_path_difference	Optical path difference phase/k_0.
phase_shifted_holograms	Record a stack of holograms with phase-shifted references.
quantitative_phase	Quantitative phase map of a reconstructed (and refocused) field.
reconstruct_gabor	Reconstruct an in-line (Gabor) hologram recorded with plane wave illumination.
reconstruct_off_axis	Reconstruct the object field from an off-axis hologram by sideband filtering.
reconstruct_phase_shifting	Reconstruct the object field from phase-shifted holograms.
reference_wave	Tilted plane reference wave on the camera.
refocus	Refocus a reconstructed field numerically.
remove_tilt	Subtract a least squares plane (tilt and optionally piston) from a phase map.
suppress_zero_order	Remove the zero order |O|**2 + |R|**2 from a hologram.

pyoptics.holography.reference_wave(x, y, wavelength, ref_angle=(0.0, 0.0), ref_amp=1.0, phase_shift=0.0, n=1.0)

Tilted plane reference wave on the camera.

Parameters:

x, yarrays

Linear or meshed coordinates of the camera.

wavelengthnumber

ref_angletuple of two numbers, optional

Tilt angles (theta_x, theta_y) in radians: the reference has the transverse wave vector k*(sin(theta_x), sin(theta_y)).

ref_ampnumber, optional

Amplitude A of the reference.

phase_shiftnumber, optional

Additional phase delta (e.g. of a phase shifter).

nnumber, optional

Refractive index.

Returns:

Rarray

A*exp(i*(k_x*x + k_y*y + delta)).

pyoptics.holography.max_object_bandwidth(x, y)

Largest object bandwidth for which an off-axis hologram on this grid can separate the orders.

The +1 order (radius B) has to stay clear of the zero order (|O|**2 has radius 2B), so the carrier must satisfy |k_c| >= 3B, and it has to stay below the Nyquist frequency, |k_c| + B <= pi/dx. Both are possible for B <= pi/(4*dx) with the carrier along one axis of the grid. (A diagonal carrier allows somewhat more, but the axis-aligned value is the simple and safe bound.)

Parameters:

x, yarrays

Camera coordinates.

Returns:

B_maxnumber

Maximum object bandwidth k*NA (radius in k space) in rad per length unit, along the coarser axis.

pyoptics.holography.off_axis_angle_range(bandwidth, x, y, wavelength, n=1.0, axis='x')

Range of reference tilt angles that separate the orders of an off-axis hologram.

For an object with the bandwidth B = k*NA (radius of its angular spectrum), the +1 order must not overlap the zero order, |k_c| >= 3B, and must be sampled by the camera, |k_c| + B <= pi/dx, where k_c = k*sin(theta) is the carrier.

Parameters:

bandwidthnumber

Object bandwidth B (radius in k space, rad per length unit), e.g. k*NA of the imaging optics in front of the camera.

x, yarrays

Camera coordinates.

wavelengthnumber

nnumber, optional

axis“x” or “y”, optional

Direction of the tilt.

Returns:

theta_min, theta_maxnumbers

Smallest and largest tilt in radians. Raises ValueError if no angle works (theta_min > theta_max); then the object has to be band limited further or the camera pixels have to be smaller.

pyoptics.holography.hologram(obj_field, x, y, wavelength, ref_angle=(0.0, 0.0), ref_amp=1.0, phase_shift=0.0, n=1.0)

Record the hologram |O + R|**2 of an object field O.

Parameters:

obj_fieldarray

Object wave on the camera (i.e. already propagated there).

x, yarrays

Camera coordinates.

wavelengthnumber

ref_angletuple of two numbers, optional

Reference tilt (theta_x, theta_y) in radians; (0, 0) gives an in-line hologram, see reference_wave().

ref_ampnumber, optional

Reference amplitude.

phase_shiftnumber, optional

Phase shift of the reference (phase-shifting holography).

nnumber, optional

Returns:

Harray

Hologram intensity. Apply processing.CameraNoiser for a noisy recording (after scaling to photon numbers with processing.Scaler). Warns if the carrier exceeds the Nyquist frequency of the camera.

pyoptics.holography.phase_shifted_holograms(obj_field, x, y, wavelength, shifts, ref_angle=(0.0, 0.0), ref_amp=1.0, n=1.0)

Record a stack of holograms with phase-shifted references.

Parameters:

obj_field, x, y, wavelength, ref_angle, ref_amp, n

As for hologram().

shiftssequence of numbers

Reference phase shifts delta_j in radians, e.g. interferometry.phase_shifts(4) for the 4-step method.

Returns:

framesarray of shape (Ny, Nx, len(shifts))

Frame index last, as in pyoptics.processing and pyoptics.interferometry.

pyoptics.holography.suppress_zero_order(H, x=None, y=None, method='mean', I_ref=None, I_obj=None, radius=None)

Remove the zero order |O|**2 + |R|**2 from a hologram.

Parameters:

Harray

Hologram.

x, yarrays, optional

Camera coordinates (needed for method=”highpass”).

methodstr, optional

• “mean”: subtract the mean of H (removes the DC peak only).

• “highpass”: subtract a Gaussian low pass of H with the 1/e radius radius in k space (rad per length unit), which removes the slowly varying |O|**2 + |R|**2 for objects with little bandwidth.

• “recorded”: subtract separately recorded intensities I_ref = |R|**2 and I_obj = |O|**2 (exact).

I_ref, I_objarrays or numbers, optional

Intensities of reference and object alone (method=”recorded”; a missing I_obj is treated as zero).

radiusnumber, optional

Low pass radius for method=”highpass”.

Returns:

H_filteredarray

H without the zero order: ideally O*conj(R) + conj(O)*R.

pyoptics.holography.find_carrier(H, x, y, exclude_radius=None)

Locate the off-axis +1 order in the spectrum of a hologram.

Parameters:

Harray

Hologram.

x, yarrays

Camera coordinates.

exclude_radiusnumber, optional

Radius around k = 0 (rad per length unit) that is ignored, so that the zero order is not picked. Defaults to a tenth of the Nyquist frequency.

Returns:

k_x, k_ynumbers

Carrier of the reference (the +1 order O*conj(R) sits at -(k_x, k_y)), refined to sub-pixel precision by a parabola fit of log|FT(H)|.

Notes

The spectrum of a real hologram is symmetric, |FT(H)(k)| = |FT(H)(-k)|, so the +1 order cannot be told from the twin. The carrier with k_x > 0 (or k_y > 0 if k_x = 0) is returned. If the true reference has the opposite sign, the reconstruction gives conj(O), which is in focus at +d instead of -d. Pass the known reference angle whenever possible.

pyoptics.holography.reconstruct_off_axis(H, x, y, wavelength, ref_angle=None, carrier=None, window='circle', radius=None, ref_amp=1.0, n=1.0, reference=None, subtract_mean=True, return_spectrum=False)

Reconstruct the object field from an off-axis hologram by sideband filtering.

The hologram is multiplied by the (numerical) reference wave R, which moves the object term O*|R|**2 to k = 0, the zero order to the carrier k_c and the twin to 2*k_c. A window of the given radius around k = 0 then selects the object.

Parameters:

Harray

Off-axis hologram.

x, yarrays

Camera coordinates.

wavelengthnumber

ref_angletuple of two numbers, optional

Reference tilt (theta_x, theta_y) in radians.

carriertuple of two numbers, optional

Carrier (k_x, k_y) instead of ref_angle. If neither is given, it is found with find_carrier() (with the sign ambiguity described there).

windowstr or array, optional

“circle” (default; hard edge: keeps the full object spectrum inside the radius, but rings if the spectrum does not vanish there), “tukey” (flat up to half the radius, cosine edge), “hann”, “gaussian” (smooth: less ringing and noise, but they attenuate the object spectrum like an apodized pupil, so use a radius well above the object bandwidth), or an array on the k grid (centred at k = 0) that multiplies the shifted spectrum.

radiusnumber, optional

Window radius in rad per length unit; the object bandwidth k*NA if it is known. Defaults to |k_c|/3, the largest radius that keeps the zero order out for an object fulfilling |k_c| >= 3B.

ref_ampnumber, optional

Reference amplitude A; the result is divided by A.

nnumber, optional

referencearray, optional

The reference wave R on the camera, if it is not a plane wave (e.g. a spherical reference). It replaces the plane wave built from ref_angle or carrier, so that its curvature is removed before filtering; ref_angle or carrier still set the default window radius. A curved reference that is treated as a plane wave leaves a lens phase on O, which shifts the focus and the magnification on refocusing.

subtract_meanbool, optional

If True, subtract the mean of H first. This removes the strong DC peak, whose spectral leakage (the reference carrier is in general not periodic in the window) would otherwise spread into the object order.

return_spectrumbool, optional

If True, also return the shifted spectrum FT(H*R) and the window.

Returns:

Oarray

Object field in the hologram plane.

spectrum, windowarrays

Only if return_spectrum is True.

pyoptics.holography.reconstruct_phase_shifting(frames, shifts=None, x=None, y=None, wavelength=None, ref_angle=(0.0, 0.0), ref_amp=1.0, n=1.0, reference=None)

Reconstruct the object field from phase-shifted holograms.

Each frame is H_j = |O + R*exp(i*delta_j)|**2 = a + b*cos(phi - delta_j) with phi = arg(O*conj(R)) and b = 2*|O*conj(R)|. The least squares phase-shifting algorithm interferometry.psi_n_step() gives phi, a and the modulation b/a per pixel, hence C = O*conj(R) = (b/2)*exp(i*phi) and O = C*R/|R|**2. This is exact for any N >= 3 distinct known shifts and removes the zero order and the twin image, also for in-line holograms.

Parameters:

framesarray of shape (Ny, Nx, N)

Phase-shifted holograms, frame index last.

shiftssequence of N numbers, optional

Reference phase shifts delta_j in radians. Defaults to 2*pi*j/N.

x, y, wavelengthoptional

Camera coordinates and wavelength, needed to build the reference from ref_angle (not needed for an in-line reference).

ref_angletuple of two numbers, optional

Reference tilt (theta_x, theta_y); (0, 0) for in-line holography.

ref_ampnumber, optional

nnumber, optional

referencearray, optional

Reference wave R (at delta = 0) on the camera, overrides ref_angle and ref_amp, e.g. for a spherical reference.

Returns:

Oarray

Object field in the hologram plane.

pyoptics.holography.reconstruct_gabor(H, x, y, z, wavelength, n_iter=0, support=None, ref_amp=1.0, n=1.0, propagator='asm', pad_factor=2, return_errors=False)

Reconstruct an in-line (Gabor) hologram recorded with plane wave illumination.

The object, a transmission t = 1 + o at the distance z in front of the camera, is illuminated by a plane wave of amplitude A. The camera records H = A**2*|1 + O|**2 with O the scattered wave on the camera. Propagating H/A**2 = 1 + O + conj(O) + |O|**2 back by -z gives 1 + o (in focus) plus the twin image, conj(O) propagated by -z, which is o* defocused by 2z, and the small |O|**2 term.

With n_iter > 0, the twin image is suppressed iteratively (Gerchberg-Saxton type, cf. Latychevskaia and Fink, PRL 98, 233901 (2007)): in the object plane the transmission is forced to |t| <= 1 (absorption, no gain) and to t = 1 outside the support; in the camera plane the amplitude is replaced by the measured sqrt(H)/A.

Parameters:

Harray

In-line hologram.

x, yarrays

Camera coordinates.

znumber

Object-camera distance (> 0).

wavelengthnumber

n_iterint, optional

Number of twin image suppression iterations (0: plain back propagation).

supportbool array, optional

Pixels that may contain the object; outside, t = 1. Without support only the absorption constraint is used.

ref_ampnumber, optional

Illumination amplitude A.

n, propagator, pad_factor

As for refocus().

return_errorsbool, optional

If True, also return the rms amplitude mismatch in the camera plane for each iteration.

Returns:

tarray

Object transmission in the object plane (not in the camera plane, as the other reconstructions).

errorsarray

Only if return_errors is True.

pyoptics.holography.refocus(field, x, y, z, wavelength, n=1.0, propagator='asm', pad_factor=2, **kwargs)

Refocus a reconstructed field numerically.

Parameters:

fieldarray

Complex field, e.g. the output of reconstruct_off_axis().

x, yarrays

Coordinates (centred grid, see pyoptics.propagators).

znumber

Propagation distance. Positive z propagates downstream (away from the object); an object at distance d in front of the camera is in focus for z = -d.

wavelengthnumber

nnumber, optional

propagator“asm”, “fresnel” or callable, optional

“asm”: band-limited angular spectrum (rayleigh_sommerfeld_I_TF with band_limit=True), exact for all angles; evanescent waves are removed, because backward propagation would amplify them. “fresnel”: paraxial transfer function (fresnel_TF). A callable is used as a pyoptics propagator.

pad_factornumber, optional

Zero padding factor (1 for none) against wrap-around.

kwargs

Passed on to the propagator.

Returns:

field_refocusedarray

pyoptics.holography.focus_metric(field, metric='tamura', mask=None)

Sharpness metric of the amplitude of a field.

Parameters:

fieldarray

Complex field (only its amplitude |field| enters).

metricstr, optional

• “variance”: variance of the normalized amplitude, var(|U|)/mean(|U|)**2

• “tamura”: Tamura coefficient sqrt(std(|U|)/mean(|U|))

• “gradient”: mean squared gradient of |U|, normalized by mean(|U|)**2

• “laplacian”: mean squared Laplacian of |U|, normalized likewise

maskbool array, optional

Evaluate the metric only where mask is True.

Returns:

valuefloat

Notes

For amplitude objects all metrics are largest in focus. Pure phase objects have a uniform amplitude in focus, the metrics are smallest there (they rise on both sides, with phase contrast of opposite sign).

pyoptics.holography.autofocus(field, x, y, z_range, wavelength, metric='tamura', object_type='amplitude', n=1.0, propagator='asm', pad_factor=2, mask=None, refine=True)

Find the best focus of a reconstructed field.

The field is refocused to all distances in z_range and the focus metric is evaluated in each plane.

Parameters:

fieldarray

Field in the hologram plane.

x, yarrays

z_rangesequence of numbers

Candidate refocus distances (e.g. np.linspace(-d_max, -d_min, 41)).

wavelengthnumber

metricstr or callable, optional

One of focus_metrics, or a function of the field returning a number.

object_type“amplitude” or “phase”, optional

Amplitude objects are in focus at the maximum of the metric, phase objects at its minimum.

n, propagator, pad_factor

As for refocus().

maskbool array, optional

Region in which the metric is evaluated.

refinebool, optional

If True, refine the best distance by a parabola through the best sample and its neighbours (for equidistant z_range).

Returns:

z_bestfloat

Best refocus distance.

valuesarray

Metric for each entry of z_range.

pyoptics.holography.remove_tilt(phase, x, y, mask=None, piston=True)

Subtract a least squares plane (tilt and optionally piston) from a phase map.

Parameters:

phasearray

Unwrapped phase.

x, yarrays

maskbool array, optional

Pixels used for the fit (e.g. the background around a sample).

pistonbool, optional

If True, also subtract the mean (over mask).

Returns:

phase_correctedarray

coeffsarray

(piston, a_x, a_y) of the fitted plane piston + a_x*x + a_y*y.

pyoptics.holography.quantitative_phase(field, x=None, y=None, reference_field=None, unwrap='quality', mask=None, tilt_mask=None, weights=None)

Quantitative phase map of a reconstructed (and refocused) field.

Parameters:

fieldarray

Complex object field.

x, yarrays, optional

Coordinates, needed for tilt removal (tilt_mask).

reference_fieldarray, optional

Field reconstructed without sample (or a fitted background). The phase of field*conj(reference_field) is free of the aberrations and the residual carrier tilt common to both (“double exposure” or background compensation).

unwrap“quality”, “ls”, None or callable, optional

Phase unwrapping: “quality” (default) uses pyoptics.unwrap.unwrap_quality_guided (with the amplitude as quality if no weights are given), which stays congruent with the wrapped phase; “ls” uses pyoptics.unwrap.unwrap_ls (weighted with weights if given), which is robust against residues but smooths steep phases. None returns the wrapped phase, a callable is called with the wrapped phase.

maskbool array, optional

Valid pixels, passed to the unwrapper.

tilt_maskbool array, optional

If given, remove_tilt() is applied with this mask (e.g. the sample free background), which also fixes the piston there.

weightsarray, optional

Weights or quality map for the unwrapper.

Returns:

phasearray

Phase in radians (unwrapped unless unwrap is None).

pyoptics.holography.optical_path_difference(phase, wavelength)

Optical path difference phase/k_0.

The thickness of a sample with index n_s in a medium n_m is OPD/(n_s - n_m).

Parameters:

phasenumber or array

Phase in radians.

wavelengthnumber

Vacuum wavelength.

Returns:

opdnumber or array

In the unit of the wavelength.