GSW Phase Mask Computation

ScanImage synthesizes the phase mask for multi-point holographic stimulation using a weighted Gerchberg–Saxton algorithm (GSW). This page describes the algorithm, the data structures and coordinate systems involved, the configurable options exposed through the Slm Pattern GUI, and the path the result takes from the CPU/GPU compute kernel through to on-screen rendering.

The Gerchberg–Saxton–Weighted Algorithm

A phase-only SLM places a target spot at deflection \((\theta_x,\theta_y)\) and vergence \(v\) by adding a carrier phase

\[\phi_m(x,y) = \frac{2\pi}{\lambda}\sin\theta_x\, x + \frac{2\pi}{\lambda}\sin\theta_y\, y + \frac{\pi}{\lambda}\, v\,(x^2+y^2)\]

For a single point, the mask \(\phi=\phi_m\) is exact. For \(M\) simultaneous targets there is no closed-form mask: GSW iterates between target intensities and a unit-modulus phase mask, weighting each target so their delivered intensities converge toward the user-specified weights \(w_m\).

The complex amplitude delivered to target \(m\) by mask \(\phi\) under source amplitude \(A(x,y)=\sqrt{I(x,y)}\) is the field-overlap integral

\[V_m \;=\; \iint A(x,y)\, e^{i\phi(x,y)}\, e^{-i\phi_m(x,y)}\,dx\,dy\]

ScanImage normalizes \(|V_m|^2\) by \(V_{\text{norm}}^2 = \big(\iint A\,dx\,dy\big)^2\), the maximum possible single-point throughput. The result reads as “fraction of the perfect single-point efficiency delivered to target \(m\).”

The GSW iteration:

  1. Initialize phase weights \(w_m^{(0)}=1\) and start the mask from a random superposition of carrier phases.

  2. Compute \(V_m^{(k)}\) from the current mask.

  3. Update weights:

    \[w_m^{(k+1)} \;=\; w_m^{(k)}\cdot \left(\frac{\langle |V|/\sqrt{w}\rangle}{|V_m|/\sqrt{w_m}}\right)^{\!\alpha}\]

    The exponent \(\alpha\) (relaxationAlpha) damps the update so the algorithm doesn’t overshoot. \(\alpha=1\) is textbook GSW.

  4. Resynthesize the mask:

    \[\phi^{(k+1)} \;=\; \arg\!\left( \sum_{m=1}^{M} e^{i\phi_m}\, w_m^{(k+1)}\,\frac{V_m^{(k)}}{|V_m^{(k)}|} \right)\]

  5. Repeat until gswIterations is reached or targetUniformityThreshold is satisfied.

The algorithm minimizes the deviation between achieved and requested intensities, not zero-order leakage. Random phase across the SLM aperture produces near-zero DC by destructive interference, so the zero-order figure printed at compute time is typically < 1% — the physical zero-order in a real system is dominated by hardware imperfections (pixel fill factor, phase quantization, calibration error) which this idealized calculation does not model.

The slmPattern Coordinate System

For each stim ROI’s StimulusField, the field slmPattern is an [Nx4] matrix of points with columns

Column

Meaning

1 (X)

Lateral X offset from scanfield.centerXY, in reference space (optical degrees scaled by the imaging pipeline).

2 (Y)

Lateral Y offset from scanfield.centerXY, same units as X.

3 (Z)

Depth offset, in reference space Z (microns).

4 (Weight)

Relative target weight \(w_m\). Larger values request more intensity at that target. Internally normalized so \(\sum_m w_m = 1\) before GSW iterates. Zero-weight rows are dropped before computation.

