diff --git a/cortex/export/headless.py b/cortex/export/headless.py
index a8eadb458..855cd2792 100644
--- a/cortex/export/headless.py
+++ b/cortex/export/headless.py
@@ -99,6 +99,60 @@ def _wait_for_viewer_loaded(handle, timeout: float = 60.0) -> None:
 # --------------------------------------------------------------------------- #
 
 
+#: Browser messages that mean WebGL itself failed, as opposed to any of the
+#: unrelated javascript a page may log. Kept deliberately narrow.
+#:
+#: Checking *all* browser errors is not usable: a healthy viewer already logs a
+#: console.error for the Leap Motion websocket it cannot reach
+#: (ws://127.0.0.1:6437), so anything failing on "any error" fails every run.
+#: That is why the existing assertions in the test suite filter to [pageerror].
+#:
+#: These two are unambiguous failures:
+#:   - "Could not initialise shader" is three.js reporting a shader that compiled
+#:     but failed to *link*. It is a console.error, not an exception, so nothing
+#:     raises and the render silently comes back as a blank canvas. This is
+#:     exactly how Vertex2D broke (gh-714).
+#:   - "Error creating WebGL context" is thrown, so it also arrives as a
+#:     [pageerror]; matching the text covers it either way.
+#:
+#: gl.getProgramInfoLog / gl.getShaderInfoLog warnings are deliberately NOT
+#: listed. Drivers emit those benignly, so matching them would reintroduce the
+#: false positives this list exists to avoid.
+# A link failure makes three.js r69 emit three console.error lines together,
+# from one unconditional block (resources/js/three.js, ~line 24578):
+#
+#     THREE.WebGLProgram: Could not initialise shader.
+#     gl.VALIDATE_STATUS false
+#     gl.getError() 0
+#
+# Matching only the first is enough -- the three are inseparable, so there is no
+# case where the later ones appear alone.
+#
+# Do not be tempted to match "gl.VALIDATE_STATUS false" instead or as well. It
+# does not mean validation failed: nothing in this codebase ever calls
+# gl.validateProgram(), so VALIDATE_STATUS is simply false because validation was
+# never run. It reads like a diagnosis and is not one. Likewise "gl.getError() 0"
+# is the *absence* of a GL error -- the failure is in the link status, which
+# three.js has already reported on the line above.
+#: How often the worker thread calls into Playwright to dispatch queued browser
+#: events. Console messages are only delivered during such a call, so this is
+#: the granularity at which ``browser_errors`` becomes current.
+EVENT_PUMP_INTERVAL = 0.25
+
+WEBGL_FAILURE_PATTERNS = (
+    "THREE.WebGLProgram: Could not initialise shader",
+    "Error creating WebGL context",
+)
+
+
+def webgl_failures(browser_errors: list[str]) -> list[str]:
+    """Return only those browser messages that indicate WebGL itself failed."""
+    return [
+        e for e in browser_errors
+        if any(pattern in e for pattern in WEBGL_FAILURE_PATTERNS)
+    ]
+
+
 class _PlaywrightThread:
     """Manages the Playwright lifecycle on a private daemon thread.
 
@@ -230,7 +284,23 @@ def _worker(self) -> None:
             return
 
         # Keep the thread (and therefore Playwright) alive until shutdown.
-        self._shutdown_event.wait()
+        #
+        # Parking in wait() is not enough. Playwright's sync API only dispatches
+        # queued events while something is calling into it, so a thread that
+        # blocks for the viewer's whole lifetime leaves console messages sitting
+        # undelivered in the driver until _cleanup() closes the page. Anything
+        # reading browser_errors before then (like save_3d_views does, once per
+        # view) would see only what arrived during page load, and would miss
+        # any OpenGL errors (e.g. shader link failures).
+        #
+        # So poll instead, making a cheap round-trip into Playwright each time,
+        # which pumps the event loop and delivers whatever has queued up.
+        while not self._shutdown_event.wait(EVENT_PUMP_INTERVAL):
+            try:
+                self._page.evaluate("0")
+            except Exception:
+                # Page or browser is gone; _cleanup() handles the rest.
+                break
         self._cleanup()
 
     # -- Playwright event handlers (called on the worker thread) ---------- #
diff --git a/cortex/export/save_views.py b/cortex/export/save_views.py
index e240afb90..2ebaceda2 100644
--- a/cortex/export/save_views.py
+++ b/cortex/export/save_views.py
@@ -196,6 +196,25 @@ def save_3d_views(
                 )
             time.sleep(1)
 
+            # A WebGL failure does not raise on its own: three.js reports a
+            # shader that failed to link via console.error, so the png is still
+            # written and simply comes back blank, which is indistinguishable
+            # from a successful render. Fail loudly instead. Only genuine WebGL
+            # failures are matched -- a healthy viewer logs unrelated console
+            # errors (the Leap Motion websocket), so checking every browser
+            # error would fail every run.
+            pw_thread = getattr(handle, "_pw_thread", None)
+            if pw_thread is not None:
+                from cortex.export.headless import webgl_failures
+
+                failures = webgl_failures(pw_thread.browser_errors)
+                if failures:
+                    raise RuntimeError(
+                        f"WebGL failed while rendering {view_name!r}/{surface!r}; "
+                        f"{file_name!r} is likely blank.\n  "
+                        + "\n  ".join(sorted(set(failures)))
+                    )
+
             # Trim transparent edges
             if trim:
                 try:
diff --git a/cortex/tests/reference_images/README.md b/cortex/tests/reference_images/README.md
new file mode 100644
index 000000000..5da2b0a55
--- /dev/null
+++ b/cortex/tests/reference_images/README.md
@@ -0,0 +1,206 @@
+# Reference images
+
+Stored renders that `cortex/tests/test_visual_regression.py` asserts against, so
+a change in rendering output fails a test instead of needing to be spotted by
+eye.
+
+## `alpha_dataviews/`
+
+Twelve images: each of the six public dataview classes (`Volume`, `Vertex`,
+`Volume2D`, `Vertex2D`, `VolumeRGB`, `VertexRGB`) rendered through both paths --
+`quickflat_*` via `cortex.quickshow` (matplotlib) and `webgl_*` via
+`cortex.export.plot_panels` (headless WebGL).
+
+Between them they exercise every way pycortex encodes alpha: `Volume`/`Vertex` are
+a no-alpha baseline, `Volume2D`/`Vertex2D` use the 2D alpha colormap
+`RdBu_r_alpha`, and `VolumeRGB`/`VertexRGB` use the native `alpha=` keyword. All
+six also composite the curvature underlay.
+
+## `nan_dataviews/`
+
+Twelve images, laid out exactly as `alpha_dataviews/`: the same six dataview
+classes, but with NaNs standing in for "missing data" over roughly half of the
+volume/surface (the primary data channel, not the alpha channel). Both
+renderers are expected to draw those elements as fully transparent, falling
+through to the curvature underlay, rather than mapping NaN through the
+colormap as if it were a real value.
+
+## `nan_alpha_dataviews/`
+
+Four images: `VolumeRGB` and `VertexRGB` -- the only dataviews taking an explicit
+`alpha=` -- with the NaNs in the **alpha map** rather than in the data. That is a
+distinct path, since alpha is not colour-mapped but used directly as a blend
+weight, so the NaN reaches the compositing arithmetic rather than a colormap
+lookup. Both renderers currently draw those elements fully transparent.
+
+This behaviour is not settled: work to unify NaN and alpha handling across
+quickflat, WebGL and the RGB dataviews (`cb976270`, not in this branch's history)
+changes how the surviving RGB is blended without changing the transparency
+itself. Expect these four to need regenerating if that lands.
+
+## `nonflat_views/`
+
+Four images, and the only ones in this directory that are **webgl-only**:
+`Volume` and `Vertex` rendered on the inflated and fiducial surfaces at
+`lateral_pivot`. There is no quickflat counterpart because `cortex.quickshow`
+renders flatmaps and nothing else, so these get a reference check and no
+cross-renderer check -- an absence forced by the renderers, not a choice.
+
+Everything else here is a flatmap, which left the 3D views covered only by
+`test_webgl_headless.py`'s smoke tests, and those assert the file is over 1000
+bytes: enough to catch a render that never happened, not one that came out
+wrong. Both shader paths appear because the flatmap turns out not to be
+representative of how they behave -- the two known webgl lighting bugs
+reproduce on flatmaps only.
+
+They are also rendered with `save_3d_views` rather than `plot_panels`, so the
+stored pixels are the browser screenshot itself rather than matplotlib's
+interpolation of it (`plot_panels` composes the screenshot into a figure,
+upsampling it). The flatmap suites cannot do the same: `save_3d_views` output is
+~44% transparent while `quickshow`'s figure is fully opaque, so their
+cross-renderer check would be diffing a white background against an empty one.
+
+These keep pycortex's **default** thresholded curvature, unlike the flatmap
+suites. Those un-threshold it to reduce cross-renderer disagreement, which is
+not a concern without a second renderer, so the default is used instead and
+recovers coverage of the default curvature path that the flatmap references
+deliberately give up.
+
+## Within- vs cross-renderer checks
+
+The three flatmap tests -- `..._alpha_dataviews`, `..._nan_dataviews` and
+`..._nan_alpha_dataviews` -- check every render two ways. (`..._nonflat_views`
+gets the within-renderer check only, for the reason given above.)
+
+- **Within-renderer**: `quickflat_*`/`webgl_*` against their own stored
+  reference above, at a tight tolerance. Four criteria, all of which must pass:
+  `MAX_MEAN_ABS_DIFF`, `MAX_FRACTION_DIFFERING`,
+  `MAX_FRACTION_GROSSLY_DIFFERING` and `MAX_SSIM_LOSS`. They are complementary,
+  not redundant -- the mean is blind to a change that moves a few pixels a long
+  way, the gross fraction is blind to a broad low-amplitude shift, and SSIM is
+  blind to a channel permutation. See the test file for the calibration. This is
+  what catches a regression in one renderer's pipeline.
+- **Cross-renderer**: quickflat vs webgl for the *same* dataview, diffed
+  directly against each other with no stored fixture, after the affine
+  correction below puts them on a common grid, at a looser tolerance
+  (`CROSS_MAX_MEAN_ABS_DIFF` / `CROSS_MAX_FRACTION_DIFFERING`). matplotlib and
+  Three.js still differ in colormap sampling and anti-aliasing, so this can't
+  use the within-renderer tolerance, but it still catches a large disagreement
+  (wrong colormap, dropped alpha, swapped channels) between the two paths even
+  if a stale or wrongly-regenerated reference would otherwise hide it.
+
+## The cross-renderer affine correction
+
+The two renderers do not share a pixel grid: each computes its own trim/extent,
+so webgl's flatmap lands in quickflat's frame off by a fixed anisotropic scale
+(x 0.9690, y 0.9365) plus a ~7 px translation. `_check_cross_renderer` undoes
+that before diffing, via the `CROSS_RENDERER_*` constants in
+`cortex/tests/test_visual_regression.py`.
+
+Measured over the fourteen stored quickflat/webgl pairs, the correction takes
+mean|diff| from 4.83-9.42 down to 1.06-2.01, and the fraction of pixels past
+`CROSS_DIFF_THRESHOLD` from 2.99-8.49% down to 0.10-1.64%. On curvature-only
+content -- what the fit itself uses, and the most exposed case, since there is
+no opaque data covering the sulcal texture -- the uncorrected figure is ~19
+rather than ~9.
+
+Note this is *not* perspective distortion, despite webgl rendering through a
+`THREE.PerspectiveCamera` (FOV 45°, `axes3d.js`). The flatmap view looks
+straight down at a planar surface, where a pinhole projection degenerates to a
+uniform scale — fitting a full homography returns a projective row of
+`[~0, ~0, 1]` and ~0° rotation, i.e. it collapses to the affine. Chasing a
+reverse perspective projection is a dead end.
+
+To re-derive the constants:
+
+```
+uv run --with opencv-python-headless \
+    python cortex/tests/reference_images/fit_cross_renderer_affine.py
+```
+
+It renders curvature-only content (a `VolumeRGB` with `alpha=0` everywhere, so
+only the curvature underlay is drawn — the one layer both paths composite
+identically, which isolates the coordinate-frame mismatch from any colormap
+difference), fits webgl → quickflat with `cv2.findTransformECC`, and prints the
+constants ready to paste. OpenCV is only needed for the fit, so it is not a
+project dependency — hence `--with`.
+
+Run it if the cross-renderer check starts failing broadly with a mean|diff|
+around 5-9 rather than 1-2, which is what losing the correction looks like —
+e.g. after a change to `plot_panels`' figure composition, to quickflat's
+`height=`/`dpi=`, or to either renderer's trim logic. Note "broadly", not
+"everywhere": dropping the correction fails 13 of the 14 pairs, and the survivor
+(`alpha_dataviews/Volume2D`, at 4.83 and 2.99%) clears both limits. A single
+passing pair is not evidence the correction is intact. The script's render
+settings mirror the test's; if the test's change, update the script to match or
+the fit will not apply. It also writes `fit_affine_residual.png`, which should
+show a faint sulcal outline only — broad structure there means the fit didn't
+take.
+
+## Provenance
+
+Generated on this branch, on top of `claude/issue-714-irf5g0` (the Vertex2D
+blank-render fix).
+
+That fix is load-bearing for four of the 32 images and no others. Regenerating
+the whole set against `main` instead -- reverting only the two JS files the fix
+touches, `mriview_surface.js` and `shaderlib.js` -- reproduces 28 of them
+**byte-for-byte**. The four it cannot produce are the `Vertex2D` pairs in
+`alpha_dataviews/` and `nan_dataviews/`: on `main` that render comes out blank,
+`save_3d_views` raises, and neither the webgl nor the quickflat image is
+written. So the rest of the suite can be regenerated on either commit, and a
+diff against `main`-generated references isolates exactly what the fix changed.
+
+They can **no longer be reproduced at `5af26a86`**, which earlier revisions of
+this file named. Regenerating there and testing against the current tree fails
+the six vertex-based cases: `#679` moved every pure-Vertex flatmap's trim box
+from 594 to 596px, and that is intended new geometry which its follow-up fix
+does not revert. The failures are a genuine content change, not a resampling
+artifact -- optimal affine alignment removes only 12% of the difference, and the
+difference is concentrated in the interior rather than at the silhouette.
+
+Every *quickflat* reference that existed then -- the twelve outside
+`nan_alpha_dataviews/`, which came later -- is byte-identical at either commit;
+only the webgl ones moved.
+
+| | |
+| --- | --- |
+| chromium | 151.0.7922.34 (headless shell, SwiftShader software rendering) |
+| matplotlib | 3.11.1 |
+| pillow | 12.3.0 |
+
+## Format
+
+Lossless WebP (`method=6`, `quality=100`, `exact=True`): bit-exact after decode,
+and 59% the size of optimized PNG (1408 KiB versus 2368 KiB for the set of 32).
+
+## Distribution
+
+These are test fixtures with no runtime use, so they are **kept out of the wheel**
+(`exclude_package_data` in `setup.py`) and **kept in the source tarball**
+(`MANIFEST.in`'s `recursive-include cortex *`). A build from source can therefore
+run the test; a `pip install` does not carry ~1.4 MiB of fixtures into
+site-packages for data no user will read.
+
+The test skips, rather than fails, when the images are absent, so a test run
+against an installed wheel degrades gracefully.
+
+## Regenerating
+
+The renders are deterministic: repeated runs on one machine produce bit-identical
+output, including the WebGL ones under software rendering. They are, however,
+coupled to the Chromium and matplotlib builds above, so a browser or matplotlib
+upgrade can shift anti-aliasing and rasterization slightly. The test's tolerances
+absorb that; if an upgrade moves output beyond them, inspect the `diff_*.png`
+files the failure writes, confirm the change is cosmetic, then:
+
+```
+REGENERATE_REFERENCE_IMAGES=1 pytest cortex/tests/test_visual_regression.py
+```
+
+(This regenerates all four directories in one run, since all the tests live in
+that file. It does not touch the affine correction, which is fitted separately —
+see above.)
+
+Review the resulting diff before committing -- regenerating is how a real
+regression gets silently blessed.
diff --git a/cortex/tests/reference_images/alpha_dataviews/quickflat_Vertex.webp b/cortex/tests/reference_images/alpha_dataviews/quickflat_Vertex.webp
new file mode 100644
index 000000000..cf0ce02d9
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/quickflat_Vertex.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/quickflat_Vertex2D.webp b/cortex/tests/reference_images/alpha_dataviews/quickflat_Vertex2D.webp
new file mode 100644
index 000000000..0ce9b0606
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/quickflat_Vertex2D.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/quickflat_VertexRGB.webp b/cortex/tests/reference_images/alpha_dataviews/quickflat_VertexRGB.webp
new file mode 100644
index 000000000..359042ceb
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/quickflat_VertexRGB.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/quickflat_Volume.webp b/cortex/tests/reference_images/alpha_dataviews/quickflat_Volume.webp
new file mode 100644
index 000000000..e69478f6a
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/quickflat_Volume.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/quickflat_Volume2D.webp b/cortex/tests/reference_images/alpha_dataviews/quickflat_Volume2D.webp
new file mode 100644
index 000000000..f892c97a4
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/quickflat_Volume2D.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/quickflat_VolumeRGB.webp b/cortex/tests/reference_images/alpha_dataviews/quickflat_VolumeRGB.webp
new file mode 100644
index 000000000..ab0e47684
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/quickflat_VolumeRGB.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/webgl_Vertex.webp b/cortex/tests/reference_images/alpha_dataviews/webgl_Vertex.webp
new file mode 100644
index 000000000..70e7896c9
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/webgl_Vertex.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/webgl_Vertex2D.webp b/cortex/tests/reference_images/alpha_dataviews/webgl_Vertex2D.webp
new file mode 100644
index 000000000..55f34e3ed
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/webgl_Vertex2D.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/webgl_VertexRGB.webp b/cortex/tests/reference_images/alpha_dataviews/webgl_VertexRGB.webp
new file mode 100644
index 000000000..abda89bfc
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/webgl_VertexRGB.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/webgl_Volume.webp b/cortex/tests/reference_images/alpha_dataviews/webgl_Volume.webp
new file mode 100644
index 000000000..2e93c7efa
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/webgl_Volume.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/webgl_Volume2D.webp b/cortex/tests/reference_images/alpha_dataviews/webgl_Volume2D.webp
new file mode 100644
index 000000000..1ae994708
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/webgl_Volume2D.webp differ
diff --git a/cortex/tests/reference_images/alpha_dataviews/webgl_VolumeRGB.webp b/cortex/tests/reference_images/alpha_dataviews/webgl_VolumeRGB.webp
new file mode 100644
index 000000000..d40cb591c
Binary files /dev/null and b/cortex/tests/reference_images/alpha_dataviews/webgl_VolumeRGB.webp differ
diff --git a/cortex/tests/reference_images/fit_cross_renderer_affine.py b/cortex/tests/reference_images/fit_cross_renderer_affine.py
new file mode 100644
index 000000000..ec96448e1
--- /dev/null
+++ b/cortex/tests/reference_images/fit_cross_renderer_affine.py
@@ -0,0 +1,137 @@
+"""Re-derive the cross-renderer affine transform used by test_visual_regression.py.
+
+quickflat and webgl do not agree on a common pixel grid: each computes its own
+trim/extent, so webgl's flatmap lands in quickflat's frame off by a fixed
+anisotropic scale plus a translation. ``_check_cross_renderer`` corrects for
+that before diffing, using the ``CROSS_RENDERER_*`` constants in
+``cortex/tests/test_visual_regression.py``. This script recomputes them.
+
+Run it if the cross-renderer check starts failing everywhere at once with a
+large mean|diff| (~19 rather than ~2), which is what losing the correction looks
+like -- e.g. after a change to ``plot_panels``' figure composition, to
+quickflat's ``height=``/``dpi=``, or to either renderer's trim logic.
+
+    uv run --with opencv-python-headless \
+        python cortex/tests/reference_images/fit_cross_renderer_affine.py
+
+OpenCV is only needed for the fit, so it is deliberately not a project
+dependency -- hence ``--with``.
+
+Method: render curvature-only content (a VolumeRGB with alpha=0 everywhere, so
+the data layer is fully transparent and only the curvature underlay is drawn)
+through both renderers. Curvature is the one layer both paths composite
+identically, so the fit measures the coordinate-frame mismatch rather than any
+dataview-specific colormap difference. Then fit webgl -> quickflat with
+``cv2.findTransformECC``, first as a full homography to confirm there is no real
+projective component, then as an affine to read off the constants.
+
+Note on conventions: ``cv2.warpAffine(..., WARP_INVERSE_MAP)`` and PIL's
+``Image.AFFINE`` use the same inverse mapping (output pixel -> source coords),
+so the fitted matrix rows drop straight into PIL's coefficient tuple in
+``_check_cross_renderer`` with no inversion or transpose.
+"""
+
+import numpy as np
+
+import cortex
+import cortex.export
+
+SUBJ = "S1"
+XFMNAME = "fullhead"
+
+# Must match what test_visual_regression.py's _render_and_check_dataview does,
+# or the fitted transform will not apply to the renders the test produces.
+QUICKFLAT_HEIGHT = 256
+QUICKFLAT_DPI = 80
+WEBGL_FIGSIZE = (6, 3)
+WEBGL_WINDOWSIZE = (512, 384)
+WEBGL_SLEEP = 10
+
+FLATMAP_PANEL = [
+    {"extent": [0.0, 0.0, 1.0, 1.0], "view": {"angle": "flatmap", "surface": "flatmap"}}
+]
+
+
+def _curvature_only_dataview():
+    """A dataview whose data layer is fully transparent, leaving only curvature."""
+    zeros = np.zeros((31, 100, 100))
+    chan = lambda: cortex.Volume(zeros, SUBJ, XFMNAME, vmin=0, vmax=1)
+    return cortex.VolumeRGB(
+        chan(), chan(), chan(), SUBJ, XFMNAME, alpha=chan(),
+    )
+
+
+def main():
+    import cv2
+    import matplotlib.pyplot as plt
+    from PIL import Image
+
+    view = _curvature_only_dataview()
+
+    qf_path = "fit_affine_quickflat.png"
+    fig = cortex.quickshow(
+        view, with_curvature=True, with_rois=False, with_labels=False,
+        with_colorbar=False, with_sulci=False, with_borders=False,
+        height=QUICKFLAT_HEIGHT,
+    )
+    fig.savefig(qf_path, bbox_inches="tight", pad_inches=0, dpi=QUICKFLAT_DPI)
+    plt.close(fig)
+
+    wg_path = "fit_affine_webgl.png"
+    fig = cortex.export.plot_panels(
+        view, panels=FLATMAP_PANEL, figsize=WEBGL_FIGSIZE,
+        windowsize=WEBGL_WINDOWSIZE, save_name=wg_path, sleep=WEBGL_SLEEP,
+        viewer_params=dict(labels_visible=[], overlays_visible=[]),
+        headless=True,
+    )
+    plt.close(fig)
+
+    qf_im = Image.open(qf_path).convert("L")
+    wg_im = Image.open(wg_path).convert("L")
+    print(f"quickflat {qf_im.size}, webgl {wg_im.size}")
+
+    # Fit in quickflat's frame, which is what _check_cross_renderer diffs in.
+    qf = np.asarray(qf_im).astype(np.float32) / 255.0
+    wg = np.asarray(wg_im.resize(qf_im.size, Image.BILINEAR)).astype(np.float32) / 255.0
+    h, w = qf.shape
+    criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 5000, 1e-6)
+
+    print(f"\nmean|diff| with a plain resize and no correction: "
+          f"{np.abs(qf - wg).mean() * 255:.3f}")
+
+    # Homography first, purely as a check that there is nothing projective to
+    # undo. webgl uses a THREE.PerspectiveCamera, but the flatmap view looks
+    # straight down at a planar surface, where a pinhole projection degenerates
+    # to a uniform scale -- so the bottom row should come out as [0, 0, 1].
+    hom, _ = np.eye(3, 3, dtype=np.float32), None
+    _, hom = cv2.findTransformECC(qf, wg, hom, cv2.MOTION_HOMOGRAPHY, criteria, None, 5)
+    print(f"\nhomography projective row: [{hom[2, 0]:.3e}, {hom[2, 1]:.3e}, {hom[2, 2]:.3f}]")
+    print("  near [0, 0, 1] => no real perspective distortion; affine is sufficient")
+
+    aff = np.eye(2, 3, dtype=np.float32)
+    cc, aff = cv2.findTransformECC(qf, wg, aff, cv2.MOTION_AFFINE, criteria, None, 5)
+    aligned = cv2.warpAffine(
+        wg, aff, (w, h), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP
+    )
+    mask = aligned != 0
+    print(f"\naffine fit: correlation={cc:.4f}")
+    print(f"mean|diff| after correction: {np.abs(qf[mask] - aligned[mask]).mean() * 255:.3f}")
+
+    (a, b, tx), (c, d, ty) = aff[0], aff[1]
+    rotation = np.degrees(np.arctan2(c, a))
+    print(f"rotation={rotation:.4f} deg (expect ~0), shear terms b={b:.2e} c={c:.2e}")
+
+    print("\nPaste into cortex/tests/test_visual_regression.py:\n")
+    print(f"CROSS_RENDERER_SCALE_X = {np.hypot(a, c):.4f}")
+    print(f"CROSS_RENDERER_SCALE_Y = {np.hypot(b, d):.4f}")
+    print(f"CROSS_RENDERER_TRANSLATE_X_FRAC = {tx:.4f} / {w}")
+    print(f"CROSS_RENDERER_TRANSLATE_Y_FRAC = {ty:.4f} / {h}")
+
+    amp = np.clip(np.abs(qf - aligned) * 255 * 4, 0, 255).astype("uint8")
+    Image.fromarray(amp).save("fit_affine_residual.png")
+    print("\nresidual (4x amplified) written to fit_affine_residual.png -- expect a")
+    print("faint sulcal outline only; broad structure means the fit did not take")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/cortex/tests/reference_images/nan_alpha_dataviews/quickflat_VertexRGB.webp b/cortex/tests/reference_images/nan_alpha_dataviews/quickflat_VertexRGB.webp
new file mode 100644
index 000000000..512f51d13
Binary files /dev/null and b/cortex/tests/reference_images/nan_alpha_dataviews/quickflat_VertexRGB.webp differ
diff --git a/cortex/tests/reference_images/nan_alpha_dataviews/quickflat_VolumeRGB.webp b/cortex/tests/reference_images/nan_alpha_dataviews/quickflat_VolumeRGB.webp
new file mode 100644
index 000000000..77f3ca1ca
Binary files /dev/null and b/cortex/tests/reference_images/nan_alpha_dataviews/quickflat_VolumeRGB.webp differ
diff --git a/cortex/tests/reference_images/nan_alpha_dataviews/webgl_VertexRGB.webp b/cortex/tests/reference_images/nan_alpha_dataviews/webgl_VertexRGB.webp
new file mode 100644
index 000000000..92a4604a5
Binary files /dev/null and b/cortex/tests/reference_images/nan_alpha_dataviews/webgl_VertexRGB.webp differ
diff --git a/cortex/tests/reference_images/nan_alpha_dataviews/webgl_VolumeRGB.webp b/cortex/tests/reference_images/nan_alpha_dataviews/webgl_VolumeRGB.webp
new file mode 100644
index 000000000..7a3950399
Binary files /dev/null and b/cortex/tests/reference_images/nan_alpha_dataviews/webgl_VolumeRGB.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/quickflat_Vertex.webp b/cortex/tests/reference_images/nan_dataviews/quickflat_Vertex.webp
new file mode 100644
index 000000000..07c55504a
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/quickflat_Vertex.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/quickflat_Vertex2D.webp b/cortex/tests/reference_images/nan_dataviews/quickflat_Vertex2D.webp
new file mode 100644
index 000000000..a71a710ba
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/quickflat_Vertex2D.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/quickflat_VertexRGB.webp b/cortex/tests/reference_images/nan_dataviews/quickflat_VertexRGB.webp
new file mode 100644
index 000000000..32754d2d6
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/quickflat_VertexRGB.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/quickflat_Volume.webp b/cortex/tests/reference_images/nan_dataviews/quickflat_Volume.webp
new file mode 100644
index 000000000..881b85ba1
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/quickflat_Volume.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/quickflat_Volume2D.webp b/cortex/tests/reference_images/nan_dataviews/quickflat_Volume2D.webp
new file mode 100644
index 000000000..640c93869
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/quickflat_Volume2D.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/quickflat_VolumeRGB.webp b/cortex/tests/reference_images/nan_dataviews/quickflat_VolumeRGB.webp
new file mode 100644
index 000000000..003ec55e7
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/quickflat_VolumeRGB.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/webgl_Vertex.webp b/cortex/tests/reference_images/nan_dataviews/webgl_Vertex.webp
new file mode 100644
index 000000000..37c4a8d16
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/webgl_Vertex.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/webgl_Vertex2D.webp b/cortex/tests/reference_images/nan_dataviews/webgl_Vertex2D.webp
new file mode 100644
index 000000000..f38837b51
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/webgl_Vertex2D.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/webgl_VertexRGB.webp b/cortex/tests/reference_images/nan_dataviews/webgl_VertexRGB.webp
new file mode 100644
index 000000000..931df43e1
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/webgl_VertexRGB.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/webgl_Volume.webp b/cortex/tests/reference_images/nan_dataviews/webgl_Volume.webp
new file mode 100644
index 000000000..8a5e94858
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/webgl_Volume.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/webgl_Volume2D.webp b/cortex/tests/reference_images/nan_dataviews/webgl_Volume2D.webp
new file mode 100644
index 000000000..49679ce3a
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/webgl_Volume2D.webp differ
diff --git a/cortex/tests/reference_images/nan_dataviews/webgl_VolumeRGB.webp b/cortex/tests/reference_images/nan_dataviews/webgl_VolumeRGB.webp
new file mode 100644
index 000000000..db0dc0f44
Binary files /dev/null and b/cortex/tests/reference_images/nan_dataviews/webgl_VolumeRGB.webp differ
diff --git a/cortex/tests/reference_images/nonflat_views/webgl_fiducial_lateral_pivot_Vertex.webp b/cortex/tests/reference_images/nonflat_views/webgl_fiducial_lateral_pivot_Vertex.webp
new file mode 100644
index 000000000..3719de400
Binary files /dev/null and b/cortex/tests/reference_images/nonflat_views/webgl_fiducial_lateral_pivot_Vertex.webp differ
diff --git a/cortex/tests/reference_images/nonflat_views/webgl_fiducial_lateral_pivot_Volume.webp b/cortex/tests/reference_images/nonflat_views/webgl_fiducial_lateral_pivot_Volume.webp
new file mode 100644
index 000000000..b73232e6b
Binary files /dev/null and b/cortex/tests/reference_images/nonflat_views/webgl_fiducial_lateral_pivot_Volume.webp differ
diff --git a/cortex/tests/reference_images/nonflat_views/webgl_inflated_lateral_pivot_Vertex.webp b/cortex/tests/reference_images/nonflat_views/webgl_inflated_lateral_pivot_Vertex.webp
new file mode 100644
index 000000000..998552e3b
Binary files /dev/null and b/cortex/tests/reference_images/nonflat_views/webgl_inflated_lateral_pivot_Vertex.webp differ
diff --git a/cortex/tests/reference_images/nonflat_views/webgl_inflated_lateral_pivot_Volume.webp b/cortex/tests/reference_images/nonflat_views/webgl_inflated_lateral_pivot_Volume.webp
new file mode 100644
index 000000000..ebf1e63d2
Binary files /dev/null and b/cortex/tests/reference_images/nonflat_views/webgl_inflated_lateral_pivot_Volume.webp differ
diff --git a/cortex/tests/test_export.py b/cortex/tests/test_export.py
index cbb0aa9e4..a80102e0f 100644
--- a/cortex/tests/test_export.py
+++ b/cortex/tests/test_export.py
@@ -81,4 +81,91 @@ def test_plot_panels_headless():
         # provided, should have written the file to disk.
         assert fig is not None
         assert os.path.isfile(save_name)
-        assert os.path.getsize(save_name) > 0
\ No newline at end of file
+        assert os.path.getsize(save_name) > 0
+
+def test_webgl_failures_filters_to_real_failures():
+    """Only genuine WebGL failures match; ordinary console noise does not.
+
+    The filter has to stay narrow. A healthy viewer always logs a console.error
+    for the Leap Motion websocket it cannot reach, and the driver emits shader
+    warnings, so anything broader would fail every run.
+    """
+    from cortex.export.headless import webgl_failures
+
+    noise = [
+        "[console.error] WebSocket connection to 'ws://127.0.0.1:6437/v6.json' "
+        "failed: Error in connection establishment: net::ERR_CONNECTION_REFUSED",
+        "[console.warning] THREE.WebGLShader: gl.getShaderInfoLog() WARNING: 0:87",
+        "[console.warning] [.WebGL-0x21b400157a00]GL Driver Message (OpenGL, Perf)",
+        # Emitted alongside a link failure, but meaningless alone: nothing ever
+        # calls gl.validateProgram(), so VALIDATE_STATUS is false because
+        # validation never ran, not because it failed.
+        "[console.error] gl.VALIDATE_STATUS false",
+        "[console.error] gl.getError() 0",
+    ]
+    assert webgl_failures(noise) == []
+
+    link_failure = "[console.error] THREE.WebGLProgram: Could not initialise shader."
+    context_failure = "[pageerror] Error creating WebGL context."
+    assert webgl_failures(noise + [link_failure]) == [link_failure]
+    assert webgl_failures(noise + [context_failure]) == [context_failure]
+
+
+def test_browser_errors_are_delivered_before_teardown():
+    """Console messages must arrive while the viewer is alive, not at teardown.
+
+    Playwright's sync API only dispatches queued events while something is
+    calling into it. The worker thread used to park in ``_shutdown_event.wait()``
+    for the viewer's whole lifetime, so messages emitted during rendering sat
+    undelivered until ``_cleanup()`` closed the page -- which made
+    ``save_3d_views``' failure check dead code, since it reads per view.
+
+    Rendering therefore has to grow ``browser_errors`` beyond whatever arrived
+    during page load. If the pump in ``_PlaywrightThread._run`` is removed, the
+    count stays at its page-load value and this fails.
+    """
+    import time
+
+    vol = cortex.Volume(np.random.randn(*volshape), subj, xfmname)
+    with cortex.export.headless_viewer(vol, viewer_params={}) as handle:
+        after_load = len(handle._pw_thread.browser_errors)
+        with tempfile.TemporaryDirectory() as tmpdir:
+            outfile = os.path.join(tmpdir, "pump.png")
+            handle.getImage(outfile, (512, 384))
+            for _ in range(300):
+                if os.path.exists(outfile) and os.path.getsize(outfile) > 0:
+                    break
+                time.sleep(0.1)
+        time.sleep(2)
+        during_render = len(handle._pw_thread.browser_errors)
+
+    assert during_render > after_load, (
+        "browser_errors did not grow while the viewer was alive (%d -> %d); "
+        "queued console events are not being dispatched"
+        % (after_load, during_render)
+    )
+
+
+def test_save_3d_views_raises_on_webgl_failure(monkeypatch):
+    """A reported WebGL failure aborts the render rather than writing a blank png.
+
+    The failure is injected rather than provoked, so this does not depend on a
+    shader being broken in the tree under test.
+    """
+    monkeypatch.setattr(
+        "cortex.export.headless.webgl_failures",
+        lambda errors: ["[console.error] THREE.WebGLProgram: Could not initialise shader."],
+    )
+    vol = cortex.Volume(np.random.randn(*volshape), subj, xfmname)
+    with tempfile.TemporaryDirectory() as tmpdir:
+        with pytest.raises(RuntimeError, match="WebGL failed while rendering"):
+            cortex.export.save_3d_views(
+                vol,
+                base_name=os.path.join(tmpdir, "boom"),
+                list_angles=["lateral_pivot"],
+                list_surfaces=["inflated"],
+                size=(512, 384),
+                trim=False,
+                sleep=10,
+                headless=True,
+            )
diff --git a/cortex/tests/test_visual_regression.py b/cortex/tests/test_visual_regression.py
new file mode 100644
index 000000000..74260644a
--- /dev/null
+++ b/cortex/tests/test_visual_regression.py
@@ -0,0 +1,897 @@
+"""Visual regression tests: quickflat and webgl renders vs stored references.
+
+Four suites. Three of them render flatmaps of the six public dataview classes
+(``Volume``, ``Vertex``, ``Volume2D``, ``Vertex2D``, ``VolumeRGB``,
+``VertexRGB``) through both matplotlib (``cortex.quickshow``) and the headless
+WebGL viewer (``cortex.export.plot_panels``), varying what the data carries:
+alpha-bearing values, NaNs in the data channels, and NaNs in the alpha map. The
+last covers only the two RGB classes, the only ones taking an explicit
+``alpha=``. Every one of those renders is checked twice -- against its own
+stored reference at a tight tolerance, and directly against the other renderer's
+render of the same dataview at a loose one.
+
+The fourth suite renders non-flatmap views, ``Volume`` and ``Vertex`` on the
+inflated and fiducial surfaces, through ``save_3d_views``. Those are
+webgl-only and get the reference check alone: ``cortex.quickshow`` produces
+flatmaps and nothing else, so there is nothing to diff them against.
+
+See ``reference_images/README.md`` for how the references were produced and how
+to regenerate them.
+
+All tests are skipped if playwright is not installed.
+"""
+
+import os
+from pathlib import Path
+from typing import Optional
+
+import numpy as np
+import numpy.typing as npt
+import pytest
+
+import cortex
+import cortex.export
+import cortex.polyutils
+from cortex.dataset import Dataview
+from cortex.tests.testing_utils import has_playwright
+
+pytestmark = pytest.mark.skipif(
+    not has_playwright, reason="playwright and chromium are required"
+)
+
+# These all xfailed until the Vertex2D webgl fix landed: the #679 lighting
+# refactor ported the HASFLAT bump-displacement block into the vertex shader,
+# which under headless/SwiftShader rendering left pure-Vertex flatmaps either
+# blank (Vertex2D, which raised in plot_panels' border-trim) or trimmed to a
+# bounding box 2px wider than the then-current references (Vertex/VertexRGB).
+# The blank render is fixed; the 2px bounding box is the shader change's intended
+# new geometry and persists, so the references were regenerated to match it.
+DATAVIEW_NAMES = ["Volume", "Vertex", "Volume2D", "Vertex2D", "VolumeRGB", "VertexRGB"]
+
+subj = "S1"
+xfmname = "fullhead"
+
+#: Stored renders this test asserts against. See that directory's README for how
+#: they were produced and how to regenerate them.
+REFERENCE_DIR = Path(__file__).parent / "reference_images" / "alpha_dataviews"
+
+#: As REFERENCE_DIR, but for dataviews whose source data contains NaNs.
+NAN_REFERENCE_DIR = Path(__file__).parent / "reference_images" / "nan_dataviews"
+
+#: As NAN_REFERENCE_DIR, but with the NaNs in the *alpha map* rather than in the
+#: data. Only the RGB dataviews take an explicit ``alpha=``, so only those two.
+NAN_ALPHA_REFERENCE_DIR = Path(__file__).parent / "reference_images" / "nan_alpha_dataviews"
+
+#: Dataviews that accept an explicit alpha map, and so can carry NaNs in it.
+NAN_ALPHA_DATAVIEW_NAMES = ["VolumeRGB", "VertexRGB"]
+
+#: Non-flatmap views, checked against a webgl reference only. quickflat renders
+#: nothing but flatmaps, so these have no counterpart to diff against and no
+#: cross-renderer leg -- see test_visual_comparison_nonflat_views.
+NONFLAT_REFERENCE_DIR = Path(__file__).parent / "reference_images" / "nonflat_views"
+
+#: (surface, angle, dataview). Volume and Vertex cover both shader paths, which
+#: matters because the flatmap suite exercises them under conditions that turn
+#: out to be a different regime: the two known webgl lighting bugs reproduce on
+#: flatmaps only.
+NONFLAT_VIEWS = [
+    ("inflated", "lateral_pivot", "Volume"),
+    ("inflated", "lateral_pivot", "Vertex"),
+    ("fiducial", "lateral_pivot", "Volume"),
+    ("fiducial", "lateral_pivot", "Vertex"),
+]
+
+#: Lossless WebP: bit-exact after decode and 59% the size of optimized PNG. AVIF
+#: is smaller still but Pillow cannot write it losslessly -- it was measured at
+#: max|difference| 27-29, larger than DIFF_THRESHOLD below, so it would corrupt
+#: the comparison it is meant to feed. Higher PNG compression levels are pointless
+#: here: level 6 and level 9 produce byte-identical output.
+REFERENCE_SUFFIX = ".webp"
+
+#: Rewrite the references from this run instead of comparing against them.
+REGENERATE_REFERENCES = bool(os.environ.get("REGENERATE_REFERENCE_IMAGES"))
+
+# Tolerances. The renders are deterministic -- repeated runs on one machine are
+# bit-identical -- so these are not absorbing noise. They exist because the
+# references are coupled to the Chromium and matplotlib builds that produced them,
+# and an upgrade can shift anti-aliasing and rasterization slightly. They are far
+# tighter than any real regression: a wrong colormap, a dropped alpha channel or
+# swapped colour channels all move large areas of the image by much more.
+MAX_MEAN_ABS_DIFF = 2.0        # mean |difference| over all pixels/channels, of 255
+DIFF_THRESHOLD = 16            # a pixel "differs" if any channel moves by more
+MAX_FRACTION_DIFFERING = 0.02  # at most this fraction of pixels may differ
+
+# The two limits above are both weak against a change that moves a *small* number
+# of pixels by a *large* amount, which is what a geometry or contour shift looks
+# like: the mean is diluted by the ~97% of pixels that did not move, and a shifted
+# contour only just clears the fraction limit. The #679 vertex-shader change was
+# caught with 1.4x margin on the fraction and would have passed the mean outright
+# (1.574 against a limit of 2.0), despite moving 2.8% of pixels by a median of
+# 67/255.
+#
+# So two further criteria, each covering what the others miss. Measured against
+# simulated cosmetic drift (a 596<->594 resample, +-1 LSB quantisation, gamma
+# 1.02) versus the real #679 change:
+#
+#            metric                cosmetic     #679    separation
+#            mean                    0.596      1.574       2.6x
+#            fraction > 16           0.48%      2.79%       5.8x
+#            fraction > 32          0.0052%     2.16%       large
+#            SSIM loss              0.0022     0.0393      18.3x
+#
+# They are complementary rather than ranked, so all four are checked:
+#   - mean and fraction>16 catch broad, low-amplitude shifts (a gamma change
+#     scores 6.68 mean but 0.00% on fraction>32).
+#   - fraction>32 catches sparse, high-amplitude ones (the #679 case).
+#   - SSIM catches structural change, but is computed on luminance and is
+#     therefore blind to a channel permutation -- an R/B swap scores 0.0000 SSIM
+#     loss while scoring 21.75 on the mean. It only adds sensitivity alongside
+#     the others; it cannot replace them.
+#
+# Caveat on the limits below: the cosmetic figures are *simulated*, not measured
+# across a real Chromium or matplotlib upgrade. A genuine toolchain bump changes
+# anti-aliasing, which can move a few pixels a long way and so does show up in
+# fraction>32. Both limits are set well above the simulated drift rather than at
+# it, to leave room for that.
+#
+# The gross threshold is 32 rather than 64 because of a second worked example.
+# cb976270 ("unify NaN and alpha handling") changes four quickflat volumetric
+# renders through premultiplied-alpha thickness averaging, and at 64 every one of
+# them scores 0.000% -- the suite missed it entirely. Its mean cannot be tightened
+# into a catch either: at 0.19 it sits *below* the simulated cosmetic floor of
+# 0.60, so a tighter mean yields false positives before it yields a catch. At 32
+# the signal is 0.256% against this 0.1% limit while the worst simulated cosmetic
+# drift is 0.0052%, i.e. 19x below it. Dropping to 32 also strengthens the #679
+# catch (1.46% -> 2.16%). The cost is headroom: at 64 no cosmetic perturbation
+# registered at all, so a real toolchain bump now has 19x of room rather than
+# effectively unlimited.
+GROSS_DIFF_THRESHOLD = 32             # a pixel differs "grossly" if any channel moves by more
+MAX_FRACTION_GROSSLY_DIFFERING = 0.001
+MAX_SSIM_LOSS = 0.01                  # 1 - mean SSIM over the luminance channel
+
+# quickflat and webgl's flatmap renders disagree by a fixed anisotropic
+# scale + translation, not a genuine content difference: fitting a homography
+# (cv2.findTransformECC, MOTION_HOMOGRAPHY) between the two on curvature-only
+# content (data alpha=0, so no dataview-specific signal) landed on essentially
+# zero rotation and zero projective terms -- i.e. no real perspective
+# distortion -- but scale_x=0.9690, scale_y=0.9365, and a ~7px translation.
+# Applying just that (affine, webgl -> quickflat's frame) dropped curvature-only
+# mean|diff| from 18.9 to 4.4. Translation is stored as a fraction of
+# quickflat's own (width, height) so it scales if that size varies slightly
+# between dataviews; the scale factors are already dimensionless.
+CROSS_RENDERER_SCALE_X = 0.9690
+CROSS_RENDERER_SCALE_Y = 0.9365
+CROSS_RENDERER_TRANSLATE_X_FRAC = 6.9008 / 392
+CROSS_RENDERER_TRANSLATE_Y_FRAC = 6.9650 / 204
+
+# Cross-renderer tolerances: quickflat vs webgl for the *same* dataview, rather
+# than each against its own reference. Looser than the within-renderer ones
+# above because matplotlib and Three.js still differ in anti-aliasing and
+# colormap sampling even after the affine correction above removes the
+# dominant coordinate-frame mismatch. What must not happen is the kind of
+# broad color/shape disagreement (wrong colormap, dropped alpha, swapped
+# channels) that a real regression causes, which is far larger than this.
+#
+# Calibrated across all fourteen flatmap pairs, which with the affine correction
+# span mean|diff| 1.06-2.01 and 0.10-1.64% of pixels differing. The limits sit
+# 3-4x above those worst cases: tight enough that losing the correction, or a
+# real colormap/alpha/channel regression, fails immediately, loose enough to
+# absorb a Chromium or matplotlib upgrade.
+CROSS_MAX_MEAN_ABS_DIFF = 8.0
+CROSS_DIFF_THRESHOLD = 64
+CROSS_MAX_FRACTION_DIFFERING = 0.05
+
+# Render settings chosen to minimise cross-renderer disagreement, from a factorial
+# sweep of both renderers' settings (curvature brightness/contrast/threshold,
+# depth, sampler, thick/layers, and webgl's three lighting controls).
+#
+# Only one setting was worth changing from the defaults: curvature thresholding.
+# Thresholded curvature puts a hard binary edge at curvature=0 whose sub-pixel
+# placement each rasteriser resolves differently, so it is maximally sensitive to
+# the residual misalignment; smooth curvature is low-frequency and resamples
+# cleanly. Averaged over the brightness/contrast grid the thresholded pairing
+# scored 8.50 against 3.66 for the smooth one, and 4.98 vs 2.66 at the default
+# brightness/contrast. quickflat's ``curvature_threshold`` and webgl's
+# ``curvature.smoothness`` are the same knob from opposite ends -- smoothness 0.0
+# *is* thresholded -- so both have to move together or they disagree by more.
+#
+# Everything else stayed at its default, on the evidence:
+#   - lighting: the default (all off) already ties the best combination at 2.304;
+#     `uniform_illumination=1` renders identically on a flatmap (0.003). Only
+#     `topleft_lighting=1` hurts (15.30), and it is off by default.
+#   - sampler: trilinear beat nearest by ~1% (2.336 vs 2.365), i.e. noise.
+#   - depth / thick / layers: total spread 0.38 over the whole grid, below the
+#     noise floor -- the polarity difference between quickflat's `depth` and
+#     webgl's `thickmix` is not resolvable and does not matter here.
+#
+# NB this deliberately renders with curvature *un*-thresholded, which is not
+# pycortex's default appearance, so these references do not cover the default
+# curvature path. That is the trade for a ~2x tighter cross-renderer floor.
+QUICKFLAT_CURVATURE_THRESHOLD = False
+WEBGL_CURVATURE_SMOOTHNESS = 1.0
+
+
+def _ssim(a: npt.NDArray, b: npt.NDArray) -> float:
+    """Mean structural similarity between two RGBA images, over luminance.
+
+    Standard SSIM with an 11x11 Gaussian window (sigma 1.5) and the usual
+    stabilising constants, implemented on scipy because scikit-image, which
+    would otherwise supply it, is not a pycortex dependency. Returns 1.0 for
+    identical input.
+
+    Computed on the channel mean, so it is invariant to a channel permutation --
+    see the note on MAX_SSIM_LOSS. It is a structural check, not a colour one.
+    """
+    from scipy.ndimage import gaussian_filter
+
+    x = a[..., :3].mean(axis=-1).astype(np.float64)
+    y = b[..., :3].mean(axis=-1).astype(np.float64)
+    c1, c2 = (0.01 * 255) ** 2, (0.03 * 255) ** 2
+
+    blur = lambda img: gaussian_filter(img, sigma=1.5, truncate=(11 - 1) / 2 / 1.5)
+    mu_x, mu_y = blur(x), blur(y)
+    var_x = blur(x * x) - mu_x**2
+    var_y = blur(y * y) - mu_y**2
+    cov = blur(x * y) - mu_x * mu_y
+
+    num = (2 * mu_x * mu_y + c1) * (2 * cov + c2)
+    den = (mu_x**2 + mu_y**2 + c1) * (var_x + var_y + c2)
+    return float((num / den).mean())
+
+
+def _check_against_reference(
+    name: str,
+    actual_path: Path,
+    debug_dir: Path,
+    reference_dir: Path = REFERENCE_DIR,
+) -> Optional[str]:
+    """Compare one render to its reference.
+
+    Checks four criteria, all of which must pass -- see ``MAX_MEAN_ABS_DIFF``
+    and the tolerances below it for why they are complementary rather than
+    redundant. Returns a description of every breach, or None if the render
+    matches.
+
+    A shape mismatch is not an immediate failure: the render is resized to the
+    reference's size before diffing, since the two are still expected to show
+    the same content at a slightly different trim/crop. That is a plain resize,
+    without the affine that ``_check_cross_renderer`` also applies -- one
+    renderer against itself has no coordinate-frame mismatch to undo.
+
+    On any breach, writes the (possibly resized) render and an amplified
+    difference image into ``debug_dir``, so the change can be inspected rather
+    than guessed at. With ``REGENERATE_REFERENCES`` set it overwrites the
+    reference instead and reports no mismatch.
+    """
+    from PIL import Image
+
+    ref_path = reference_dir / f"{name}{REFERENCE_SUFFIX}"
+
+    if REGENERATE_REFERENCES:
+        reference_dir.mkdir(parents=True, exist_ok=True)
+        Image.open(actual_path).convert("RGBA").save(
+            ref_path, format="WEBP", lossless=True, method=6, quality=100, exact=True
+        )
+        return None
+
+    actual_im = Image.open(actual_path).convert("RGBA")
+    ref_im = Image.open(ref_path).convert("RGBA")
+    shape_note = ""
+    if actual_im.size != ref_im.size:
+        shape_note = f" (aligned: render was {actual_im.size}, reference is {ref_im.size})"
+        actual_im = actual_im.resize(ref_im.size, Image.BILINEAR)
+
+    actual = np.asarray(actual_im).astype(np.int16)
+    ref = np.asarray(ref_im).astype(np.int16)
+
+    diff = np.abs(actual - ref)
+    per_pixel = diff.max(axis=-1)
+    mean_abs = float(diff.mean())
+    fraction = float((per_pixel > DIFF_THRESHOLD).mean())
+    gross = float((per_pixel > GROSS_DIFF_THRESHOLD).mean())
+    ssim_loss = 1.0 - _ssim(actual, ref)
+
+    breaches = []
+    if mean_abs > MAX_MEAN_ABS_DIFF:
+        breaches.append(f"mean|diff|={mean_abs:.3f} (limit {MAX_MEAN_ABS_DIFF})")
+    if fraction > MAX_FRACTION_DIFFERING:
+        breaches.append(
+            f"{fraction:.2%} of pixels differ by more than {DIFF_THRESHOLD} "
+            f"(limit {MAX_FRACTION_DIFFERING:.0%})"
+        )
+    if gross > MAX_FRACTION_GROSSLY_DIFFERING:
+        breaches.append(
+            f"{gross:.3%} of pixels differ by more than {GROSS_DIFF_THRESHOLD} "
+            f"(limit {MAX_FRACTION_GROSSLY_DIFFERING:.1%})"
+        )
+    if ssim_loss > MAX_SSIM_LOSS:
+        breaches.append(f"SSIM loss={ssim_loss:.4f} (limit {MAX_SSIM_LOSS})")
+    if not breaches:
+        return None
+
+    actual_im.save(debug_dir / f"actual_{name}.png")
+    amplified = np.clip(diff[..., :3] * 8, 0, 255).astype("uint8")
+    Image.fromarray(amplified).save(debug_dir / f"diff_{name}.png")
+    return f"{name}: " + "; ".join(breaches) + shape_note
+
+
+def _check_cross_renderer(
+    name: str,
+    quickflat_path: Path,
+    webgl_path: Path,
+    debug_dir: Path,
+) -> Optional[str]:
+    """Compare a quickflat render directly against its webgl counterpart.
+
+    Unlike ``_check_against_reference`` this has no stored fixture: it diffs
+    the two renders produced by *this* test run against each other. webgl's
+    render is resized to quickflat's own size, then the fixed affine
+    correction above (scale + translation) is applied to align it onto
+    quickflat's coordinate frame before diffing -- see
+    ``CROSS_RENDERER_SCALE_X`` for how that was derived. Returns a mismatch
+    description, or None if they agree within the (loose) cross-renderer
+    tolerance.
+    """
+    from PIL import Image
+
+    qf_im = Image.open(quickflat_path).convert("RGBA")
+    wg_im = Image.open(webgl_path).convert("RGBA")
+    shape_note = f" (aligned: quickflat was {qf_im.size}, webgl was {wg_im.size})"
+
+    size = qf_im.size
+    w, h = size
+    tx = CROSS_RENDERER_TRANSLATE_X_FRAC * w
+    ty = CROSS_RENDERER_TRANSLATE_Y_FRAC * h
+    wg_aligned = wg_im.resize(size, Image.BILINEAR).transform(
+        size, Image.AFFINE,
+        (CROSS_RENDERER_SCALE_X, 0.0, tx, 0.0, CROSS_RENDERER_SCALE_Y, ty),
+        resample=Image.BILINEAR,
+    )
+
+    qf = np.asarray(qf_im).astype(np.int16)
+    wg = np.asarray(wg_aligned).astype(np.int16)
+
+    diff = np.abs(qf - wg)
+    mean_abs = float(diff.mean())
+    fraction = float((diff.max(axis=-1) > CROSS_DIFF_THRESHOLD).mean())
+    if mean_abs <= CROSS_MAX_MEAN_ABS_DIFF and fraction <= CROSS_MAX_FRACTION_DIFFERING:
+        return None
+
+    amplified = np.clip(diff[..., :3] * 4, 0, 255).astype("uint8")
+    Image.fromarray(amplified).save(debug_dir / f"cross_diff_{name}.png")
+    return (
+        f"cross_{name}: quickflat vs webgl mean|diff|={mean_abs:.3f} "
+        f"(limit {CROSS_MAX_MEAN_ABS_DIFF}), {fraction:.2%} of pixels differ by "
+        f"more than {CROSS_DIFF_THRESHOLD} (limit {CROSS_MAX_FRACTION_DIFFERING:.0%})"
+        f"{shape_note}"
+    )
+
+
+def _build_alpha_dataview(name: str) -> Dataview:
+    """Build a single alpha-bearing dataview by name.
+
+    Returns the dataview object. All data synthesis (volumetric and vertex) is
+    done here, shared across all six dataview types.
+    """
+    # Volumetric
+    zz, yy, xx = np.mgrid[0:31, 0:100, 0:100]
+    data_vol = (xx - 50) / 50.0  # ~ [-1, 1]
+    center = np.array([15, 50, 50])
+    sigma_v = 25.0
+    dist2 = (
+        (zz - center[0]) ** 2 + (yy - center[1]) ** 2 + (xx - center[2]) ** 2
+    )
+    accuracy_vol = np.exp(-dist2 / (2 * sigma_v**2))  # [0, 1] bump
+    red_vol = np.clip(xx / 99.0, 0, 1)
+    green_vol = np.clip(yy / 99.0, 0, 1)
+    blue_vol = np.clip(zz / 30.0, 0, 1)
+
+    # Surface (vertex) — encode by spatial coordinate, not vertex index
+    surfs = [
+        cortex.polyutils.Surface(*d)
+        for d in cortex.db.get_surf(subj, "fiducial")
+    ]
+    num_verts = [s.pts.shape[0] for s in surfs]
+    pts = np.vstack([surfs[0].pts, surfs[1].pts])
+    y_centered = pts[:, 1] - pts[:, 1].mean()
+    data_vtx = y_centered / np.abs(y_centered).max()  # [-1, 1]
+    xyz_norm = (pts - pts.min(axis=0)) / (pts.max(axis=0) - pts.min(axis=0))
+
+    # Gaussian falloff from `seed`, used as the accuracy/alpha channel.
+    def _bump(surf: cortex.polyutils.Surface, seed: int, sigma: float) -> npt.NDArray[np.floating]:
+        d = np.linalg.norm(surf.pts - surf.pts[seed], axis=1)
+        return np.exp(-(d**2) / (2 * sigma**2))
+
+    accuracy_vtx = np.hstack(
+        [
+            _bump(surfs[0], num_verts[0] // 2, sigma=40.0),
+            _bump(surfs[1], num_verts[1] // 2, sigma=40.0),
+        ]
+    )
+
+    cmap_plain = "viridis"
+    cmap_2d = "RdBu_r_alpha"
+
+    if name == "Volume":
+        return cortex.Volume(data_vol, subj, xfmname, cmap=cmap_plain, vmin=-1, vmax=1)
+    elif name == "Vertex":
+        return cortex.Vertex(data_vtx, subj, cmap=cmap_plain, vmin=-1, vmax=1)
+    elif name == "Volume2D":
+        return cortex.Volume2D(
+            data_vol, accuracy_vol, subj, xfmname, cmap=cmap_2d,
+            vmin=-1, vmax=1, vmin2=0, vmax2=1,
+        )
+    elif name == "Vertex2D":
+        return cortex.Vertex2D(
+            data_vtx, accuracy_vtx, subj, cmap=cmap_2d,
+            vmin=-1, vmax=1, vmin2=0, vmax2=1,
+        )
+    elif name == "VolumeRGB":
+        return cortex.VolumeRGB(
+            cortex.Volume(red_vol, subj, xfmname, vmin=0, vmax=1),
+            cortex.Volume(green_vol, subj, xfmname, vmin=0, vmax=1),
+            cortex.Volume(blue_vol, subj, xfmname, vmin=0, vmax=1),
+            subj, xfmname,
+            alpha=cortex.Volume(accuracy_vol, subj, xfmname, vmin=0, vmax=1),
+        )
+    elif name == "VertexRGB":
+        return cortex.VertexRGB(
+            cortex.Vertex(xyz_norm[:, 0], subj, vmin=0, vmax=1),
+            cortex.Vertex(xyz_norm[:, 1], subj, vmin=0, vmax=1),
+            cortex.Vertex(xyz_norm[:, 2], subj, vmin=0, vmax=1),
+            subj,
+            alpha=cortex.Vertex(accuracy_vtx, subj, vmin=0, vmax=1),
+        )
+    else:
+        raise ValueError(f"Unknown dataview: {name}")
+
+
+def _build_nan_dataview(name: str) -> Dataview:
+    """Build a single NaN-bearing dataview by name.
+
+    Returns the dataview object. All data synthesis is done here, with NaNs
+    masking roughly half of each volume/surface in the primary data channel.
+    """
+    zz, yy, xx = np.mgrid[0:31, 0:100, 0:100]
+    data_vol = (xx - 50) / 50.0  # ~ [-1, 1]
+    center = np.array([15, 50, 50])
+    sigma_v = 25.0
+    dist2 = (
+        (zz - center[0]) ** 2 + (yy - center[1]) ** 2 + (xx - center[2]) ** 2
+    )
+    accuracy_vol = np.exp(-dist2 / (2 * sigma_v**2))  # [0, 1] bump
+    red_vol = np.clip(xx / 99.0, 0, 1)
+    green_vol = np.clip(yy / 99.0, 0, 1)
+    blue_vol = np.clip(zz / 30.0, 0, 1)
+
+    # cb976270's rule is that a NaN *anywhere* at a voxel -- the data, either 2D
+    # dimension, any RGB channel, or the alpha map -- renders fully transparent.
+    # Each of those gets NaN'd over its own region, so a single render exercises
+    # several branches of the rule at once and a failure still says which one
+    # moved. The vertex regions are disjoint; the volume ones (x>=50, y>=50,
+    # z>=15) overlap, which is harmless and additionally covers voxels carrying
+    # more than one NaN at once. Blue is deliberately left clean, as a control
+    # that not everything has simply gone transparent.
+    #
+    # The alpha map is NaN'd here too, on a third axis. That is not a duplicate
+    # of the nan_alpha suite: this covers alpha NaNs superposed on colour NaNs,
+    # where nan_alpha isolates them with every colour channel clean.
+    #
+    # Expect the alpha map's surviving NaN fraction to look smaller than its mask.
+    # Where a colour channel is already NaN the pipeline writes alpha's vmin over
+    # it, so of the z>=15 region only the quarter with red and green both clean
+    # stays NaN; the rest is resolved to a hard 0. Both halves of that are worth
+    # rendering, which is why the regions are allowed to overlap.
+    vol_nan_mask = xx >= 50          # primary channels: data, and red for RGB
+    vol_nan_mask2 = yy >= 50         # secondary: 2D dimension 2, and green for RGB
+    vol_nan_mask3 = zz >= 15         # the alpha map, on a third independent axis
+    data_vol_nan = data_vol.copy()
+    data_vol_nan[vol_nan_mask] = np.nan
+    red_vol_nan = red_vol.copy()
+    red_vol_nan[vol_nan_mask] = np.nan
+    accuracy_vol_nan = accuracy_vol.copy()
+    accuracy_vol_nan[vol_nan_mask2] = np.nan
+    green_vol_nan = green_vol.copy()
+    green_vol_nan[vol_nan_mask2] = np.nan
+    alpha_vol_nan = accuracy_vol.copy()
+    alpha_vol_nan[vol_nan_mask3] = np.nan
+
+    surfs = [
+        cortex.polyutils.Surface(*d)
+        for d in cortex.db.get_surf(subj, "fiducial")
+    ]
+    num_verts = [s.pts.shape[0] for s in surfs]
+    pts = np.vstack([surfs[0].pts, surfs[1].pts])
+    y_centered = pts[:, 1] - pts[:, 1].mean()
+    data_vtx = y_centered / np.abs(y_centered).max()  # [-1, 1]
+    xyz_norm = (pts - pts.min(axis=0)) / (pts.max(axis=0) - pts.min(axis=0))
+
+    # Gaussian falloff from `seed`, used as the accuracy/alpha channel.
+    def _bump(surf: cortex.polyutils.Surface, seed: int, sigma: float) -> npt.NDArray[np.floating]:
+        d = np.linalg.norm(surf.pts - surf.pts[seed], axis=1)
+        return np.exp(-(d**2) / (2 * sigma**2))
+
+    accuracy_vtx = np.hstack(
+        [
+            _bump(surfs[0], num_verts[0] // 2, sigma=40.0),
+            _bump(surfs[1], num_verts[1] // 2, sigma=40.0),
+        ]
+    )
+
+    # As above, in two disjoint index ranges rather than spatial ones.
+    total_verts = sum(num_verts)
+    vtx_nan_mask = np.arange(total_verts) >= total_verts // 2
+    vtx_nan_mask2 = np.arange(total_verts) < total_verts // 4
+    vtx_nan_mask3 = (np.arange(total_verts) >= total_verts // 4) & (
+        np.arange(total_verts) < total_verts // 2
+    )
+    data_vtx_nan = data_vtx.copy()
+    data_vtx_nan[vtx_nan_mask] = np.nan
+    x_norm_nan = xyz_norm[:, 0].copy()
+    x_norm_nan[vtx_nan_mask] = np.nan
+    accuracy_vtx_nan = accuracy_vtx.copy()
+    accuracy_vtx_nan[vtx_nan_mask2] = np.nan
+    y_norm_nan = xyz_norm[:, 1].copy()
+    y_norm_nan[vtx_nan_mask2] = np.nan
+    alpha_vtx_nan = accuracy_vtx.copy()
+    alpha_vtx_nan[vtx_nan_mask3] = np.nan
+
+    cmap_plain = "viridis"
+    cmap_2d = "RdBu_r_alpha"
+
+    if name == "Volume":
+        return cortex.Volume(data_vol_nan, subj, xfmname, cmap=cmap_plain, vmin=-1, vmax=1)
+    elif name == "Vertex":
+        return cortex.Vertex(data_vtx_nan, subj, cmap=cmap_plain, vmin=-1, vmax=1)
+    elif name == "Volume2D":
+        return cortex.Volume2D(
+            data_vol_nan, accuracy_vol_nan, subj, xfmname, cmap=cmap_2d,
+            vmin=-1, vmax=1, vmin2=0, vmax2=1,
+        )
+    elif name == "Vertex2D":
+        return cortex.Vertex2D(
+            data_vtx_nan, accuracy_vtx_nan, subj, cmap=cmap_2d,
+            vmin=-1, vmax=1, vmin2=0, vmax2=1,
+        )
+    elif name == "VolumeRGB":
+        return cortex.VolumeRGB(
+            cortex.Volume(red_vol_nan, subj, xfmname, vmin=0, vmax=1),
+            cortex.Volume(green_vol_nan, subj, xfmname, vmin=0, vmax=1),
+            cortex.Volume(blue_vol, subj, xfmname, vmin=0, vmax=1),
+            subj, xfmname,
+            alpha=cortex.Volume(alpha_vol_nan, subj, xfmname, vmin=0, vmax=1),
+        )
+    elif name == "VertexRGB":
+        return cortex.VertexRGB(
+            cortex.Vertex(x_norm_nan, subj, vmin=0, vmax=1),
+            cortex.Vertex(y_norm_nan, subj, vmin=0, vmax=1),
+            cortex.Vertex(xyz_norm[:, 2], subj, vmin=0, vmax=1),
+            subj,
+            alpha=cortex.Vertex(alpha_vtx_nan, subj, vmin=0, vmax=1),
+        )
+    else:
+        raise ValueError(f"Unknown dataview: {name}")
+
+
+def _build_nan_alpha_dataview(name: str) -> Dataview:
+    """Build an RGB dataview whose *alpha map* carries NaNs, colour channels clean.
+
+    The other NaN suite puts NaNs in the data; this puts them in the alpha map,
+    which is a separate code path -- alpha is not colour-mapped, it is used
+    directly as a blend weight, so a NaN reaches the compositing arithmetic
+    rather than a colormap lookup.
+
+    Current behaviour is that those elements render fully transparent, i.e. the
+    curvature underlay shows through, which is what the other NaN cases do too.
+    """
+    zz, yy, xx = np.mgrid[0:31, 0:100, 0:100]
+    center = np.array([15, 50, 50])
+    dist2 = (zz - center[0]) ** 2 + (yy - center[1]) ** 2 + (xx - center[2]) ** 2
+    accuracy_vol = np.exp(-dist2 / (2 * 25.0**2))
+
+    if name == "VolumeRGB":
+        alpha_nan = accuracy_vol.copy()
+        alpha_nan[xx >= 50] = np.nan
+        return cortex.VolumeRGB(
+            cortex.Volume(np.clip(xx / 99.0, 0, 1), subj, xfmname, vmin=0, vmax=1),
+            cortex.Volume(np.clip(yy / 99.0, 0, 1), subj, xfmname, vmin=0, vmax=1),
+            cortex.Volume(np.clip(zz / 30.0, 0, 1), subj, xfmname, vmin=0, vmax=1),
+            subj, xfmname,
+            alpha=cortex.Volume(alpha_nan, subj, xfmname, vmin=0, vmax=1),
+        )
+    elif name == "VertexRGB":
+        surfs = [cortex.polyutils.Surface(*d) for d in cortex.db.get_surf(subj, "fiducial")]
+        num_verts = [s.pts.shape[0] for s in surfs]
+        pts = np.vstack([surfs[0].pts, surfs[1].pts])
+        xyz_norm = (pts - pts.min(axis=0)) / (pts.max(axis=0) - pts.min(axis=0))
+
+        # Gaussian falloff from `seed`, used as the accuracy/alpha channel.
+        def _bump(
+            surf: cortex.polyutils.Surface, seed: int, sigma: float
+        ) -> npt.NDArray[np.floating]:
+            d = np.linalg.norm(surf.pts - surf.pts[seed], axis=1)
+            return np.exp(-(d**2) / (2 * sigma**2))
+
+        accuracy_vtx = np.hstack([
+            _bump(surfs[0], num_verts[0] // 2, sigma=40.0),
+            _bump(surfs[1], num_verts[1] // 2, sigma=40.0),
+        ])
+        total = sum(num_verts)
+        alpha_nan = accuracy_vtx.copy()
+        alpha_nan[np.arange(total) >= total // 2] = np.nan
+        return cortex.VertexRGB(
+            cortex.Vertex(xyz_norm[:, 0], subj, vmin=0, vmax=1),
+            cortex.Vertex(xyz_norm[:, 1], subj, vmin=0, vmax=1),
+            cortex.Vertex(xyz_norm[:, 2], subj, vmin=0, vmax=1),
+            subj,
+            alpha=cortex.Vertex(alpha_nan, subj, vmin=0, vmax=1),
+        )
+    else:
+        raise ValueError(f"Unknown dataview: {name}")
+
+
+def _render_and_check_dataview(
+    name: str,
+    view: Dataview,
+    reference_dir: Path,
+    tmp_path: Path,
+) -> list[str]:
+    """Render a single dataview via quickshow + webgl and check against reference.
+
+    Each render is checked three ways: quickflat vs its own reference, webgl vs
+    its own reference (both tight tolerances, see ``_check_against_reference``),
+    and quickflat vs webgl directly (loose tolerance, see
+    ``_check_cross_renderer``).
+
+    Returns a list of failure messages (empty if no failures). Skips the test
+    if reference images are missing, and regenerates them if ``REGENERATE_REFERENCES``
+    is set.
+    """
+    import matplotlib.pyplot as plt
+
+    from cortex.export.save_views import angle_view_params
+
+    # webgl settings have to ride along with the angle params, as an
+    # ``(name, params)`` tuple: ``plot_panels``' ``viewer_params`` is forwarded to
+    # ``cortex.webshow`` -> ``show()``, whose ``**kwargs`` silently swallows
+    # anything that is not one of its named arguments, so viewer state passed that
+    # way is a no-op. The tuple form reaches ``_set_view``, which is the channel
+    # that actually applies it. Keeping ``surface`` as the plain string "flatmap"
+    # leaves the output filename clean.
+    flatmap_angle = (
+        "flatmap",
+        {
+            **angle_view_params["flatmap"],
+            "surface.{subject}.curvature.smoothness": WEBGL_CURVATURE_SMOOTHNESS,
+        },
+    )
+
+    if not REGENERATE_REFERENCES:
+        for prefix in ("quickflat", "webgl"):
+            ref_path = reference_dir / f"{prefix}_{name}{REFERENCE_SUFFIX}"
+            if not ref_path.exists():
+                pytest.skip(
+                    f"No reference for {prefix}_{name} in {reference_dir}. "
+                    "See that directory's README for regeneration."
+                )
+
+    # quickshow → low-res PNG
+    qf_path = tmp_path / f"quickflat_{name}.png"
+    qf_fig = cortex.quickshow(
+        view,
+        with_curvature=True,
+        with_rois=False,
+        with_labels=False,
+        with_colorbar=False,
+        with_sulci=False,
+        with_borders=False,
+        height=256,
+        curvature_threshold=QUICKFLAT_CURVATURE_THRESHOLD,
+    )
+    qf_fig.savefig(qf_path, bbox_inches="tight", pad_inches=0, dpi=80)
+    plt.close(qf_fig)
+
+    # webgl → trimmed flatmap PNG via plot_panels (single flatmap panel).
+    # Not save_3d_views, whose output is transparent where quickshow's figure is
+    # opaque; going through matplotlib keeps both sides comparable for the
+    # cross-renderer check.
+    flatmap_panel = [
+        {
+            "extent": [0.0, 0.0, 1.0, 1.0],
+            "view": {"angle": flatmap_angle, "surface": "flatmap"},
+        }
+    ]
+    wg_path = str(tmp_path / f"webgl_{name}.png")
+    wg_fig = cortex.export.plot_panels(
+        view,
+        panels=flatmap_panel,
+        figsize=(6, 3),
+        windowsize=(512, 384),
+        save_name=wg_path,
+        sleep=10,
+        viewer_params=dict(labels_visible=[], overlays_visible=[]),
+        headless=True,
+    )
+    plt.close(wg_fig)
+    wg_path = Path(wg_path)
+
+    failures = []
+
+    # Within-renderer checks (always run, but skip assertion if regenerating)
+    for prefix, path in [("quickflat", qf_path), ("webgl", wg_path)]:
+        msg = _check_against_reference(f"{prefix}_{name}", path, tmp_path, reference_dir)
+        if msg is not None and not REGENERATE_REFERENCES:
+            failures.append(msg)
+
+    if REGENERATE_REFERENCES:
+        pytest.skip(f"Regenerated {name} references in {reference_dir}")
+
+    # Cross-renderer check (never regenerates, always compares)
+    msg = _check_cross_renderer(name, qf_path, wg_path, tmp_path)
+    if msg is not None:
+        failures.append(msg)
+
+    return failures
+
+
+def _render_and_check_webgl_only(
+    tag: str,
+    view: Dataview,
+    surface: str,
+    angle: str,
+    reference_dir: Path,
+    tmp_path: Path,
+) -> list[str]:
+    """Render one non-flatmap view through webgl and check it against a reference.
+
+    A cut-down ``_render_and_check_dataview``: there is no quickflat render to
+    compare against, because ``cortex.quickshow`` produces flatmaps and nothing
+    else, so both the quickflat reference check and the cross-renderer check are
+    absent by necessity rather than by choice.
+
+    It therefore calls ``save_3d_views`` directly rather than ``plot_panels``,
+    which would compose the screenshot into a matplotlib figure and store an
+    upsampled interpolation of the render instead of the render. The flatmap
+    path pays that cost to keep both renderers comparable; without a second
+    renderer there is nothing to buy with it.
+
+    Curvature is left at pycortex's default (thresholded) here, unlike the
+    flatmap suites. Those un-threshold it to reduce cross-renderer
+    disagreement -- a reason that does not apply when there is no second
+    renderer -- so using the default recovers coverage of the default curvature
+    path, which the flatmap references explicitly do not provide.
+    """
+    from cortex.export.save_views import save_3d_views
+
+    if not REGENERATE_REFERENCES:
+        ref_path = reference_dir / f"webgl_{tag}{REFERENCE_SUFFIX}"
+        if not ref_path.exists():
+            pytest.skip(
+                f"No reference for webgl_{tag} in {reference_dir}. "
+                "See that directory's README for regeneration."
+            )
+
+    wg_path = save_3d_views(
+        view,
+        base_name=str(tmp_path / f"webgl_{tag}"),
+        list_angles=[angle],
+        list_surfaces=[surface],
+        trim=True,
+        size=(512, 384),
+        sleep=10,
+        viewer_params=dict(labels_visible=[], overlays_visible=[]),
+        headless=True,
+    )[0]
+
+    msg = _check_against_reference(
+        f"webgl_{tag}", Path(wg_path), tmp_path, reference_dir
+    )
+    if REGENERATE_REFERENCES:
+        pytest.skip(f"Regenerated webgl_{tag} in {reference_dir}")
+    return [msg] if msg is not None else []
+
+
+@pytest.mark.parametrize("name", DATAVIEW_NAMES)
+def test_visual_comparison_alpha_dataviews(tmp_path, name):
+    """Render an alpha-bearing dataview via quickshow + webgl, and assert it matches.
+
+    Plain Volume / Vertex have no native per-element alpha (pycortex's
+    bundled ``*_alpha`` colormaps are all 2D and only apply to the 2D
+    dataview types), so those two act as a no-alpha baseline. The other four
+    exercise alpha: Volume2D / Vertex2D via the 2D-alpha cmap ``RdBu_r_alpha``,
+    VolumeRGB / VertexRGB via the ``alpha=`` kwarg.
+
+    Compared both within-renderer (against a stored reference) and
+    cross-renderer (quickflat vs webgl); see ``_render_and_check_dataview``. A
+    mismatch leaves ``actual_*.png`` and an amplified ``diff_*.png`` in the
+    test's ``tmp_path``.
+    """
+    view = _build_alpha_dataview(name)
+    failures = _render_and_check_dataview(name, view, REFERENCE_DIR, tmp_path)
+    assert not failures, (
+        f"Renders differ from expectations:\n  "
+        + "\n  ".join(failures)
+        + f"\n\nFor details, see {tmp_path} and {REFERENCE_DIR.parent}/README.md"
+    )
+
+
+@pytest.mark.parametrize("name", DATAVIEW_NAMES)
+def test_visual_comparison_nan_dataviews(tmp_path, name):
+    """Render a NaN-bearing dataview via quickshow + webgl, and assert it matches.
+
+    NaN is pycortex's convention for "no data at this voxel/vertex" -- both
+    renderers are expected to draw those elements as fully transparent (falling
+    through to the curvature underlay) rather than mapping NaN through the
+    colormap as if it were a real value. This test renders the six dataview
+    classes with the *primary* data channel (not alpha) containing NaNs over
+    roughly half of each volume/surface.
+
+    Compared both within-renderer (against a stored reference) and
+    cross-renderer (quickflat vs webgl); see ``_render_and_check_dataview``. A
+    mismatch leaves ``actual_*.png`` and an amplified ``diff_*.png`` in the
+    test's ``tmp_path``.
+    """
+    view = _build_nan_dataview(name)
+    failures = _render_and_check_dataview(name, view, NAN_REFERENCE_DIR, tmp_path)
+    assert not failures, (
+        f"Renders differ from expectations:\n  "
+        + "\n  ".join(failures)
+        + f"\n\nFor details, see {tmp_path} and {NAN_REFERENCE_DIR.parent}/README.md"
+    )
+
+
+@pytest.mark.parametrize("name", NAN_ALPHA_DATAVIEW_NAMES)
+def test_visual_comparison_nan_alpha_dataviews(tmp_path, name):
+    """Render an RGB dataview whose alpha map carries NaNs, and assert it matches.
+
+    The other NaN suite puts NaNs in the data channels; this one puts them in the
+    alpha map. That is a distinct path -- alpha is not colour-mapped, it is used
+    directly as a blend weight, so the NaN lands in the compositing arithmetic
+    rather than in a colormap lookup. Only ``VolumeRGB``/``VertexRGB`` take an
+    explicit ``alpha=``, so only those two are covered.
+
+    Current behaviour, which these references encode, is that NaN-alpha elements
+    render fully transparent and the curvature underlay shows through -- the same
+    outcome as a NaN in the data.
+
+    Be aware that this behaviour is not settled. Work to unify NaN and alpha
+    handling across quickflat, WebGL and the RGB dataviews (cb976270, not in this
+    branch's history) changes how the surviving RGB is blended, without changing
+    the transparency itself. If that lands, expect these four references to need
+    regenerating; the transparency assertion should survive, the exact blend will
+    not.
+    """
+    view = _build_nan_alpha_dataview(name)
+    failures = _render_and_check_dataview(name, view, NAN_ALPHA_REFERENCE_DIR, tmp_path)
+    assert not failures, (
+        f"Renders differ from expectations:\n  "
+        + "\n  ".join(failures)
+        + f"\n\nFor details, see {tmp_path} and {NAN_ALPHA_REFERENCE_DIR.parent}/README.md"
+    )
+
+
+@pytest.mark.parametrize("surface,angle,name", NONFLAT_VIEWS)
+def test_visual_comparison_nonflat_views(tmp_path, surface, angle, name):
+    """Render a non-flatmap view through webgl and assert it matches its reference.
+
+    Everything else in this file renders the flatmap. That leaves the 3D views
+    covered only by test_webgl_headless.py's smoke tests, which assert the file
+    is larger than 1000 bytes -- enough to catch a render that never happened,
+    not one that came out wrong.
+
+    webgl only, necessarily: quickshow renders flatmaps and nothing else, so
+    these views have no quickflat counterpart and therefore no cross-renderer
+    check. Volume and Vertex both appear because they take different shader
+    paths, and the flatmap is demonstrably not representative of how those
+    behave -- both known webgl lighting bugs reproduce on flatmaps only.
+    """
+    view = _build_alpha_dataview(name)
+    tag = f"{surface}_{angle}_{name}"
+    failures = _render_and_check_webgl_only(
+        tag, view, surface, angle, NONFLAT_REFERENCE_DIR, tmp_path
+    )
+    assert not failures, (
+        f"Render differs from expectations:\n  "
+        + "\n  ".join(failures)
+        + f"\n\nFor details, see {tmp_path} and {NONFLAT_REFERENCE_DIR.parent}/README.md"
+    )
diff --git a/cortex/tests/test_webgl_headless.py b/cortex/tests/test_webgl_headless.py
index f97340ebc..446cec3a9 100644
--- a/cortex/tests/test_webgl_headless.py
+++ b/cortex/tests/test_webgl_headless.py
@@ -72,6 +72,79 @@ def _wait_for_file(path, timeout=30):
     raise RuntimeError(f"File {path!r} not written within {timeout}s")
 
 
+def _assert_no_browser_failures(handle):
+    """Fail on uncaught JS exceptions, or on WebGL reporting its own failure.
+
+    Two separate things, deliberately kept narrow:
+
+    - ``[pageerror]`` is an uncaught JS exception, which includes three.js
+      throwing on a WebGL context it could not create.
+    - ``webgl_failures`` additionally catches a shader that compiled but failed
+      to *link*, which three.js reports only via ``console.error`` -- nothing
+      raises, and the render silently comes back blank.
+
+    Everything else on the console is ignored on purpose. A healthy viewer
+    already logs a console.error for the Leap Motion websocket it cannot reach
+    (ws://127.0.0.1:6437), so asserting on every browser message fails every run.
+
+    **Call this after the viewer has been torn down.** Playwright's sync API only
+    dispatches queued events while something is calling into it, and the worker
+    thread parks in ``_shutdown_event.wait()`` for the viewer's whole lifetime
+    (``headless.py``). So console messages emitted during rendering sit in the
+    driver, undelivered, until ``_cleanup()`` closes the page. Called inside the
+    ``with`` block this sees only what arrived during page load -- measured: an
+    8-second wait mid-test delivered nothing, while closing the page delivered
+    everything at once. ``handle._pw_thread`` outlives the context manager, so
+    reading it afterwards is safe and is the only way to see the full list.
+
+    An earlier version of this note claimed the gh-714 link failure does not
+    reproduce under SwiftShader, on the evidence that capturing every console
+    message yielded nothing about shaders. That was wrong, and wrong for the
+    reason above: those captures all read *before* teardown, so they were
+    measuring the delivery bug rather than the driver. Read afterwards,
+    SwiftShader does report it -- "THREE.WebGLProgram: Could not initialise
+    shader." arrives in the teardown flush and this assertion fires on it.
+
+    Still, do not treat this as the only defence. It can only catch what the
+    driver chooses to report, and a driver that silently links an over-allocating
+    shader and draws nothing would say nothing at all. ``_assert_not_blank``
+    covers that by inspecting the rendered pixels instead.
+
+    Neither, though, covers a shader *variant* that no test ever renders. The
+    viewer generates 24 surface-shader variants (vertex/pixel x cmap/rgb/2d x
+    hasflat x equivolume) plus the pick and depth shaders; rendering the six
+    dataview types links only the handful of combinations those happen to use.
+    ``equivolume``, and the pick and depth shaders, are linked by nothing here.
+    ``test_webgl_shaders.py`` on the gh-714 branch covers that, linking all 26
+    variants directly; it is deliberately not duplicated into this branch.
+    """
+    from cortex.export.headless import webgl_failures
+
+    errors = handle._pw_thread.browser_errors
+    pageerrors = [e for e in errors if "[pageerror]" in e]
+    assert not pageerrors, f"JS errors: {pageerrors}"
+    failures = webgl_failures(errors)
+    assert not failures, f"WebGL reported a failure: {failures}"
+
+
+def _assert_not_blank(path):
+    """Fail if the render came out as a single flat color.
+
+    A shader that fails to compile or link (or geometry that never made it to
+    the GPU) leaves the canvas showing nothing but the background, which is
+    otherwise indistinguishable from a successful render: the png is written,
+    and no javascript exception is raised.
+    """
+    from PIL import Image
+
+    rgb = np.asarray(Image.open(path).convert("RGB")).reshape(-1, 3)
+    ncolors = len(np.unique(rgb, axis=0))
+    assert ncolors > 10, (
+        f"{path} has only {ncolors} distinct color(s); the brain was probably "
+        "never drawn."
+    )
+
+
 # ---------------------------------------------------------------------------
 # Group 1: Data type smoke tests
 # ---------------------------------------------------------------------------
@@ -90,9 +163,12 @@ def test_datatype_renders(dtype_name, tmp_path):
         _wait_for_file(outfile)
         assert os.path.isfile(outfile)
         assert os.path.getsize(outfile) > 0
-        # No uncaught JS errors
-        pageerrors = [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e]
-        assert len(pageerrors) == 0, f"JS errors: {pageerrors}"
+
+    # Both checks run outside the with, because teardown is what flushes
+    # Playwright's queued console events (see _assert_no_browser_failures).
+    # Browser errors are checked first, since one may explain a blank file.
+    _assert_no_browser_failures(handle)
+    _assert_not_blank(outfile)
 
 
 # ---------------------------------------------------------------------------
@@ -687,6 +763,9 @@ def _addData_viewer():
     vol1 = cortex.Volume(np.random.randn(*volshape), subj, xfmname)
     with cortex.export.headless_viewer(vol1, viewer_params={}) as handle:
         yield handle
+    # As above: only now are the queued console events delivered. The per-test
+    # calls inside the class see the page-load messages only.
+    _assert_no_browser_failures(handle)
 
 
 def _served_metadata(handle):
@@ -744,8 +823,7 @@ def test_adds_dataview(self, _addData_viewer):
         handle.addData(second=vol2)
         time.sleep(2)
 
-        pageerrors = [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e]
-        assert len(pageerrors) == 0, f"JS errors after addData: {pageerrors}"
+        _assert_no_browser_failures(handle)
 
         # "data" is the name webshow gives to a bare Dataview.
         assert set(handle.dataviews.attrs) == {"data", "second"}
@@ -800,8 +878,7 @@ def test_replaces_existing_name(self, _addData_viewer):
         assert len(metadata["images"]) == 2
         assert metadata["views"][1]["data"][0] not in previous_brains
 
-        pageerrors = [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e]
-        assert len(pageerrors) == 0, f"JS errors after addData: {pageerrors}"
+        _assert_no_browser_failures(handle)
 
     def test_rejects_unknown_subject(self, _addData_viewer):
         """Surfaces cannot be added to a running viewer, so neither can subjects."""
@@ -847,210 +924,5 @@ def test_addData_vertex_data(tmp_path):
         assert "mosaic" not in metadata["data"][vertex_name]
         assert _fetch(handle, metadata["images"][vertex_name][0])[1:6] == b"NUMPY"
 
-        pageerrors = [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e]
-        assert len(pageerrors) == 0, f"JS errors after addData: {pageerrors}"
-
-
-# ---------------------------------------------------------------------------
-# Group 10: Manual visual A/B comparison across all alpha-bearing dataviews
-# ---------------------------------------------------------------------------
-
-
-@pytest.mark.skipif(
-    not os.environ.get("RUN_VISUAL_COMPARISON"),
-    reason="Manual visual comparison; set RUN_VISUAL_COMPARISON=1 to run.",
-)
-def test_visual_comparison_alpha_dataviews(tmp_path):
-    """Render all 6 dataview types via quickshow + webgl, side-by-side.
-
-    Skipped by default — set ``RUN_VISUAL_COMPARISON=1`` to run. Builds a
-    grid where each row is one dataview type (Volume, Vertex, Volume2D,
-    Vertex2D, VolumeRGB, VertexRGB) and the two columns are the matplotlib
-    (``cortex.quickshow``) reference vs the headless WebGL flatmap render.
-    Used as a manual smoke check that the alpha-blend fix
-    (``Package``-side premultiply for VertexRGB + cmap-LUT
-    ``premultiplyAlpha=true`` for the 2D-cmap path) keeps both viewers in
-    visual agreement across every alpha-encoding pattern.
-
-    Plain Volume / Vertex have no native per-element alpha (pycortex's
-    bundled ``*_alpha`` colormaps are all 2D and only apply to the 2D
-    dataview types), so those two rows act as a no-alpha baseline. The
-    other four rows exercise alpha: Volume2D / Vertex2D via the 2D-alpha
-    cmap ``RdBu_r_alpha``, VolumeRGB / VertexRGB via the ``alpha=`` kwarg.
-
-    Renders are intentionally low-resolution (quickshow ``height=256``,
-    webgl ``size=(512, 384)``) so the final composite PNG stays small.
-    Both viewers run with no labels, no ROIs, and curvature underlay on.
-
-    The composite PNG is written under ``tmp_path`` and the absolute path
-    is printed at the end of the test so the file is easy to open.
-    """
-    import matplotlib.pyplot as plt
-
-    import cortex.polyutils
-
-    # ------- Synthesize data and alpha maps (mirrors plot_data_with_alpha.py) -
-
-    # Volumetric
-    zz, yy, xx = np.mgrid[0:31, 0:100, 0:100]
-    data_vol = (xx - 50) / 50.0  # ~ [-1, 1]
-    center = np.array([15, 50, 50])
-    sigma_v = 25.0
-    dist2 = (
-        (zz - center[0]) ** 2 + (yy - center[1]) ** 2 + (xx - center[2]) ** 2
-    )
-    accuracy_vol = np.exp(-dist2 / (2 * sigma_v**2))  # [0, 1] bump
-    red_vol = np.clip(xx / 99.0, 0, 1)
-    green_vol = np.clip(yy / 99.0, 0, 1)
-    blue_vol = np.clip(zz / 30.0, 0, 1)
-
-    # Surface (vertex) — encode by spatial coordinate, not vertex index
-    surfs = [
-        cortex.polyutils.Surface(*d)
-        for d in cortex.db.get_surf(subj, "fiducial")
-    ]
-    num_verts = [s.pts.shape[0] for s in surfs]
-    pts = np.vstack([surfs[0].pts, surfs[1].pts])
-    y_centered = pts[:, 1] - pts[:, 1].mean()
-    data_vtx = y_centered / np.abs(y_centered).max()  # [-1, 1]
-    xyz_norm = (pts - pts.min(axis=0)) / (pts.max(axis=0) - pts.min(axis=0))
-
-    def _bump(surf, seed, sigma):
-        d = np.linalg.norm(surf.pts - surf.pts[seed], axis=1)
-        return np.exp(-(d**2) / (2 * sigma**2))
-
-    accuracy_vtx = np.hstack(
-        [
-            _bump(surfs[0], num_verts[0] // 2, sigma=40.0),
-            _bump(surfs[1], num_verts[1] // 2, sigma=40.0),
-        ]
-    )
-
-    # ------- Build the six dataviews ----------------------------------------
-    # Volume / Vertex have no native per-element alpha — pycortex's bundled
-    # `*_alpha` colormaps are all 2D LUTs and only apply to Volume2D /
-    # Vertex2D. So plain Volume / Vertex use a non-alpha cmap (`viridis`)
-    # and serve as the no-alpha baseline; Volume2D / Vertex2D pair data
-    # against accuracy via the 2D-alpha cmap `RdBu_r_alpha`; VolumeRGB /
-    # VertexRGB use the native `alpha=` kwarg.
-
-    cmap_plain = "viridis"
-    cmap_2d = "RdBu_r_alpha"
-
-    dataviews = [
-        (
-            "Volume",
-            cortex.Volume(
-                data_vol, subj, xfmname,
-                cmap=cmap_plain, vmin=-1, vmax=1,
-            ),
-        ),
-        (
-            "Vertex",
-            cortex.Vertex(
-                data_vtx, subj,
-                cmap=cmap_plain, vmin=-1, vmax=1,
-            ),
-        ),
-        (
-            "Volume2D",
-            cortex.Volume2D(
-                data_vol, accuracy_vol, subj, xfmname,
-                cmap=cmap_2d,
-                vmin=-1, vmax=1, vmin2=0, vmax2=1,
-            ),
-        ),
-        (
-            "Vertex2D",
-            cortex.Vertex2D(
-                data_vtx, accuracy_vtx, subj,
-                cmap=cmap_2d,
-                vmin=-1, vmax=1, vmin2=0, vmax2=1,
-            ),
-        ),
-        (
-            "VolumeRGB",
-            cortex.VolumeRGB(
-                cortex.Volume(red_vol, subj, xfmname, vmin=0, vmax=1),
-                cortex.Volume(green_vol, subj, xfmname, vmin=0, vmax=1),
-                cortex.Volume(blue_vol, subj, xfmname, vmin=0, vmax=1),
-                subj, xfmname,
-                alpha=cortex.Volume(accuracy_vol, subj, xfmname, vmin=0, vmax=1),
-            ),
-        ),
-        (
-            "VertexRGB",
-            cortex.VertexRGB(
-                cortex.Vertex(xyz_norm[:, 0], subj, vmin=0, vmax=1),
-                cortex.Vertex(xyz_norm[:, 1], subj, vmin=0, vmax=1),
-                cortex.Vertex(xyz_norm[:, 2], subj, vmin=0, vmax=1),
-                subj,
-                alpha=cortex.Vertex(accuracy_vtx, subj, vmin=0, vmax=1),
-            ),
-        ),
-    ]
-
-    # ------- Render each dataview through both paths ------------------------
-    # Each WebGL render spins up its own headless browser via plot_panels;
-    # six sequential launches × ~15s sleep = ~90s+ end to end. That's fine
-    # for a manual A/B and avoids the broken `addData` path on headless.
-
-    n = len(dataviews)
-    fig, axes = plt.subplots(n, 2, figsize=(7, 2.2 * n))
-
-    flatmap_panel = [
-        {
-            "extent": [0.0, 0.0, 1.0, 1.0],
-            "view": {"angle": "flatmap", "surface": "flatmap"},
-        }
-    ]
-
-    for row, (name, view) in enumerate(dataviews):
-        # quickshow → low-res PNG
-        qs_path = tmp_path / f"qs_{name}.png"
-        qs_fig = cortex.quickshow(
-            view,
-            with_curvature=True,
-            with_rois=False,
-            with_labels=False,
-            with_colorbar=False,
-            with_sulci=False,
-            with_borders=False,
-            height=256,
-        )
-        qs_fig.savefig(qs_path, bbox_inches="tight", pad_inches=0, dpi=80)
-        plt.close(qs_fig)
-
-        # webgl → trimmed flatmap PNG via plot_panels (single flatmap panel)
-        wg_path = str(tmp_path / f"wg_{name}.png")
-        wg_fig = cortex.export.plot_panels(
-            view,
-            panels=flatmap_panel,
-            figsize=(6, 3),
-            windowsize=(512, 384),
-            save_name=wg_path,
-            sleep=10,
-            viewer_params=dict(labels_visible=[], overlays_visible=[]),
-            headless=True,
-        )
-        plt.close(wg_fig)
-
-        ax_qs, ax_wg = axes[row]
-        ax_qs.imshow(plt.imread(qs_path))
-        ax_qs.set_title(f"{name} — quickshow", fontsize=9)
-        ax_qs.axis("off")
-        ax_wg.imshow(plt.imread(wg_path))
-        ax_wg.set_title(f"{name} — webgl (flatmap)", fontsize=9)
-        ax_wg.axis("off")
-
-    fig.suptitle(
-        "Alpha-bearing dataviews: quickshow vs WebGL", fontsize=11,
-    )
-    fig.tight_layout()
-    out_path = tmp_path / "alpha_dataview_comparison.png"
-    fig.savefig(out_path, dpi=100, bbox_inches="tight")
-    plt.close(fig)
+        _assert_no_browser_failures(handle)
 
-    print(f"\nVisual comparison saved to:\n  {out_path}\n")
-    assert out_path.exists()
-    assert out_path.stat().st_size > 0
diff --git a/cortex/tests/test_webgl_shaders.py b/cortex/tests/test_webgl_shaders.py
new file mode 100644
index 000000000..9d4d8df25
--- /dev/null
+++ b/cortex/tests/test_webgl_shaders.py
@@ -0,0 +1,195 @@
+"""Tests that every webgl shader variant compiles and links.
+
+WebGL only guarantees 16 vertex attribute slots (``MAX_VERTEX_ATTRIBS``), and
+the surface shaders use nearly all of them. A shader that asks for one slot too
+many still compiles: it fails at *link* time, which three.js only reports on
+the browser console and which shows up in the viewer as an unexplained black
+screen. That is how ``Vertex2D`` data broke (gh-714), so every combination of
+options the viewer generates shaders with is linked here.
+
+These tests only need Chromium; no subject database or viewer is involved.
+"""
+
+import json
+import os
+
+import pytest
+
+import cortex.webgl
+from cortex.tests.testing_utils import has_playwright
+
+pytestmark = pytest.mark.skipif(
+    not has_playwright, reason="playwright and chromium are required"
+)
+
+JS_PATH = os.path.join(os.path.dirname(cortex.webgl.__file__), "resources", "js")
+
+# The declarations THREE.WebGLProgram prepends to every shader it builds
+# (three.js r69, resources/js/three.js). Only the ones the surface shaders
+# actually rely on are listed; a missing one shows up as a compile error rather
+# than as a silently passing test.
+VERTEX_PREFIX = """
+precision highp float;
+precision highp int;
+#define MAX_DIR_LIGHTS 3
+#define MAX_POINT_LIGHTS 0
+#define MAX_SPOT_LIGHTS 0
+#define MAX_HEMI_LIGHTS 0
+#define MAX_SHADOWS 0
+uniform mat4 modelMatrix;
+uniform mat4 modelViewMatrix;
+uniform mat4 projectionMatrix;
+uniform mat4 viewMatrix;
+uniform mat3 normalMatrix;
+uniform vec3 cameraPosition;
+attribute vec3 position;
+attribute vec3 normal;
+attribute vec2 uv;
+attribute vec2 uv2;
+"""
+
+FRAGMENT_PREFIX = """
+precision highp float;
+precision highp int;
+#define MAX_DIR_LIGHTS 3
+#define MAX_POINT_LIGHTS 0
+#define MAX_SPOT_LIGHTS 0
+#define MAX_HEMI_LIGHTS 0
+#define MAX_SHADOWS 0
+uniform mat4 viewMatrix;
+uniform vec3 cameraPosition;
+"""
+
+# Loads the viewer's shader library into a page and exposes a hook that builds
+# one shader variant and links it, the way THREE.WebGLProgram does.
+PAGE = """
+<html><body><canvas id="c" width="32" height="32"></canvas>
+<script src="file://__JSDIR__/three.js"></script>
+<script src="file://__JSDIR__/shaderlib.js"></script>
+<script>
+var gl = document.getElementById('c').getContext('webgl');
+// THREE.WebGLRenderer asks for these too; the fragment shaders use fwidth
+// (derivatives) and float textures.
+gl.getExtension('OES_standard_derivatives');
+gl.getExtension('OES_texture_float');
+var VERTEX_PREFIX = __VERTEX_PREFIX__;
+var FRAGMENT_PREFIX = __FRAGMENT_PREFIX__;
+
+function compile(type, source) {
+    var shader = gl.createShader(type);
+    gl.shaderSource(shader, source);
+    gl.compileShader(shader);
+    return shader;
+}
+
+window.linkShader = function(shadername, opts) {
+    var code = Shaders[shadername](opts);
+    // The pick shader returns one fragment shader per axis; any of them will
+    // do, they all go with the vertex shader that holds the attributes.
+    var frag = code.fragment instanceof Array ? code.fragment[0] : code.fragment;
+    var vs = compile(gl.VERTEX_SHADER, VERTEX_PREFIX + code.vertex);
+    var fs = compile(gl.FRAGMENT_SHADER, FRAGMENT_PREFIX + frag);
+    var program = gl.createProgram();
+    gl.attachShader(program, vs);
+    gl.attachShader(program, fs);
+    gl.linkProgram(program);
+
+    var result = {
+        compiled: !!(gl.getShaderParameter(vs, gl.COMPILE_STATUS) &&
+                     gl.getShaderParameter(fs, gl.COMPILE_STATUS)),
+        linked: !!gl.getProgramParameter(program, gl.LINK_STATUS),
+        log: [gl.getShaderInfoLog(vs), gl.getShaderInfoLog(fs),
+              gl.getProgramInfoLog(program)].join("\\n"),
+        max_attributes: gl.getParameter(gl.MAX_VERTEX_ATTRIBS),
+        attributes: [],
+    };
+    var nattr = gl.getProgramParameter(program, gl.ACTIVE_ATTRIBUTES) || 0;
+    for (var i = 0; i < nattr; i++)
+        result.attributes.push(gl.getActiveAttrib(program, i).name);
+    return result;
+};
+</script></body></html>
+"""
+
+# The options the viewer generates surface shaders with. ``morphs`` is the
+# number of surfaces to mix between (anatomical, inflated and flat), ``volume``
+# says the subject has a white matter surface; the rest come from the dataview
+# and from the surface menu.
+SURFACE_OPTS = dict(morphs=3, volume=1, layers=1, rois=True, extratex=False,
+                    halo=False, dither=False, voxline=False, sampler="nearest")
+
+
+def _surface_variants():
+    """Every (shader, opts) pair the viewer can ask for a surface shader."""
+    for shader in ("surface_vertex", "surface_pixel"):
+        for rgb in (False, True):
+            for twod in (False, True):
+                if rgb and twod:
+                    continue  # RGB data has no second dimension
+                for hasflat in (False, True):
+                    for equivolume in (False, True):
+                        opts = dict(SURFACE_OPTS, rgb=rgb, twod=twod,
+                                    hasflat=hasflat, equivolume=equivolume)
+                        name = "%s-%s%s%s%s" % (
+                            shader,
+                            "rgb" if rgb else "cmap",
+                            "-2d" if twod else "",
+                            "-flat" if hasflat else "",
+                            "-equivolume" if equivolume else "",
+                        )
+                        yield pytest.param(shader, opts, id=name)
+
+
+def _variants():
+    yield from _surface_variants()
+    # The shaders the picker renders with; they morph the same geometry but
+    # carry no data.
+    yield pytest.param("pick", dict(morphs=3, volume=1), id="pick")
+    yield pytest.param("depth", dict(morphs=3, volume=1), id="depth")
+
+
+@pytest.fixture(scope="module")
+def link_shader(tmp_path_factory):
+    """Return a function linking one shader variant in a real GL context."""
+    from playwright.sync_api import sync_playwright
+
+    page_path = tmp_path_factory.mktemp("shaders") / "shaders.html"
+    page_path.write_text(
+        PAGE.replace("__JSDIR__", JS_PATH)
+        .replace("__VERTEX_PREFIX__", json.dumps(VERTEX_PREFIX))
+        .replace("__FRAGMENT_PREFIX__", json.dumps(FRAGMENT_PREFIX))
+    )
+
+    with sync_playwright() as playwright:
+        browser = playwright.chromium.launch(
+            headless=True,
+            args=["--enable-webgl", "--use-gl=swiftshader", "--no-sandbox",
+                  "--disable-dev-shm-usage"],
+        )
+        page = browser.new_page()
+        page.goto("file://%s" % page_path, wait_until="load", timeout=60000)
+        if not page.evaluate("() => !!window.linkShader"):
+            browser.close()
+            pytest.skip("no WebGL context available in this browser")
+        yield lambda shader, opts: page.evaluate(
+            "args => window.linkShader(args[0], args[1])", [shader, opts]
+        )
+        browser.close()
+
+
+@pytest.mark.parametrize("shader,opts", list(_variants()))
+def test_shader_links(shader, opts, link_shader):
+    """Each shader variant has to compile *and* link.
+
+    A variant that uses more vertex attributes than the driver has slots for
+    compiles fine and fails to link, which leaves the viewer showing nothing at
+    all.
+    """
+    result = link_shader(shader, opts)
+    assert result["compiled"], "%s did not compile:\n%s" % (shader, result["log"])
+    assert result["linked"], (
+        "%s compiled but did not link, using %d of the %d available vertex "
+        "attributes:\n%s" % (shader, len(result["attributes"]),
+                             result["max_attributes"], result["log"])
+    )
+    assert len(result["attributes"]) <= result["max_attributes"]
diff --git a/cortex/webgl/resources/js/mriview_surface.js b/cortex/webgl/resources/js/mriview_surface.js
index fdbea89b2..23c748a9b 100644
--- a/cortex/webgl/resources/js/mriview_surface.js
+++ b/cortex/webgl/resources/js/mriview_surface.js
@@ -210,13 +210,18 @@ var mriview = (function(module) {
                 pialareas = module.iterativelySmoothVertexData(hemi.attributes.position, hemi.attributes.index, hemi.offsets, pialareas, areasmoothfactor, areasmoothiter);
                 hemi.pialareas = pialareas;
 
-		var pialarea_attr = new THREE.BufferAttribute(pialareas, 1);
-                pialarea_attr.needsUpdate = true;
-                var wmarea_attr = new THREE.BufferAttribute(wmareas, 1);
-                wmarea_attr.needsUpdate = true;
-		
-                hemi.addAttribute('pialarea', pialarea_attr);
-                hemi.addAttribute('wmarea', wmarea_attr);
+                //The vertex areas are only read by the equivolume depth
+                //sampling in the vertex shaders, and WebGL only guarantees 16
+                //vertex attribute slots -- which these shaders already use up
+                //-- so they ride in the two unused components of auxdat
+                //(x is the medial wall mask, y the curvature) instead of
+                //taking two slots of their own.
+                var auxdat = hemi.attributes.auxdat;
+                for (var v = 0; v < wmareas.length; v++) {
+                    auxdat.array[v*4+2] = wmareas[v];
+                    auxdat.array[v*4+3] = pialareas[v];
+                }
+                auxdat.needsUpdate = true;
 
                 if (this.flatlims !== undefined) {
                     var flats = this._makeFlat(hemi.attributes.uv.array, json.flatlims, names[name]);
@@ -257,9 +262,20 @@ var mriview = (function(module) {
                     // console.log(flatoff_geom);
                     // this.flatoff = flatoff_geom;
 
-                    // hemi.addAttribute('flatBumpNorms', flatoff_geom.attributes.normal);
-                    hemi.addAttribute('flatheight', flatheights);
-                    hemi.addAttribute('flatBumpNorms', module.computeNormal(flat_offset_verts, hemi.attributes.index, hemi.offsets) );
+                    //Same story as the areas above: the bump height rides in
+                    //the fourth component of the bumped normals rather than in
+                    //an attribute of its own.
+                    var bumpnorms = module.computeNormal(flat_offset_verts, hemi.attributes.index, hemi.offsets);
+                    var flatbump = new Float32Array(flatheights.array.length * 4);
+                    for (var v = 0; v < flatheights.array.length; v++) {
+                        flatbump[v*4]   = bumpnorms.array[v*3];
+                        flatbump[v*4+1] = bumpnorms.array[v*3+1];
+                        flatbump[v*4+2] = bumpnorms.array[v*3+2];
+                        flatbump[v*4+3] = flatheights.array[v];
+                    }
+                    var flatbump_attr = new THREE.BufferAttribute(flatbump, 4);
+                    flatbump_attr.needsUpdate = true;
+                    hemi.addAttribute('flatbump', flatbump_attr);
                 } else {
                     // Fill these attributes so the shader doesn't choke, even though
                     // there's no flatmap
diff --git a/cortex/webgl/resources/js/shaderlib.js b/cortex/webgl/resources/js/shaderlib.js
index dc4deb82e..8f0b0cc63 100644
--- a/cortex/webgl/resources/js/shaderlib.js
+++ b/cortex/webgl/resources/js/shaderlib.js
@@ -174,19 +174,24 @@ var Shaderlib = (function() {
             return glsl;
         },
 
-        // thickmixer: header code that loads the uniforms and attributes needed to
-        // do equivolume sampling, for vertex shaders
+        // thickmixer: header code that loads the uniforms needed to do
+        // equivolume sampling, for vertex shaders. The white matter and pial
+        // vertex areas it needs ride along in auxdat.zw, which mriview_surface
+        // fills in when the surfaces load, rather than in attributes of their
+        // own: WebGL only guarantees 16 vertex attribute slots and these
+        // shaders are right up against that limit, so whatever fits in the
+        // spare components of an attribute that is already there goes there.
         thickmixer: [
             "uniform float thickmix;",
             "uniform int equivolume;",
-            "attribute float wmarea;",
-            "attribute float pialarea;",
         ].join("\n"),
 
         // thickmixer_main: translates a desired volume fraction into linear mixing
-        // parameter.
+        // parameter. Requires auxdat to be declared by the including shader.
         thickmixer_main: [
             "#ifdef EQUIVOLUME",
+                "float wmarea = auxdat.z;",
+                "float pialarea = auxdat.w;",
                 "float use_thickmix = 1. - (1. / (pialarea - wmarea) * (-1. * wmarea + sqrt((1. - thickmix) * pialarea * pialarea + thickmix * wmarea * wmarea)));",
             "#else",
                 "float use_thickmix = thickmix;",
@@ -368,8 +373,9 @@ var Shaderlib = (function() {
             "attribute vec4 auxdat;",
 
             "#ifdef HASFLAT",
-                "attribute vec3 flatBumpNorms;",
-                "attribute float flatheight;",
+                //xyz: normal of the bump-displaced flatmap, w: bump height.
+                //Packed into one attribute to stay under the 16 slot limit.
+                "attribute vec4 flatbump;",
             "#endif",
             // "attribute float dropout;",
             
@@ -425,20 +431,17 @@ var Shaderlib = (function() {
                 "vec3 pos, norm;",
                 "mixfunc(mpos, mnorm, pos, norm);",
 
-                // "norm = mix(flatBumpNorms, normalize(onorm), thickmix);",
-                // "norm = normalize(flatBumpNorms);",
-
             "#ifdef CORTSHEET",
                 // 
                 "#ifdef HASFLAT",
-                    "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * mix(1., 0., use_thickmix) * flatheight * f_bumpyflat;",
+                    "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * mix(1., 0., use_thickmix) * flatbump.w * f_bumpyflat;",
                 "#else",
                     "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * .62 * distance(position, wm.xyz) * mix(1., 0., use_thickmix);",
                 "#endif",
             "#endif",
 
                 "#ifdef HASFLAT",
-                    "vNormal = normalMatrix * mix(norm, flatBumpNorms, (1.0 - use_thickmix) * clamp(surfmix*"+(morphs-1)+". - 1.0, 0., 1.) * f_bumpyflat);",
+                    "vNormal = normalMatrix * mix(norm, flatbump.xyz, (1.0 - use_thickmix) * clamp(surfmix*"+(morphs-1)+". - 1.0, 0., 1.) * f_bumpyflat);",
                 "#else",
                     "vNormal = normalMatrix * norm;",
                 "#endif",
@@ -665,14 +668,9 @@ var Shaderlib = (function() {
                 wm: { type: 'v4', value:null },
                 wmnorm: { type: 'v3', value:null },
                 auxdat: { type: 'v4', value:null },
-                wmarea: { type: 'f', value:null },
-                pialarea: { type: 'f', value:null },
-                // flatBumpNorms: { type: 'v3', value:null },
-                // flatheight: { type: 'f', value:null },
             };
             if (opts.hasflat) {
-                attributes.flatBumpNorms = { type: 'v3', value:null };
-                attributes.flatheight = { type: 'f', value:null };
+                attributes.flatbump = { type: 'v4', value:null };
             }
             for (var i = 0; i < morphs-1; i++) {
                 attributes['mixSurfs'+i] = { type:'v4', value:null};
@@ -730,8 +728,9 @@ var Shaderlib = (function() {
             "attribute vec4 auxdat;",
 
             "#ifdef HASFLAT",
-                "attribute vec3 flatBumpNorms;",
-                "attribute float flatheight;",
+                //xyz: normal of the bump-displaced flatmap, w: bump height.
+                //Packed into one attribute to stay under the 16 slot limit.
+                "attribute vec4 flatbump;",
             "#endif",
             // "attribute float dropout;",
             
@@ -788,14 +787,14 @@ var Shaderlib = (function() {
 
             "#ifdef CORTSHEET",
                 "#ifdef HASFLAT",
-                    "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * mix(1., 0., use_thickmix) * flatheight * f_bumpyflat;",
+                    "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * mix(1., 0., use_thickmix) * flatbump.w * f_bumpyflat;",
                 "#else",
                     "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * .62 * distance(position, wm.xyz) * mix(1., 0., use_thickmix);",
                 "#endif",
             "#endif",
 
                 "#ifdef HASFLAT",
-                    "vNormal = normalMatrix * mix(norm, flatBumpNorms, (1.0 - use_thickmix) * clamp(surfmix*"+(morphs-1)+". - 1.0, 0., 1.) * f_bumpyflat);",
+                    "vNormal = normalMatrix * mix(norm, flatbump.xyz, (1.0 - use_thickmix) * clamp(surfmix*"+(morphs-1)+". - 1.0, 0., 1.) * f_bumpyflat);",
                 "#else",
                     "vNormal = normalMatrix * norm;",
                 "#endif",
@@ -888,13 +887,10 @@ var Shaderlib = (function() {
                 wm: { type: 'v4', value:null },
                 wmnorm: { type: 'v3', value:null },
                 auxdat: { type: 'v4', value:null },
-                wmarea: { type: 'f', value:null },
-                pialarea: { type: 'f', value:null },
             };
 
             if (opts.hasflat) {
-                attributes.flatBumpNorms = { type: 'v3', value:null };
-                attributes.flatheight = { type: 'f', value:null };
+                attributes.flatbump = { type: 'v4', value:null };
             }
 
             for (var i = 0; i < 4; i++)
diff --git a/setup.py b/setup.py
index 91f934e95..cb50e7aa0 100644
--- a/setup.py
+++ b/setup.py
@@ -141,6 +141,13 @@ def run(self):
       # Don't use `extras_require` here. Put them in pyproject.toml .
       cmdclass=dict(install=my_install),
       include_package_data=True,
+      # Reference renders for the visual-regression test are fixtures, of no use
+      # at runtime. MANIFEST.in's recursive-include keeps them in the source
+      # tarball, so a build from source can still run the test; this keeps them
+      # out of the wheel, where they would just sit in every user's
+      # site-packages. The test skips when they are absent.
+      exclude_package_data={'cortex.tests': ['reference_images/*',
+                                             'reference_images/*/*']},
       classifiers=[
           'Development Status :: 6 - Mature',
           'Intended Audience :: Science/Research',
diff --git a/temp.txt b/temp.txt
new file mode 100644
index 000000000..e69de29bb