Zernike Holography

Optical aberration is a pervasive roadblock to good imaging. Aberration can be modeled as a phase distortion upon an optical wavefront: while a flat wavefront is free of aberration, non-ideal optical elements impart a sort of distortion hologram on the beam leading to distorted spots. Zernike polynomials act as a good (orthonormal under some conditions) basis to describe optical aberration, with terms corresponding directly to common forms of aberration.

We have already investigated simple wavefront calibration in a previous example. The action of this calibration is (in part) to measure the phase distortion and then compensate for the known distortion by subtracting the phase from the phase of desired holograms. However, a major problem with this technique is that the measured phase is only valid for the \(k\)-vector (corresponding to a point on the camera) at which the calibration was measured. Other \(k\)-vectors (points across the field of view of the camera) experience a different optical path through the optical train and thus experience a different aberration. In systems with strong distortion, this can lead to good sharp wavefront correction only in the neighborhood surrounding the calibration point such as seen in slmsuite issue #30.

This phenomenon is an invariant for additive correction; one cannot find a single additive correction mask to correct all \(k\)-vectors for a given hologram. The solution is to instead “bake” the correction into the hologram itself. In the case of spot holograms, instead of thinking of mapping power to points in 2D space, it is instead more practical to think of mapping power to points in Zernike aberration space, where each spot has a unique Zernike calibration.

Initialization

To start, we initialize our system:

