#!/usr/bin/env python3 """ ftp_walkthrough.py ================== Step-by-step demonstration of Fourier Transform Profilometry (Takeda & Mutoh 1983) and its main limits. Figure 1 ftp_1_pipeline.png One measurement, every stage: fringe image -> 2-D spectrum -> 1-D spectrum with the band-pass filter -> wrapped phase -> unwrapped height -> error map and profile. The projector has an uncorrected gamma so the harmonic lobes are visible, and the object carries a dark reflectance mark so you can see the zero order broaden. Figure 2 ftp_2_slope_limit.png A ridge whose steepest slope is swept through Takeda's limit |dh/dx| < L/(3d). Spectra below and above the limit, and RMS error vs slope for FTP and pi-shift FTP, each with the best band-pass width from a small set (an "oracle" choice that shows the intrinsic limit, not a practical filter-selection rule). Figure 3 ftp_3_edges_unwrapping.png (a) Ringing at the image border and at a height step, plain FFT vs mirror padding. (b) A step taller than half a fringe: spatial unwrapping assigns the wrong fringe order, so the step comes out 2*pi (one fringe height) too small. Model: small-height (linear) phase-height relation dphi = k*h with k = 2*pi*f0*d/L. The exact Takeda conversion h = L*dphi/(dphi - 2*pi*f0*d) is in phase_to_height_exact() for use with real data. Requires numpy, scipy, matplotlib. Uses scikit-image's 2-D unwrapper if it is installed, otherwise a simple row/column unwrap. """ from __future__ import annotations import argparse import os from dataclasses import dataclass import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from scipy.ndimage import gaussian_filter try: from skimage.restoration import unwrap_phase as _sk_unwrap except Exception: # scikit-image not installed _sk_unwrap = None TWO_PI = 2 * np.pi # ---------------------------------------------------------------------------------- # Setup # ---------------------------------------------------------------------------------- @dataclass class Setup: n: int = 512 # image is n x n pixels period_px: float = 16.0 # fringe period on the reference plane, in camera pixels pixel_mm: float = 0.0625 # object-space pixel size -> period = 1 mm L_mm: float = 500.0 # pupils to reference plane d_mm: float = 150.0 # projector-camera baseline gamma: float = 1.8 # uncorrected projector gamma (makes harmonics) noise: float = 0.004 # camera noise, fraction of full scale seed: int = 3 @property def f0(self): # carrier, cycles per pixel return 1.0 / self.period_px @property def f0_mm(self): # carrier, cycles per mm return 1.0 / (self.period_px * self.pixel_mm) @property def k(self): # rad of phase per mm of height (small-height) return TWO_PI * self.f0_mm * self.d_mm / self.L_mm @property def slope_limit(self): # Takeda: |dh/dx| < L/(3d) return self.L_mm / (3 * self.d_mm) def xy_mm(self): y, x = np.mgrid[0:self.n, 0:self.n] return x * self.pixel_mm, y * self.pixel_mm def phase_to_height_exact(dphi, s: Setup): """Takeda's crossed-axes formula. Sign convention: projector on the side that makes dphi positive for positive h with this code's carrier direction.""" return s.L_mm * dphi / (dphi + TWO_PI * s.f0_mm * s.d_mm) def capture(h_mm, R, s: Setup, rng, delta=0.0): """Camera image of fringes projected on a surface of height h and reflectance R.""" x = np.arange(s.n)[None, :] t = 0.5 + 0.5 * np.cos(TWO_PI * s.f0 * x + s.k * h_mm + delta) I = 0.03 + 0.85 * R * t ** s.gamma I = I + rng.normal(0, s.noise, I.shape) return np.clip(I, 0, 1) # ---------------------------------------------------------------------------------- # FTP core # ---------------------------------------------------------------------------------- def bandpass(n, f0, half_width_frac, kind="hann"): """Filter in FFT order (np.fft.fftfreq), centred on +f0.""" f = np.fft.fftfreq(n) u = (f - f0) / (half_width_frac * f0) if kind == "rect": return (np.abs(u) < 1).astype(float) return np.where(np.abs(u) < 1, 0.5 * (1 + np.cos(np.pi * u)), 0.0) def analytic(img, f0, half_width_frac=0.5, pad="none", kind="hann"): """Steps 2-4: FFT along x, keep the +1 lobe, inverse FFT. pad='mirror' reflects each row so the FFT sees no jump at the borders.""" n = img.shape[1] if pad == "mirror": work = np.concatenate([img[:, ::-1], img, img[:, ::-1]], axis=1) else: work = img H = bandpass(work.shape[1], f0, half_width_frac, kind) a = np.fft.ifft(np.fft.fft(work, axis=1) * H[None, :], axis=1) return a[:, n:2 * n] if pad == "mirror" else a def unwrap2d(p): if _sk_unwrap is not None: return _sk_unwrap(p) u = np.unwrap(p, axis=1) col = np.unwrap(u[:, u.shape[1] // 2]) return u + (col - u[:, u.shape[1] // 2])[:, None] def ftp_height(img_obj, img_ref, s: Setup, half_width_frac=0.5, pad="none", img_obj_pi=None, img_ref_pi=None): """Full FTP (or pi-shift FTP if the pi-shifted frames are given). Returns wrapped phase, unwrapped height (mm), and the analytic signals.""" go = img_obj - img_obj_pi if img_obj_pi is not None else img_obj gr = img_ref - img_ref_pi if img_ref_pi is not None else img_ref ao = analytic(go, s.f0, half_width_frac, pad) ar = analytic(gr, s.f0, half_width_frac, pad) wrapped = np.angle(ao * np.conj(ar)) # step 5-6: carrier removed by the reference un = unwrap2d(wrapped) # step 7 return wrapped, un / s.k, ao, ar # step 8 (linear); see phase_to_height_exact def remove_piston(h_meas, h_true, m): """Spatial unwrapping leaves an unknown 2*pi offset; fix it from the interior mean.""" core = (h_meas - h_true)[m:-m, m:-m] return h_meas - np.median(core) def mean_row_spectrum(img): F = np.fft.fft(img - img.mean(), axis=1) return np.fft.fftshift(np.mean(np.abs(F), axis=0)), np.fft.fftshift(np.fft.fftfreq(img.shape[1])) # ---------------------------------------------------------------------------------- # Scenes # ---------------------------------------------------------------------------------- def scene_bumps(s: Setup): x, y = s.xy_mm() c = s.n * s.pixel_mm / 2 h = (2.5 * np.exp(-((x - c + 6) ** 2 + (y - c) ** 2) / (2 * 4.5 ** 2)) + 1.2 * np.exp(-((x - c - 7) ** 2 + (y - c + 5) ** 2) / (2 * 3.0 ** 2)) - 0.8 * np.exp(-((x - c - 4) ** 2 + (y - c - 8) ** 2) / (2 * 3.0 ** 2))) R = np.full_like(h, 0.9) mark = ((x - c - 7) ** 2 / 3.0 ** 2 + (y - c - 9) ** 2 / 1.2 ** 2) < 1 # dark mark R = np.where(mark, 0.3, R) R = gaussian_filter(R, 1.0) return h, R def scene_ridge(s: Setup, slope): """Ridge along y, Gaussian cross-section in x with peak |dh/dx| = slope.""" x, _ = s.xy_mm() c = s.n * s.pixel_mm / 2 sigma = 3.0 A = slope * sigma * np.sqrt(np.e) # max |d/dx| of A*exp(-x^2/2s^2) = A/(s*sqrt(e)) return A * np.exp(-(x - c) ** 2 / (2 * sigma ** 2)), np.full((s.n, s.n), 0.9) def scene_step(s: Setup, step_mm, tilt=0.0): """Raised box on a plane. A tilt makes the phase at the left and right image borders differ, which is what a real object/scene usually does; the FFT treats each row as periodic, so that difference becomes a jump it has to ring around.""" x, y = s.xy_mm() c = s.n * s.pixel_mm / 2 box = (np.abs(x - c) < 7) & (np.abs(y - c) < 9) return np.where(box, step_mm, 0.0) + tilt * x, np.full((s.n, s.n), 0.9) # ---------------------------------------------------------------------------------- # Figures # ---------------------------------------------------------------------------------- def fig_pipeline(s: Setup, out): rng = np.random.default_rng(s.seed) h, R = scene_bumps(s) I_obj = capture(h, R, s, rng) I_ref = capture(np.zeros_like(h), np.full_like(h, 0.9), s, rng) hw = 0.5 wrapped, h_meas, ao, _ = ftp_height(I_obj, I_ref, s, hw) m = 24 h_meas = remove_piston(h_meas, h, m) err_um = (h_meas - h) * 1000 rms = np.sqrt(np.mean(err_um[m:-m, m:-m] ** 2)) ext = [0, s.n * s.pixel_mm, s.n * s.pixel_mm, 0] fig, axs = plt.subplots(2, 4, figsize=(18, 8.6)) ax = axs[0, 0] ax.imshow(I_obj, cmap="gray", extent=ext) ax.set_title("① captured fringe image g(x,y)\n(gamma 1.8, dark mark lower right)") ax = axs[0, 1] F2 = np.fft.fftshift(np.abs(np.fft.fft2(I_obj - I_obj.mean()))) fx = np.fft.fftshift(np.fft.fftfreq(s.n)) / s.f0 ax.imshow(np.log10(F2 + 1e-3), cmap="magma", extent=[fx[0], fx[-1], fx[-1], fx[0]], aspect="auto") ax.set_xlim(-3.2, 3.2) ax.set_ylim(1.5, -1.5) ax.set_xlabel("f$_x$ / f$_0$") ax.set_ylabel("f$_y$ / f$_0$") ax.set_title("② 2-D spectrum |G| (log)\nlobes at 0, ±f₀, ±2f₀, ±3f₀") ax = axs[0, 2] S_obj, f = mean_row_spectrum(I_obj) S_ref, _ = mean_row_spectrum(I_ref) H = np.fft.fftshift(bandpass(s.n, s.f0, hw)) ax.semilogy(f / s.f0, S_ref, color="0.6", lw=1, label="reference plane") ax.semilogy(f / s.f0, S_obj, color="C0", lw=1.2, label="object") ax2 = ax.twinx() ax2.fill_between(f / s.f0, 0, H, color="C3", alpha=0.18) ax2.plot(f / s.f0, H, color="C3", lw=1.5, label="band-pass H(f)") ax2.set_ylim(0, 1.6) ax2.set_yticks([0, 1]) ax2.set_ylabel("filter", color="C3") for xv, lab in [(0, "zero\norder"), (1, "+1"), (-1, "−1"), (2, "2nd"), (3, "3rd")]: ax.annotate(lab, (xv, S_obj.max() * 2), ha="center", fontsize=8, color="k") ax.set_xlim(-3.5, 3.5) ax.set_ylim(S_obj.max() * 1e-4, S_obj.max() * 8) ax.set_xlabel("f / f₀ (along x)") ax.set_ylabel("mean |FFT| of rows") ax.set_title("③ 1-D spectrum + Hann band-pass (½-width 0.5 f₀)\nobject lobes broaden where slopes are steep") h1, l1 = ax.get_legend_handles_labels() h2, l2 = ax2.get_legend_handles_labels() ax.legend(h1 + h2, l1 + l2, fontsize=8, loc="lower left") ax = axs[0, 3] im = ax.imshow(np.abs(ao) / np.abs(ao).max(), cmap="viridis", extent=ext) ax.set_title("④ |analytic signal| after inverse FFT\n= fringe modulation b/2 (quality map)") fig.colorbar(im, ax=ax, shrink=0.8) ax = axs[1, 0] im = ax.imshow(wrapped, cmap="twilight", extent=ext, vmin=-np.pi, vmax=np.pi) ax.set_title("⑤⑥ wrapped Δφ = arg(ĝ·ĝ₀*)\ncarrier removed by the reference") fig.colorbar(im, ax=ax, shrink=0.8, label="rad") ax = axs[1, 1] im = ax.imshow(h_meas, cmap="viridis", extent=ext) ax.set_title(f"⑦⑧ unwrapped → height (mm)\nk = {s.k:.3f} rad/mm") fig.colorbar(im, ax=ax, shrink=0.8, label="mm") ax = axs[1, 2] v = np.percentile(np.abs(err_um[m:-m, m:-m]), 99) im = ax.imshow(err_um, cmap="RdBu_r", vmin=-v, vmax=v, extent=ext) ax.set_title(f"height error (µm), interior RMS {rms:.1f} µm\nborders, steep flanks, dark mark") fig.colorbar(im, ax=ax, shrink=0.8, label="µm") ax = axs[1, 3] row = int(s.n / 2) xm = np.arange(s.n) * s.pixel_mm ax.plot(xm, h[row], "k", lw=2, label="true") ax.plot(xm, h_meas[row], "C1", lw=1, label="FTP") ax3 = ax.twinx() ax3.plot(xm, err_um[row], "C3", lw=0.8, alpha=0.7) ax3.set_ylabel("error (µm)", color="C3") ax.set_xlabel("x (mm)") ax.set_ylabel("height (mm)") ax.set_title(f"profile through row y = {row * s.pixel_mm:.0f} mm") ax.legend(fontsize=8, loc="upper right") for a in [axs[0, 0], axs[0, 3], axs[1, 0], axs[1, 1], axs[1, 2]]: a.set_xlabel("x (mm)") a.set_ylabel("y (mm)") fig.suptitle(f"Fourier Transform Profilometry, step by step — period {s.period_px:.0f} px = " f"{s.period_px * s.pixel_mm:.1f} mm, L = {s.L_mm:.0f} mm, d = {s.d_mm:.0f} mm", fontsize=13) fig.tight_layout() fig.savefig(out, dpi=130) plt.close(fig) return rms def fig_slope(s: Setup, out): lim = s.slope_limit ratios = np.array([0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.5, 1.8, 2.2, 2.6, 3.0]) widths = [0.3, 0.5, 0.7, 0.9] m = 24 res = {"FTP": [], "π-FTP": []} for r in ratios: h, R = scene_ridge(s, r * lim) for name, gam in [("FTP", s.gamma), ("π-FTP", s.gamma)]: sp = Setup(**{**s.__dict__, "gamma": gam}) best = np.inf for w in widths: rng = np.random.default_rng(s.seed) Io = capture(h, R, sp, rng) Ir = capture(np.zeros_like(h), R, sp, rng) if name.startswith("π"): Iop = capture(h, R, sp, rng, np.pi) Irp = capture(np.zeros_like(h), R, sp, rng, np.pi) _, hm, _, _ = ftp_height(Io, Ir, sp, w, img_obj_pi=Iop, img_ref_pi=Irp) else: _, hm, _, _ = ftp_height(Io, Ir, sp, w) hm = remove_piston(hm, h, m) e = np.sqrt(np.mean(((hm - h) * 1000)[m:-m, m:-m] ** 2)) best = min(best, e) res[name].append(best) fig, axs = plt.subplots(1, 3, figsize=(17, 4.8)) for ax, r in zip(axs[:2], [0.6, 1.8]): h, R = scene_ridge(s, r * lim) rng = np.random.default_rng(s.seed) Io = capture(h, R, s, rng) S, f = mean_row_spectrum(Io) ax.semilogy(f / s.f0, S, "C0", lw=1.2, label="object row spectrum") lo, hi = 1 - (s.d_mm / s.L_mm) * r * lim, 1 + (s.d_mm / s.L_mm) * r * lim for k_ in [1, 2]: ax.axvspan(k_ * lo, k_ * hi, color=f"C{k_}", alpha=0.15, label=f"{'+1' if k_ == 1 else '2nd'} lobe, predicted span") ax.axvspan(-(hi - 1), hi - 1, color="0.5", alpha=0.15, label="zero-order span (same width)") ax.set_xlim(-0.5, 3.2) ax.set_ylim(S.max() * 1e-4, S.max() * 3) ax.set_xlabel("f / f₀") ax.set_title(f"max slope = {r:.1f} × L/(3d)" + (" — lobes separate" if r < 1 else " — lobes overlap")) ax.legend(fontsize=7, loc="upper right") ax = axs[2] for name, sty in zip(res, ["o-", "s-"]): ax.semilogy(ratios, res[name], sty, label=name) ax.axvline(1.0, color="k", ls=":", lw=1) ax.axvline(3.0, color="k", ls=":", lw=1) top = max(max(v) for v in res.values()) ax.text(1.03, top * 0.6, "L/(3d)", fontsize=8) ax.text(2.75, top * 0.6, "L/d", fontsize=8) ax.set_xlabel(f"max |dh/dx| / (L/3d) [L/(3d) = {lim:.2f}]") ax.set_ylabel("RMS height error (µm), best filter width") ax.set_title("Slope limit: FTP vs π-shift FTP") ax.grid(alpha=0.3, which="both") ax.legend(fontsize=8) fig.tight_layout() fig.savefig(out, dpi=130) plt.close(fig) return ratios, res def fig_edges(s: Setup, out): m = 24 s_lin = Setup(**{**s.__dict__, "gamma": 1.0}) fig, axs = plt.subplots(1, 3, figsize=(17, 4.8)) xm = np.arange(s.n) * s.pixel_mm row = s.n // 2 # (a) small step + borders: plain vs mirror padding h, R = scene_step(s_lin, 0.25, tilt=0.05) rng = np.random.default_rng(s.seed) Io = capture(h, R, s_lin, rng) Ir = capture(np.zeros_like(h), R, s_lin, rng) outs = {} for pad in ["none", "mirror"]: _, hm, _, _ = ftp_height(Io, Ir, s_lin, 0.5, pad=pad) outs[pad] = remove_piston(hm, h, m) ax = axs[0] ax.plot(xm, h[row] * 1000, "k", lw=2, label="true (0.25 mm box on a 0.05 tilt)") ax.plot(xm, outs["none"][row] * 1000, "C3", lw=1, label="FTP, plain FFT") ax.plot(xm, outs["mirror"][row] * 1000, "C0", lw=1, label="FTP, mirror-padded rows") ax.set_xlabel("x (mm)") ax.set_ylabel("height (µm)") ax.set_title("Ringing at steps and image borders\n(the band-pass acts as a low-pass on the phase)") ax.legend(fontsize=8) ax.grid(alpha=0.3) ax = axs[1] for pad, c in [("none", "C3"), ("mirror", "C0")]: e = np.abs(outs[pad] - h) * 1000 ax.semilogy(xm, np.mean(e, axis=0) + 1e-3, c, lw=1, label=f"{'plain' if pad == 'none' else 'mirror-padded'}") ax.axvspan(xm[s.n // 2] - 7, xm[s.n // 2] + 7, color="0.85", zorder=-1, label="raised box") ax.set_xlabel("x (mm)") ax.set_ylabel("|error| averaged over rows (µm)") ax.set_title("Where the error lives: borders and step edges\n(≈ one fringe period wide)") ax.legend(fontsize=8) ax.grid(alpha=0.3, which="both") # (b) step taller than half a fringe -> unwrapping ambiguity step = 2.2 # mm; phase step = k*2.2 rad h, R = scene_step(s_lin, step) rng = np.random.default_rng(s.seed) Io = capture(h, R, s_lin, rng) Ir = capture(np.zeros_like(h), R, s_lin, rng) wr, hm, _, _ = ftp_height(Io, Ir, s_lin, 0.5, pad="mirror") hm = hm - np.median(hm[m:-m, m:m + 40]) # zero on the base plane ax = axs[2] ax.plot(xm, h[row], "k", lw=2, label=f"true ({step} mm step = {s.k * step:.2f} rad)") ax.plot(xm, hm[row], "C3", lw=1.2, label="FTP + spatial unwrap") fringe_h = TWO_PI / s.k ax.axhline(step - fringe_h, color="C3", ls=":", lw=1) ax.text(xm[20], step - fringe_h + 0.05, f"true − one fringe height ({fringe_h:.2f} mm)", fontsize=8, color="C3") ax.set_xlabel("x (mm)") ax.set_ylabel("height (mm)") ax.set_title("Step > half a fringe: the fringe order is lost\n(needs a 2nd frequency or temporal unwrapping)") ax.legend(fontsize=8, loc="upper left") ax.grid(alpha=0.3) fig.tight_layout() fig.savefig(out, dpi=130) plt.close(fig) def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--outdir", default="ftp_walkthrough_results") ap.add_argument("--gamma", type=float, default=1.8) ap.add_argument("--period", type=float, default=16.0, help="fringe period in pixels") args = ap.parse_args() os.makedirs(args.outdir, exist_ok=True) s = Setup(gamma=args.gamma, period_px=args.period, pixel_mm=1.0 / args.period) print(f"carrier f0 = {s.f0:.4f} cyc/px = {s.f0_mm:.2f} cyc/mm, k = {s.k:.3f} rad/mm, " f"slope limit L/(3d) = {s.slope_limit:.2f}, one fringe of height = {TWO_PI / s.k:.2f} mm") print("unwrapper:", "scikit-image" if _sk_unwrap else "row/column fallback") rms = fig_pipeline(s, os.path.join(args.outdir, "ftp_1_pipeline.png")) print(f"pipeline figure: interior RMS error {rms:.1f} um") ratios, res = fig_slope(s, os.path.join(args.outdir, "ftp_2_slope_limit.png")) print("slope sweep (RMS um):") print(" slope/(L/3d) " + " ".join(f"{r:6.1f}" for r in ratios)) for k, v in res.items(): print(f" {k:<22s}" + " ".join(f"{x:6.1f}" for x in v)) fig_edges(s, os.path.join(args.outdir, "ftp_3_edges_unwrapping.png")) print(f"figures written to {os.path.abspath(args.outdir)}") if __name__ == "__main__": main()