diff --git a/.github/workflows/build_docs.yml b/.github/workflows/build_docs.yml
index 004d8e542..fcf547083 100644
--- a/.github/workflows/build_docs.yml
+++ b/.github/workflows/build_docs.yml
@@ -11,11 +11,15 @@ on:
- main
workflow_dispatch: # Manual trigger for publishing docs
+concurrency:
+ group: ${{ github.workflow }}-${{ github.head_ref || github.ref }}
+ cancel-in-progress: ${{ github.event_name == 'pull_request' }}
+
jobs:
build-docs:
runs-on: ubuntu-latest
steps:
- - uses: actions/checkout@v6
+ - uses: actions/checkout@v7
with:
fetch-depth: 0 # Required for setuptools-scm to get version from tags
@@ -48,6 +52,8 @@ jobs:
with:
path: ~/.cache/ms-playwright
key: ${{ runner.os }}-playwright-${{ hashFiles('**/pyproject.toml') }}
+ restore-keys: |
+ ${{ runner.os }}-playwright-
- name: Install Playwright Chromium
run: playwright install --with-deps --only-shell chromium
diff --git a/.github/workflows/codespell.yml b/.github/workflows/codespell.yml
index c2416cea7..c40d06fd1 100644
--- a/.github/workflows/codespell.yml
+++ b/.github/workflows/codespell.yml
@@ -7,6 +7,10 @@ on:
pull_request:
branches: [main]
+concurrency:
+ group: ${{ github.workflow }}-${{ github.head_ref || github.ref }}
+ cancel-in-progress: ${{ github.event_name == 'pull_request' }}
+
permissions:
contents: read
@@ -17,6 +21,6 @@ jobs:
steps:
- name: Checkout
- uses: actions/checkout@v6
+ uses: actions/checkout@v7
- name: Codespell
uses: codespell-project/actions-codespell@v2
diff --git a/.github/workflows/install_from_wheel.yml b/.github/workflows/install_from_wheel.yml
index 8656db430..81a5596fb 100644
--- a/.github/workflows/install_from_wheel.yml
+++ b/.github/workflows/install_from_wheel.yml
@@ -8,6 +8,10 @@ on:
branches:
- main
+concurrency:
+ group: ${{ github.workflow }}-${{ github.head_ref || github.ref }}
+ cancel-in-progress: ${{ github.event_name == 'pull_request' }}
+
jobs:
install-from-wheel:
runs-on: ubuntu-latest
@@ -17,7 +21,7 @@ jobs:
max-parallel: 5
steps:
- - uses: actions/checkout@v6
+ - uses: actions/checkout@v7
- name: Set up Python
uses: actions/setup-python@v6
with:
diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
index 8d4a410c9..13d5c183d 100644
--- a/.github/workflows/publish.yml
+++ b/.github/workflows/publish.yml
@@ -12,7 +12,7 @@ jobs:
contents: read # Required to checkout the repository
steps:
- - uses: actions/checkout@v6
+ - uses: actions/checkout@v7
with:
fetch-depth: 0 # Required for setuptools-scm to get version from tags
diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml
index c66bef78c..6a4a8a763 100644
--- a/.github/workflows/run_tests.yml
+++ b/.github/workflows/run_tests.yml
@@ -8,6 +8,10 @@ on:
branches:
- main
+concurrency:
+ group: ${{ github.workflow }}-${{ github.head_ref || github.ref }}
+ cancel-in-progress: ${{ github.event_name == 'pull_request' }}
+
jobs:
run-tests:
runs-on: ubuntu-latest
@@ -17,7 +21,7 @@ jobs:
max-parallel: 5
steps:
- - uses: actions/checkout@v6
+ - uses: actions/checkout@v7
- name: Set up Python
uses: actions/setup-python@v6
with:
@@ -47,16 +51,18 @@ jobs:
with:
path: ~/.cache/ms-playwright
key: ${{ runner.os }}-playwright-${{ hashFiles('**/pyproject.toml') }}
+ restore-keys: |
+ ${{ runner.os }}-playwright-
- name: Install Playwright Chromium
run: playwright install --with-deps --only-shell chromium
- name: Test with pytest
- timeout-minutes: 20 # in case webviewer tests hang
+ timeout-minutes: 25 # in case webviewer tests hang
run: pytest --cov=./
- name: Upload coverage to Codecov
- uses: codecov/codecov-action@v6
+ uses: codecov/codecov-action@v7
with:
env_vars: OS,PYTHON
fail_ci_if_error: true
diff --git a/.gitignore b/.gitignore
index b54f158ba..613aa48c5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -69,3 +69,7 @@ docs/colormaps.rst
# Git worktrees
.worktrees
+
+# Claude Code working directory (per-worktree scratch: launchers, verify
+# scripts, plans). Not part of the project.
+.claude
diff --git a/cortex/dataset/braindata.py b/cortex/dataset/braindata.py
index 7cbe94efb..817c024bb 100644
--- a/cortex/dataset/braindata.py
+++ b/cortex/dataset/braindata.py
@@ -1,4 +1,5 @@
import hashlib
+import warnings
from copy import deepcopy
import sys
from typing import Generic, Optional, TypeVar, Union, cast
@@ -552,7 +553,33 @@ def right(self):
def blend_curvature(self, alpha, threshold=0, brightness=0.5,
contrast=0.25, smooth=20):
"""Blend the data with a curvature map depending on a transparency map.
-
+
+ .. deprecated::
+ Per-vertex/voxel alpha is now honored directly by both the WebGL
+ viewer and ``cortex.quickshow``, so this curvature-blending hack
+ is no longer needed. The recommended replacement for scalar data
+ with a transparency map is :class:`Vertex2D` (or
+ :class:`Volume2D`) with a 2D colormap whose second axis encodes
+ alpha (e.g. ``"fire_alpha"``, ``"PU_RdBu_covar_alpha"``)::
+
+ # Was:
+ # blended = vtx.blend_curvature(alpha)
+ # cortex.quickshow(blended)
+ # Now:
+ v2d = cortex.Vertex2D(vtx.data, alpha, subject,
+ cmap="fire_alpha",
+ vmin=vtx.vmin, vmax=vtx.vmax,
+ vmin2=0, vmax2=1)
+ cortex.quickshow(v2d) # or cortex.webgl.show(v2d)
+
+ The 2D colormap path keeps colormap parameters (``cmap``,
+ ``vmin``, ``vmax``) editable on the resulting object, and the
+ curvature underlay is composited through automatically by both
+ the matplotlib and WebGL renderers.
+
+ For data that is already RGB, pass ``alpha=`` to
+ :class:`VertexRGB` / :class:`VolumeRGB` directly instead.
+
Vertex objects cannot use transparency as Volume objects. This method
is a hack to mimic the transparency of Volume objects, blending the
Vertex data with a curvature map. This method returns a VertexRGB
@@ -577,6 +604,19 @@ def blend_curvature(self, alpha, threshold=0, brightness=0.5,
blended : VertexRGB object
The original map blended with a curvature map.
"""
+ warnings.warn(
+ "blend_curvature is deprecated and will be removed in a future "
+ "release. Per-vertex/voxel alpha is now honored directly by both "
+ "the WebGL viewer and quickshow, so this curvature-blending hack "
+ "is no longer needed. For scalar data with a transparency map, "
+ "use Vertex2D / Volume2D with a 2D colormap whose second axis "
+ "encodes alpha (e.g. 'fire_alpha', 'PU_RdBu_covar_alpha'), e.g. "
+ "`Vertex2D(data, alpha, subject, cmap='fire_alpha', vmin=..., "
+ "vmax=..., vmin2=0, vmax2=1)`. For data that is already RGB, "
+ "pass `alpha=` to VertexRGB / VolumeRGB directly.",
+ DeprecationWarning,
+ stacklevel=2,
+ )
from .views import Vertex
from .viewRGB import VertexRGB
# prepare curvature map
diff --git a/cortex/export/headless.py b/cortex/export/headless.py
index 9145e2c20..d3e8f2337 100644
--- a/cortex/export/headless.py
+++ b/cortex/export/headless.py
@@ -42,6 +42,7 @@
import contextlib
import logging
import threading
+import time
from typing import Any, Mapping, Optional
import cortex
@@ -50,6 +51,49 @@
logger = logging.getLogger(__name__)
+def _wait_for_viewer_loaded(handle, timeout: float = 60.0) -> None:
+ """Block until ``window.viewer.loaded`` resolves in the browser.
+
+ The WebSocket "connect" message fires as soon as the page DOM is ready,
+ but the WebGL viewer's CTM mesh download, parse, and initial render
+ typically take a few more seconds. ``viewer.loaded`` is a jQuery
+ Deferred that resolves at the end of ``setData`` (mriview.js), so polling
+ its ``.state()`` is a faithful "the viewer is ready to be driven"
+ signal that replaces brittle fixed sleeps in callers and tests.
+
+ We bypass the JSProxy attribute machinery and call ``send`` directly so
+ each poll is exactly one WebSocket roundtrip (the JSProxy ``__getattr__``
+ path would issue an extra ``query`` call per poll).
+ """
+ deadline = time.monotonic() + timeout
+ poll_interval = 0.1
+ last_err: Optional[str] = None
+ while time.monotonic() < deadline:
+ try:
+ result = handle.send(
+ method="run", params=["window.viewer.loaded.state", []]
+ )
+ except Exception as exc:
+ last_err = repr(exc)
+ result = None
+ # WebApp.send always wraps its single per-client response in a list;
+ # unpack the leaf so the "resolved"/"pending"/error-dict check below
+ # is straightforward and last_err carries the actual value seen.
+ val = result[0] if isinstance(result, list) and result else result
+ if val == "resolved":
+ return
+ if val is not None:
+ last_err = (
+ str(val.get("error", val)) if isinstance(val, dict) else str(val)
+ )
+ time.sleep(poll_interval)
+ raise RuntimeError(
+ f"Viewer's .loaded deferred did not resolve within {timeout:.0f}s "
+ f"(last response: {last_err!r}). The CTM mesh may have failed to "
+ "download or parse, or mriview.js failed to initialise."
+ )
+
+
# --------------------------------------------------------------------------- #
# Helper: run Playwright in a dedicated thread to avoid asyncio conflicts #
# --------------------------------------------------------------------------- #
@@ -337,6 +381,12 @@ def headless_viewer(
# any point during the session (each call returns a fresh snapshot).
handle._pw_thread = pw_thread
+ # Block until the WebGL viewer has finished initialising (CTM mesh
+ # download + parse + first setData). Replaces ad-hoc time.sleep(10)
+ # calls in tests and callers, and shortens the wait when the
+ # browser is faster than the worst-case timeout.
+ _wait_for_viewer_loaded(handle, timeout=timeout)
+
yield handle
finally:
diff --git a/cortex/export/save_views.py b/cortex/export/save_views.py
index 2ea757510..c9a1e7710 100644
--- a/cortex/export/save_views.py
+++ b/cortex/export/save_views.py
@@ -114,8 +114,12 @@ def save_3d_views(
cm = contextlib.nullcontext(cortex.webshow(volume, **viewer_params))
with cm as handle:
- # Wait for the viewer to be loaded
- time.sleep(sleep)
+ # Wait for the viewer to be loaded. The headless context manager
+ # already blocks on ``viewer.loaded`` before yielding, so we only
+ # need this fixed sleep for the interactive (real-browser) path
+ # where the user is opening the page manually.
+ if not headless:
+ time.sleep(sleep)
# Add interpolation and layers params only if we have a volume
if isinstance(volume, (cortex.Volume, cortex.Volume2D, cortex.VolumeRGB)):
diff --git a/cortex/freesurfer.py b/cortex/freesurfer.py
index f30c81c9f..b90d7546b 100644
--- a/cortex/freesurfer.py
+++ b/cortex/freesurfer.py
@@ -16,8 +16,7 @@
import nibabel
import numpy as np
from nibabel import gifti
-from scipy.linalg import lstsq
-from scipy.sparse import coo_matrix
+from scipy.sparse import coo_matrix, diags
from scipy.spatial import KDTree
from . import anat, database
@@ -705,128 +704,168 @@ def mri_surf2surf(data, source_subj, target_subj, hemi, subjects_dir=None):
return output_data
-def get_mri_surf2surf_matrix(source_subj, hemi, surface_type,
- target_subj='fsaverage', subjects_dir=None,
- n_neighbors=20, random_state=0,
- n_test_images=40, coef_threshold=None,
- renormalize=True):
+def _read_sphere_reg(subject, hemi, subjects_dir=None):
+ """Read the registered sphere (``?h.sphere.reg``) vertex coordinates.
+
+ These are the coordinates on which freesurfer's spherical registration
+ defines the cross-subject vertex correspondence used by ``mri_surf2surf``.
+ """
+ surf_file = get_paths(subject, hemi, 'surf',
+ freesurfer_subject_dir=subjects_dir).format(
+ name='sphere.reg')
+ pts, _ = parse_surf(surf_file)
+ return pts
+
+
+def _surf2surf_nnfr_matrix(src_sphere, trg_sphere):
+ """Build the freesurfer ``nnfr`` surf2surf matrix from sphere coordinates.
+
+ Implements freesurfer's nearest-neighbor, forward-and-reverse (``nnfr``)
+ mapping directly from the registered-sphere geometry:
+
+ * **Forward**: every target vertex draws its value from its single nearest
+ source vertex on the sphere. This guarantees every target gets a value.
+ * **Reverse**: any source vertex that was *not* selected by some target in
+ the forward pass (an "orphan") is folded into its nearest target vertex,
+ so that no source vertex's data is silently dropped. This is the
+ "forward and reverse" part of ``nnfr`` and is what makes downsampling
+ (e.g. to a coarser subject) an average rather than a subsampling.
+
+ Each target row is then renormalized so that it is the *mean* of its
+ contributing source vertices.
- """Creates a matrix implementing freesurfer mri_surf2surf command.
-
- A surface-to-surface transform is a linear transform between vertex spaces.
- Such a transform must be highly localized in the sense that a vertex in the
- target surface only draws its values from very few source vertices.
- This function exploits the localization to create an inverse problem for
- each vertex.
- The source neighborhoods for each target vertex are found by using
- mri_surf2surf to transform the three coordinate maps from the source
- surface to the target surface, yielding three coordinate values for each
- target vertex, for which we find the nearest neighbors in the source space.
- A small number of test images is transformed from source surface to
- target surface.
- For each target vertex in the transformed test images, a regression is
- performed using only the corresponding source image neighborhood, yielding
- the entries for a sparse matrix encoding the transform.
-
Parameters
- ==========
-
- source_subj: str
- Freesurfer name of source subject
-
- hemi: str in ("lh", "rh")
- Indicator for hemisphere
-
- surface_type: str in ("white", "pial", ...)
- Indicator for surface layer
-
- target_subj: str, default "fsaverage"
- Freesurfer name of target subject
-
- subjects_dir: str, default os.environ["SUBJECTS_DIR"]
- The freesurfer subjects directory
-
- n_neighbors: int, default 20
- The size of the neighborhood to take into account when estimating
- the source support of a vertex
-
- random_state: int, default 0
- Random number generator or seed for generating test images
-
- n_test_images: int, default 40
- Number of test images transformed to compute inverse problem. This
- should be greater than n_neighbors or equal.
-
- coef_treshold: float, default 1 / (10 * n_neighbors)
- Value under which to set a weight to zero in the inverse problem.
-
- renormalize: boolean, default True
- Determines whether the rows of the output matrix should add to 1,
- implementing what is sensible: a weighted averaging
-
+ ----------
+ src_sphere : ndarray, shape (n_src, 3)
+ Source ``sphere.reg`` vertex coordinates.
+ trg_sphere : ndarray, shape (n_trg, 3)
+ Target ``sphere.reg`` vertex coordinates.
+
+ Returns
+ -------
+ matrix : scipy.sparse.csr_matrix, shape (n_trg, n_src)
+ Linear operator mapping source vertex data to target vertex data via
+ ``target_data = matrix.dot(source_data)``.
+
Notes
- =====
- It turns out that freesurfer seems to do the following: For each target
- vertex, find, on the sphere, the nearest source vertices, and average their
- values. Try to be as one-to-one as possible.
+ -----
+ For ordinary subjects (and full-resolution ``fsaverage``), the resulting
+ matrix reproduces ``mri_surf2surf`` bit-for-bit, because the spheres are
+ irregular enough that nearest-neighbor assignment is unambiguous.
+
+ The icosahedrally-subsampled targets (``fsaverage6``/``5``/``4``/``3``)
+ are the one exception: their meshes are perfectly regular, so a sizeable
+ number of fine vertices land *exactly* equidistant between two coarse
+ target vertices. Freesurfer breaks these exact ties using its internal
+ vertex ordering, which this implementation does not reproduce. The result
+ is that a small fraction (~0.3-2%) of coarse vertices average a different,
+ equidistant neighbor, giving correlations of ~0.997-0.9999 rather than an
+ exact match. The numerical difference is tiny (the alternative neighbor is
+ the same distance away) and not worth chasing freesurfer's tie-breaking
+ bookkeeping to eliminate.
"""
+ src_sphere = np.asarray(src_sphere)
+ trg_sphere = np.asarray(trg_sphere)
+ n_src = len(src_sphere)
+ n_trg = len(trg_sphere)
+
+ # Forward: each target vertex -> its nearest source vertex.
+ _, fwd_src = KDTree(src_sphere).query(trg_sphere, k=1)
+ # Reverse: each source vertex -> its nearest target vertex.
+ _, rev_trg = KDTree(trg_sphere).query(src_sphere, k=1)
+
+ # Orphan source vertices are those not chosen by any target's forward
+ # mapping; fold them into their nearest target so their data is preserved.
+ used = np.zeros(n_src, dtype=bool)
+ used[fwd_src] = True
+ orphan = np.flatnonzero(~used)
+
+ rows = np.concatenate([np.arange(n_trg), rev_trg[orphan]])
+ cols = np.concatenate([fwd_src, orphan])
+ matrix = coo_matrix((np.ones(len(rows)), (rows, cols)),
+ shape=(n_trg, n_src)).tocsr()
+ # A target/source pair can be added by both passes; collapse duplicates.
+ matrix.sum_duplicates()
+
+ # Renormalize each row so the target value is the mean of its sources.
+ row_sums = np.asarray(matrix.sum(axis=1)).ravel()
+ row_sums[row_sums == 0] = 1.0
+ matrix = (diags(1.0 / row_sums) @ matrix).tocsr()
+ return matrix
- source_verts, _, _ = get_surf(source_subj, hemi, surface_type,
- freesurfer_subject_dir=subjects_dir)
-
- transformed_coords = mri_surf2surf(source_verts.T,
- source_subj, target_subj, hemi,
- subjects_dir=subjects_dir)
-
- kdt = KDTree(source_verts)
- print("Getting nearest neighbors")
- distances, indices = kdt.query(transformed_coords.T, k=n_neighbors)
- print("Done")
-
- rng = (np.random.RandomState(random_state)
- if isinstance(random_state, int) else random_state)
- test_images = rng.randn(n_test_images, len(source_verts))
- transformed_test_images = mri_surf2surf(test_images, source_subj,
- target_subj, hemi,
- subjects_dir=subjects_dir)
-
- # Solve linear problems to get coefficients
- all_coefs = []
- residuals = []
- print("Computing coefficients")
- i = 0
- for target_activation, source_inds in zip(
- transformed_test_images.T, indices):
- i += 1
- print("{i}".format(i=i), end="\r")
- source_values = test_images[:, source_inds]
- r = lstsq(source_values, target_activation,
- overwrite_a=True, overwrite_b=True)
- all_coefs.append(r[0])
- residuals.append(r[1])
- print("Done")
-
- all_coefs = np.array(all_coefs)
-
- if coef_threshold is None: # we know now that coefs are doing averages
- coef_threshold = (1 / 10. / n_neighbors )
- all_coefs[np.abs(all_coefs) < coef_threshold] = 0
- if renormalize:
- all_coefs /= np.abs(all_coefs).sum(axis=1)[:, np.newaxis] + 1e-10
-
- # there seem to be like 7 vertices that don't constitute an average over
- # 20 vertices or less, but all the others are such an average.
-
- # Let's make a matrix that does the transform:
- col_indices = indices.ravel()
- row_indices = (np.arange(indices.shape[0])[:, np.newaxis] *
- np.ones(indices.shape[1], dtype='int')).ravel()
- data = all_coefs.ravel()
- shape = (transformed_coords.shape[1], source_verts.shape[0])
-
- matrix = coo_matrix((data, (row_indices, col_indices)), shape=shape)
- return matrix
+# Legacy keyword arguments accepted by the old (regression-based)
+# implementation of ``get_mri_surf2surf_matrix``; kept only so existing
+# callers do not break.
+_SURF2SURF_LEGACY_KWARGS = frozenset(
+ ('n_neighbors', 'random_state', 'n_test_images', 'coef_threshold',
+ 'renormalize'))
+
+
+def get_mri_surf2surf_matrix(source_subj, hemi, surface_type=None,
+ target_subj='fsaverage', subjects_dir=None,
+ **kwargs):
+ """Create a sparse matrix implementing freesurfer's ``mri_surf2surf``.
+
+ A surface-to-surface resampling is a linear, highly localized transform
+ between two subjects' vertex spaces. Freesurfer defines this mapping
+ entirely from the spherical registration (``?h.sphere.reg``): for each
+ target vertex it takes the value of its nearest source vertex, and it
+ additionally averages in any source vertex that no target selected, so
+ that no source data is lost (the ``nnfr`` -- nearest-neighbor,
+ forward-and-reverse -- method).
+
+ This implementation builds that matrix directly from the registered-sphere
+ geometry with two nearest-neighbor queries (see
+ :func:`_surf2surf_nnfr_matrix`). It is exact for ordinary subjects and
+ full-resolution ``fsaverage``, requires no calls to the ``mri_surf2surf``
+ binary, and is deterministic. (The previous implementation reverse
+ -engineered the matrix by probing ``mri_surf2surf`` with random test images
+ and solving a per-vertex least-squares problem.)
+
+ Parameters
+ ----------
+ source_subj : str
+ Freesurfer name of the source subject.
+ hemi : str in ("lh", "rh")
+ Hemisphere.
+ surface_type : str, optional
+ Ignored. Retained for backwards compatibility with the previous
+ signature. The surf2surf correspondence depends only on the spherical
+ registration and is therefore identical for every surface
+ (``white``/``pial``/``inflated``/...) of a given subject.
+ target_subj : str, default "fsaverage"
+ Freesurfer name of the target subject.
+ subjects_dir : str, default os.environ["SUBJECTS_DIR"]
+ The freesurfer subjects directory.
+
+ Returns
+ -------
+ matrix : scipy.sparse.csr_matrix, shape (n_target_verts, n_source_verts)
+ Apply with ``target_data = matrix.dot(source_data)``.
+
+ Notes
+ -----
+ See :func:`_surf2surf_nnfr_matrix` for the algorithm and for the one known
+ caveat (exact tie-breaking on the regular icosahedral targets such as
+ ``fsaverage6``).
+ """
+ legacy = _SURF2SURF_LEGACY_KWARGS.intersection(kwargs)
+ if legacy:
+ warnings.warn(
+ "get_mri_surf2surf_matrix no longer uses the regression-based "
+ "parameters {}; they are ignored. The matrix is now built "
+ "directly from the spherical registration.".format(sorted(legacy)),
+ DeprecationWarning, stacklevel=2)
+ unexpected = set(kwargs) - _SURF2SURF_LEGACY_KWARGS
+ if unexpected:
+ raise TypeError(
+ "get_mri_surf2surf_matrix got unexpected keyword argument(s) "
+ "{}".format(sorted(unexpected)))
+
+ src_sphere = _read_sphere_reg(source_subj, hemi, subjects_dir)
+ trg_sphere = _read_sphere_reg(target_subj, hemi, subjects_dir)
+ return _surf2surf_nnfr_matrix(src_sphere, trg_sphere)
def get_curv(fs_subject, hemi, type='wm', freesurfer_subject_dir=None):
diff --git a/cortex/rois.py b/cortex/rois.py
index d5d9bbace..e9553c90b 100644
--- a/cortex/rois.py
+++ b/cortex/rois.py
@@ -117,8 +117,8 @@ def to_svg(self, open_inkscape=False, filename=None):
# Find polys
allbpolys = np.unique(surf.connected[inbound+exbound].indices)
selbpolys = surf.polys[allbpolys]
- inpolys = np.in1d(selbpolys, inbound).reshape(selbpolys.shape)
- expolys = np.in1d(selbpolys, exbound).reshape(selbpolys.shape)
+ inpolys = np.isin(selbpolys, inbound)
+ expolys = np.isin(selbpolys, exbound)
badpolys = np.logical_or(inpolys.all(1), expolys.all(1))
boundpolys = np.logical_and(np.logical_or(inpolys, expolys).all(1), ~badpolys)
diff --git a/cortex/segment.py b/cortex/segment.py
index f274b6ab6..3d5c24808 100644
--- a/cortex/segment.py
+++ b/cortex/segment.py
@@ -384,7 +384,7 @@ def flatten_slim(
)
# Cull pts that are not in manifold
pi = np.arange(len(pts))
- pii = np.in1d(pi, polys.flatten())
+ pii = np.isin(pi, polys.flatten())
idx = np.nonzero(pii)[0]
pts_new = pts[idx]
# Match indices in polys to new index for pts
diff --git a/cortex/surfinfo.py b/cortex/surfinfo.py
index 798403335..f59525648 100644
--- a/cortex/surfinfo.py
+++ b/cortex/surfinfo.py
@@ -186,7 +186,7 @@ def make_surface_graph(tris):
mwallset = set.union(*(set(g[v]) for v in fog.nodes())) & set(allbounds)
#cutset = (set(g.nodes()) - mwallset) & set(allbounds)
- mwallbounds = [np.in1d(b, mwallset) for b in bounds]
+ mwallbounds = [np.isin(b, mwallset) for b in bounds]
changes = [np.nonzero(np.diff(b.astype(float))!=0)[0]+1 for b in mwallbounds]
#splitbounds = [np.split(b, c) for b,c in zip(bounds, changes)]
diff --git a/cortex/testing_utils.py b/cortex/testing_utils.py
index 24a3221e0..85c5f2e76 100644
--- a/cortex/testing_utils.py
+++ b/cortex/testing_utils.py
@@ -11,12 +11,19 @@ def inkscape_version():
if not has_installed('inkscape'):
return None
cmd = 'inkscape --version'
- output = sp.check_output(cmd.split(), stderr=sp.PIPE)
- # b'Inkscape 1.0 (4035a4f, 2020-05-01)\n'
- version = output.split()[1]
- if isinstance(version, bytes):
- version = version.decode('utf-8')
- return version
+ result = sp.run(cmd.split(), stdout=sp.PIPE, stderr=sp.PIPE, check=True)
+ # Combine stdout and stderr; some systems print diagnostic messages
+ # (e.g. "Setting _INKSCAPE_GC=disable …") before the version line.
+ combined = result.stdout + result.stderr
+ if isinstance(combined, bytes):
+ combined = combined.decode('utf-8')
+ # Find the line that starts with 'Inkscape' to get the real version,
+ # e.g. 'Inkscape 1.2.2 (b0a8486, 2022-12-01)'
+ for line in combined.splitlines():
+ if line.strip().startswith('Inkscape'):
+ version = line.split()[1]
+ return version
+ return None
INKSCAPE_VERSION = inkscape_version()
diff --git a/cortex/tests/test_dataset.py b/cortex/tests/test_dataset.py
index e45ee082d..d3224bd4f 100644
--- a/cortex/tests/test_dataset.py
+++ b/cortex/tests/test_dataset.py
@@ -312,23 +312,27 @@ def test_blend_curvature():
view = cortex.Vertex.empty(subj)
alpha = np.linspace(0, 1, view.data.size).reshape(view.data.shape)
- # test alpha with float
- view_rgb = view.blend_curvature(alpha)
- # test alpha with bool
- view_rgb = view.blend_curvature(alpha > 0.3)
+ # blend_curvature is deprecated; the warning should fire on every call.
+ with pytest.warns(DeprecationWarning, match="blend_curvature is deprecated"):
+ view_rgb = view.blend_curvature(alpha)
+ with pytest.warns(DeprecationWarning):
+ view_rgb = view.blend_curvature(alpha > 0.3)
# test that it returns a VertexRGB
assert isinstance(view_rgb, cortex.VertexRGB)
# test on Vertex2D
view_2d = cortex.Vertex2D(view_rgb.red.data, view_rgb.green.data, subj)
- view_rgb = view_2d.blend_curvature(alpha)
+ with pytest.warns(DeprecationWarning):
+ view_rgb = view_2d.blend_curvature(alpha)
# test on VertexRGB
- view_rgb_new = view_rgb.blend_curvature(alpha)
+ with pytest.warns(DeprecationWarning):
+ view_rgb_new = view_rgb.blend_curvature(alpha)
# test that it returns a different VertexRGB
assert not np.allclose(view_rgb.red.data, view_rgb_new.red.data)
# test that it returns a VertexRGB with same values when alpha is ones
- view_rgb_new = view_rgb.blend_curvature(np.ones_like(alpha))
+ with pytest.warns(DeprecationWarning):
+ view_rgb_new = view_rgb.blend_curvature(np.ones_like(alpha))
assert np.allclose(view_rgb.red.data, view_rgb_new.red.data)
diff --git a/cortex/tests/test_freesurfer.py b/cortex/tests/test_freesurfer.py
index e7668fdfb..2a16e3ade 100644
--- a/cortex/tests/test_freesurfer.py
+++ b/cortex/tests/test_freesurfer.py
@@ -1,6 +1,15 @@
+import os
+import shutil
+
import numpy as np
+import pytest
-from cortex.freesurfer import _remove_disconnected_polys
+import cortex.freesurfer as fs
+from cortex.freesurfer import (
+ _remove_disconnected_polys,
+ _surf2surf_nnfr_matrix,
+ get_mri_surf2surf_matrix,
+)
def test_remove_disconnected_polys_examples():
@@ -27,3 +36,145 @@ def test_remove_disconnected_polys_idempotence():
# make sure calling the function does not change anything
polys_2 = _remove_disconnected_polys(polys_1)
np.testing.assert_array_equal(polys_1, polys_2)
+
+
+# ---------------------------------------------------------------------------
+# surf2surf (nnfr) matrix construction
+# ---------------------------------------------------------------------------
+
+def _random_sphere(n, seed):
+ """n points on the unit sphere (so KDTree distances behave like the
+ real ?h.sphere.reg geometry)."""
+ rng = np.random.RandomState(seed)
+ pts = rng.randn(n, 3)
+ return pts / np.linalg.norm(pts, axis=1, keepdims=True)
+
+
+def test_surf2surf_identity_when_source_equals_target():
+ # If source and target spheres are identical, every target maps to the
+ # co-located source vertex and there are no orphans -> identity matrix.
+ sphere = _random_sphere(60, seed=0)
+ m = _surf2surf_nnfr_matrix(sphere, sphere)
+ assert m.shape == (60, 60)
+ np.testing.assert_allclose(m.toarray(), np.eye(60))
+
+
+def test_surf2surf_rows_sum_to_one():
+ src = _random_sphere(80, seed=1)
+ trg = _random_sphere(50, seed=2)
+ m = _surf2surf_nnfr_matrix(src, trg)
+ assert m.shape == (50, 80)
+ row_sums = np.asarray(m.sum(axis=1)).ravel()
+ np.testing.assert_allclose(row_sums, np.ones(50))
+
+
+def test_surf2surf_preserves_constants():
+ # A row-normalized averaging matrix maps a constant map to itself.
+ src = _random_sphere(80, seed=3)
+ trg = _random_sphere(50, seed=4)
+ m = _surf2surf_nnfr_matrix(src, trg)
+ const = np.full(80, 7.0)
+ np.testing.assert_allclose(m.dot(const), np.full(50, 7.0))
+
+
+def test_surf2surf_no_source_vertex_is_dropped():
+ # The reverse pass guarantees every source vertex contributes to at least
+ # one target vertex (this is the whole point of "forward and reverse").
+ src = _random_sphere(120, seed=5)
+ trg = _random_sphere(40, seed=6) # heavy downsampling
+ m = _surf2surf_nnfr_matrix(src, trg)
+ col_nnz = m.getnnz(axis=0)
+ assert (col_nnz > 0).all()
+
+
+def test_surf2surf_known_averaging_example():
+ # Four collinear source vertices, two target vertices placed so that each
+ # target's nearest source is an endpoint, leaving the two middle source
+ # vertices as orphans that get folded into their nearest target.
+ src = np.array([[0., 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0]])
+ trg = np.array([[0.4, 0, 0], [2.6, 0, 0]])
+ m = _surf2surf_nnfr_matrix(src, trg)
+
+ expected = np.array([[0.5, 0.5, 0.0, 0.0], # mean of src 0 and 1
+ [0.0, 0.0, 0.5, 0.5]]) # mean of src 2 and 3
+ np.testing.assert_allclose(m.toarray(), expected)
+
+ data = np.array([10., 20., 30., 40.])
+ np.testing.assert_allclose(m.dot(data), [15., 35.])
+
+
+def test_get_mri_surf2surf_matrix_ignores_legacy_kwargs(monkeypatch):
+ # Legacy regression-based kwargs are accepted but ignored, with a warning.
+ src = _random_sphere(20, seed=7)
+ trg = _random_sphere(15, seed=8)
+ monkeypatch.setattr(
+ fs, "_read_sphere_reg",
+ lambda subj, hemi, subjects_dir=None: src if subj == "A" else trg)
+
+ with pytest.warns(DeprecationWarning):
+ m = get_mri_surf2surf_matrix("A", "lh", target_subj="B",
+ n_neighbors=20, n_test_images=40)
+ assert m.shape == (15, 20)
+
+
+def test_get_mri_surf2surf_matrix_rejects_unknown_kwargs():
+ with pytest.raises(TypeError):
+ get_mri_surf2surf_matrix("A", "lh", target_subj="B", bogus_kwarg=1)
+
+
+def _have_template(subjects_dir, name, hemi="lh"):
+ return bool(subjects_dir) and os.path.exists(
+ os.path.join(subjects_dir, name, "surf", hemi + ".sphere.reg"))
+
+
+@pytest.mark.skipif(shutil.which("mri_surf2surf") is None,
+ reason="freesurfer mri_surf2surf not available")
+def test_surf2surf_matches_freesurfer_identity():
+ """fsaverage -> fsaverage must be the identity, matching mri_surf2surf.
+
+ Uses only the standard fsaverage template (no individual subjects), so it
+ is reproducible anywhere freesurfer is installed and skipped otherwise.
+ """
+ subjects_dir = os.environ.get("SUBJECTS_DIR")
+ hemi = "lh"
+ if not _have_template(subjects_dir, "fsaverage", hemi):
+ pytest.skip("fsaverage template with sphere.reg not found in SUBJECTS_DIR")
+
+ m = get_mri_surf2surf_matrix("fsaverage", hemi, target_subj="fsaverage",
+ subjects_dir=subjects_dir)
+ rng = np.random.RandomState(0)
+ data = rng.randn(4, m.shape[1]).astype(np.float32)
+ reference = fs.mri_surf2surf(data, "fsaverage", "fsaverage", hemi,
+ subjects_dir=subjects_dir)
+ got = np.stack([m.dot(data[i]) for i in range(data.shape[0])])
+
+ np.testing.assert_allclose(got, data, atol=1e-4) # identity
+ np.testing.assert_allclose(got, reference, atol=1e-4) # matches freesurfer
+
+
+@pytest.mark.skipif(shutil.which("mri_surf2surf") is None,
+ reason="freesurfer mri_surf2surf not available")
+def test_surf2surf_matches_freesurfer_downsample():
+ """fsaverage -> fsaverage6 (icosahedral downsample) against mri_surf2surf.
+
+ Uses only standard fsaverage templates. This is the documented inexact
+ case: freesurfer's tie-breaking on exactly-equidistant vertices of the
+ regular mesh differs, so we only require a high correlation rather than a
+ bit-exact match (see _surf2surf_nnfr_matrix notes).
+ """
+ subjects_dir = os.environ.get("SUBJECTS_DIR")
+ hemi = "lh"
+ if not (_have_template(subjects_dir, "fsaverage", hemi)
+ and _have_template(subjects_dir, "fsaverage6", hemi)):
+ pytest.skip("fsaverage/fsaverage6 templates not found in SUBJECTS_DIR")
+
+ m = get_mri_surf2surf_matrix("fsaverage", hemi, target_subj="fsaverage6",
+ subjects_dir=subjects_dir)
+ rng = np.random.RandomState(0)
+ data = rng.randn(4, m.shape[1]).astype(np.float32)
+ reference = fs.mri_surf2surf(data, "fsaverage", "fsaverage6", hemi,
+ subjects_dir=subjects_dir)
+ got = np.stack([m.dot(data[i]) for i in range(data.shape[0])])
+
+ corr = np.corrcoef(reference.ravel(), got.ravel())[0, 1]
+ assert corr > 0.99
diff --git a/cortex/tests/test_quickflat.py b/cortex/tests/test_quickflat.py
index 11782f1cc..e0ec115bf 100644
--- a/cortex/tests/test_quickflat.py
+++ b/cortex/tests/test_quickflat.py
@@ -3,7 +3,9 @@
import tempfile
import pytest
+from cortex import dataset
from cortex.testing_utils import has_installed
+from cortex.webgl.data import Package
no_inkscape = not has_installed('inkscape')
@@ -40,3 +42,80 @@ def test_make_flatmap_image_nanmean(type_, nanmean):
vol, nanmean=nanmean)
# assert that the nanmean only returns NaNs and 1s
assert np.nanmin(img) == 1
+
+
+def test_quickshow_webgl_alpha_equivalence():
+ """quickshow (matplotlib) and WebGL must render the same VertexRGB+α identically.
+
+ Issue #631: the WebGL shader uses a premultiplied "over" composite, while
+ matplotlib's imshow layering uses straight alpha. The fix premultiplies α
+ into RGB at the WebGL serialization step only, so both paths converge on
+ the same composite formula out = α·rgb + (1-α)·bg for any background.
+ This test asserts that equivalence at the per-vertex level for an
+ arbitrary curvature gray.
+ """
+ subj = "S1"
+ nverts = cortex.db.get_surf(subj, "fiducial", merge=True)[0].shape[0]
+ rng = np.random.default_rng(631)
+ r = rng.uniform(0, 1, nverts).astype(np.float32)
+ g = rng.uniform(0, 1, nverts).astype(np.float32)
+ b = rng.uniform(0, 1, nverts).astype(np.float32)
+ alpha = rng.uniform(0, 1, nverts).astype(np.float32)
+
+ vrgb = cortex.VertexRGB(
+ r, g, b, subj,
+ alpha=cortex.Vertex(alpha, subj, vmin=0, vmax=1),
+ )
+
+ raw = vrgb.vertices # what quickshow/matplotlib will composite (non-premult)
+ pkg = Package(dataset.Dataset(view=vrgb))
+ packaged = pkg.images[vrgb.name][0] # what the shader will composite (premult)
+
+ # Sanity: alpha is shared between the two paths.
+ assert np.array_equal(raw[..., 3], packaged[..., 3])
+
+ # Composite both against an arbitrary curvature gray. matplotlib's
+ # imshow with two layered images uses straight alpha; the GLSL shader at
+ # shaderlib.js line 851 uses gl_FragColor = vColor + (1-α)·bg.
+ a_norm = raw[..., 3:4].astype(np.float32) / 255.0
+ rgb_raw = raw[..., :3].astype(np.float32) / 255.0
+ rgb_pkg = packaged[..., :3].astype(np.float32) / 255.0
+ for curv in (0.0, 0.25, 0.5, 0.75, 1.0):
+ bg = np.full_like(rgb_raw, curv)
+ matplotlib_out = a_norm * rgb_raw + (1.0 - a_norm) * bg
+ webgl_out = rgb_pkg + (1.0 - a_norm) * bg
+ # 1 LSB of uint8 rounding on each side -> 2/255 worst case.
+ np.testing.assert_allclose(matplotlib_out, webgl_out, atol=2.0 / 255.0)
+
+
+def test_make_flatmap_image_vertexrgb_alpha_unchanged():
+ """The matplotlib path must keep using NON-premultiplied RGBA bytes.
+
+ Premultiplying inside .vertices would silently double-attenuate the
+ quickshow output. Pin that .vertices stays straight-alpha by checking
+ a uniform bright-red, half-transparent VertexRGB survives
+ make_flatmap_image without losing red intensity.
+ """
+ subj = "S1"
+ nverts = cortex.db.get_surf(subj, "fiducial", merge=True)[0].shape[0]
+ # Uniform bright red, half transparent everywhere. Pass explicit Vertex
+ # objects with vmin/vmax to avoid auto-range degeneracy on the flat
+ # green/blue channels.
+ r = cortex.Vertex(np.ones(nverts, dtype=np.float32), subj, vmin=0, vmax=1)
+ g = cortex.Vertex(np.zeros(nverts, dtype=np.float32), subj, vmin=0, vmax=1)
+ b = cortex.Vertex(np.zeros(nverts, dtype=np.float32), subj, vmin=0, vmax=1)
+ alpha = cortex.Vertex(np.full(nverts, 0.5, dtype=np.float32), subj,
+ vmin=0, vmax=1)
+ vrgb = cortex.VertexRGB(r, g, b, subj, alpha=alpha)
+ img, _ = cortex.quickflat.utils.make_flatmap_image(vrgb)
+ # img is the rasterized RGBA flatmap. The data layer's red channel (where
+ # mask is filled and pixmap is non-degenerate) must be ~255, not ~127 --
+ # if we ever start premultiplying inside .vertices, this drops to ~127.
+ rgba_in_mask = img[img[..., 3] > 0]
+ assert rgba_in_mask.size > 0
+ # Filled pixels should have red close to 255 (bright red, with alpha=128).
+ bright_red_pixels = rgba_in_mask[rgba_in_mask[..., 0] > 200]
+ assert bright_red_pixels.size > 0, (
+ "VertexRGB.vertices appears to be premultiplied -- the matplotlib "
+ "path will double-attenuate. The fix should live in webgl/data.py."
+ )
diff --git a/cortex/tests/test_testing_utils.py b/cortex/tests/test_testing_utils.py
new file mode 100644
index 000000000..3abe4be00
--- /dev/null
+++ b/cortex/tests/test_testing_utils.py
@@ -0,0 +1,88 @@
+"""Tests for cortex/testing_utils.py"""
+import subprocess
+from unittest import mock
+
+import pytest
+
+from cortex.testing_utils import inkscape_version
+
+
+def _make_run_result(stdout=b'', stderr=b'', returncode=0):
+ """Helper to build a mock subprocess.CompletedProcess."""
+ result = mock.MagicMock()
+ result.stdout = stdout
+ result.stderr = stderr
+ result.returncode = returncode
+ return result
+
+
+@mock.patch('cortex.testing_utils.has_installed', return_value=False)
+def test_inkscape_version_not_installed(mock_has):
+ """Returns None when inkscape is not on PATH."""
+ assert inkscape_version() is None
+
+
+@mock.patch('cortex.testing_utils.has_installed', return_value=True)
+@mock.patch('cortex.testing_utils.sp.run')
+def test_inkscape_version_clean_stdout(mock_run, mock_has):
+ """Parses version correctly from clean stdout output."""
+ mock_run.return_value = _make_run_result(
+ stdout=b'Inkscape 1.0 (4035a4f, 2020-05-01)\n'
+ )
+ assert inkscape_version() == '1.0'
+
+
+@mock.patch('cortex.testing_utils.has_installed', return_value=True)
+@mock.patch('cortex.testing_utils.sp.run')
+def test_inkscape_version_newer(mock_run, mock_has):
+ """Parses a newer version number correctly."""
+ mock_run.return_value = _make_run_result(
+ stdout=b'Inkscape 1.2.2 (b0a8486, 2022-12-01)\n'
+ )
+ assert inkscape_version() == '1.2.2'
+
+
+@mock.patch('cortex.testing_utils.has_installed', return_value=True)
+@mock.patch('cortex.testing_utils.sp.run')
+def test_inkscape_version_with_diagnostic_noise(mock_run, mock_has):
+ """Returns correct version even when a diagnostic message precedes it.
+
+ This is the regression test for the bug where systems that print
+ "Setting _INKSCAPE_GC=disable as a workaround for broken libgc"
+ caused INKSCAPE_VERSION to be set to '_INKSCAPE_GC=disable'.
+ """
+ mock_run.return_value = _make_run_result(
+ stdout=(
+ b'Setting _INKSCAPE_GC=disable as a workaround for broken libgc\n'
+ b'Inkscape 1.2.2 (b0a8486, 2022-12-01)\n'
+ )
+ )
+ assert inkscape_version() == '1.2.2'
+
+
+@mock.patch('cortex.testing_utils.has_installed', return_value=True)
+@mock.patch('cortex.testing_utils.sp.run')
+def test_inkscape_version_noise_in_stderr(mock_run, mock_has):
+ """Returns correct version when the diagnostic message is on stderr."""
+ mock_run.return_value = _make_run_result(
+ stderr=b'Setting _INKSCAPE_GC=disable as a workaround for broken libgc\n',
+ stdout=b'Inkscape 1.2.2 (b0a8486, 2022-12-01)\n',
+ )
+ assert inkscape_version() == '1.2.2'
+
+
+@mock.patch('cortex.testing_utils.has_installed', return_value=True)
+@mock.patch('cortex.testing_utils.sp.run')
+def test_inkscape_version_no_inkscape_line(mock_run, mock_has):
+ """Returns None when no 'Inkscape …' line is found in output."""
+ mock_run.return_value = _make_run_result(stdout=b'some unexpected output\n')
+ assert inkscape_version() is None
+
+
+@mock.patch('cortex.testing_utils.has_installed', return_value=True)
+@mock.patch('cortex.testing_utils.sp.run')
+def test_inkscape_version_subprocess_error(mock_run, mock_has):
+ """Propagates CalledProcessError when inkscape exits non-zero."""
+ mock_run.side_effect = subprocess.CalledProcessError(1, 'inkscape')
+ with pytest.raises(subprocess.CalledProcessError):
+ inkscape_version()
diff --git a/cortex/tests/test_utils.py b/cortex/tests/test_utils.py
index 65c1c813c..22ee0c937 100644
--- a/cortex/tests/test_utils.py
+++ b/cortex/tests/test_utils.py
@@ -1,14 +1,76 @@
+import shutil
+import tarfile
+from unittest import mock
+
+import pytest
+
import cortex
-def test_download_subject():
- # Test that newly downloaded subjects are added to the current database.
- # remove fsaverage from the list of available subjects if present.
- if "fsaverage" in cortex.db.subjects:
- cortex.db._subjects.pop("fsaverage")
+@pytest.fixture
+def fake_fsaverage_tarball(tmp_path):
+ """Build a minimal fsaverage.tar.gz that pycortex will recognize as a subject."""
+ subj_src = tmp_path / "src" / "fsaverage"
+ subj_src.mkdir(parents=True)
+ (subj_src / "marker").write_text("ok")
+ tarball = tmp_path / "fsaverage.tar.gz"
+ with tarfile.open(tarball, "w:gz") as tar:
+ tar.add(subj_src, arcname="fsaverage")
+ return tarball
+
+
+@pytest.fixture
+def isolated_filestore(tmp_path, monkeypatch):
+ """Redirect cortex.db.filestore to an empty tmp dir and reset the subject cache."""
+ store = tmp_path / "store"
+ store.mkdir()
+ monkeypatch.setattr(cortex.db, "filestore", str(store))
+ original_subjects = cortex.db._subjects
+ cortex.db._subjects = None
+ yield store
+ cortex.db._subjects = original_subjects
+
+
+def test_download_subject(isolated_filestore, fake_fsaverage_tarball, monkeypatch):
+ # Newly downloaded subjects are added to the current database.
+ def fake_retrieve(url, dest):
+ shutil.copy(fake_fsaverage_tarball, dest)
+ return dest, None
+
+ mock_retrieve = mock.Mock(side_effect=fake_retrieve)
+ monkeypatch.setattr(cortex.utils.urllib.request, "urlretrieve", mock_retrieve)
assert "fsaverage" not in cortex.db.subjects
cortex.utils.download_subject(subject_id='fsaverage')
assert "fsaverage" in cortex.db.subjects
- # test that downloading it again works
+ assert mock_retrieve.call_count == 1
+
+
+def test_download_subject_skips_when_present(isolated_filestore, monkeypatch):
+ # If the subject is already in the database and download_again is False,
+ # download_subject warns and returns without touching the network.
+ cortex.db._subjects = {"fsaverage": mock.MagicMock()}
+
+ mock_retrieve = mock.Mock()
+ monkeypatch.setattr(cortex.utils.urllib.request, "urlretrieve", mock_retrieve)
+
+ with pytest.warns(UserWarning, match="already present"):
+ cortex.utils.download_subject(subject_id='fsaverage')
+ mock_retrieve.assert_not_called()
+
+
+def test_download_subject_download_again(
+ isolated_filestore, fake_fsaverage_tarball, monkeypatch
+):
+ # download_again=True re-downloads even when the subject is already present.
+ def fake_retrieve(url, dest):
+ shutil.copy(fake_fsaverage_tarball, dest)
+ return dest, None
+
+ mock_retrieve = mock.Mock(side_effect=fake_retrieve)
+ monkeypatch.setattr(cortex.utils.urllib.request, "urlretrieve", mock_retrieve)
+
+ cortex.utils.download_subject(subject_id='fsaverage')
+ assert "fsaverage" in cortex.db.subjects
cortex.utils.download_subject(subject_id='fsaverage', download_again=True)
+ assert mock_retrieve.call_count == 2
diff --git a/cortex/tests/test_webgl_data.py b/cortex/tests/test_webgl_data.py
new file mode 100644
index 000000000..35302e1af
--- /dev/null
+++ b/cortex/tests/test_webgl_data.py
@@ -0,0 +1,194 @@
+"""Tests for the WebGL serialization layer in cortex.webgl.data."""
+
+import numpy as np
+import pytest
+
+import cortex
+from cortex import dataset
+from cortex.webgl.data import Package
+
+
+subj, xfmname, nverts, volshape = "S1", "fullhead", 304380, (31, 100, 100)
+
+
+def _packaged_rgba(brain):
+ """Run a dataview through Package and recover the pre-mosaic uint8 RGBA bytes."""
+ pkg = Package(dataset.Dataset(view=brain))
+ images = pkg.images[brain.name]
+ # Vertex* paths store the raw uint8 array directly; Volume* paths mosaic+PNG.
+ if isinstance(brain, dataset.VertexRGB):
+ return images[0]
+ raise NotImplementedError("Use _packaged_rgba_volume for VolumeRGB")
+
+
+def _expected_premultiplied(raw_uint8):
+ """Compute alpha-premultiplied RGB bytes the same way Package should."""
+ a = raw_uint8[..., 3:4].astype(np.float32) / 255.0
+ out = raw_uint8.copy()
+ out[..., :3] = np.round(raw_uint8[..., :3].astype(np.float32) * a).astype(np.uint8)
+ return out
+
+
+def test_vertexrgb_alpha_is_premultiplied_in_package():
+ """WebGL Package should ship alpha-premultiplied RGB bytes (issue #631).
+
+ The shader formula is gl_FragColor = vColor + (1-α)·bg, which only yields
+ correct "over" compositing when vColor is premultiplied. The fix lives in
+ Package.__init__; this test pins the contract.
+ """
+ rng = np.random.default_rng(0)
+ r = rng.uniform(0, 1, nverts).astype(np.float32)
+ g = rng.uniform(0, 1, nverts).astype(np.float32)
+ b = rng.uniform(0, 1, nverts).astype(np.float32)
+ alpha = rng.uniform(0, 1, nverts).astype(np.float32)
+
+ vrgb = cortex.VertexRGB(
+ r,
+ g,
+ b,
+ subj,
+ alpha=cortex.Vertex(alpha, subj, vmin=0, vmax=1),
+ )
+
+ # Sanity: the .vertices property itself stays NON-premultiplied so the
+ # quickshow (matplotlib) path keeps working.
+ raw = vrgb.vertices
+ assert raw.dtype == np.uint8
+ assert raw.shape == (1, nverts, 4)
+ # Raw alpha must round-trip the input alpha to within 1 LSB.
+ expected_a = (np.clip(alpha, 0, 1) * 255).astype(np.uint8)
+ assert np.allclose(raw[0, :, 3].astype(int), expected_a.astype(int), atol=1)
+ # Raw RGB bytes must NOT already be premultiplied -- if they were, the
+ # fix is in the wrong layer and quickshow would double-attenuate.
+ nontrivial = expected_a < 200 # avoid α≈1 where premult ≈ raw
+ naive_premult = np.round(
+ raw[0, :, 0].astype(np.float32) * raw[0, :, 3].astype(np.float32) / 255.0
+ ).astype(np.uint8)
+ assert (
+ np.mean(
+ np.abs(
+ raw[0, nontrivial, 0].astype(int)
+ - naive_premult[nontrivial].astype(int)
+ )
+ )
+ > 5
+ ), "vertices property looks already-premultiplied; quickshow would break"
+
+ # Now check what Package serializes for the WebGL viewer.
+ packaged = _packaged_rgba(vrgb)
+ expected = _expected_premultiplied(raw)
+ assert packaged.shape == raw.shape
+ assert packaged.dtype == np.uint8
+ # Alpha channel must be passed through unchanged (shader needs it for 1-α).
+ assert np.array_equal(packaged[..., 3], expected[..., 3])
+ # RGB channels must be alpha-premultiplied.
+ assert np.array_equal(packaged[..., :3], expected[..., :3])
+
+ # And the un-packaged property must NOT have been mutated by Package
+ # (Package should defensive-copy).
+ raw_after = vrgb.vertices
+ assert np.array_equal(raw_after, raw)
+
+
+def test_vertexrgb_alpha_one_is_passthrough():
+ """When α=1 everywhere, premultiplication is a no-op (bug was invisible at α=1)."""
+ rng = np.random.default_rng(1)
+ r = rng.uniform(0, 1, nverts).astype(np.float32)
+ g = rng.uniform(0, 1, nverts).astype(np.float32)
+ b = rng.uniform(0, 1, nverts).astype(np.float32)
+
+ vrgb = cortex.VertexRGB(r, g, b, subj) # default alpha = 1
+ raw = vrgb.vertices
+ packaged = _packaged_rgba(vrgb)
+ assert np.array_equal(packaged[..., 3], raw[..., 3]) # alpha all 255
+ # RGB unchanged because α=255 → premultiply by 1.
+ assert np.array_equal(packaged[..., :3], raw[..., :3])
+
+
+def test_vertexrgb_alpha_zero_zeros_rgb():
+ """α=0 must drive packaged RGB to 0 -- the shader then renders pure curvature."""
+ rng = np.random.default_rng(2)
+ r = rng.uniform(0, 1, nverts).astype(np.float32)
+ g = rng.uniform(0, 1, nverts).astype(np.float32)
+ b = rng.uniform(0, 1, nverts).astype(np.float32)
+ alpha = np.zeros(nverts, dtype=np.float32)
+
+ vrgb = cortex.VertexRGB(
+ r,
+ g,
+ b,
+ subj,
+ alpha=cortex.Vertex(alpha, subj, vmin=0, vmax=1),
+ )
+ packaged = _packaged_rgba(vrgb)
+ assert np.array_equal(packaged[..., 3], np.zeros_like(packaged[..., 3]))
+ assert np.array_equal(packaged[..., :3], np.zeros_like(packaged[..., :3]))
+
+
+def test_volumergb_alpha_is_NOT_premultiplied_in_package():
+ """VolumeRGB must ship straight-alpha bytes -- Three.js premultiplies on upload.
+
+ Three.js sets ``tex.premultiplyAlpha = true`` for raw VolumeRGB textures
+ (cortex/webgl/resources/js/dataset.js:335-338), which makes WebGL apply
+ UNPACK_PREMULTIPLY_ALPHA_WEBGL on texture upload. So the shader sees
+ premultiplied RGB by the time vColor is sampled, but ONLY because the
+ texture-upload pipeline does it for us. If Package also premultiplied
+ here, partial-alpha VolumeRGB would render double-attenuated (too dark).
+ """
+ rng = np.random.default_rng(3)
+ shape = volshape
+ r = rng.uniform(0, 1, shape).astype(np.float32)
+ g = rng.uniform(0, 1, shape).astype(np.float32)
+ b = rng.uniform(0, 1, shape).astype(np.float32)
+ alpha = rng.uniform(0, 1, shape).astype(np.float32)
+
+ vrgb = cortex.VolumeRGB(
+ r,
+ g,
+ b,
+ subj,
+ xfmname,
+ alpha=cortex.Volume(alpha, subj, xfmname, vmin=0, vmax=1),
+ )
+
+ raw = vrgb.volume
+ assert raw.dtype == np.uint8
+
+ # Spy on the Package internals: wrap volume.mosaic to capture the array
+ # Package actually ships before it gets PNG-encoded. mock.patch handles
+ # restoration even if the patched call raises before reaching the
+ # assertions below.
+ from unittest import mock
+
+ from cortex.webgl import data as webgl_data
+
+ captured = []
+ real_mosaic = webgl_data.volume.mosaic
+
+ def spy_mosaic(arr, show=False):
+ captured.append(arr.copy())
+ return real_mosaic(arr, show=show)
+
+ with mock.patch.object(webgl_data.volume, "mosaic", side_effect=spy_mosaic):
+ Package(dataset.Dataset(view=vrgb))
+
+ assert len(captured) == 1
+ packaged_frame = captured[0]
+
+ # Bytes shipped to the browser must equal the raw .volume bytes verbatim
+ # (NOT premultiplied) -- Three.js will premultiply once on texture upload.
+ assert packaged_frame.shape == raw[0].shape
+ assert np.array_equal(packaged_frame, raw[0])
+
+ # And specifically, RGB should NOT have been premultiplied by alpha.
+ naive_premult = _expected_premultiplied(raw[0])
+ nontrivial = (raw[0, ..., 3] < 200) & (raw[0, ..., 0] > 50)
+ assert (
+ np.mean(
+ np.abs(
+ packaged_frame[nontrivial][..., 0].astype(int)
+ - naive_premult[nontrivial][..., 0].astype(int)
+ )
+ )
+ > 5
+ ), "VolumeRGB Package output looks premultiplied; Three.js will then double-attenuate"
diff --git a/cortex/tests/test_webgl_headless.py b/cortex/tests/test_webgl_headless.py
index f998ba73c..908084124 100644
--- a/cortex/tests/test_webgl_headless.py
+++ b/cortex/tests/test_webgl_headless.py
@@ -83,7 +83,6 @@ def test_datatype_renders(dtype_name, tmp_path):
"""Each data type should render in the headless viewer without errors."""
vol = make_dataview(dtype_name)
with cortex.export.headless_viewer(vol, viewer_params={}) as handle:
- time.sleep(10)
outfile = str(tmp_path / "test.png")
handle.getImage(outfile, (512, 384))
_wait_for_file(outfile)
@@ -111,7 +110,6 @@ def _setup_viewer(self, tmp_path_factory):
cls = type(self)
cls.tmp_dir = tmp_path_factory.mktemp("angles")
with cortex.export.headless_viewer(vol, viewer_params={}) as handle:
- time.sleep(10)
cls.handle = handle
yield
@@ -149,7 +147,6 @@ def _setup_viewer(self, tmp_path_factory):
cls = type(self)
cls.tmp_dir = tmp_path_factory.mktemp("surfaces")
with cortex.export.headless_viewer(vol, viewer_params={}) as handle:
- time.sleep(10)
cls.handle = handle
yield
@@ -208,7 +205,6 @@ def test_capture_view_roundtrip():
"""Setting view parameters and capturing them back should match."""
vol = cortex.Volume(np.random.randn(*volshape), subj, xfmname)
with cortex.export.headless_viewer(vol, viewer_params={}) as handle:
- time.sleep(10)
target_params = {
"camera.azimuth": 90,
"camera.altitude": 90,
@@ -243,7 +239,6 @@ def test_overlay_visibility_changes_image(tmp_path):
with cortex.export.headless_viewer(
vol, viewer_params=dict(overlays_visible=["rois"])
) as handle:
- time.sleep(10)
handle._set_view(**view)
time.sleep(1)
handle.getImage(f1, (512, 384))
@@ -254,7 +249,6 @@ def test_overlay_visibility_changes_image(tmp_path):
with cortex.export.headless_viewer(
vol, viewer_params=dict(overlays_visible=[])
) as handle:
- time.sleep(10)
handle._set_view(**view)
time.sleep(1)
handle.getImage(f2, (512, 384))
@@ -266,7 +260,343 @@ def test_overlay_visibility_changes_image(tmp_path):
# ---------------------------------------------------------------------------
-# Group 7: addData dataset switching
+# Group 7: Vertex NaN-mask regression tests (#612, #626)
+# ---------------------------------------------------------------------------
+
+
+def _count_red_pixels(png_path):
+ """Count strongly red-dominant pixels (R - max(G, B) > 50)."""
+ from PIL import Image
+
+ rgb = np.array(Image.open(png_path))[..., :3].astype(int)
+ return int((rgb[..., 0] - np.maximum(rgb[..., 1], rgb[..., 2]) > 50).sum())
+
+
+def test_vertex_no_nan_renders_data(tmp_path):
+ """A NaN-free Vertex must render visibly, not fall through to transparent.
+
+ Regression test for #626: prior to the fix, the surface_vertex shader's
+ nanmask attribute defaulted to zeros when the Python data had no NaNs,
+ causing every vertex to be discarded and the brain to render with only
+ the grayscale curvature underlay.
+ """
+ np.random.seed(0)
+ # Constant high values + chromatic colormap so colored pixels are
+ # easily distinguishable from the grayscale curvature underlay.
+ data = np.full(nverts, 5.0)
+ vtx = cortex.Vertex(data, subj, vmin=0, vmax=1, cmap="Reds")
+
+ view = {
+ **default_view_params,
+ **angle_view_params["lateral_pivot"],
+ **unfold_view_params["inflated"],
+ }
+
+ with cortex.export.headless_viewer(vtx, viewer_params={}) as handle:
+ handle._set_view(**view)
+ time.sleep(1)
+ outfile = str(tmp_path / "vtx.png")
+ handle.getImage(outfile, (512, 384))
+ _wait_for_file(outfile)
+
+ n_red = _count_red_pixels(outfile)
+ assert n_red > 1000, (
+ f"Vertex data does not appear to be rendering "
+ f"(only {n_red} red-dominant pixels). "
+ "Surface may be falling through to grayscale curvature (#626)."
+ )
+
+
+def test_vertex_with_nan_renders_partial(tmp_path):
+ """A Vertex with some NaN values still renders the non-NaN portion (#612).
+
+ Sanity check that the per-vertex NaN mask path keeps working: half-NaN
+ data should render strictly fewer red pixels than fully-valid data, but
+ still meaningfully more than zero.
+ """
+ np.random.seed(0)
+
+ full = np.full(nverts, 5.0)
+ half_nan = full.copy()
+ half_nan[: nverts // 2] = np.nan
+
+ view = {
+ **default_view_params,
+ **angle_view_params["lateral_pivot"],
+ **unfold_view_params["inflated"],
+ }
+
+ def render(data, name):
+ vtx = cortex.Vertex(data, subj, vmin=0, vmax=1, cmap="Reds")
+ with cortex.export.headless_viewer(vtx, viewer_params={}) as handle:
+ handle._set_view(**view)
+ time.sleep(1)
+ outfile = str(tmp_path / f"{name}.png")
+ handle.getImage(outfile, (512, 384))
+ _wait_for_file(outfile)
+ return _count_red_pixels(outfile)
+
+ n_full = render(full, "full")
+ n_half = render(half_nan, "half_nan")
+
+ assert n_full > 1000, "Fully-valid Vertex should render visibly"
+ assert (
+ n_half > 100
+ ), "Half-NaN Vertex should still render the non-NaN half (#612 regression)"
+ assert n_half < n_full, (
+ f"Expected half-NaN render ({n_half} red px) to have fewer red "
+ f"pixels than fully-valid render ({n_full} red px)"
+ )
+
+
+# ---------------------------------------------------------------------------
+# Group 8: VertexRGB alpha attenuation regression test (#631)
+# ---------------------------------------------------------------------------
+
+
+def test_vertexrgb_alpha_zero_renders_curvature_only(tmp_path):
+ """VertexRGB with α=0 must render the curvature underlay, not bright color.
+
+ Regression test for #631: prior to the fix, the WebGL fragment shader's
+ premultiplied-alpha composite formula (gl_FragColor = vColor + (1-α)·bg)
+ consumed un-premultiplied RGB bytes from VertexRGB.vertices, so α=0 left
+ the foreground color fully opaque and clipped toward white instead of
+ falling through to the gray curvature.
+
+ With the fix, RGB is premultiplied at the WebGL serialization step
+ (cortex/webgl/data.py), so packaged vColor.rgb=0 when α=0, and the
+ shader produces pure curvature gray.
+ """
+ from PIL import Image
+
+ rng = np.random.default_rng(631)
+ # Bright, saturated colors -- if the bug returns these will leak through
+ # as red/green/blue pixels. With the fix and α=0, only neutral (curvature)
+ # gray pixels should remain in the brain region.
+ r = rng.uniform(0.7, 1.0, nverts).astype(np.float32)
+ g = rng.uniform(0.0, 0.3, nverts).astype(np.float32)
+ b = rng.uniform(0.0, 0.3, nverts).astype(np.float32)
+ alpha = np.zeros(nverts, dtype=np.float32)
+
+ vrgb = cortex.VertexRGB(
+ r,
+ g,
+ b,
+ subj,
+ alpha=cortex.Vertex(alpha, subj, vmin=0, vmax=1),
+ )
+
+ view = {
+ **default_view_params,
+ **angle_view_params["lateral_pivot"],
+ **unfold_view_params["inflated"],
+ }
+ with cortex.export.headless_viewer(vrgb, viewer_params={}) as handle:
+ handle._set_view(**view)
+ time.sleep(1)
+ outfile = str(tmp_path / "alpha_zero.png")
+ handle.getImage(outfile, (512, 384))
+ _wait_for_file(outfile)
+
+ rgb = np.array(Image.open(outfile))[..., :3].astype(int)
+ # Count strongly red-dominant pixels: with the bug, α=0 lets the
+ # bright reds through and we'd see thousands of them. With the fix,
+ # the brain renders curvature gray (R≈G≈B) and red-dominant pixels
+ # fall to near zero (a handful from anti-aliased ROI overlays).
+ n_red = int((rgb[..., 0] - np.maximum(rgb[..., 1], rgb[..., 2]) > 50).sum())
+ assert n_red < 500, (
+ f"VertexRGB with α=0 produced {n_red} red-dominant pixels; "
+ "expected near-zero. The shader composite is consuming "
+ "un-premultiplied RGB (issue #631)."
+ )
+
+
+def test_volumergb_alpha_half_renders_correct_blend(tmp_path):
+ """VolumeRGB with α=0.5 must blend halfway, not double-attenuate.
+
+ Companion regression to test_vertexrgb_alpha_zero_renders_curvature_only
+ (#631). VolumeRGB ships through the PNG texture path: Three.js sets
+ ``tex.premultiplyAlpha = true`` on upload, so the texture is premultiplied
+ once by WebGL itself. Package therefore must NOT premultiply on the
+ Python side -- if it does, the shader sees double-attenuated RGB.
+
+ α=0 won't catch that bug (0·anything = 0), so we use α=0.5 with bright
+ uniform red. With curvature contribution included, observed shader
+ output for the brain region is:
+ - correct (single premult by JS): median R ≈ 145-160
+ - bug (double premult: Py + JS): median R ≈ 90-105
+ Threshold at 125 sits in the middle of the gap and tolerates 20+ LSB
+ of boundary/interpolation noise on either side.
+ """
+ from PIL import Image
+
+ # Uniform saturated red over the whole volume, half transparent. Wrap in
+ # explicit Volume(vmin=0, vmax=1) so the .volume property doesn't
+ # auto-normalize a constant array to NaN.
+ r = cortex.Volume(
+ np.full(volshape, 1.0, dtype=np.float32), subj, xfmname, vmin=0, vmax=1
+ )
+ g = cortex.Volume(
+ np.full(volshape, 0.0, dtype=np.float32), subj, xfmname, vmin=0, vmax=1
+ )
+ b = cortex.Volume(
+ np.full(volshape, 0.0, dtype=np.float32), subj, xfmname, vmin=0, vmax=1
+ )
+ alpha = cortex.Volume(
+ np.full(volshape, 0.5, dtype=np.float32), subj, xfmname, vmin=0, vmax=1
+ )
+ vrgb = cortex.VolumeRGB(r, g, b, subj, xfmname, alpha=alpha)
+
+ view = {
+ **default_view_params,
+ **angle_view_params["lateral_pivot"],
+ **unfold_view_params["inflated"],
+ }
+ with cortex.export.headless_viewer(vrgb, viewer_params={}) as handle:
+ handle._set_view(**view)
+ time.sleep(1)
+ outfile = str(tmp_path / "volumergb_alpha_half.png")
+ handle.getImage(outfile, (512, 384))
+ _wait_for_file(outfile)
+
+ rgb = np.array(Image.open(outfile))[..., :3].astype(int)
+ # Brain-region pixels are red-dominant under both correct and buggy
+ # paths, but their R intensity differs. Pick the strongly red
+ # pixels (R clearly > G,B) and check median R.
+ red_mask = rgb[..., 0] - np.maximum(rgb[..., 1], rgb[..., 2]) > 30
+ assert red_mask.sum() > 1000, (
+ "Expected a large red-dominant region for half-transparent red "
+ f"VolumeRGB; got only {red_mask.sum()} pixels. Did the brain render?"
+ )
+ median_r = float(np.median(rgb[red_mask, 0]))
+ # Discriminator: correct path produces ~145-160, double-premult
+ # produces ~90-105. Threshold at 125 sits in the middle.
+ assert median_r > 125, (
+ f"VolumeRGB α=0.5 brain pixels have median R={median_r:.0f}; "
+ "expected ~150. R<125 indicates Package is double-premultiplying "
+ "VolumeRGB (issue #631 regression)."
+ )
+
+
+def test_vertex2d_alpha_half_renders_correct_blend(tmp_path):
+ """Vertex2D with α=0.5 must blend halfway, not over-attenuate the bg.
+
+ Companion regression to the issue #631 fix on the colormap-texture
+ path. The 2D dataview ships dim1 / dim2 as separate scalar maps and
+ the LUT lookup happens on the GPU via
+ ``texture2D(colormap, vec2(dim1_norm, dim2_norm))``. The shader's
+ composite (shaderlib.js:851) uses the premultiplied-over formula
+ ``vColor + (1-α)·bg``, so the colormap texture itself must be
+ premultiplied on upload (``tex.premultiplyAlpha = true`` in
+ mriview.js). Without that, alpha-bearing colormaps like
+ ``RdBu_r_alpha`` produce ``R + (1-α)·bg`` -- where the foreground
+ is added on top of a partially-attenuated curvature -- instead of
+ the correct ``α·R + (1-α)·bg``.
+
+ α=0 doesn't catch this bug because most alpha colormaps store
+ ``(0, 0, 0, 0)`` at the transparent end of the LUT (so neither the
+ buggy nor the correct shader produces foreground there). At α=0.5
+ the LUT stores its full RGB with α=127, and the difference between
+ buggy and correct composites is maximal in the brain region.
+
+ Empirical pixel stats for RdBu_r_alpha at data=+1, alpha=0.5,
+ inflated lateral_pivot view, default viewer params (S1):
+
+ - correct (premultiplied): red_dom median R ≈ 93 (25/50/75 = 80/93/110)
+ - buggy (un-premultiplied): red_dom median R ≈ 129 (25/50/75 = 105/129/149)
+
+ Threshold at 115 sits between the two distributions.
+ """
+ from PIL import Image
+
+ # data=+1 puts every vertex at the deep red end of RdBu_r_alpha,
+ # alpha=0.5 puts every vertex at mid-α (LUT row ~128).
+ data = np.full(nverts, 1.0, dtype=np.float32)
+ alpha = np.full(nverts, 0.5, dtype=np.float32)
+
+ vtx2d = cortex.Vertex2D(
+ data, alpha, subj,
+ cmap="RdBu_r_alpha",
+ vmin=-1, vmax=1,
+ vmin2=0, vmax2=1,
+ )
+
+ view = {
+ **default_view_params,
+ **angle_view_params["lateral_pivot"],
+ **unfold_view_params["inflated"],
+ }
+ # The cmap ![]()
elements decode asynchronously in Chromium. If the first
+ # render frame happens before the LUT image has decoded, three.js skips
+ # the texImage2D upload and the data layer renders against a 1×1 black
+ # texture (R==G==B everywhere -- looks like the curvature underlay).
+ # Retry _set_view + getImage until we observe a colored data layer
+ # (some pixels with R clearly > G or B, or vice-versa). We then run the
+ # premultiplication discriminator on that frame.
+ with cortex.export.headless_viewer(vtx2d, viewer_params={}) as handle:
+ # viewer.loaded already resolved by the context manager; a short
+ # extra pause covers the gap before the cmap ![]()
decodes. The
+ # retry loop below is the real guard for slow decodes.
+ time.sleep(2)
+ rgb = None
+ outfile = None
+ for attempt in range(6):
+ handle._set_view(**view)
+ time.sleep(3)
+ # Use a fresh filename each retry so we never read a partial PNG
+ # left over from a prior iteration (getImage writes async).
+ outfile = str(tmp_path / f"vertex2d_alpha_half_{attempt}.png")
+ handle.getImage(outfile, (512, 384))
+ _wait_for_file(outfile)
+ # Give the PNG writer a moment to finish flushing.
+ time.sleep(1)
+ try:
+ rgb = np.array(Image.open(outfile))[..., :3].astype(int)
+ except Exception:
+ continue
+ # "Colored" = at least some pixels deviate strongly from R==G==B.
+ # In a curvature-only (cmap-unbound) frame all brain pixels have
+ # R==G==B exactly; any non-zero count of channel-divergent pixels
+ # means the cmap texture is bound.
+ channel_spread = np.abs(rgb[..., 0] - rgb[..., 1]) + np.abs(
+ rgb[..., 1] - rgb[..., 2]
+ )
+ if (channel_spread > 5).sum() > 1000:
+ break
+ else:
+ pytest.skip(
+ "Cmap texture never bound in headless Chromium across 6 "
+ "render retries; can't discriminate fix vs bug."
+ )
+
+ # Both fix and bug produce a red-dominant brain region (the bug
+ # doesn't zero the foreground, just over-brightens it). The
+ # discriminator is the *median R intensity* of those red-dominant
+ # pixels: the buggy un-premultiplied path adds the full R on top of
+ # half the curvature, biasing R upward; the correct premultiplied
+ # path attenuates R by α before adding curvature.
+ #
+ # First, the brain must render as a red-dominant region (this is
+ # also satisfied by the bug, but if even this fails the cmap is
+ # unbound and we can't discriminate).
+ red_mask = rgb[..., 0] - np.maximum(rgb[..., 1], rgb[..., 2]) > 20
+ assert red_mask.sum() > 1000, (
+ f"Vertex2D α=0.5 deep-red rendered only {red_mask.sum()} "
+ "red-dominant pixels. Check that the cmap LUT bound and the "
+ "data layer rendered at all."
+ )
+ median_r = float(np.median(rgb[red_mask, 0]))
+ assert median_r < 115, (
+ f"Vertex2D α=0.5 brain pixels have median R={median_r:.0f}; "
+ "expected ≈93 (correct), saw ≥115 which is in the buggy range "
+ "(~129). The colormap texture is being sampled straight-alpha "
+ "while the shader applies a premultiplied composite -- check "
+ "mriview.js cmap texture premultiplyAlpha."
+ )
+
+
+# ---------------------------------------------------------------------------
+# Group 9: addData dataset switching
# ---------------------------------------------------------------------------
@@ -280,8 +610,212 @@ def test_addData_no_crash():
vol1 = cortex.Volume(np.random.randn(*volshape), subj, xfmname)
vol2 = cortex.Volume(np.random.randn(*volshape), subj, xfmname)
with cortex.export.headless_viewer(vol1, viewer_params={}) as handle:
- time.sleep(10)
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}"
+
+
+# ---------------------------------------------------------------------------
+# 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)
+
+ 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/utils.py b/cortex/utils.py
index 8cfaf201d..6ed7ab520 100644
--- a/cortex/utils.py
+++ b/cortex/utils.py
@@ -811,12 +811,12 @@ def get_roi_masks(subject, xfmname, roi_list=None, gm_sampler='cortical', split_
vert_in_scan = np.hstack([np.array((m>0).sum(1)).flatten() for m in mapper.masks])
vert_in_scan = vert_in_scan[roi_verts[roi]]
elif use_cortex_mask:
- vox_in_roi = np.in1d(vox_idx.flatten(), roi_verts[roi]).reshape(vox_idx.shape)
+ vox_in_roi = np.isin(vox_idx, roi_verts[roi])
roi_voxels[roi] = vox_in_roi & cortex_mask
# This is not accurate... because vox_idx only contains the indices of the *nearest*
# vertex to each voxel, it excludes many vertices. I can't think of a way to compute
# this accurately for non-mapper gm_samplers for now... ML 2017.07.14
- vert_in_scan = np.in1d(roi_verts[roi], vox_idx[cortex_mask])
+ vert_in_scan = np.isin(roi_verts[roi], vox_idx[cortex_mask])
# Compute ROI coverage
pct_coverage[roi] = vert_in_scan.mean() * 100
if use_mapper:
diff --git a/cortex/webgl/data.py b/cortex/webgl/data.py
index dc15eccb8..0dfed7d4c 100644
--- a/cortex/webgl/data.py
+++ b/cortex/webgl/data.py
@@ -7,6 +7,7 @@
images=(__braindata_name=["img1.png", "img2.png"]),
)
"""
+
import os
import json
from io import BytesIO
@@ -16,14 +17,15 @@
from .. import volume
-#TODO: How to package multiviews?
+# TODO: How to package multiviews?
class Package(object):
"""Package the data into a form usable by javascript"""
+
def __init__(self, data):
self.dataset = dataset.normalize(data)
self.uniques = list(data.uniques(collapse=True))
self.subjects = set()
-
+
self.brains = dict()
self.images = dict()
for brain in self.uniques:
@@ -36,19 +38,41 @@ def __init__(self, data):
encdata = brain.volume
if isinstance(brain, (dataset.VolumeRGB, dataset.VertexRGB)):
encdata = encdata.astype(np.uint8)
- self.brains[name]['raw'] = True
+ # The WebGL fragment shader (shaderlib.js) composites with a
+ # premultiplied-alpha "over" formula
+ # (gl_FragColor = vColor + (1-α)·bg). We only need to pre-
+ # multiply on the Python side for VertexRGB, where the bytes
+ # are uploaded as raw vertex attributes and Three.js does NOT
+ # premultiply (see dataset.js VertexData path). VolumeRGB ships
+ # through the PNG texture path (dataset.js:335-338, raw=true),
+ # where Three.js sets `tex.premultiplyAlpha = true` and the
+ # WebGL UNPACK_PREMULTIPLY_ALPHA_WEBGL hook premultiplies the
+ # texture once on upload -- premultiplying here would double-
+ # attenuate it. The .vertices/.volume properties stay
+ # non-premultiplied so the matplotlib (quickshow) path keeps
+ # using matplotlib's straight-alpha imshow compositor.
+ if isinstance(brain, dataset.VertexRGB):
+ # Note: encdata is already a fresh uint8 copy from the
+ # .astype(np.uint8) call above, so we can write into it
+ # in place. The assignment to a uint8 slice handles the
+ # float→uint8 cast for us.
+ a = encdata[..., 3:4].astype(np.float32) / 255.0
+ encdata[..., :3] = np.round(
+ encdata[..., :3].astype(np.float32) * a
+ )
+ self.brains[name]["raw"] = True
else:
encdata = encdata.astype(np.float32)
- self.brains[name]['raw'] = False
+ self.brains[name]["raw"] = False
- #VertexData requires reordering, only save normalized version for now
+ # VertexData requires reordering, only save normalized version for now
if isinstance(brain, (dataset.Vertex, dataset.VertexRGB)):
self.images[name] = [encdata]
else:
self.images[name] = [volume.mosaic(vol, show=False) for vol in encdata]
if len(set([shape for m, shape in self.images[name]])) != 1:
- raise ValueError('Internal error in mosaic')
- self.brains[name]['mosaic'] = self.images[name][0][1]
+ raise ValueError("Internal error in mosaic")
+ self.brains[name]["mosaic"] = self.images[name][0][1]
self.images[name] = [_pack_png(m) for m, shape in self.images[name]]
@property
@@ -56,22 +80,24 @@ def views(self):
metadata = []
for name, view in self.dataset:
meta = view.to_json(simple=False)
- meta['name'] = name
- if 'stim' in meta['attrs']:
- meta['attrs']['stim'] = os.path.split(meta['attrs']['stim'])[1]
+ meta["name"] = name
+ if "stim" in meta["attrs"]:
+ meta["attrs"]["stim"] = os.path.split(meta["attrs"]["stim"])[1]
metadata.append(meta)
return metadata
def reorder(self, subjects):
- indices = dict((k, np.load(os.path.splitext(v)[0]+".npz")) for k, v in subjects.items())
+ indices = dict(
+ (k, np.load(os.path.splitext(v)[0] + ".npz")) for k, v in subjects.items()
+ )
for brain in self.uniques:
if isinstance(brain, (dataset.Vertex, dataset.VertexRGB)):
data = np.array(self.images[brain.name])[0]
npyform = BytesIO()
- if self.brains[brain.name]['raw']:
- data = data[..., indices[brain.subject]['index'], :]
+ if self.brains[brain.name]["raw"]:
+ data = data[..., indices[brain.subject]["index"], :]
else:
- data = data[..., indices[brain.subject]['index']]
+ data = data[..., indices[brain.subject]["index"]]
np.save(npyform, np.ascontiguousarray(data))
npyform.seek(0)
self.images[brain.name] = [npyform.read()]
@@ -81,8 +107,10 @@ def reorder(self, subjects):
def metadata(self, submap=None, **kwargs):
if submap is not None:
for data in self.brains.values():
- data['subject'] = submap[data['subject']]
- return dict(views=self.views, data=self.brains, images=self.image_names(**kwargs))
+ data["subject"] = submap[data["subject"]]
+ return dict(
+ views=self.views, data=self.brains, images=self.image_names(**kwargs)
+ )
def image_names(self, fmt="/data/{name}/{frame}/"):
names = dict()
@@ -90,14 +118,16 @@ def image_names(self, fmt="/data/{name}/{frame}/"):
names[name] = [fmt.format(name=name, frame=i) for i in range(len(imgs))]
return names
+
def _pack_png(mosaic):
from PIL import Image
+
buf = BytesIO()
if mosaic.dtype not in (np.float32, np.uint8):
raise TypeError
y, x = mosaic.shape[:2]
- im = Image.frombuffer('RGBA', (x,y), mosaic.data, 'raw', 'RGBA', 0, 1)
- im.save(buf, format='PNG')
+ im = Image.frombuffer("RGBA", (x, y), mosaic.data, "raw", "RGBA", 0, 1)
+ im.save(buf, format="PNG")
buf.seek(0)
return buf.read()
diff --git a/cortex/webgl/resources/css/mriview.css b/cortex/webgl/resources/css/mriview.css
index f365dc8bd..47962cbbc 100644
--- a/cortex/webgl/resources/css/mriview.css
+++ b/cortex/webgl/resources/css/mriview.css
@@ -116,25 +116,10 @@ a:visited { color:white; }
margin:20px;
border-radius:10px;
background:rgba(255,255,255,.2);
- max-width:400px;
+ width:max-content;
+ min-width:200px;
+ max-width:90vw;
display:none;
- /*width:320px;
- height:58px;
-*/
- transition:all .3s;
- transition-timing-function:ease;
- transition-delay:1s;
- -moz-transition: all .3s;
- -moz-transition-timing-function: ease;
- -moz-transition-delay:1s;
- -webkit-transition: all .3s;
- -webkit-transition-timing-function: ease;
- -webkit-transition-delay:1s;
-}
-#dataopts:hover, #dataopts:active {
- transition-delay:.1s;
- -moz-transition-delay:.1s;
- -webkit-transition-delay:.1s;
}
#dataname {
color:white;
@@ -486,12 +471,14 @@ input.vlim {
/* datasets */
ul#datasets {
- max-width: 400px;
list-style: none;
margin: 0;
padding: 0;
- overflow:auto;
+ overflow-y: auto;
+ overflow-x: hidden;
max-height:70vh;
+ width: max-content;
+ max-width: 100%;
}
ul#datasets li {
background: white;
@@ -501,6 +488,9 @@ ul#datasets li {
border: 1px solid black;
list-style: none;
padding-left: 42px;
+ white-space: nowrap;
+ overflow: hidden;
+ text-overflow: ellipsis;
}
ul#datasets li .handle {
background: #CCC;
diff --git a/cortex/webgl/resources/js/dataset.js b/cortex/webgl/resources/js/dataset.js
index 762210eb8..9a0a71643 100644
--- a/cortex/webgl/resources/js/dataset.js
+++ b/cortex/webgl/resources/js/dataset.js
@@ -109,34 +109,40 @@ var dataset = (function(module) {
this._dispatch = this.dispatchEvent.bind(this);
- //Aggregate all Volume/VertexData deferreds to determine when to resolve
+ // Aggregate all child Volume/VertexData deferreds to determine when
+ // to resolve. Each child is tracked separately by registering on its
+ // own .loaded — using $.when().progress and looping over data.length
+ // would mark every child ready on a single child's notify, since
+ // $.when's combined progress event doesn't say which source fired.
+ // That ordering bug caused setData → active.set() to dispatch
+ // verts/textures for a sibling that hadn't pushed yet.
var allready = [];
for (var i = 0; i < this.data.length; i++) {
allready.push(false);
}
-
- var deferred = this.data.length == 1 ?
- $.when(this.data[0].loaded) :
- $.when(this.data[0].loaded, this.data[1].loaded);
- deferred
- .progress(function(available) {
- for (var i = 0; i < this.data.length; i++) {
- //TODO: fix this load order
- if (available > this.delay && !allready[i]) {
- allready[i] = true;
-
- //Resolve this deferred if ALL the BrainData objects are loaded (for multiviews)
- var test = true;
- for (var i = 0; i < allready.length; i++)
- test = test && allready[i];
- if (test)
- this.loaded.resolve();
- }
- }
- }.bind(this))
- .done(function() {
+ var checkResolve = function() {
+ for (var j = 0; j < allready.length; j++)
+ if (!allready[j]) return;
this.loaded.resolve();
- }.bind(this));
+ }.bind(this);
+ var markReady = function(idx) {
+ if (!allready[idx]) {
+ allready[idx] = true;
+ checkResolve();
+ }
+ };
+ for (var i = 0; i < this.data.length; i++) {
+ (function(idx) {
+ this.data[idx].loaded
+ .progress(function(available) {
+ if (available > this.delay) markReady(idx);
+ }.bind(this))
+ .done(function() { markReady(idx); });
+ }.bind(this))(i);
+ }
+ // Handle the empty-data edge case so this.loaded still resolves
+ // (matches the pre-refactor $.when() semantics for no arguments).
+ checkResolve();
this.ui = new jsplot.Menu();
@@ -275,6 +281,31 @@ var dataset = (function(module) {
for (var i = 0; i < this.data.length; i++) {
this.data[i].set(this.uniforms, i, fframe, this._dispatch);
}
+ // Combine per-dim NaN masks into the single shared nanmask
+ // attribute. Vertex2D dispatches each dim's data separately
+ // (data0/1 vs data2/3) but shares one nanmask attribute in the
+ // shader; if either dim's value is NaN at a vertex, that vertex
+ // must be discarded.
+ if (this.vertex && !this.data[0].raw && this.data[0].nanmasks.length > 0) {
+ var dim0 = this.data[0].nanmasks[fframe];
+ var combined;
+ if (this.data.length === 1) {
+ combined = dim0;
+ } else {
+ combined = [0, 1].map(function(side) {
+ var a = dim0[side].array;
+ var b = this.data[1].nanmasks[fframe][side].array;
+ var out = new Float32Array(a.length);
+ for (var i = 0; i < a.length; i++) {
+ out[i] = (a[i] < 0.5 || b[i] < 0.5) ? 0.0 : 1.0;
+ }
+ var attr = new THREE.BufferAttribute(out, 1);
+ attr.needsUpdate = true;
+ return attr;
+ }.bind(this));
+ }
+ this._dispatch({type:"attribute", name:"nanmask", value:combined});
+ }
}
module.DataView.prototype.setFilter = function(interp) {
this.filter = interp;
@@ -419,31 +450,34 @@ var dataset = (function(module) {
}
} else {
- // Remap indices and detect NaN in a single pass.
- // WebGL drivers may sanitize NaN in vertex attributes,
- // so we build a mask and replace NaN with 0 here.
- var hasNaN = false;
- for (var i = 0; i < sleft.length; i++) {
- sleft[i] = left[hemis.left.reverseIndexMap[i]];
- if (isNaN(sleft[i])) hasNaN = true;
- }
- for (var i = 0; i < sright.length; i++) {
- sright[i] = right[hemis.right.reverseIndexMap[i]];
- if (isNaN(sright[i])) hasNaN = true;
- }
- if (hasNaN) {
- var masks = [sleft, sright].map(function(arr) {
- var mask = new Float32Array(arr.length);
- for (var i = 0; i < arr.length; i++) {
- if (isNaN(arr[i])) { mask[i] = 0.0; arr[i] = 0.0; }
- else { mask[i] = 1.0; }
+ // Remap indices and build the NaN mask in a single
+ // pass. WebGL drivers may sanitize NaN in vertex
+ // attributes, so we always build a mask attribute
+ // (1=valid, 0=NaN) and replace NaN with 0 in the
+ // data. The shader uses the mask to discard NaN
+ // vertices. We always push a mask (all 1s when
+ // there are no NaNs) so per-frame indexing in
+ // VertexData.set stays aligned with this.verts.
+ var masks = [
+ {dest: sleft, src: left, map: hemis.left.reverseIndexMap},
+ {dest: sright, src: right, map: hemis.right.reverseIndexMap}
+ ].map(function(h) {
+ var mask = new Float32Array(h.dest.length);
+ for (var i = 0; i < h.dest.length; i++) {
+ var val = h.src[h.map[i]];
+ if (isNaN(val)) {
+ h.dest[i] = 0.0;
+ mask[i] = 0.0;
+ } else {
+ h.dest[i] = val;
+ mask[i] = 1.0;
}
- var attr = new THREE.BufferAttribute(mask, 1);
- attr.needsUpdate = true;
- return attr;
- });
- this.nanmasks.push(masks);
- }
+ }
+ var attr = new THREE.BufferAttribute(mask, 1);
+ attr.needsUpdate = true;
+ return attr;
+ });
+ this.nanmasks.push(masks);
}
var lattr = new THREE.BufferAttribute(sleft, this.raw?4:1);
var rattr = new THREE.BufferAttribute(sright, this.raw?4:1);
@@ -467,9 +501,8 @@ var dataset = (function(module) {
var name = dim == 0 ? "data0":"data2";
dispatch({type:"attribute", name:"data"+(2*dim), value:this.verts[fframe]});
dispatch({type:"attribute", name:"data"+(2*dim+1), value:this.verts[(fframe+1).mod(this.verts.length)]});
- if (this.nanmasks.length > 0) {
- dispatch({type:"attribute", name:"nanmask", value:this.nanmasks[fframe]});
- }
+ // The combined nanmask is dispatched by DataView.setFrame after
+ // every dim's data has been set, so we don't dispatch it here.
}
return module;
diff --git a/cortex/webgl/resources/js/facepick.js b/cortex/webgl/resources/js/facepick.js
index a6cedf708..0747ff576 100644
--- a/cortex/webgl/resources/js/facepick.js
+++ b/cortex/webgl/resources/js/facepick.js
@@ -82,8 +82,27 @@ function PickPosition(surf, posdata) {
PickPosition.prototype = {
//Set the transform for any voxels
apply: function(dataview) {
- if (dataview)
- this.xfm.set.apply(this.xfm, dataview.xfm);
+ if (dataview) {
+ if (dataview.xfm) {
+ // For 2D volume dataviews, dataview.xfm is [xfm_dim1, xfm_dim2];
+ // both share xfmname (enforced server-side), so dim1's transform
+ // is correct for picker -> voxel conversion. For 1D volume,
+ // dataview.xfm is a flat 16-element array. Mirror the dataset.js
+ // volxfm length check (see VolumeData.init).
+ var xfm = (dataview.xfm.length === 16) ? dataview.xfm : dataview.xfm[0];
+ this.xfm.set.apply(this.xfm, xfm);
+ } else {
+ // Vertex dataviews have no xfm. Setting the matrix to NaN
+ // suppresses voxel-axis marker rendering (Three.js culls NaN
+ // positions). This mirrors the implicit pre-existing
+ // behavior where `set.apply(xfm, undefined)` invoked
+ // `Matrix4.set()` with no args, leaving the matrix
+ // un-numerable — voxel markers don't belong on a vertex
+ // view (you're picking a vertex, not a voxel).
+ this.xfm.set(NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
+ NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN);
+ }
+ }
for (var i = 0; i < this.axes.length; i++) {
var ax = this.axes[i];
diff --git a/cortex/webgl/resources/js/mriview.js b/cortex/webgl/resources/js/mriview.js
index 8c220a604..a407b03c7 100644
--- a/cortex/webgl/resources/js/mriview.js
+++ b/cortex/webgl/resources/js/mriview.js
@@ -1,5 +1,61 @@
var mriview = (function(module) {
var grid_shapes = [null, [1,1], [2, 1], [3, 1], [2, 2], [2, 2], [3, 2], [3, 2]];
+
+ // Returns true iff every VolumeData entry in `dataview.data` shares the
+ // same shape, mosaic, and numslices, AND every dim's transform matches
+ // dim0's. The hover/click value readout for 2D volume views computes
+ // mouse_index from dim0's voxel coordinate (via picker.xfm = dim0 xfm)
+ // and reuses it across all dims; that is only correct when both shape/
+ // mosaic AND the transforms match. Python-side Volume2D enforces matching
+ // xfmname (so xfms agree); dataset.makeFrom can pair volumes with
+ // matching shape but different transforms — bail out in that case.
+ module.volumeGeomsMatch = function (dataview) {
+ var data = dataview && dataview.data;
+ if (!data || data.length < 2) return true;
+ var d0 = data[0];
+ for (var i = 1; i < data.length; i++) {
+ var di = data[i];
+ if (!d0.shape || !di.shape || !d0.mosaic || !di.mosaic) return false;
+ if (d0.shape[0] !== di.shape[0] || d0.shape[1] !== di.shape[1]) return false;
+ if (d0.mosaic[0] !== di.mosaic[0] || d0.mosaic[1] !== di.mosaic[1]) return false;
+ if (d0.numslices !== di.numslices) return false;
+ }
+ // For 2D volume dataviews, dataview.xfm is [xfm_dim0, xfm_dim1, ...];
+ // each entry is a flat 16-element array. (For 1D the xfm is itself a
+ // flat 16-element array, no per-dim entries to compare — handled by
+ // the early return above.)
+ var xfm = dataview.xfm;
+ if (Array.isArray(xfm) && xfm.length === data.length && Array.isArray(xfm[0])) {
+ var x0 = xfm[0];
+ for (var j = 1; j < xfm.length; j++) {
+ var xj = xfm[j];
+ if (!Array.isArray(xj) || xj.length !== x0.length) return false;
+ for (var k = 0; k < x0.length; k++) {
+ if (x0[k] !== xj[k]) return false;
+ }
+ }
+ }
+ return true;
+ };
+
+ // Returns true iff every entry in `data` has the per-frame buffers the
+ // hover/click handlers need (verts for vertex data, textures for
+ // volume data). Loading is async — the dataview's `loaded` deferred
+ // can resolve a tick before all child VertexData/VolumeData buffers
+ // are populated, and the setData refire fires inside that window.
+ module.dataBuffersReady = function (data) {
+ if (!data || data.length === 0) return false;
+ for (var i = 0; i < data.length; i++) {
+ var d = data[i];
+ if (d.mosaic !== undefined) {
+ if (!d.textures || d.textures.length === 0) return false;
+ } else {
+ if (!d.verts || d.verts.length === 0) return false;
+ }
+ }
+ return true;
+ };
+
module.Viewer = function(figure) {
jsplot.Axes.call(this, figure);
@@ -21,7 +77,16 @@ var mriview = (function(module) {
var tex = new THREE.Texture(this);
tex.minFilter = THREE.LinearFilter;
tex.magFilter = THREE.LinearFilter;
- tex.premultiplyAlpha = false;
+ // Premultiply alpha on upload so the shader's premultiplied-over
+ // composite (gl_FragColor = vColor + (1-α)·cColor in shaderlib.js)
+ // produces the correct α·rgb + (1-α)·bg result when sampling
+ // alpha-bearing colormaps (e.g. RdBu_r_alpha, fire_alpha). For
+ // non-alpha cmaps every row has α=255, so premultiplication is a
+ // no-op. Without this, alpha-cmap-rendered Vertex2D/Volume2D
+ // foregrounds clip toward white instead of fading to the
+ // curvature underlay (companion fix to issue #631 for the LUT
+ // texture path).
+ tex.premultiplyAlpha = true;
tex.flipY = true;
tex.needsUpdate = true;
colormaps[this.parentNode.id] = tex;
@@ -229,8 +294,28 @@ var mriview = (function(module) {
this.setData(data[0].name);
};
+ module.Viewer.prototype.fitDataname = function() {
+ var dn = document.getElementById("dataname");
+ if (!dn || !dn.offsetWidth) return;
+ dn.style.fontSize = "";
+ var available = Math.floor(window.innerWidth * 0.9) - 60;
+ var width = dn.scrollWidth;
+ if (width <= available) return;
+ var maxPx = 32, minPx = 14;
+ var fitted = Math.max(minPx, Math.floor(maxPx * available / width));
+ dn.style.fontSize = fitted + "px";
+ };
+
module.Viewer.prototype.setData = function(name) {
+ // Hide any previously displayed hover/click values immediately. They
+ // refer to the OLD dataview; we will re-fire them against the NEW
+ // dataview once it has loaded (see this.active.loaded.done() below).
+ // Without this transient hide, a stale number is briefly visible
+ // during the (sometimes slow) load of the new dataset.
+ $('#mouseover_value').css('display', 'none')
+ $('#picked_value').css('display', 'none')
+
// blur any selected input elements
let ids = [
['#vmin', '#vmin-input'],
@@ -291,7 +376,28 @@ var mriview = (function(module) {
this.active.loaded.done(function() {
this.active.set();
// Register event for dataset switching
- this.dispatchEvent({type:"setData", name:this.active.name});
+ this.dispatchEvent({type:"setData", name:this.active.name});
+ // Push the new dataview into each surface (shaders + picker xfm)
+ // before re-firing hover/click. drawView normally does this on
+ // the next render frame, but pick() needs the picker's xfm to
+ // already match the new dataview — otherwise a vertex→volume
+ // switch picks with the prior (vertex-view, NaN) transform.
+ for (var i = 0; i < this.surfs.length; i++) {
+ if (this.surfs[i].apply !== undefined)
+ this.surfs[i].apply(this.active);
+ }
+ // Re-fire hover and click indicators so the displayed values
+ // reflect the newly-active dataview at the same screen positions
+ // (see _lastHoverEvent / _lastPickEvt cached by the handlers).
+ if (this._lastHoverEvent) {
+ $("#brain").trigger($.Event('mousemove', {
+ clientX: this._lastHoverEvent.clientX,
+ clientY: this._lastHoverEvent.clientY,
+ }));
+ }
+ if (this._lastPickEvt) {
+ this.pick({x: this._lastPickEvt.x, y: this._lastPickEvt.y});
+ }
}.bind(this));
// var surf, scene, grid = grid_shapes[this.active.data.length];
@@ -590,6 +696,7 @@ var mriview = (function(module) {
$("#dataname").text(name);
$("#dataopts").show();
}
+ this.fitDataname();
this.schedule();
this.loaded.resolve();
@@ -645,40 +752,64 @@ var mriview = (function(module) {
$("#brain").on('mousemove',
function (event) {
- // only implemented for 1d volume datasets or vertex datasets
- if (this.active.data.length != 1 || this.active.data[0].raw) {
+ // Cache last mouse position so setData() can recompute the
+ // hover indicator for the newly-active dataset at the same
+ // screen coordinate.
+ this._lastHoverEvent = event;
+ // Length check first so the short-circuit guards us against
+ // an empty data array before we read data[0]. Then skip RGB
+ // (no underlying scalars), then ensure all child buffers
+ // have populated (the setData refire can briefly precede
+ // buffer fill).
+ if ((this.active.data.length !== 1 && this.active.data.length !== 2) ||
+ this.active.data[0].raw ||
+ !module.dataBuffersReady(this.active.data)) {
$('#mouseover_value').css('display', 'none')
return
}
// We need to use a different logic if we have a VolumeData or a VertexData object
- let value = null;
+ let values = null;
if (this.active.vertex) {
- coords = this.getCoords(event)
+ let coords = this.getCoords(event)
if (coords !== -1) {
- hemiIdx = (coords.hemi == 'left') ? 0 : 1
- // console.log('hemiIdx: ', hemiIdx)
- vertex = coords.vertex
- // console.log('vertex: ', vertex)
- // Now we need to map back with the index map
- // First figure out the subject, then get the index map
- subject = this.active.data[0].subject
- indexMap = subjects[subject].hemis[coords.hemi].indexMap
- // console.log("vertex before: " + vertex);
- vertex = indexMap[vertex]
- // console.log("vertex after: " + vertex);
- // Now access the data
- value = this.active.data[0].verts[0][hemiIdx].array[vertex]
+ let hemiIdx = (coords.hemi == 'left') ? 0 : 1
+ // Map the picked vertex through the subject's indexMap.
+ // dim1 and dim2 share subject for 2D views (server-side).
+ let subject = this.active.data[0].subject
+ let indexMap = subjects[subject].hemis[coords.hemi].indexMap
+ let vertex = indexMap[coords.vertex]
+ // Now access the data for each channel (1 for 1D, 2 for 2D)
+ values = this.active.data.map(function (d) {
+ return d.verts[0][hemiIdx].array[vertex]
+ })
}
} else {
- // Get index of the mosaic to get the value
- mouse_index = this.getMouseIndex(event)
+ // Volume branch. For 2D views we reuse data[0]'s mouse_index
+ // across all dims; that is safe iff every dim shares the
+ // same shape/mosaic/numslices. Python-side Volume2D enforces
+ // matching xfmname → matching geometry. Client-side
+ // dataset.makeFrom (mriview.js:283) can pair arbitrary 1D
+ // volumes; if their geometries diverge, hide the readout
+ // rather than display a wrong-but-plausible value.
+ if (!module.volumeGeomsMatch(this.active)) {
+ $('#mouseover_value').css('display', 'none')
+ return
+ }
+ let mouse_index = this.getMouseIndex(event)
if (mouse_index !== -1) {
- value = this.active.data[0].textures[0].image.data[mouse_index]
+ values = this.active.data.map(function (d) {
+ return d.textures[0].image.data[mouse_index]
+ })
}
}
- // console.log("Value on mouseover: " + value);
- if (value !== null) {
- $('#mouseover_value').text(parseFloat(value).toPrecision(3))
+ if (values !== null) {
+ let formatted = values.map(function (v) {
+ return parseFloat(v).toPrecision(3)
+ }).join(', ')
+ if (values.length > 1) {
+ formatted = '(' + formatted + ')'
+ }
+ $('#mouseover_value').text(formatted)
$('#mouseover_value').css('display', 'block')
} else {
$('#mouseover_value').css('display', 'none')
@@ -817,44 +948,63 @@ var mriview = (function(module) {
}
module.Viewer.prototype.pick = function(evt) {
+ // Cache last pick position so setData() can refresh the picked
+ // indicator for the newly-active dataset at the same screen point.
+ this._lastPickEvt = {x: evt.x, y: evt.y};
let coords
for (var i = 0; i < this.surfs.length; i++) {
if (this.surfs[i].pick)
coords = this.surfs[i].pick(this.renderer, this.camera, evt.x, evt.y);
}
// set the picked value display
- // only implemented for 1d volume datasets or vertex datasets
- if (this.active.data.length != 1 || this.active.data[0].raw) {
+ // Length check first so we don't index data[0] on an empty array.
+ // Skip RGB, then ensure all child buffers have populated.
+ if ((this.active.data.length !== 1 && this.active.data.length !== 2) ||
+ this.active.data[0].raw ||
+ !module.dataBuffersReady(this.active.data)) {
$('#picked_value').css('display', 'none')
return
}
// We need to use a different logic if we have a VolumeData or a VertexData object
- let value = null;
+ let values = null;
if (this.active.vertex) {
if (coords !== -1) {
- hemiIdx = (coords.hemi == 'left') ? 0 : 1
- vertex = coords.vertex
- // Now we need to map back with the index map
- // First figure out the subject, then get the index map
- subject = this.active.data[0].subject
- indexMap = subjects[subject].hemis[coords.hemi].indexMap
- vertex = indexMap[vertex]
- // Now access the data
- value = this.active.data[0].verts[0][hemiIdx].array[vertex]
+ let hemiIdx = (coords.hemi == 'left') ? 0 : 1
+ // Map the picked vertex through the subject's indexMap.
+ // dim1 and dim2 share subject for 2D views (server-side).
+ let subject = this.active.data[0].subject
+ let indexMap = subjects[subject].hemis[coords.hemi].indexMap
+ let vertex = indexMap[coords.vertex]
+ // Now access the data for each channel (1 for 1D, 2 for 2D)
+ values = this.active.data.map(function (d) {
+ return d.verts[0][hemiIdx].array[vertex]
+ })
}
} else {
if (coords !== -1) {
- // Get index of the mosaic to get the value
+ // Volume branch. See the matching note in the mousemove handler
+ // for why we hide on geometry mismatch.
+ if (!module.volumeGeomsMatch(this.active)) {
+ $('#picked_value').css('display', 'none')
+ return
+ }
let mouse_index = this.xyxToI(coords.voxel.x, coords.voxel.y, coords.voxel.z)
if (mouse_index !== -1) {
- value = this.active.data[0].textures[0].image.data[mouse_index]
+ values = this.active.data.map(function (d) {
+ return d.textures[0].image.data[mouse_index]
+ })
}
}
}
- console.log("Value on click: " + value);
- if (value !== null) {
- $('#picked_value').text(parseFloat(value).toPrecision(3))
+ if (values !== null) {
+ let formatted = values.map(function (v) {
+ return parseFloat(v).toPrecision(3)
+ }).join(', ')
+ if (values.length > 1) {
+ formatted = '(' + formatted + ')'
+ }
+ $('#picked_value').text(formatted)
$('#picked_value').css('display', 'block')
} else {
$('#picked_value').css('display', 'none')
@@ -1000,7 +1150,7 @@ var mriview = (function(module) {
var _bound = false;
module.Viewer.prototype._bindUI = function() {
$(window).scrollTop(0);
- $(window).resize(function() { this.resize(); }.bind(this));
+ $(window).resize(function() { this.resize(); this.fitDataname(); }.bind(this));
this.canvas.resize(function() { this.resize(); }.bind(this));
var cam_ui = this.ui.addFolder("camera", true);
diff --git a/cortex/webgl/resources/js/shaderlib.js b/cortex/webgl/resources/js/shaderlib.js
index 0ec91c8e5..1c40ed603 100644
--- a/cortex/webgl/resources/js/shaderlib.js
+++ b/cortex/webgl/resources/js/shaderlib.js
@@ -755,6 +755,8 @@ var Shaderlib = (function() {
"vColor = texture2D(colormap, cuv);",
// NaN mask: WebGL drivers sanitize NaN in vertex attributes,
// so we detect NaN in JavaScript and pass a mask (0=NaN, 1=valid).
+ // For 2D vertex views the JS layer combines per-dim masks
+ // before dispatch, so a single shared attribute is enough.
"if (nanmask < 0.5) vColor = vec4(0.);",
"#endif",
diff --git a/examples/datasets/plot_data_with_alpha.py b/examples/datasets/plot_data_with_alpha.py
new file mode 100644
index 000000000..b4df461e0
--- /dev/null
+++ b/examples/datasets/plot_data_with_alpha.py
@@ -0,0 +1,195 @@
+"""
+==========================
+Plot Data with Alpha Values
+==========================
+
+It is often useful to plot a primary map (the "data" you are interested in)
+masked or attenuated by a secondary map (a "confidence" or "weight"
+map). For example, an encoding model's tuning maps are
+interpretable where the model fits well, so one could plot
+tuning maps with opacity proportional to the per-voxel/per-vertex prediction
+accuracy. Voxels/vertices where the model fits poorly fade into the
+gray curvature underlay; voxels/vertices where the model fits well are
+fully opaque.
+
+pycortex supports two patterns for this:
+
+1. **Scalar data with an alpha map** -- use :class:`Volume2D` /
+ :class:`Vertex2D` with a 2D colormap whose second axis encodes alpha
+ (the colormap LUT itself goes from transparent to opaque along
+ ``dim2``). No extra arithmetic is needed; the alpha channel is
+ composited correctly by both the matplotlib (``quickshow``) and the
+ WebGL renderers.
+
+2. **RGB data with an alpha map** -- pass ``alpha=`` directly to
+ :class:`VolumeRGB` / :class:`VertexRGB`. The alpha can be any
+ per-voxel/per-vertex array (or a :class:`Volume`/:class:`Vertex`)
+ in ``[0, 1]``.
+
+Below, we illustrate both patterns with a synthetic "model accuracy"
+mask -- a 3D Gaussian bump for the volume case and a vertex-distance
+falloff for the surface case -- so cortex near the bump centre stays
+opaque while the periphery fades into the curvature.
+"""
+
+import cortex
+import cortex.polyutils
+import numpy as np
+import matplotlib.pyplot as plt
+
+subject = "S1"
+xfm = "fullhead"
+
+# %%
+# Synthesize the data and alpha maps
+# ----------------------------------
+#
+# All four patterns below reuse the same synthetic inputs, so we set
+# everything up once here and only show the *plotting* call in each
+# pattern's cell. In a real analysis these would come from your model
+# fits (e.g. ``data`` = regression coefficients, ``accuracy`` =
+# cross-validated prediction r^2).
+
+# --- Volumetric data + alpha -------------------------------------------------
+# Signed gradient across the brain stands in for a regression coefficient
+# or tuning preference.
+zz, yy, xx = np.mgrid[0:31, 0:100, 0:100]
+data_vol = (xx - 50) / 50.0 # range ~ [-1, 1]
+
+# A 3D Gaussian bump centered in the volume stands in for a per-voxel
+# model accuracy / prediction r in [0, 1].
+center = np.array([15, 50, 50])
+sigma = 25.0
+dist2 = (zz - center[0]) ** 2 + (yy - center[1]) ** 2 + (xx - center[2]) ** 2
+accuracy_vol = np.exp(-dist2 / (2 * sigma**2)) # in [0, 1]
+
+# RGB volumetric channels: anatomical x/y/z normalized to [0, 1]. Three
+# smoothly-varying volumetric channels stand in for, e.g., three latent
+# RGB tuning axes from a model.
+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) data + alpha -------------------------------------------
+# Encode by *spatial coordinate*, not vertex index: vertex indices on the
+# cortical surface are not arranged by spatial neighborhood, so a
+# vertex-index ramp would render as visual noise.
+surfs = [cortex.polyutils.Surface(*d) for d in cortex.db.get_surf(subject, "fiducial")]
+num_verts = [s.pts.shape[0] for s in surfs]
+total_verts = sum(num_verts)
+pts = np.vstack([surfs[0].pts, surfs[1].pts]) # (total_verts, 3)
+
+# Scalar surface data: anterior-posterior gradient (anatomical y),
+# centered at zero so a diverging colormap reads naturally.
+y_centered = pts[:, 1] - pts[:, 1].mean()
+data_vtx = y_centered / np.abs(y_centered).max() # in [-1, 1]
+
+# RGB surface channels: anatomical x/y/z normalized to [0, 1].
+xyz_norm = (pts - pts.min(axis=0)) / (pts.max(axis=0) - pts.min(axis=0))
+
+
+# Surface alpha: a soft bump centered at a particular vertex in each hemi.
+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),
+ ]
+)
+
+# %%
+# Pattern 1a: scalar Volume + alpha via Volume2D + 2D alpha colormap
+# ------------------------------------------------------------------
+#
+# The 2D colormap ``"RdBu_r_alpha"`` maps ``(data, alpha) -> RGBA``: along
+# the first axis, blue-white-red diverging; along the second axis,
+# transparent-to-opaque. So passing the data as ``dim1`` and the accuracy
+# as ``dim2`` yields exactly "diverging colormap, opacity = accuracy".
+#
+# Other 2D alpha colormaps shipped with pycortex include ``"fire_alpha"``
+# (sequential, perceptually uniform), ``"PU_RdBu_covar_alpha"`` (diverging,
+# perceptually uniform), ``"plasma_alpha"``, and ``"autumn_alpha"``.
+v2d = cortex.Volume2D(
+ data_vol,
+ accuracy_vol,
+ subject,
+ xfm,
+ cmap="RdBu_r_alpha",
+ vmin=-1,
+ vmax=1, # range for the data (dim1)
+ vmin2=0,
+ vmax2=1, # range for the alpha (dim2)
+)
+cortex.quickshow(v2d, with_colorbar=True, with_curvature=True)
+plt.suptitle("Volume2D + RdBu_r_alpha: data masked by 'accuracy'")
+plt.show()
+
+# %%
+# Pattern 1b: scalar Vertex + alpha via Vertex2D + 2D alpha colormap
+# ------------------------------------------------------------------
+#
+# Same idea on the surface.
+vtx2d = cortex.Vertex2D(
+ data_vtx,
+ accuracy_vtx,
+ subject,
+ cmap="RdBu_r_alpha",
+ vmin=-1,
+ vmax=1,
+ vmin2=0,
+ vmax2=1,
+)
+cortex.quickshow(vtx2d, with_colorbar=True, with_curvature=True)
+plt.suptitle("Vertex2D + RdBu_r_alpha: data masked by 'accuracy'")
+plt.show()
+
+# %%
+# Pattern 2a: RGB Volume + alpha via VolumeRGB(alpha=...)
+# -------------------------------------------------------
+#
+# When the "data" is itself three independent channels, use
+# :class:`VolumeRGB` and pass the accuracy as the ``alpha=`` argument.
+red = cortex.Volume(red_vol, subject, xfm, vmin=0, vmax=1)
+green = cortex.Volume(green_vol, subject, xfm, vmin=0, vmax=1)
+blue = cortex.Volume(blue_vol, subject, xfm, vmin=0, vmax=1)
+alpha_vol = cortex.Volume(accuracy_vol, subject, xfm, vmin=0, vmax=1)
+
+vrgb = cortex.VolumeRGB(red, green, blue, subject, alpha=alpha_vol)
+cortex.quickshow(vrgb, with_colorbar=False, with_curvature=True)
+plt.suptitle("VolumeRGB(alpha=accuracy): RGB tuning masked by 'accuracy'")
+plt.show()
+
+# %%
+# Pattern 2b: RGB Vertex + alpha via VertexRGB(alpha=...)
+# -------------------------------------------------------
+#
+# Same idea on the surface.
+red_v = cortex.Vertex(xyz_norm[:, 0], subject, vmin=0, vmax=1)
+green_v = cortex.Vertex(xyz_norm[:, 1], subject, vmin=0, vmax=1)
+blue_v = cortex.Vertex(xyz_norm[:, 2], subject, vmin=0, vmax=1)
+alpha_v = cortex.Vertex(accuracy_vtx, subject, vmin=0, vmax=1)
+
+vrgb_vtx = cortex.VertexRGB(red_v, green_v, blue_v, subject, alpha=alpha_v)
+cortex.quickshow(vrgb_vtx, with_colorbar=False, with_curvature=True)
+plt.suptitle("VertexRGB(alpha=accuracy): RGB channels masked by 'accuracy'")
+plt.show()
+
+# %%
+# Notes
+# -----
+#
+# * Both patterns produce the same composite formula at the pixel level:
+# ``out = alpha * data + (1 - alpha) * curvature_underlay``. Choose
+# based on what the "data" is: scalar (use Pattern 1) or RGB (use
+# Pattern 2).
+# * The same objects work in the WebGL viewer:
+# ``cortex.webgl.show(v2d)`` etc.; opacity is honored identically.
+# * The deprecated ``Vertex.blend_curvature(alpha)`` helper produced a
+# pre-blended :class:`VertexRGB` that lost ``cmap``/``vmin``/``vmax``
+# editability. The Pattern 1 :class:`Vertex2D` route above is the
+# recommended replacement: it keeps the colormap parameters editable
+# on the resulting object and renders identically.
diff --git a/filestore/colormaps/BuWtRd.png b/filestore/colormaps/BuWtRd.png
index 676f80734..7ae596af6 100644
Binary files a/filestore/colormaps/BuWtRd.png and b/filestore/colormaps/BuWtRd.png differ
diff --git a/filestore/colormaps/BuWtRd_alpha.png b/filestore/colormaps/BuWtRd_alpha.png
index 08d2f8d05..2101a2db5 100644
Binary files a/filestore/colormaps/BuWtRd_alpha.png and b/filestore/colormaps/BuWtRd_alpha.png differ
diff --git a/pyproject.toml b/pyproject.toml
index e65c04616..6acd93aad 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,7 +1,7 @@
[build-system]
# Minimum requirements for the build system to execute, according to PEP518
# specification.
-requires = ["setuptools>=64", "setuptools-scm>=8", "build", "numpy", "cython", "wheel"]
+requires = ["setuptools>=64", "setuptools-scm>=8", "build", "numpy>=1.13.0", "cython", "wheel"]
build-backend = "setuptools.build_meta"
[project]
@@ -20,6 +20,7 @@ test = [
"nbformat>=5.10.4",
"pytest",
"pytest-cov",
+ "pytest-timeout",
]
[project.optional-dependencies]
diff --git a/pytest.ini b/pytest.ini
index 84083f96f..25191cf6d 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -4,3 +4,8 @@ testpaths =
addopts =
-r a
-v
+# Per-test timeout (in seconds) so a single hung headless browser session
+# does not consume the entire CI budget. Individual tests can override with
+# @pytest.mark.timeout(N). Requires the optional ``pytest-timeout``
+# plugin; if not installed, this key is silently ignored.
+timeout = 240
diff --git a/requirements.txt b/requirements.txt
index b05b6dd0b..56fa497c0 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,6 +1,6 @@
setuptools
future
-numpy
+numpy>=1.13.0
scipy
tornado>=4.3
shapely