pyoptics.interferometry

Interferometry: simulation of interferograms, phase-shifting interferometry (PSI) and Fourier-transform fringe analysis (Takeda).

Conventions

Fields vary as exp(i(k.r - omega*t)) (see the README). The test beam carries the phase phi = k*W of its wavefront error (or optical path difference) W, so a wavefront that runs ahead has phi > 0. Superposing it with a reference beam E_ref gives the two-beam interferogram

I = |E_test + E_ref*exp(i*delta)|**2
  = I_0*(1 + gamma*cos(phi - delta)),

where phi = arg(E_test*conj(E_ref)) is the phase of the test beam relative to the reference and gamma the visibility (fringe contrast). A phase shift delta > 0 delays the reference by delta (its optical path gets longer by delta/k), so the fringes move towards larger phi. All algorithms here return phi, i.e. the phase of the test beam relative to the reference, with the sign convention of pyoptics: phase = +k*OPD. Note that many textbooks write cos(phi + delta) instead; their formulas then give -phi.

In reflection (Twyman-Green, Fizeau), the light passes the surface twice, so a surface height h gives phi = 2*k*h.

A tilt between the beams adds a linear carrier phase k_c.r = k_cx*x + k_cy*y (given here as a wavenumber in rad per length unit, like the k of utils.freq_grid()).

Image stacks follow the convention of pyoptics.processing: the last index numbers the images, i.e. frames[:, :, j] is the j-th phase-shifted frame (any number of leading dimensions is fine). Stacks can thus be passed through processing.Scaler and processing.CameraNoiser directly.

The phase-shifting algorithms return the wrapped phase in (-pi, pi] and the modulation gamma = b/a, where I = a + b*cos(phi - delta), i.e. the local fringe visibility (between 0 and 1 for physical intensities). Pixels without light (a = 0) get gamma = 0. Unwrap the phase with pyoptics.unwrap.

Contents

Simulation:

• interferogram(): interferogram (or stack) from a phase map

• interfere(): interferogram (or stack) from two complex fields

• zernike_phase(): phase map from Fringe Zernike coefficients in waves

• phase_shifts(): nominal or miscalibrated phase shifts

Phase-shifting algorithms:

• psi_n_step(): least squares for N frames with arbitrary known shifts

• psi_4_step(): four-bucket algorithm (shifts 0, pi/2, pi, 3*pi/2)

• psi_carre(): Carré, four frames with an unknown but constant step

• psi_hariharan(): Hariharan’s five-frame algorithm, insensitive to linear shift errors

• shift_error_sensitivity(): rms phase error vs. phase shifter miscalibration

Single-frame analysis:

• fourier_fringe_analysis(): Takeda’s Fourier-transform method for interferograms with a carrier

• find_carrier(): locate the carrier sideband in the spectrum

Spherical surfaces (radius of curvature tests):

• aplanatic_directions(): ray directions of a converging test wave

• sphere_displacement_phase(): phase of a displaced test sphere

• fit_sphere_displacement(): displacement of a test sphere from the null position (confocal or cat’s eye) from its measured phase

A sphere tested with a converging spherical wave (Fizeau transmission sphere, Twyman-Green with a diverger lens) gives a null interferogram in two positions: confocal, with its centre of curvature in the focus, and cat’s eye, with its vertex in the focus. The radius of curvature is the distance between the two, R = z_cat - z_confocal (signed as in pyoptics.raytrace, R > 0 for a centre of curvature in +z, i.e. downstream). Near either position, a small displacement d of the sphere changes the optical path of the ray with the direction n by d.n per pass, i.e. the phase by

phi = passes*k*(d_x*n_x + d_y*n_y + d_z*n_z),

tilt fringes for a lateral and power (defocus) fringes for an axial displacement. Fitting this model to the measured phase gives the null positions from slightly defocused interferograms. In cat’s eye, lateral displacements cancel to first order (the beam returns on the opposite side of the pupil).

References

J. H. Bruning, D. R. Herriott, J. E. Gallagher, D. P. Rosenfeld, A. D. White and D. J. Brangaccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses”, Appl. Opt. 13, 2693-2703 (1974).

