#!/usr/bin/env python3 """ lcos_fpp_error_sim.py ===================== Error-budget simulation for LCoS-based fringe-projection profilometry. Compares Fourier-transform profilometry (FTP) and phase-shift profilometry (PSP) for three ways of using an LCoS device: A. Amplitude-mode LCoS projector (gray-level sinusoids between crossed polarizers) error sources: display gamma, 8-bit quantization, fringing-field blur, pixel fill factor, finite contrast, LC settling (frame-to-frame crosstalk), digital-drive (PWM) intensity ripple vs. camera exposure, camera noise. B. Phase-only LCoS showing a grating; the +1/-1 orders interfere to form the fringes. Phase steps come from shifting the grating laterally by whole SLM pixels (exact). error sources: speckle, SLM phase flicker (common-mode -> cancels), camera noise. C. Phase-only LCoS used as a uniform phase shifter in one arm of an interferometric fringe projector. error sources: pixel-to-pixel PRETILT nonuniformity (e.g. from oblique/GLAD SiO2 alignment-layer column-tilt variation) -> phase-vs-voltage curve differs from the calibration LUT -> spatially varying phase-step error; SLM phase flicker. Decoders: FTP, pi-shift FTP, N-step PSP (N = 3, 4, 8), Hariharan 5-step (first-order insensitive to phase-step scale error). Every measurement is done like a real system: a flat reference plane and the object are both measured with the same simulated hardware, and the height phase is their difference. Errors are reported as RMS phase error (rad) and RMS height error (um) after removing piston, inside an interior region (FTP edge effects excluded). Sweeps of systematic errors (gamma, settling, ripple, pretilt, phase flicker) use a noise-free camera so each source is seen on its own; a dotted line marks the camera-noise floor. The overview, period and speckle studies use the noisy camera. The LC electro-optic model for approach C is a deliberately simple mean-tilt model (documented in PhaseLCoS); swap in your own director-profile solver or measured phase-vs-gray curves for quantitative work. To use a pretilt map from the GLAD Monte Carlo, save it as a 2-D .npy array in degrees and pass --pretilt-map. Outputs (in --outdir, default ./lcos_fpp_results): fig1_overview.png true height and error maps for representative cases fig2_gamma.png RMS height error vs display gamma, all decoders fig3_timing.png LC settling time and PWM ripple vs exposure fig4_pretilt.png phase-vs-voltage curves and error vs pretilt spread fig5_period_speckle_flicker.png summary.csv every simulated point Requires numpy, scipy, matplotlib. """ from __future__ import annotations import argparse import csv import os import sys from dataclasses import dataclass, replace from typing import Callable, Optional import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from scipy.ndimage import gaussian_filter, gaussian_filter1d, zoom TWO_PI = 2 * np.pi def wrap(p): """Wrap phase to (-pi, pi].""" return np.angle(np.exp(1j * p)) # ---------------------------------------------------------------------------------- # Geometry, scene, camera # ---------------------------------------------------------------------------------- @dataclass class Geometry: n: int = 512 # camera image is n x n pixels period_px: float = 16.0 # fringe period seen on camera (projector:camera = 1:1) pixel_mm: float = 0.0625 # object-space pixel size -> 1 fringe = 1 mm here L_mm: float = 500.0 # camera/projector pupil to reference plane d_mm: float = 150.0 # projector-camera baseline margin: int = 32 # border excluded from statistics @property def k_rad_per_mm(self) -> float: """Phase per unit height, small-height approximation of h = L*dphi/(dphi - 2*pi*f0*d) -> dphi = 2*pi*f0*d*h/L.""" period_mm = self.period_px * self.pixel_mm return TWO_PI * self.d_mm / (self.L_mm * period_mm) @dataclass class Scene: h_mm: np.ndarray # object height map R_obj: np.ndarray # object reflectance R_ref: np.ndarray # reference-plane reflectance def make_scene(g: Geometry, seed: int = 1) -> Scene: n = g.n y, x = np.mgrid[0:n, 0:n] * g.pixel_mm c = n * g.pixel_mm / 2 s = g.n / 512 # keep features the same fraction of the field for other n h = (3.0 * np.exp(-((x - c + 6 * s) ** 2 + (y - c) ** 2) / (2 * (4.0 * s) ** 2)) + 1.5 * np.exp(-((x - c - 7 * s) ** 2 + (y - c + 5 * s) ** 2) / (2 * (2.5 * s) ** 2)) - 1.0 * np.exp(-((x - c - 4 * s) ** 2 + (y - c - 8 * s) ** 2) / (2 * (3.0 * s) ** 2))) rng = np.random.default_rng(seed) r = gaussian_filter(rng.standard_normal((n, n)), 20 * s) r = (r - r.min()) / (r.max() - r.min()) R_obj = 0.45 + 0.45 * r R_ref = np.full((n, n), 0.85) return Scene(h, R_obj, R_ref) @dataclass class Camera: gain: float = 0.85 # projector full white on R=1 -> fraction of full scale ambient: float = 0.03 full_well_e: float = 12000.0 read_sigma: float = 0.003 # fraction of full scale bits: int = 10 noise: bool = True # False: only quantization (isolates systematic errors) def capture(self, I_proj, R, rng): I = R * I_proj * self.gain + self.ambient if self.noise: I = rng.poisson(np.clip(I, 0, None) * self.full_well_e) / self.full_well_e I = I + rng.normal(0.0, self.read_sigma, I.shape) q = 2 ** self.bits - 1 return np.round(np.clip(I, 0, 1) * q) / q # ---------------------------------------------------------------------------------- # Approach A: amplitude LCoS projector # ---------------------------------------------------------------------------------- @dataclass class AmpLCoS: gamma: float = 2.2 # T = drive**gamma (normalized) gamma_est: Optional[float] = None # if set, patterns are pre-distorted with 1/gamma_est bits: int = 8 contrast: float = 1000.0 # on/off fill_factor: float = 0.93 fringing_sigma_px: float = 0.6 # director smoothing between pixels (fringing field) defocus_sigma_px: float = 0.8 # projection-lens blur at the object tau_ms: float = 2.5 # LC 1/e response time frame_ms: float = 1000 / 120 settle_ms: float = 5.0 # delay from pattern update to exposure start exposure_ms: Optional[float] = None # default: frame - settle first_prev: str = "black" # what preceded frame 0: 'black' or 'cyclic' ripple_amp: float = 0.0 # fractional intensity ripple from digital (PWM) drive ripple_hz: float = 1440.0 oversample: int = 8 def exposure(self): return self.exposure_ms if self.exposure_ms is not None else self.frame_ms - self.settle_ms def settle_mix(self): """Fraction of the previous pattern still present, averaged over the exposure, for a first-order LC response with time constant tau.""" te, ts, tau = self.exposure(), self.settle_ms, self.tau_ms if tau <= 0: return 0.0 return tau / te * (np.exp(-ts / tau) - np.exp(-(ts + te) / tau)) def ripple_gain(self, rng): """Exposure-averaged gain of a sinusoidal PWM ripple with random trigger phase. Zero error when exposure is an integer number of ripple periods.""" if self.ripple_amp == 0: return 1.0 w = TWO_PI * self.ripple_hz / 1000.0 T = self.exposure() psi = rng.uniform(0, TWO_PI) return 1.0 + self.ripple_amp * (np.cos(psi) - np.cos(w * T + psi)) / (w * T) def profile(self, delta, period_px, n): """1-D projected intensity vs projector column u (fine grid).""" M, os_ = 64, self.oversample j = np.arange(-M, n + M) t = 0.5 + 0.5 * np.cos(TWO_PI * j / period_px + delta) drive = t ** (1.0 / self.gamma_est) if self.gamma_est else t q = 2 ** self.bits - 1 drive = np.round(drive * q) / q frac = (np.tile(np.arange(os_), len(j)) + 0.5) / os_ u = np.repeat(j, os_) + frac - 0.5 fine = gaussian_filter1d(np.repeat(drive, os_), max(self.fringing_sigma_px * os_, 1e-6), mode="nearest") Tb = 1.0 / self.contrast T = Tb + (1 - Tb) * np.clip(fine, 0, 1) ** self.gamma gap = (1 - self.fill_factor) / 2 lo, hi = frac - 0.5 / os_, frac + 0.5 / os_ mask = np.clip((np.minimum(hi, 1 - gap) - np.maximum(lo, gap)) * os_, 0, 1) T = T * mask T = gaussian_filter1d(T, max(self.defocus_sigma_px * os_, 1e-6), mode="nearest") return u, T class AmplitudeProjector: name = "A: amplitude LCoS" def __init__(self, p: AmpLCoS, cam: Camera): self.p, self.cam = p, cam def frames(self, dphi, deltas, R, g: Geometry, rng): n = dphi.shape[0] profs = [self.p.profile(d, g.period_px, n) for d in deltas] u = profs[0][0] Ts = [pr[1] for pr in profs] m = self.p.settle_mix() ucam = np.arange(n)[None, :] + dphi / TWO_PI * g.period_px out = [] for i, T in enumerate(Ts): if i == 0 and len(Ts) > 1 and self.p.first_prev == "black": prev = np.full_like(T, 1.0 / self.p.contrast) else: prev = Ts[i - 1] D = (T + m * (prev - T)) * self.p.ripple_gain(rng) out.append(self.cam.capture(np.interp(ucam, u, D), R, rng)) return out # ---------------------------------------------------------------------------------- # Approach B: phase-only LCoS grating, +/-1 orders interfere # ---------------------------------------------------------------------------------- @dataclass class GratingSLM: visibility: float = 0.95 # partial coherence / order imbalance speckle_contrast: float = 0.0 # 0 = LED/SLD-like, 1 = fully developed laser speckle speckle_grain_px: float = 2.0 phase_flicker_rad: float = 0.0 # per-frame SLM phase flicker (std) diff_flicker_fraction: float = 0.0 # part of flicker NOT common to both orders class GratingInterferometer: """Fringes from two diffraction orders of an SLM grating. A lateral grating shift of dx gives the +/-1 orders phases +/-2*pi*dx/p, so the fringe phase steps by 4*pi*dx/p. With p = 32 SLM px, dx = 4 px gives exactly pi/2, so steps are exact; flicker that is common to the whole SLM cancels in |E1+E2|^2.""" name = "B: SLM grating (+/-1 orders)" def __init__(self, p: GratingSLM, cam: Camera): self.p, self.cam = p, cam def _speckle(self, n, rng): if self.p.speckle_contrast == 0: return 1.0 f = gaussian_filter(rng.standard_normal((n, n)), self.p.speckle_grain_px) + \ 1j * gaussian_filter(rng.standard_normal((n, n)), self.p.speckle_grain_px) S = np.abs(f) ** 2 S /= S.mean() c = self.p.speckle_contrast return (1 - c) + c * S def frames(self, dphi, deltas, R, g: Geometry, rng): n = dphi.shape[0] arg = TWO_PI * np.arange(n)[None, :] / g.period_px + dphi S = self._speckle(n, rng) # static during one measurement out = [] for d in deltas: f = rng.normal(0, self.p.phase_flicker_rad) fd = self.p.diff_flicker_fraction * rng.normal(0, self.p.phase_flicker_rad) E1 = np.exp(1j * (arg / 2 + d / 2 + f + fd / 2)) E2 = np.exp(-1j * (arg / 2 + d / 2) + 1j * (f - fd / 2)) I = 0.5 * (1 + self.p.visibility * np.real(E1 * np.conj(E2))) out.append(self.cam.capture(I * S, R, rng)) return out # ---------------------------------------------------------------------------------- # Approach C: phase-only LCoS as a uniform phase shifter, with pretilt nonuniformity # ---------------------------------------------------------------------------------- @dataclass class PhaseLCoS: """Mean-tilt ECB (parallel-aligned) reflective cell. theta(V, theta_p) = theta_p + (pi/2 - theta_p) * (1 - exp(-(V - Vth)/V0)), V > Vth Vth(theta_p) = Vth0 * (1 - 2*theta_p/pi) (pretilt lowers the threshold) n_eff(theta) = ne*no / sqrt(ne^2 sin^2 theta + no^2 cos^2 theta) phi(V) = 4*pi*d*(n_eff - no)/lambda (double pass) This is a toy model, adequate for showing how pretilt variation turns into phase-step error. Replace phase() with a director-profile solver or measured curves for quantitative work.""" ne: float = 1.65 no: float = 1.50 d_um: float = 3.0 lam_um: float = 0.532 pretilt0_deg: float = 3.0 # pretilt the LUT is calibrated at Vth0: float = 1.0 V0: float = 1.2 def tilt(self, V, tp): Vth = self.Vth0 * (1 - 2 * tp / np.pi) s = 1 - np.exp(-np.clip(V - Vth, 0, None) / self.V0) return tp + (np.pi / 2 - tp) * s def phase(self, V, tp_deg): tp = np.deg2rad(tp_deg) th = self.tilt(V, tp) neff = self.ne * self.no / np.sqrt(self.ne ** 2 * np.sin(th) ** 2 + self.no ** 2 * np.cos(th) ** 2) return 4 * np.pi * self.d_um * (neff - self.no) / self.lam_um def calibrate(self, targets): """Voltages giving phase shifts `targets` (rad, >= 0) at the nominal pretilt.""" V = np.linspace(0, 12, 24001) shift = self.phase(0.0, self.pretilt0_deg) - self.phase(V, self.pretilt0_deg) if np.max(targets) > shift.max(): raise ValueError("cell cannot reach required phase shift; increase d or dn") return np.interp(targets, shift, V) def make_pretilt_map(n, lc: PhaseLCoS, spread_pv_deg, seed=3, path=None): """Pretilt map in degrees: nominal + peak-to-valley spread made of a linear gradient (deposition-angle variation across the substrate) plus a smooth random component. If `path` is given, load a 2-D .npy map (degrees) and resample it.""" if path: m = np.load(path).astype(float) return zoom(m, (n / m.shape[0], n / m.shape[1]), order=1) rng = np.random.default_rng(seed) y, x = np.mgrid[0:n, 0:n] / (n - 1) lin = 0.7 * (2 * x - 1) + 0.3 * (2 * y - 1) rnd = gaussian_filter(rng.standard_normal((n, n)), n / 10) base = 0.6 * lin / np.abs(lin).max() + 0.4 * rnd / np.abs(rnd).max() base = 2 * (base - base.min()) / (base.max() - base.min()) - 1 # -> [-1, 1] return lc.pretilt0_deg + 0.5 * spread_pv_deg * base class ShifterInterferometer: name = "C: LCoS phase shifter" def __init__(self, lc: PhaseLCoS, pretilt_map_deg, cam: Camera, visibility=0.95, phase_flicker_rad=0.0): self.lc, self.tp, self.cam = lc, pretilt_map_deg, cam self.visibility, self.phase_flicker_rad = visibility, phase_flicker_rad def actual_steps(self, deltas): d = np.asarray(deltas, float) off = -d.min() V = self.lc.calibrate(d + off) phi = [self.lc.phase(v, self.tp) for v in V] return [phi[0] - p - off + d[0] for p in phi] # per-pixel actual steps def frames(self, dphi, deltas, R, g: Geometry, rng): n = dphi.shape[0] arg = TWO_PI * np.arange(n)[None, :] / g.period_px + dphi out = [] for da in self.actual_steps(deltas): f = rng.normal(0, self.phase_flicker_rad) I = 0.5 * (1 + self.visibility * np.cos(arg + da + f)) out.append(self.cam.capture(I, R, rng)) return out # ---------------------------------------------------------------------------------- # Decoders # ---------------------------------------------------------------------------------- def dec_psp(frames, deltas, g): s = sum(I * np.sin(d) for I, d in zip(frames, deltas)) c = sum(I * np.cos(d) for I, d in zip(frames, deltas)) return np.arctan2(-s, c) def dec_hariharan(frames, deltas, g): I1, I2, I3, I4, I5 = frames return np.arctan2(2 * (I2 - I4), 2 * I3 - I1 - I5) def _ftp_core(I, g, bw=0.5): n = I.shape[1] f = np.fft.fftfreq(n) f0 = 1.0 / g.period_px x = (f - f0) / (bw * f0) H = np.where(np.abs(x) < 1, 0.5 * (1 + np.cos(np.pi * x)), 0.0) a = np.fft.ifft(np.fft.fft(I, axis=1) * H[None, :], axis=1) return np.angle(a) def dec_ftp(frames, deltas, g): return _ftp_core(frames[0], g) def dec_piftp(frames, deltas, g): return _ftp_core(frames[0] - frames[1], g) def psp_deltas(N): return [TWO_PI * k / N for k in range(N)] METHODS: dict[str, tuple[list, Callable]] = { "FTP": ([0.0], dec_ftp), "pi-FTP": ([0.0, np.pi], dec_piftp), "PSP3": (psp_deltas(3), dec_psp), "PSP4": (psp_deltas(4), dec_psp), "PSP8": (psp_deltas(8), dec_psp), "Hariharan5": ([-np.pi, -np.pi / 2, 0.0, np.pi / 2, np.pi], dec_hariharan), } # ---------------------------------------------------------------------------------- # Measurement # ---------------------------------------------------------------------------------- def measure(system, method, scene: Scene, g: Geometry, rng): deltas, dec = METHODS[method] dphi_true = g.k_rad_per_mm * scene.h_mm ph_ref = dec(system.frames(np.zeros_like(dphi_true), deltas, scene.R_ref, g, rng), deltas, g) ph_obj = dec(system.frames(dphi_true, deltas, scene.R_obj, g, rng), deltas, g) dphi = wrap(ph_obj - ph_ref) err = wrap(dphi - dphi_true) m = g.margin core = err[m:-m, m:-m] piston = np.angle(np.mean(np.exp(1j * core))) err = wrap(err - piston) core = err[m:-m, m:-m] rms = float(np.sqrt(np.mean(core ** 2))) # unwrapped height for display: row-wise unwrap, fringe order fixed from column 0 un = np.unwrap(dphi, axis=1) un -= TWO_PI * np.round((un[:, :1] - dphi_true[:, :1]) / TWO_PI) return dict(rms_rad=rms, rms_um=rms / g.k_rad_per_mm * 1000.0, err_um=err / g.k_rad_per_mm * 1000.0, h_mm=un / g.k_rad_per_mm) class Recorder: def __init__(self): self.rows = [] def add(self, study, approach, method, param, value, res): self.rows.append(dict(study=study, approach=approach, method=method, param=param, value=value, rms_phase_rad=res["rms_rad"], rms_height_um=res["rms_um"])) def save(self, path): # utf-8-sig: Greek letters (γ, τ) in labels; the BOM lets Excel on Windows read them correctly with open(path, "w", newline="", encoding="utf-8-sig") as f: w = csv.DictWriter(f, fieldnames=list(self.rows[0].keys())) w.writeheader() w.writerows(self.rows) def avg_measure(system, method, scene, g, seed, reps): rs = [measure(system, method, scene, g, np.random.default_rng(seed + 101 * r)) for r in range(reps)] out = dict(rs[0]) out["rms_rad"] = float(np.mean([r["rms_rad"] for r in rs])) out["rms_um"] = float(np.mean([r["rms_um"] for r in rs])) return out # ---------------------------------------------------------------------------------- # Studies # ---------------------------------------------------------------------------------- def study_gamma(g, scene, cam, rec, seed, floor): gammas = [1.0, 1.4, 1.8, 2.2, 2.6] curves = {} fast = dict(tau_ms=0.3, ripple_amp=0.0) # isolate gamma from timing effects for meth in ["FTP", "pi-FTP", "PSP3", "PSP4", "PSP8"]: curves[meth] = [] for gm in gammas: sysA = AmplitudeProjector(AmpLCoS(gamma=gm, **fast), cam) r = measure(sysA, meth, scene, g, np.random.default_rng(seed)) curves[meth].append(r["rms_um"]) rec.add("gamma", "A", meth, "gamma", gm, r) for meth in ["PSP3", "PSP4"]: key = meth + " +gamma LUT (5% err)" curves[key] = [] for gm in gammas: sysA = AmplitudeProjector(AmpLCoS(gamma=gm, gamma_est=gm * 1.05, **fast), cam) r = measure(sysA, meth, scene, g, np.random.default_rng(seed)) curves[key].append(r["rms_um"]) rec.add("gamma", "A", key, "gamma", gm, r) fig, ax = plt.subplots(figsize=(9.5, 4.5)) for k, v in curves.items(): ax.plot(gammas, v, "o-" if "LUT" not in k else "s--", label=k) ax.axhline(floor, color="gray", ls=":", lw=1.5, label="camera-noise floor (PSP4, ideal fringes)") ax.set_yscale("log") ax.set_xlabel("display gamma (T = drive$^\\gamma$)") ax.set_ylabel("RMS height error (µm)") ax.set_title("A: amplitude LCoS — gamma nonlinearity (noise-free camera)") ax.grid(alpha=0.3, which="both") ax.legend(fontsize=8, loc="center left", bbox_to_anchor=(1.01, 0.5)) ax.text(0.02, 0.02, "FTP and pi-FTP overlap: their floor is filter leakage, not gamma", transform=ax.transAxes, fontsize=7, color="gray") fig.tight_layout() return fig def study_timing(g, scene, cam, rec, seed, floor): fig, axs = plt.subplots(1, 2, figsize=(11, 4.3)) # (a) LC settling delay settles = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5] for meth, prev, sty in [("PSP4", "black", "o-"), ("PSP4", "cyclic", "o--"), ("PSP8", "black", "s-"), ("pi-FTP", "black", "^-")]: vals = [] for ts in settles: p = AmpLCoS(gamma=2.2, gamma_est=2.2, settle_ms=ts, first_prev=prev) r = measure(AmplitudeProjector(p, cam), meth, scene, g, np.random.default_rng(seed)) vals.append(r["rms_um"]) rec.add("settle", "A", f"{meth} ({prev} before frame 0)", "settle_ms", ts, r) axs[0].plot(settles, vals, sty, label=f"{meth}, {prev} before frame 0") p = AmpLCoS() axs[0].axhline(floor, color="gray", ls=":", lw=1.5, label="camera-noise floor (PSP4)") axs[0].set_yscale("log") axs[0].set_xlabel(f"delay before exposure (ms) [120 Hz frame, τ = {p.tau_ms} ms]") axs[0].set_ylabel("RMS height error (µm)") axs[0].set_title("LC settling (frame-to-frame crosstalk)") axs[0].grid(alpha=0.3, which="both") axs[0].legend(fontsize=8) # (b) PWM ripple vs exposure length Tr = 1000.0 / 1440.0 ratios = np.linspace(1.0, 6.0, 41) for meth, sty in [("PSP4", "-"), ("PSP8", "--")]: vals = [] for rr in ratios: p = AmpLCoS(gamma=2.2, gamma_est=2.2, tau_ms=0.01, settle_ms=0.5, exposure_ms=rr * Tr, ripple_amp=0.05) r = avg_measure(AmplitudeProjector(p, cam), meth, scene, g, seed, reps=3) vals.append(r["rms_um"]) rec.add("ripple", "A", meth, "exposure/ripple_period", float(rr), r) axs[1].plot(ratios, vals, sty, label=meth) axs[1].axhline(floor, color="gray", ls=":", lw=1.5, label="camera-noise floor (PSP4)") axs[1].set_yscale("log") axs[1].set_xlabel("camera exposure / PWM ripple period (5% ripple, 1440 Hz)") axs[1].set_ylabel("RMS height error (µm)") axs[1].set_title("Digital-drive ripple: lock exposure to integer periods") axs[1].grid(alpha=0.3) axs[1].legend(fontsize=8) fig.tight_layout() return fig def study_pretilt(g, scene, cam, rec, seed, lc, pretilt_path, floor): fig, axs = plt.subplots(1, 2, figsize=(11, 4.3)) V = np.linspace(0, 6, 600) for tp in [1, 3, 5, 7]: axs[0].plot(V, lc.phase(0, lc.pretilt0_deg) - lc.phase(V, tp), label=f"pretilt {tp}°" + (" (LUT)" if tp == lc.pretilt0_deg else "")) Vn = lc.calibrate(np.array(psp_deltas(4))) for v in Vn: axs[0].axvline(v, color="k", lw=0.6, ls=":") axs[0].set_xlabel("drive voltage (V)") axs[0].set_ylabel("phase shift vs 0 V at LUT pretilt (rad)") axs[0].set_title("C: toy ECB model — LUT voltages for PSP4 (dotted)") axs[0].grid(alpha=0.3) axs[0].legend(fontsize=8) spreads = [0, 0.5, 1, 2, 3, 4, 6] for meth, sty in [("PSP4", "o-"), ("PSP8", "s-"), ("Hariharan5", "^-")]: vals = [] for sp in spreads: tp = make_pretilt_map(g.n, lc, sp) r = measure(ShifterInterferometer(lc, tp, cam), meth, scene, g, np.random.default_rng(seed)) vals.append(r["rms_um"]) rec.add("pretilt", "C", meth, "pretilt_PV_deg", sp, r) axs[1].plot(spreads, vals, sty, label=meth) axs[1].axhline(floor, color="gray", ls=":", lw=1.5, label="camera-noise floor (PSP4)") axs[1].set_yscale("log") if pretilt_path: tp = make_pretilt_map(g.n, lc, 0, path=pretilt_path) for meth in ["PSP4", "Hariharan5"]: r = measure(ShifterInterferometer(lc, tp, cam), meth, scene, g, np.random.default_rng(seed)) rec.add("pretilt_user_map", "C", meth, "map", os.path.basename(pretilt_path), r) print(f" user pretilt map, {meth}: {r['rms_um']:.2f} um RMS") axs[1].set_xlabel(f"pretilt peak-to-valley across aperture (°), LUT at {lc.pretilt0_deg}°") axs[1].set_ylabel("RMS height error (µm)") axs[1].set_title("C: pretilt nonuniformity → phase-step error") axs[1].grid(alpha=0.3) axs[1].legend(fontsize=8) fig.tight_layout() return fig def study_period_speckle_flicker(g, scene, cam, cam_noisy, rec, seed, lc, floor): fig, axs = plt.subplots(1, 3, figsize=(15, 4.3)) # (a) fringe period vs fringing-field blur (amplitude LCoS, PSP4, gamma corrected) periods = [6, 8, 10, 12, 16, 24, 32] for fs in [0.3, 1.0, 1.8]: vals = [] for P in periods: gP = replace(g, period_px=float(P)) p = AmpLCoS(gamma=2.2, gamma_est=2.2, tau_ms=0.3, fringing_sigma_px=fs) r = measure(AmplitudeProjector(p, cam_noisy), "PSP4", scene, gP, np.random.default_rng(seed)) vals.append(r["rms_um"]) rec.add("period", "A", "PSP4", f"period_px (fringing sigma {fs} px)", P, r) axs[0].plot(periods, vals, "o-", label=f"fringing σ = {fs} px") axs[0].set_xlabel("fringe period (projector px)") axs[0].set_ylabel("RMS height error (µm)") axs[0].set_title("A: period vs fringing-field contrast loss (PSP4, noisy camera)") axs[0].grid(alpha=0.3) axs[0].legend(fontsize=8) # (b) speckle contrast, approach B cs = [0, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0] for meth in ["PSP4", "PSP8"]: vals = [] for c in cs: r = measure(GratingInterferometer(GratingSLM(speckle_contrast=c), cam_noisy), meth, scene, g, np.random.default_rng(seed)) vals.append(r["rms_um"]) rec.add("speckle", "B", meth, "speckle_contrast", c, r) axs[1].plot(cs, vals, "o-", label=meth) axs[1].set_xlabel("speckle contrast") axs[1].set_ylabel("RMS height error (µm)") axs[1].set_title("B: speckle (static during sequence, noisy camera)") axs[1].grid(alpha=0.3) axs[1].legend(fontsize=8) # (c) SLM phase flicker: B (common-mode) vs C (one arm) sig = [0, 0.02, 0.05, 0.1, 0.2, 0.3] tp = make_pretilt_map(g.n, lc, 0) for label, mk in [ ("B grating, PSP4", lambda s: GratingInterferometer(GratingSLM(phase_flicker_rad=s), cam)), ("B grating, 20% differential", lambda s: GratingInterferometer( GratingSLM(phase_flicker_rad=s, diff_flicker_fraction=0.2), cam)), ("C shifter, PSP4 (uniform pretilt)", lambda s: ShifterInterferometer(lc, tp, cam, phase_flicker_rad=s)), ]: vals = [] for s in sig: r = avg_measure(mk(s), "PSP4", scene, g, seed, reps=3) vals.append(r["rms_um"]) rec.add("phase_flicker", label[0], label, "flicker_rad", s, r) axs[2].plot(sig, vals, "o-", label=label) axs[2].axhline(floor, color="gray", ls=":", lw=1.5, label="camera-noise floor (PSP4)") axs[2].set_yscale("log") axs[2].set_xlabel("SLM phase flicker, per-frame std (rad)") axs[2].set_ylabel("RMS height error (µm)") axs[2].set_title("Phase flicker: common-mode cancels in B (noise-free)") axs[2].grid(alpha=0.3) axs[2].legend(fontsize=8) fig.tight_layout() return fig def overview(g, scene, cam, rec, seed, lc, pretilt_path): rng = lambda: np.random.default_rng(seed) cases = [] a_raw = AmplitudeProjector(AmpLCoS(gamma=2.2, tau_ms=0.3), cam) a_slow = AmplitudeProjector(AmpLCoS(gamma=2.2, gamma_est=2.2, settle_ms=1.5), cam) a_cor = AmplitudeProjector(AmpLCoS(gamma=2.2, gamma_est=2.2), cam) cases.append(("A PSP3, γ=2.2 uncorrected\n(fast LC)", measure(a_raw, "PSP3", scene, g, rng()))) cases.append(("A PSP4, gamma LUT, 120 Hz,\nτ=2.5 ms, 1.5 ms settle", measure(a_slow, "PSP4", scene, g, rng()))) cases.append(("A FTP single shot,\ngamma LUT", measure(a_cor, "FTP", scene, g, rng()))) tp = make_pretilt_map(g.n, lc, 3.0, path=pretilt_path) sysC = ShifterInterferometer(lc, tp, cam) cases.append(("C PSP4, pretilt map", measure(sysC, "PSP4", scene, g, rng()))) cases.append(("C Hariharan 5-step,\npretilt map", measure(sysC, "Hariharan5", scene, g, rng()))) for name, r in cases: rec.add("overview", name[0], name.replace("\n", " "), "-", 0, r) steps = sysC.actual_steps(METHODS["PSP4"][0]) step_err = steps[1] - np.pi / 2 fig, axs = plt.subplots(2, 4, figsize=(17, 8)) ext = [0, g.n * g.pixel_mm, g.n * g.pixel_mm, 0] im = axs[0, 0].imshow(scene.h_mm, cmap="viridis", extent=ext) axs[0, 0].set_title("true height (mm)") fig.colorbar(im, ax=axs[0, 0], shrink=0.8) for ax, (name, r) in zip([axs[0, 1], axs[0, 2], axs[0, 3]], cases[:3]): e = r["err_um"] v = np.percentile(np.abs(e[g.margin:-g.margin, g.margin:-g.margin]), 99) im = ax.imshow(e, cmap="RdBu_r", vmin=-v, vmax=v, extent=ext) ax.set_title(f"{name}\nRMS {r['rms_um']:.1f} µm") fig.colorbar(im, ax=ax, shrink=0.8, label="height error (µm)") im = axs[1, 0].imshow(tp, cmap="magma", extent=ext) axs[1, 0].set_title("C: pretilt map (°)") fig.colorbar(im, ax=axs[1, 0], shrink=0.8) im = axs[1, 1].imshow(np.rad2deg(step_err), cmap="RdBu_r", extent=ext) axs[1, 1].set_title("C: error of the π/2 step (deg)") fig.colorbar(im, ax=axs[1, 1], shrink=0.8) vC = max(np.percentile(np.abs(r["err_um"][g.margin:-g.margin, g.margin:-g.margin]), 99) for _, r in cases[3:]) for ax, (name, r) in zip([axs[1, 2], axs[1, 3]], cases[3:]): im = ax.imshow(r["err_um"], cmap="RdBu_r", vmin=-vC, vmax=vC, extent=ext) ax.set_title(f"{name}\nRMS {r['rms_um']:.1f} µm") fig.colorbar(im, ax=ax, shrink=0.8, label="height error (µm)") for ax in axs.ravel(): ax.set_xlabel("x (mm)") ax.set_ylabel("y (mm)") fig.suptitle(f"LCoS fringe projection — {g.n}×{g.n} px, 1 mm fringe period, " f"sensitivity {g.k_rad_per_mm:.2f} rad/mm", fontsize=12) fig.tight_layout() return fig, cases # ---------------------------------------------------------------------------------- def main(): # Windows consoles may default to cp1252; never crash on γ / τ / µ in printed labels for stream in (sys.stdout, sys.stderr): try: stream.reconfigure(errors="replace") except (AttributeError, ValueError): pass ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--outdir", default="lcos_fpp_results") ap.add_argument("--n", type=int, default=512, help="image size (px)") ap.add_argument("--seed", type=int, default=7) ap.add_argument("--pretilt-map", default=None, help=".npy pretilt map in degrees (approach C)") ap.add_argument("--pretilt0", type=float, default=3.0, help="pretilt the phase LUT is calibrated at (deg)") ap.add_argument("--quick", action="store_true", help="256 px for a fast run") args = ap.parse_args() n = 256 if args.quick else args.n g = Geometry(n=n, margin=max(16, n // 16)) # 1 mm fringe period; field = n * 62.5 um os.makedirs(args.outdir, exist_ok=True) scene = make_scene(g, seed=1) cam_noisy = Camera() # realistic camera, used for overview / speckle / period cam = Camera(noise=False) # isolates each systematic error source lc = PhaseLCoS(pretilt0_deg=args.pretilt0) rec = Recorder() print(f"Geometry: {g.n}x{g.n} px, period {g.period_px} px = {g.period_px * g.pixel_mm:.2f} mm, " f"k = {g.k_rad_per_mm:.3f} rad/mm (1 mrad = {1 / g.k_rad_per_mm:.2f} um)") fl = measure(GratingInterferometer(GratingSLM(), cam_noisy), "PSP4", scene, g, np.random.default_rng(args.seed)) floor = fl["rms_um"] rec.add("noise_floor", "B", "PSP4 ideal fringes, noisy camera", "-", 0, fl) print(f"camera-noise floor (PSP4, ideal fringes): {floor:.2f} um RMS") jobs = [ ("fig1_overview.png", lambda: overview(g, scene, cam_noisy, rec, args.seed, lc, args.pretilt_map)[0]), ("fig2_gamma.png", lambda: study_gamma(g, scene, cam, rec, args.seed, floor)), ("fig3_timing.png", lambda: study_timing(g, scene, cam, rec, args.seed, floor)), ("fig4_pretilt.png", lambda: study_pretilt(g, scene, cam, rec, args.seed, lc, args.pretilt_map, floor)), ("fig5_period_speckle_flicker.png", lambda: study_period_speckle_flicker(g, scene, cam, cam_noisy, rec, args.seed, lc, floor)), ] for fname, job in jobs: print(f"running {fname} ...", flush=True) fig = job() fig.savefig(os.path.join(args.outdir, fname), dpi=140) plt.close(fig) rec.save(os.path.join(args.outdir, "summary.csv")) print("\nOverview cases:") for r in rec.rows: if r["study"] == "overview": print(f" {r['method']:<45s} {r['rms_phase_rad'] * 1000:8.1f} mrad {r['rms_height_um']:8.2f} um") print(f"\nresults written to {os.path.abspath(args.outdir)}") if __name__ == "__main__": main()