# Try to load physical hardware.
try:
    from slmsuite.hardware.slms.texasinstruments import PLM
    from slmsuite.hardware.cameras.flir import FLIR

    # Load the hardware.
    slm = PLM('p67', 1, wav_um=0.488, settle_time_s=0.05, configure_usb=True); print()
    cam = FLIR(serial="22562470", pitch_um=2.74)

    # For this camera, we also want to crop the WOI to the center 1024x1024 area.
    h_max, w_max = cam.default_shape
    cam.set_woi(((w_max - 1024) // 2, 1024, (h_max - 1024) // 2, 1024))

    # Load the composite CameraSLM.
    fs = FourierSLM(cam, slm)
    fs.load_calibration("wavefront_superpixel");
    fs.wavefront_calibration_superpixel_process(plot=True, r2_threshold=.8)
# Fallback to virtual hardware.
except Exception as e:
    print(f"Loading physical hardware failed; loading virtual hardware. Error:\n{str(e)}")
    from slmsuite.hardware.slms.simulated import SimulatedSLM
    from slmsuite.hardware.cameras.simulated import SimulatedCamera

    # Make the SLM and camera.
    slm = SimulatedSLM((1600, 1200), pitch_um=(8,8))
    # Program the simulated and calibrated sources to be flat phase.
    slm.set_source_analytic(sim=True)
    slm.set_source_analytic()

    cam = SimulatedCamera(slm, (1440, 1100), pitch_um=(4,4), gain=200)

    # Tie the camera and SLM together with an analytic calibration.
    fs = FourierSLM(cam, slm)
    M, b = fs.fourier_calibration_build(
        f_eff=80000.,                              # f_eff of 80000 wavelengths or 80 mm with the default 1 um wavelength
        theta=5 * np.pi / 180,
    )
    cam.set_affine(M, b)
    fs.fourier_calibrate_analytic(M, b)

PLM using GPU (cupy) backend
DLPC900 connected: firmware {'app_version': 95, 'api_version': 72, 'sw_patch': 111, 'sw_minor': 108, 'sw_major': 111}
DLPC900 pre-configured (video mode, display detected)
Initializing pyglet... success
Searching for window with display_number=1... success
Creating window... Window creation successful
DLPC900 source locked, switching to video-pattern mode...
DLPC900 configured successfully - pattern sequence running

PySpin initializing... Looking for cameras... PySpin sn '22562470' initializing... PixelFormat set to Mono16 (12-bit)... Frame rate set to 33.7 Hz... Successfully initialized FLIR cam 22562470.

C:\Users\holodyne\Documents\GitHub\slmsuite\slmsuite\hardware\cameraslms.py:3720: UserWarning: remove_background is enabled and a noise floor was detected; removing this background (0.07164789736270905% of the average normalization).
  warnings.warn(

cam.set_exposure(.002)
fs.fourier_calibrate(20, 20, plot=True);

C:\Users\holodyne\Documents\GitHub\slmsuite\slmsuite\holography\toolbox\__init__.py:258: UserWarning: CameraSLM must be passed as slm for conversion 'kxy' to 'ij'
  warnings.warn(

Zernike Indexing

There are many ways to index Zernike polynomials. slmsuite supports a number of conventions, but uses the 1D ANSI "ansi" indexing \(i\) as the default. We can find the mapping \(Z_i = Z_n^l\) between ANSI and the canonical definition "radial" consisting of the radial \(n\) and azimuthal \(l\) orders:

from slmsuite.holography.toolbox.phase import zernike_convert_index

ansi = np.arange(10)
radial = zernike_convert_index(ansi, from_index="ansi", to_index="radial")
print(radial)

[[ 0  0]
 [ 1 -1]
 [ 1  1]
 [ 2 -2]
 [ 2  0]
 [ 2  2]
 [ 3 -3]
 [ 3 -1]
 [ 3  1]
 [ 3  3]]

The nature of the radial \(n\) and azimuthal \(l\) orders make it often convenient to arange the polynomials in a pyramid:

from slmsuite.holography.toolbox.phase import zernike_pyramid_plot

plt.figure(figsize=(20,20))
zernike_pyramid_plot(slm, order=3)

Notice how the polynomials are translated and scaled to correspond with the measured amplitude distribution of the SLM. This is because slmsuite is automatically computing the best center point and scaling radius for the polynomials based upon fitting to the measured amplitude distribution. This allows the Zernike basis to act accurately upon the spot, without laterally steering the beam while focusing, for instance.

slmsuite also supports derivates (computed analyticaly via the power rule before rendering, for speed):

plt.figure(figsize=(20,20))
zernike_pyramid_plot(slm, order=3, derivative=(1, 0), scale=5)    # Adjust the scale because the derivative exceeds [-1, 1]

We can do this for very high order \(n\):

from slmsuite.holography.toolbox.phase import zernike_pyramid_plot

plt.figure(figsize=(20,20))
zernike_pyramid_plot(slm, order=12, titles=["ansi"])

To actually grab the phase of a single polynomial, we can use zernike(). This is of the same style as other functions in toolbox.phase, where the SLM is passed as grid= to give the function information about the scaling of the polynomial.

from slmsuite.holography.toolbox.phase import zernike

zernke_42 = zernike(slm, 42)

slm.plot(zernke_42);

The core function zernike_sum() supports generating stacks of sums of Zernike polynomials, with derivatives are applied to each. The weights= term is here a (4, 6) matrix, 4 terms for the requested ANSI indices, 6 terms for the 6 produced sums. These weights \(w_i^{(j)}\) are multiplied directly with \(Z_i\) (normalized to \(\pm 1\)), to produce \(\phi_j = \sum_j w_i^{(j)}Z_i\).

from slmsuite.holography.toolbox.phase import zernike_sum

for derivative in [(0,0), (1,0), (0,1)]:
    images = np.pi * zernike_sum(
        slm,
        indices=[3,4,5,6],
        weights=[
            [1,0,0,1,0,1],      # 3
            [0,1,2,1,0,1],      # 4
            [0,0,0,0,1,1],      # 5
            [0,0,2,0,1,1],      # 6
        ],
        derivative=derivative
    )

    fig, axs = plt.subplots(1, 6, figsize=(20,3))

    for i, image in enumerate(images):
        plt.sca(axs[i])
        slm.plot(image, title=str(i), cbar=False)

    plt.suptitle(f"derivative = {derivative}")
    plt.tight_layout()
    plt.show()

3D “Compressed Spot” Holography

3D spots is a first simple application of the math which “bakes” aberration into a hologram. Indeed, the \(Z_4 = Z_2^0\) Zernike polynomial is exactly defocus, and it is this parabolic term which we will use for translating (aberrating) in \(z\). Our previous examples effectively only used the \(x = Z_2 = Z_1^1\) and \(y = Z_1 = Z_1^{-1}\) terms.

Let’s start by making a desired pointcloud in the camera’s pixel "ij" coordinates.

from slmsuite.holography.toolbox import convert_vector

N = 100
x = np.arange(N) / 20

zeroth_order_ij = convert_vector((0,0), from_units="norm", to_units="ij", hardware=fs)

vectors_ij = np.vstack((
    50 * (x + 1) * np.cos(2 * np.pi * x) + zeroth_order_ij[0],
    50 * (x + 1) * np.sin(2 * np.pi * x) + zeroth_order_ij[1],
    200 * np.cos(1.5 * np.pi * x)
))

Unit conversions are a major focus of slmsuite. The core method toolbox.convert_vector() described previously supports 3D points and converts between the various units. For most units, “focal power” \(1/f\) is used for \(z\), but with Fourier calibration and knowledge of the camera’s pixel pitch in microns, this extends to true cartesian vectors. If you don’t have the .pitch_um parameter of the camera set, you will probably see nan for all the \(z\) coordinates except for "ij".

from slmsuite.holography.toolbox import print_blaze_conversions

print_blaze_conversions(
    vectors_ij[:, 0],       # Print the first vector.
    from_units="ij",
    hardware=fs,
    shape=(2048, 2048)
)

'rad' : [-8.92277289e-04 -1.56471752e-05  1.16157864e-08]
'mrad' : [-8.92277289e-01 -1.56471752e-02  1.16157864e-08]
'deg' : [-5.11237228e-02 -8.96517102e-04  1.16157864e-08]
'norm' : [-8.92277289e-04 -1.56471752e-05  1.16157864e-08]
'kxy' : [-8.92277289e-04 -1.56471752e-05  1.16157864e-08]
'knm' : [9.83557898e+02 1.02329080e+03 1.16157864e-08]
'freq' : [-1.97471203e-02 -3.46289943e-04  1.16157864e-08]
'lpmm' : [-1.82843707e+00 -3.20638837e-02  1.16157864e-08]
'zernike' : [-74.49629301  -1.30638375   3.22165036]
'ij' : [631.84307109 414.44837327 200.        ]
'm' : [0.00173125 0.00113559 0.000548  ]
'cm' : [0.173125   0.11355885 0.0548    ]
'mm' : [1.73125001 1.13558854 0.548     ]
'um' : [1731.2500148  1135.58854277  548.        ]
'nm' : [1731250.01479967 1135588.54276882  548000.        ]
'mag_m' : [0.00173125 0.00113559 0.000548  ]
'mag_cm' : [0.173125   0.11355885 0.0548    ]
'mag_mm' : [1.73125001 1.13558854 0.548     ]
'mag_um' : [1731.2500148  1135.58854277  548.        ]
'mag_nm' : [1731250.01479967 1135588.54276882  548000.        ]

The "zernike" conversion is listed as the coefficients (in radians) of the \(Z_2, Z_1, Z_4\) terms, which correspond to \(x,y,z\), the tilt and focus terms (see the pyramid).

We can now make the hologram with these 3D points.

from slmsuite.holography.algorithms import CompressedSpotHologram

hologram = CompressedSpotHologram(
    spot_vectors=vectors_ij,
    basis="ij",                 # For 3D spots, the basis defaults to Zernike terms [2,1,4] = [x,y,z].
    cameraslm=fs
)
hologram.optimize()

In the future, this example will be updated to describe how this process for “compressed” holography actually works. In the meantime, let’s take a look at what we made:

hologram.plot_farfield(limit_padding=-.4);

Under the hood, slmsuite is using the cached conversion and convert_vector() to translate the "ij" coordinates into to "zernike" space and make a basis of “compressed kernels” which map to each spot in this three-dimensional aberration space. We can look at this basis of Zernike polynomial kernels (only plotting the first ten for speed) by feeding the vector in aberration space as weights= of zernike_sum().

# Get the coordinates in "Zernike space" using the [2,1,4] = [x,y,z] terms.
vectors_zernike = convert_vector(
    vectors_ij,
    from_units="ij",
    to_units="zernike",
    hardware=fs
)

# Plot only the first 10.
images = np.pi * zernike_sum(
    slm,
    indices=None,       # For 3D spots, the basis defaults to Zernike terms [2,1,4].
    weights=vectors_zernike[:, :10]
)

fig, axs = plt.subplots(2, 5, figsize=(20,6))
axs = axs.ravel()

for i, image in enumerate(images):
    plt.sca(axs[i])
    slm.plot(image, title=str(i), cbar=False)

plt.tight_layout()
plt.show()

Zernike Spot Holography

However, this technique is extensible to far more than the three \(Z_2, Z_1, Z_4\) terms. Let’s make a realized Zernike pyramid in a single hologram.

First, let’s find the "zernike" coordinates of the pyramid corresponding to desired "ij" positions for arrangement according to the radial and azimuthal orders \(n\) and \(l\).

from slmsuite.holography.toolbox import format_vectors

center = format_vectors((200, 500))
cx, cy = center.ravel()

res = (320, 320)
((i0, i1), (i2, i3)) = limits = (
    (int(cx - res[0]/2), int(cx + res[0]/2)),
    (int(cy - res[1]/2), int(cy + res[1]/2))
)

N = 28
i = np.arange(N)
nl = zernike_convert_index(i, from_index="ansi", to_index="radial")
ij = 20 * np.vstack((nl[:, 1], (3-nl[:, 0]) * np.sqrt(3))) + center
sz = convert_vector(ij, from_units="ij", to_units="zernike", hardware=fs)

plt.title("Camera ij Coordinates")
plt.scatter(ij[0], ij[1])
plt.gca().set_aspect("equal")
plt.show()

plt.title("Zernike Basis Coordinates")
plt.scatter(sz[0], sz[1])
plt.gca().set_aspect("equal")
plt.show()

base = np.zeros((N,N))

base[2, :] = sz[0]      # Tilt x term
base[1, :] = sz[1]      # Tilt y term

We impart the Zernike term for each spot by adding a perturbation in Zernike-space, which we have requested to be \(N\)-dimensional; one order for each spot. The second and third rows correspond to the steering to position the spot at the correct \(n,l\) coordinate in the pyramid in 2D space, according to the \(y = Z_1\) and \(x = Z_2\) terms. Imparting the perturbation for \(Z_i\) on the correct spot at the correct \(n,l\) coordinate is as simple as adding an identity matrix. Then, the \(i\)th spot is given the perturbation for \(Z_i\):

perturbation = np.diag(np.ones(N))
print((base + 10*perturbation)[:10, :10].astype(int))

[[ 10   0   0   0   0   0   0   0   0   0]
 [298 256 245 194 193 192 142 140 139 138]
 [564 594 545 625 565 506 656 596 536 477]
 [  0   0   0  10   0   0   0   0   0   0]
 [  0   0   0   0  10   0   0   0   0   0]
 [  0   0   0   0   0  10   0   0   0   0]
 [  0   0   0   0   0   0  10   0   0   0]
 [  0   0   0   0   0   0   0  10   0   0]
 [  0   0   0   0   0   0   0   0  10   0]
 [  0   0   0   0   0   0   0   0   0  10]]

Then, we send these desired coordinates to optimization:

hologram = CompressedSpotHologram(
    spot_vectors=base + 5*perturbation,
    basis=i,
    cameraslm=fs
)

hologram.optimize()

C:\Users\holodyne\Documents\GitHub\slmsuite\slmsuite\holography\algorithms\_spots.py:360: UserWarning: Found ANSI index '0' (Zernike piston) in the zernike_basis; this is not necessary as spot phase is controlled externally.
  warnings.warn(

Let’s see what we got:

slm.set_phase(hologram.get_phase())
cam.flush()
cam.plot(limits=limits, title="Realized Zernike Pyramid");

This result looks even cleaner with virtual holography, but our real source distribution isn’t a perfect Gaussian on our physical SLM.

To show how programmable this technique of targeting points in Zernike-space is, we make a fun .gif, modulating the level of perturbation:

import tqdm.auto as tqdm

# Make the hologram
hologram = CompressedSpotHologram(
    spot_vectors=base,
    basis=i,
    cameraslm=fs
)
hologram.optimize("WGS-Leonardo", maxiter=50, verbose=False)

# Build the term osillation movie.
imgs = []
X = 5 * np.sin(np.linspace(0, 2*np.pi, 25, endpoint=False))

for x in tqdm.tqdm(X):
    # Optimize to the new position.
    hologram.spot_zernike = base + x*perturbation
    hologram.optimize("WGS-Leonardo", maxiter=20, verbose=False)

    # Project the pattern.
    slm.set_phase(hologram.get_phase())
    cam.flush()
    imgs.append(cam.get_image()[i3:i2:-1, i0:i1])

C:\Users\holodyne\Documents\GitHub\slmsuite\slmsuite\holography\algorithms\_spots.py:360: UserWarning: Found ANSI index '0' (Zernike piston) in the zernike_basis; this is not necessary as spot phase is controlled externally.
  warnings.warn(

(We use a helper function to render stored raw data with colormaps, saving the stack of images as a .gif):

from slmsuite.holography.analysis.files import save_image

save_image("ex-zernike-spots.gif", imgs, cmap="Blues", border=127, normalize=False, loop=0)
save_image("ex-zernike-spots-dark.gif", imgs, cmap="inferno", border=127, normalize=False, loop=0)

This feature is be applied to wavefront calibration in the next example, to remove the spatial dependence of aberration across a field of view.

cam.close()
slm.close()