J. E. Greivenkamp, “Generalized data reduction for heterodyne interferometry”, Opt. Eng. 23, 350-352 (1984).

P. Carré, “Installation et utilisation du comparateur photoélectrique et interférentiel du Bureau International des Poids et Mesures”, Metrologia 2, 13-23 (1966).

P. Hariharan, B. F. Oreb and T. Eiju, “Digital phase-shifting interferometry: a simple error-compensating phase calculation algorithm”, Appl. Opt. 26, 2504-2506 (1987).

M. Takeda, H. Ina and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry”, J. Opt. Soc. Am. 72, 156-160 (1982).

K. Creath, “Phase-measurement interferometry techniques”, Prog. Opt. 26, 349-393 (1988).

4. Malacara (ed.), “Optical Shop Testing”, 3rd ed., Wiley (2007).

Functions

aplanatic_directions	Directions of the rays of a converging spherical wave from an aplanatic transmission sphere (sine condition: sin(theta) = NA*rho).
find_carrier	Locate the carrier, i.e. the +1 sideband of a fringe pattern.
fit_sphere_displacement	Least squares fit of the displacement of a test sphere from its null position to its unwrapped phase.
fourier_fringe_analysis	Takeda's Fourier-transform method (Takeda, Ina, Kobayashi 1982) for a single interferogram with a carrier.
interfere	Interferogram |E_test + E_ref*exp(i*shift)|**2 of two fields.
interferogram	Two-beam interferogram I = I0*(1 + visibility*cos(phase + k_c.r - shift)) of a test phase.
phase_shifts	Phase shifts of a phase-shifting sequence, optionally with a miscalibrated phase shifter.
psi_4_step	Four-bucket algorithm for the shifts 0, pi/2, pi and 3*pi/2.
psi_carre	Carré's algorithm: four frames with the shifts -3*alpha/2, -alpha/2, alpha/2, 3*alpha/2 and an unknown step alpha.
psi_hariharan	Hariharan's five-frame algorithm (Hariharan, Oreb, Eiju 1987) for the shifts -2*alpha, -alpha, 0, alpha, 2*alpha with alpha = pi/2.
psi_n_step	Least-squares phase-shifting algorithm for N frames with known shifts.
shift_error_sensitivity	Rms phase error of a phase-shifting algorithm vs. the calibration error of the phase shifter.
sphere_displacement_phase	Phase passes*k*(d.n) of a test sphere displaced by d from its null position (first order in d), see the module docstring.
zernike_phase	Phase map 2*pi*passes*W of a wavefront W given by Fringe Zernike coefficients in waves.

pyoptics.interferometry.zernike_phase(x, y, coeffs, R_norm, indices=None, x0=0.0, y0=0.0, passes=1)

Phase map 2*pi*passes*W of a wavefront W given by Fringe Zernike coefficients in waves.

Parameters:

x, yarrays

Linear coordinates of the grid.

coeffssequence of floats

Fringe Zernike coefficients in waves (units of the wavelength).

R_normfloat

Radius of the pupil (normalization radius of the Zernikes).

indicessequence of ints, optional

Fringe indices that belong to coeffs. Defaults to 1, 2, …, len(coeffs).

x0, y0floats, optional

Center of the pupil.

passesint, optional

Number of passes through the test object, e.g. 2 for a surface measured in reflection (the coefficients then describe the surface height in waves).

Returns:

phasearray

Phase in radians, zero outside the pupil.

maskarray of bool

Pupil mask.

pyoptics.interferometry.phase_shifts(N, step=None, start=None, linear_error=0.0, quadratic_error=0.0)

Phase shifts of a phase-shifting sequence, optionally with a miscalibrated phase shifter.

The nominal shifts are delta_j = start + j*step, j = 0, …, N - 1. A phase shifter with the relative linear calibration error epsilon_1 and a quadratic nonlinearity epsilon_2 produces

delta_j' = delta_j*(1 + epsilon_1) + epsilon_2*delta_j**2/pi,

the error being taken relative to the shifter’s zero position (delta = 0).

Parameters:

Nint

Number of frames.

stepfloat, optional

Nominal phase step in radians. Defaults to 2*pi/N.

startfloat, optional

