pyoptics.shack_hartmann

Shack-Hartmann wavefront sensor: spot images, centroids and reconstruction.

A Shack-Hartmann sensor is a square array of n x n lenslets of pitch p and focal length f in front of a camera in their focal plane. Each lenslet focuses its part of the beam to a spot; a local tilt of the wavefront moves the spot by f times the mean slope of the wavefront over the lenslet. From the spot displacements, the wavefront is reconstructed either modally (least squares fit of the gradients of Fringe Zernikes) or zonally (Southwell’s least squares integration on the lenslet grid).

ShackHartmann covers the whole chain:

Noise is not part of the model; apply processing.CameraNoiser to the image, scaled to photon numbers.

Conventions

Lengths (pitch, focal length, wavelength, pixel size, pupil radius) share one unit. The lenslet array is centred on the optical axis; lenslet (i, j) is row i (y) and column j (x), as for images. The beam footprint on the lenslets is a disk of radius pupil_radius, which maps to the unit disk of normalized pupil coordinates (P_x, P_y) = (x, y)/pupil_radius. Phases are in radians, with the field exp(i*phase); the OPD in length units is W = phase*wavelength/(2*pi), as in raytrace. A wavefront with a positive slope dW/dx (the beam turns towards +x) moves the spots by +f*dW/dx. Slopes are dimensionless (dW/dx in length per length). Fringe Zernike coefficients are in waves on the unit pupil, as in phase_retrieval and raytrace.

The centroid of the Fraunhofer pattern of a lenslet equals f times the intensity weighted mean gradient of W over the lenslet (exactly, for an infinite detector). zernike_matrix() uses this as the model, so partially illuminated lenslets at the rim are described correctly. Deviations come from the finite camera window of each lenslet, the threshold and the discrete pixels.

Limitations: no crosstalk between lenslets (each spot is computed in the window of its own lenslet only; light leaving the window is lost, which sets the dynamic range), scalar Fraunhofer diffraction with ideal thin lenslets, square lenslets without dead zones.

References

B. C. Platt and R. Shack, “History and principles of Shack-Hartmann wavefront sensing”, J. Refract. Surg. 17, S573-S577 (2001).

W. H. Southwell, “Wave-front estimation from wave-front slope measurements”, J. Opt. Soc. Am. 70, 998-1006 (1980).

R. J. Noll, “Zernike polynomials and atmospheric turbulence”, J. Opt. Soc. Am. 66, 207-211 (1976).

Classes

ShackHartmann

Shack-Hartmann sensor with n x n square lenslets.

Slopes

Measured wavefront slopes of a Shack-Hartmann sensor.
class pyoptics.shack_hartmann.Slopes(sx, sy, valid, flux)

Bases: tuple

Measured wavefront slopes of a Shack-Hartmann sensor.

sx, sy: slopes dW/dx and dW/dy per lenslet, shape (n, n), NaN where not valid. valid: lenslets with enough light. flux: light per lenslet window, relative to the brightest one.

flux

Alias for field number 3

sx

Alias for field number 0

sy

Alias for field number 1

valid

Alias for field number 2

class pyoptics.shack_hartmann.ShackHartmann(n_lenslets, pitch, focal_length, wavelength, pixels=16, samples=None, pupil_radius=None)[source]

Bases: object

Shack-Hartmann sensor with n x n square lenslets.

Parameters:
n_lensletsint

Number of lenslets across the array.

pitchnumber

Lenslet pitch (width of a lenslet).

focal_lengthnumber

Focal length of the lenslets; the camera sits in their focal plane.

wavelengthnumber

Wavelength.

pixelsint, optional

Camera pixels per lenslet across (default 16). The pixel size is pitch/pixels.

samplesint, optional

Samples of the field per lenslet across. The default is the smallest even number >= max(16, 2*N_F) with the lenslet Fresnel number N_F = pitch**2/(wavelength*focal_length), so that the Fourier sum of a lenslet does not alias over the camera window.

pupil_radiusnumber, optional

Radius of the beam footprint on the array, default n_lenslets*pitch/2 (the circle inscribed in the array).

Attributes:
x, yarrays

Field sample coordinates on the array (length n_lenslets*samples).

P_x, P_yarrays

Normalized pupil coordinates of the field samples, shape (N, N).

centersarray

Lenslet centre coordinates along x and y (length n_lenslets).

uarray

Camera pixel coordinates relative to the lenslet centre (length pixels).

property fresnel_number

Fresnel number pitch**2/(wavelength*focal_length) of a lenslet.