Before the GSW kernel runs, the points pass through three coordinate systems:

  1. Reference space — the per-ROI representation seen by the user and stored on disk.

  2. Lateral diffraction efficiency space — hCSDiffractionEfficiency on the SLM. By default this is the analytical sinc² pixel-aperture envelope from SLM.analyticalDiffractionEfficiency; on a calibrated rig it is replaced by the measured CSLut from the diffraction-efficiency calibration. The transform produces a per-point lateral throughput \(\eta_{D,m}\in[0,1]\).

  3. SLM native deflection — hCoordinateSystem on the SLM. \((\theta_x,\theta_y)\) are angular deflections in optical degrees; \(Z\) is the wavefront vergence in \(\mu m^{-1}\). The origin of this system yields a blank phase mask. GSW receives points in this space and synthesizes the carrier phases \(\phi_m\) directly from the deflection/vergence numbers.

The composite per-point efficiency the GUI reports is

\[\eta_m \;=\; \eta_{D,m}\cdot\frac{|V_m|^2}{V_{\text{norm}}^2}\]

— lateral envelope times the GSW delivery fraction. The physical total efficiency is \(\sum_m \eta_m\); the phase-mask total efficiency reported by GSW alone is \(\sum_m |V_m|^2/V_{\text{norm}}^2\). The GUI stats label and the SLM phase mask console summary expose both numbers; they always satisfy physical efficiency ≤ phase mask efficiency.

The slmOptions Struct

Every StimulusField carries an slmOptions struct, normalized through scanimage.mroi.scanners.cghFunctions.defaultSlmOptions. The fields:

Field

Default

Meaning

plainGS

false

Skip the weight-reweighting step. The algorithm becomes plain Gerchberg–Saxton: highest total efficiency, no control over per-point uniformity. Disables relaxationAlpha.

gswIterations

10

Maximum number of GSW iterations. More iterations improve uniformity at linear compute cost.

numRandomSeeds

1

Number of random initial phase trials (MT-GSW). Each trial converges from a different basin; the highest-efficiency trial is kept. Cost scales linearly. Most useful when uniformity is poor with few targets.

relaxationAlpha

1.0

Damping exponent on the per-iteration weight update. 1.0 is textbook GSW (full step). Lower (0.5–0.8) damps oscillation when the algorithm overshoots — typically symptomatic with two well-separated targets.

targetUniformityThreshold

0

Early-stop when the Michelson uniformity \(1 - \frac{\max I_w-\min I_w}{\max I_w+\min I_w}\) (over weight-normalized intensities \(I_w=|V_m|^2/w_m\)) reaches this value. 0 disables and runs the full count.

Each setter on StimulusField.slmOptions runs the input through defaultSlmOptions so that partial structs back-fill the missing fields with defaults.

CPU Implementation

scanimage.mroi.scanners.cghFunctions.GSW is the reference implementation. It operates on the full SLM grid in MATLAB doubles (or computationDatatype) and is the path used when GPU is unavailable.

Outline:

  1. Copy the per-target carrier phases \(\phi_m\) once (cell array of MxN real arrays), and pre-compute \(A(x,y)\,e^{-i\phi_m(x,y)}\) for each target. This is the inner factor used by the overlap integral every iteration; computing it once turns each \(V_m\) evaluation into a single elementwise multiply-and-sum.

  2. For trial = 1 : numRandomSeeds (seeded by a hash of the point set so the result is reproducible per pattern), run a single GSW inner loop. Track the per-iteration uniformity in a history vector.

  3. After all trials, keep the trial with the highest summed efficiency along with its phase mask, point efficiencies, uniformity, and uniformity history.

  4. Compute the zero-order leakage from the winning unit-modulus mask: \(|\sum A\,e^{i\phi}|^2/V_{\text{norm}}^2\). This is the inner product the algorithm would have computed for a hypothetical target at the origin.

  5. Examine the kept trial’s uniformity trajectory. If it shows two significant downturns or finishes more than 0.02 below its in-trajectory peak, emit a stderr hint that relaxationAlpha is too aggressive.

The CPU kernel does not print a per-call summary itself — scanimage.mroi.scanners.SLM.computeMultiPointPhaseMask does that once with the combined physical efficiency, so users see a single line that matches the GUI display.

GPU Implementation