Nominal first shift. Defaults to 0, i.e. the sequence 0, step, … Use start=-(N - 1)*step/2 for a sequence centered at 0 (as for Carré’s and Hariharan’s algorithms).

linear_errorfloat, optional

Relative linear calibration error epsilon_1, e.g. 0.05 for a shifter that moves 5 % too far.

quadratic_errorfloat, optional

Quadratic nonlinearity epsilon_2 (relative error of a shift of pi).

Returns:

shiftsarray of length N

pyoptics.interferometry.interferogram(phase, x=None, y=None, visibility=1.0, carrier=(0.0, 0.0), shift=0.0, I0=1.0, noiser=None)

Two-beam interferogram I = I0*(1 + visibility*cos(phase + k_c.r - shift)) of a test phase.

Parameters:

phasearray

Phase of the test beam relative to the reference in radians (e.g. from zernike_phase()).

x, yarrays, optional

Linear coordinates of the grid; only needed for a carrier.

visibilityfloat or array, optional

Fringe visibility gamma (may vary across the field).

carriertuple of floats, optional

Carrier wavenumbers (k_cx, k_cy) in rad per length unit of a tilt between the beams, e.g. (2*pi/period, 0) for vertical fringes.

shiftfloat or sequence of floats, optional

Phase shift(s) of the reference in radians. For a sequence, the result is a stack with the frame index last.

I0float or array, optional

Mean intensity (bias).

noisercallable, optional

Applied to the result, e.g. a processing.CameraNoiser (in that case, I0 is the mean number of photons per pixel).

Returns:

Iarray

Interferogram, or stack of interferograms of shape phase.shape + (len(shift),).

pyoptics.interferometry.interfere(E_test, E_ref, shift=0.0)

Interferogram |E_test + E_ref*exp(i*shift)|**2 of two fields.

The algorithms of this module recover arg(E_test*conj(E_ref)) from it.

Parameters:

E_test, E_refarrays (or scalar for E_ref)

Complex fields of the test and reference beams.

shiftfloat or sequence of floats, optional

Phase shift(s) of the reference beam. For a sequence, the result is a stack with the frame index last.

Returns:

Iarray

pyoptics.interferometry.psi_n_step(frames, shifts=None, return_bias=False)

Least-squares phase-shifting algorithm for N frames with known shifts.

Fits I_j = a + c*cos(delta_j) + s*sin(delta_j) per pixel (Bruning et al. 1974, Greivenkamp 1984), which gives phi = atan2(s, c) and b = sqrt(c**2 + s**2). For equally spaced shifts over one period, this is the classical N-bucket (synchronous detection) algorithm.

Parameters:

framesarray

Image stack, frame index last. At least 3 frames.

shiftssequence of floats, optional

Phase shifts of the frames in radians. Defaults to 2*pi*j/N.

return_biasbool, optional

If True, also return the bias a.

Returns:

phasearray

Wrapped phase phi in (-pi, pi].

modulationarray

gamma = b/a.

biasarray

Only if return_bias: the mean intensity a.

pyoptics.interferometry.psi_4_step(frames)

Four-bucket algorithm for the shifts 0, pi/2, pi and 3*pi/2.

The phase is

phi = atan2(I_2 - I_4, I_1 - I_3).

Parameters:

framesarray

Stack of 4 frames, frame index last.

Returns:

phase, modulationarrays

Wrapped phase and modulation gamma = b/a.

pyoptics.interferometry.psi_carre(frames, return_step=False)

Carré’s algorithm: four frames with the shifts -3*alpha/2, -alpha/2, alpha/2, 3*alpha/2 and an unknown step alpha.

With D_1 = I_1 - I_4 and D_2 = I_2 - I_3,

tan(alpha/2)**2 = (3*D_2 - D_1)/(D_1 + D_2),
tan(phi) = sgn(-D_2)*sqrt(|(3*D_2 - D_1)*(D_1 + D_2)|)
           / (I_2 + I_3 - I_1 - I_4).

The step is estimated per pixel; the modulation uses the median step. The estimate of alpha (and with it the phase) is ill-conditioned where sin(phi) or cos(phi) is close to 0, so noise hurts more than with algorithms for known shifts. Works for 0 < alpha < pi, best near alpha = 110 degrees.