property max_slope

Largest slope that keeps a spot centre in its lenslet window.

property slope_per_pixel

Slope that moves a spot by one camera pixel.

amplitude(amplitude=None)[source]

Field amplitude on the samples; default: the unit disk.

field(phase=None, amplitude=None)[source]

Complex field on the lenslet array.

Parameters:
phasecallable, array or number, optional

Pupil phase in radians: phase(P_x, P_y) on normalized pupil coordinates (e.g. raytrace.Pupil.phase), or an array on the sample grid of shape (N, N). Default: a plane wave.

amplitudecallable, array or number, optional

Field amplitude, as phase (e.g. raytrace.Pupil.mask). Default: 1 on the unit disk, 0 outside.

Returns:
Earray

Field of shape (N, N) with N = n_lenslets*samples.

spots(field)[source]

Camera image of the spots.

Each lenslet window of pixels x pixels camera pixels receives the Fraunhofer pattern of the field over that lenslet,

E(u, v) ~ sum_xy E(x, y)*exp(-i*k*(u*x + v*y)/f),

with the local coordinates x, y and u, v relative to the lenslet centre.

Parameters:
fieldarray

Field from field(), shape (N, N).

Returns:
imagearray

Power per pixel, shape (n_lenslets*pixels,)*2, normalized to the total power of the field (a sum of 1 if no light leaves the windows).

centroids(image, threshold=0.1)[source]

Spot centroids and flux per lenslet window.

Parameters:
imagearray

Camera image from spots() (possibly with noise).

thresholdnumber, optional

Fraction of the peak of each window that is subtracted before the centre of gravity is computed (values below are set to zero). This suppresses the noise of the dark pixels.

Returns:
cx, cyarrays

Centroids relative to the lenslet centres, shape (n, n), NaN for empty windows.

fluxarray

Total of each window, shape (n, n).

slopes(image, reference=None, threshold=0.1, min_flux=0.5)[source]

Wavefront slopes from a camera image.

Parameters:
imagearray

Camera image from spots().

referencearray, optional

Camera image of a reference (plane) wave, e.g. from a calibration; its centroids are subtracted. Default: the lenslet centres.

thresholdnumber, optional

See centroids().

min_fluxnumber, optional

Lenslets whose window holds less than min_flux times the flux of the brightest window are not valid (partially illuminated lenslets at the rim).

Returns:
Slopes
model_slopes(phase=None, amplitude=None)[source]

Slopes that an ideal sensor measures.

The intensity weighted mean phase gradient over each lenslet, converted to dW/dx and dW/dy. The gradient is taken by central differences on the sample grid, so the phase must be continuous (unwrapped). Samples whose difference stencil reaches a dark sample are left out, so that a phase that is only defined on the pupil (e.g. zero outside) does not produce spurious slopes at the rim.

Parameters:
phase, amplitudecallable, array or number, optional

As for field().

Returns:
sx, syarrays

Shape (n, n); NaN for dark lenslets.

zernike_matrix(indices, amplitude=None)[source]

Slopes of Fringe Zernikes of one wave each.

Parameters:
indicessequence of int

Fringe indices j.

amplitudecallable, array or number, optional

Pupil amplitude (weights of the gradients), as for field().

Returns:
Gx, Gyarrays

Shape (len(indices), n, n): model slopes dW/dx, dW/dy of 1 wave of Z_j (see model_slopes()).

reconstruct_zernikes(slopes, n_zernike=15, amplitude=None, indices=None)[source]

Modal reconstruction: least squares fit of Zernike gradients.

Parameters:
slopesSlopes

Measured slopes from slopes().

n_zernikeint, optional

Fit the Fringe Zernikes 2, …, n_zernike (piston is invisible).

amplitudecallable, array or number, optional

Pupil amplitude for the model, as for field().

indicessequence of int, optional

Fringe indices to fit instead of 2…n_zernike.

Returns:
coeffsdict

Fringe Zernike coefficients {j: c_j} in waves.

reconstruct_zonal(slopes)[source]

Zonal reconstruction (Southwell 1980) on the lenslet grid.

Solves (W_{i,j+1} - W_{i,j})/pitch = (s_x[i,j+1] + s_x[i,j])/2 and the same along y for all pairs of neighbouring valid lenslets in the least squares sense (the mean of W is zero).

Parameters:
slopesSlopes

Measured slopes from slopes().

Returns:
Warray

Wavefront at the lenslet centres in waves, shape (n, n), NaN where not valid.