scanimage.mroi.scanners.cghFunctions.GSW_GPU is the GPU port. It is mathematically identical and gated by SLM.useGPUIfAvailable plus the runtime check most.util.gpuComputingAvailable. On hardware where the SLM has hundreds of thousands of pixels, the GPU path is roughly an order of magnitude faster.

Differences from the CPU kernel:

  1. Carrier phases are not cached as full-size arrays. Instead, each target contributes only its scaled \((\theta_x,\theta_y,v)\) triple. The phase at each pixel is recomputed inside the GPU stencil via arrayfun(@calculatePixelPhase, …) so the kernel touches only the small per-target state vector and the static pixel-position grids — far better register utilization than replaying M precomputed MxN arrays.

  2. The mask is stored as a real-valued angle on the GPU, so the zero-order calculation wraps with exp(1i*phi) before the inner product. The CPU kernel keeps the mask in unit-modulus complex form throughout and skips that wrap.

  3. Both kernels return iters (actual count after early-stop) and the per-iteration uniformity history so the SLM-level summary and the alpha-too-high heuristic are consistent.

GPU memory: the geometry buffer, source amplitude, and per-trial mask are pushed to GPU once per computeMultiPointPhaseMask call. Random phase initialization happens on the host (RandStream seeded by a hash of the point set) and is cast to single complex on upload.

Rendering and User Feedback

A single computeMultiPointPhaseMask call drives several pieces of the user interface. The flow:

  1. The SlmStimulusDisplay listens on the selected scanfield’s slmPattern and slmOptions PostSet events. When either changes, recomputeEfficiencies calls stimScannerset.multiPointFovToAo, which pipes through to SLM.computeMultiPointPhaseMask.

  2. The result populates two cached, non-firing properties on the StimulusField: slmPerPointEfficiency (one entry per point, physical) and slmZeroOrderEfficiency (scalar). Their setters intentionally do not call fireChangedEvent — the recompute is triggered by content changes, and writing back the cache must not retrigger compute.

  3. SlmPatternGui listens on slmPerPointEfficiency and refreshes the table. The fifth column is

    • Efficiency — fractional \(\eta_m\) — when the photostim’s first beam modulator lacks a power calibration, or

    • Power (mW) — \(\eta_m \cdot P_{\text{beam}}\cdot 1000\), where \(P_{\text{beam}}\) comes from hBeam.convertPowerFraction2PowerWatt(scanfield.powerFractions(1)), when both powerFraction2PowerWattLut and powerFraction2ModulationVoltLut are populated on the beam.

    The stats label above the table shows total physical efficiency and Michelson uniformity:

    Efficiency: 0.832    Uniformity: 0.998

  4. DiffractionEfficiencyDisplay paints the lateral diffraction-efficiency surface in the viewport. Its position is driven by hCSScannerOffset on the SLM scanner: when a stim ROI is selected the offset is parked at the ROI’s centerXY (so the surface follows the planned target); when the user moves the linear scanner target the SLM’s own listeners on the galvo targetPosition properties snap the surface to the new physical pointing. The surface stays put on ROI deselection until one of those events fires.

  5. The compute kernel itself is silent. SLM.computeMultiPointPhaseMask prints one summary line per call:

    SLM phase mask: physical efficiency=0.832 (phase mask=0.910), uniformity=0.998, zero-order=0.003, 10 iter(s), 1 seed(s), CPU compute time: 0.124s

    physical efficiency matches the GUI stats label exactly. phase mask is the GSW-internal figure (no diffraction envelope applied) — useful for diagnosing whether throughput is lost to point placement near the edge of the field (large gap) or to GSW configuration (small gap, low total).

    If the kept trial’s uniformity oscillated during convergence, a second stderr line follows:

    GSW: uniformity oscillated during convergence — relaxationAlpha=1.00 may be too high. Try lowering it.

Editing any control in the Slm Pattern GUI — sliders for gswIterations, numRandomSeeds, relaxationAlpha, targetUniformityThreshold, or the plainGS checkbox — mutates the scanfield’s slmOptions on slider release (MATLAB’s ValueChangedFcn does not fire during drag), which triggers a single recompute through the same chain.