Parameters:

framesarray

Stack of 4 frames, frame index last.

return_stepbool, optional

If True, also return the median of the estimated step alpha.

Returns:

phase, modulationarrays

Wrapped phase and modulation gamma = b/a.

alphafloat

Only if return_step: the estimated phase step in radians.

pyoptics.interferometry.psi_hariharan(frames, return_step=False)

Hariharan’s five-frame algorithm (Hariharan, Oreb, Eiju 1987) for the shifts -2*alpha, -alpha, 0, alpha, 2*alpha with alpha = pi/2:

phi = atan2(2*(I_4 - I_2), 2*I_3 - I_1 - I_5).

A linear miscalibration alpha = pi/2*(1 + epsilon) only causes a phase error of second order in epsilon (compared to first order for the four bucket algorithm).

Parameters:

framesarray

Stack of 5 frames, frame index last.

return_stepbool, optional

If True, also return the median of the step estimated from cos(alpha) = (I_5 - I_1)/(2*(I_4 - I_2)) (Hariharan’s calibration formula).

Returns:

phase, modulationarrays

Wrapped phase and modulation gamma = b/a.

alphafloat

Only if return_step: the estimated phase step in radians.

pyoptics.interferometry.shift_error_sensitivity(algorithm, linear_errors, quadratic_error=0.0, N=None, step=None, N_phase=360)

Rms phase error of a phase-shifting algorithm vs. the calibration error of the phase shifter.

Frames with miscalibrated shifts (see phase_shifts()) of noise-free fringes with visibility 1 are evaluated for N_phase test phases evenly spread over one period; the rms of the wrapped phase error is returned for each linear error.

Parameters:

algorithmstr or callable

“4-step”, “carre”, “hariharan”, “n-step” (least squares with the nominal shifts, use N and step) or a callable f(frames) -> (phase, modulation) (then N and step are needed).

linear_errorssequence of floats

Relative linear calibration errors epsilon_1.

quadratic_errorfloat, optional

Quadratic nonlinearity epsilon_2 (see phase_shifts()).

Nint, optional

Number of frames (“n-step” and callables).

stepfloat, optional

Nominal phase step; defaults to 2*pi/N for “n-step” and callables.

N_phaseint, optional

Number of test phases.

Returns:

rms_errorarray

Rms phase error in radians for each linear error.

pyoptics.interferometry.find_carrier(I, x, y, exclude=None, mask=None)

Locate the carrier, i.e. the +1 sideband of a fringe pattern.

Searches the spectral peak in the half plane k_x > 0 (or k_x = 0, k_y > 0) outside a disk around k = 0 and returns the centroid of the power spectrum |S|**2 in a disk of radius |k|/2 around the current estimate k, starting at the peak and iterated until it settles. For a test phase without net tilt (e.g. defocus or astigmatism, which spread the sideband symmetrically), the centroid is the carrier; a tilt of the test phase is included in it.

Parameters:

Iarray

Interferogram.

x, yarrays

Linear coordinates.

excludefloat, optional

Radius around k = 0 (in rad per length unit) that is ignored. Defaults to 3 frequency samples.

maskarray of bool, optional

Region with fringes (e.g. the pupil). Without it, the edge of a pupil can outshine the sideband.

Returns:

carriertuple of floats

(k_cx, k_cy).

pyoptics.interferometry.fourier_fringe_analysis(I, x, y, carrier=None, window='tophat', radius=None, mask=None, apodize=0.0, remove_carrier=True, return_sideband=False)

Takeda’s Fourier-transform method (Takeda, Ina, Kobayashi 1982) for a single interferogram with a carrier.

The interferogram I = a + b*cos(phi + k_c.r) = a + c*exp(i k_c.r) + c.c., with c = b/2*exp(i*phi), has a sideband around k_c in its spectrum. A window around k_c isolates it; its inverse transform, multiplied with exp(-i k_c.r), gives c and hence phi and b. The bias a is the low-pass filtered interferogram (same window around k = 0).

The method needs the sideband to be separated from the zero order and the -1 order: the local fringe frequency |grad phi| must stay below the window radius (by default |k_c|/2). A warning is issued if the window holds less than 90 % of the power in the half plane of the carrier, or if it reaches beyond the Nyquist frequency.

Parameters:

Iarray

Interferogram.

x, yarrays

Linear coordinates.

carriertuple of floats, optional

Carrier (k_cx, k_cy) in rad per length unit. Found with find_carrier() if not given; then the tilt of the test phase is removed along with the carrier (remove a residual tilt with a fit if needed).

windowstr or array, optional

“tophat” (default), “hann” or “gaussian” window around k_c, or a window array on the centered frequency grid of I. The tapered windows ring less but attenuate the outer parts of a broad sideband, which biases the modulation low; the bias is filtered with the same window around k = 0 (for a window array, the global mean is used).

radiusfloat, optional

Window radius in rad per length unit. Defaults to |k_c|/2.

maskarray of bool, optional

Region with fringes (e.g. the pupil). Outside, I is replaced by its mean inside before the transform, which reduces ringing at the edge, and the outputs are set to 0.

apodizefloat, optional

Taper fraction of a separable Tukey window applied to I - mean(I) before the transform (0: no apodization, 1: Hann window). It reduces the ringing that the field edges cause in the interior, at the cost of larger errors in the tapered border. The modulation is corrected for it.

remove_carrierbool, optional

If False, the carrier phase k_c.r is kept in the result.

return_sidebandbool, optional

If True, also return the complex sideband c (after carrier removal).

Returns:

phasearray

Wrapped phase in (-pi, pi].

modulationarray

gamma = b/a.

sidebandarray

Only if return_sideband.

pyoptics.interferometry.aplanatic_directions(x, y, R_pupil, NA, x0=0.0, y0=0.0)

Directions of the rays of a converging spherical wave from an aplanatic transmission sphere (sine condition: sin(theta) = NA*rho).

The pupil point (x, y) (rho = r/R_pupil) belongs to the ray that converges towards the focus on the axis with the angle theta, i.e. it travels in +z and towards the axis.

Parameters:

x, yarrays

Linear coordinates of the grid (e.g. of the detector, on which the pupil is imaged).

R_pupilfloat

Radius of the pupil on the grid.

NAfloat

Numerical aperture of the test wave.

x0, y0floats, optional

Centre of the pupil.

Returns:

directionsarray of shape (len(y), len(x), 3)

Unit vectors (n_x, n_y, n_z), NaN outside the pupil.

pyoptics.interferometry.sphere_displacement_phase(directions, displacement, wavelength, passes=2)

Phase passes*k*(d.n) of a test sphere displaced by d from its null position (first order in d), see the module docstring.

Parameters:

directionsarray of shape (…, 3)

Unit directions n of the test wave (e.g. from aplanatic_directions() or raytrace.System.ray_map()).

displacement(d_x, d_y, d_z)

Displacement of the sphere in length units.

wavelengthfloat

passesint, optional

2 for a surface in reflection.

Returns:

phasearray

pyoptics.interferometry.fit_sphere_displacement(phase, directions, wavelength, mask=None, lateral=True, passes=2)

Least squares fit of the displacement of a test sphere from its null position to its unwrapped phase.

Fits phi = c_0 + passes*k*(d.n) (see sphere_displacement_phase()). The axial component d_z from the confocal and the cat’s eye interferograms gives the null positions and the radius of curvature: with the stage positions z_1 (near confocal) and z_2 (near cat’s eye), R = (z_2 - d_z2) - (z_1 - d_z1).

Parameters:

phasearray

Unwrapped phase of the test beam relative to the reference.

directionsarray of shape phase.shape + (3,)

Unit directions of the test wave at each pixel.

wavelengthfloat

maskarray of bool, optional

Pixels to use (default: all with finite phase and directions).

lateralbool, optional

Fit the lateral displacement (d_x, d_y) as well. In cat’s eye, set it to False and fit tilt terms instead (lateral displacements do not show there; a tilt of the reference does).

passesint, optional

Returns:

displacementarray (d_x, d_y, d_z)

(d_x, d_y are 0 for lateral=False.)

residualarray

phase minus the fit (0 outside the mask).