diff --git a/cortex/dataset/view2D.py b/cortex/dataset/view2D.py index 089deff95..f70b6e813 100644 --- a/cortex/dataset/view2D.py +++ b/cortex/dataset/view2D.py @@ -8,6 +8,7 @@ from .. import options from .views import Dataview, Volume, Vertex, VolumeRGB, VertexRGB +from .viewRGB import _warn_alpha_range from .braindata import BrainData, VolumeData, VertexData default_cmap2D = options.config.get("basic", "default_cmap2D") @@ -20,7 +21,8 @@ class Dataview2D(Dataview): def __init__(self, description: str="", cmap: Optional[str]=None, vmin: Optional[float]=None, vmax: Optional[float]=None, - vmin2: Optional[float]=None, vmax2: Optional[float]=None, state=None, **kwargs): + vmin2: Optional[float]=None, vmax2: Optional[float]=None, state=None, + alpha=None, **kwargs): self.cmap = cmap or default_cmap2D self.vmin = vmin self.vmax = vmax @@ -32,17 +34,78 @@ def __init__(self, description: str="", cmap: Optional[str]=None, if 'priority' not in self.attrs: self.attrs['priority'] = 1 self.description = description + # Optional per-voxel/vertex alpha map. Kept as a Volume/Vertex (never + # inside ``attrs``, which is JSON-serialized for the WebGL viewer) so + # that both renderers can ship it like any other brain. + self.alpha = alpha + + @property + def alpha(self): + """Optional alpha map (Volume/Vertex in [vmin, vmax]) multiplied into + the colormap alpha. NaN anywhere (dim1, dim2 or alpha) renders as + alpha 0 regardless of this map.""" + return self._alpha + + @alpha.setter + def alpha(self, alpha): + if alpha is not None and not isinstance(alpha, self._cls): + alpha = np.asarray(alpha) + _warn_alpha_range(alpha) + alpha = self._wrap_alpha(alpha) + if alpha is not None and alpha.subject != self.dim1.subject: + raise ValueError("alpha must belong to the same subject as dim1") + if alpha is not None and getattr(alpha, "xfmname", None) != getattr(self.dim1, "xfmname", None): + raise ValueError("alpha must use the same transform as dim1") + self._alpha = alpha + self._alpha_brain_cache = None + + @property + def _alpha_brain(self): + """The alpha map normalized to [0, 1] (vmin=0, vmax=1), as shipped to + the WebGL viewer and stored in HDF files. NaN is preserved.""" + if self.alpha is None: + return None + if self._alpha_brain_cache is None: + self._alpha_brain_cache = self._wrap_alpha(self._normalized_alpha()) + return self._alpha_brain_cache + + def _wrap_alpha(self, alpha): + raise NotImplementedError + + def _normalized_alpha(self, full_volume=False): + """User alpha as float array in [0, 1] (NaN preserved), in the space of + ``.data`` or, for volumes with ``full_volume=True``, of ``.volume``.""" + alpha = self.alpha + if alpha is None: + return None + raw = alpha.volume if full_volume else alpha.data + arr = np.asarray(raw, dtype=float) + if np.asarray(raw).dtype == np.uint8: + return arr / 255. + vmin = 0. if alpha.vmin is None else float(alpha.vmin) + vmax = 1. if alpha.vmax is None else float(alpha.vmax) + if vmax == vmin: + return np.where(arr >= vmax, 1., 0.) + return (arr - vmin) / (vmax - vmin) def uniques(self, collapse=False): yield self.dim1 yield self.dim2 + if self.alpha is not None: + yield self._alpha_brain def _write_hdf(self, h5, name="data"): self._cls._write_hdf(self.dim1, h5) self._cls._write_hdf(self.dim2, h5) + names = [self.dim1.name, self.dim2.name] + if self.alpha is not None: + # Stored normalized to [0, 1] so that it can be restored with + # vmin=0, vmax=1 (BrainData nodes carry no range). + self._cls._write_hdf(self._alpha_brain, h5) + names.append(self._alpha_brain.name) viewnode = Dataview._write_hdf(self, h5, name=name) - viewnode[0] = json.dumps([[self.dim1.name, self.dim2.name]]) + viewnode[0] = json.dumps([names]) viewnode[3] = json.dumps([[self.vmin, self.vmin2]]) viewnode[4] = json.dumps([[self.vmax, self.vmax2]]) return viewnode @@ -56,17 +119,30 @@ def to_json(self, simple=False): d1js = self.dim1.to_json() d2js = self.dim2.to_json() + # ``is None`` checks: a legitimate vmin/vmax of 0 must not fall back + # to the auto range (same class of bug as 5482c8bf for Volume). sdict.update(dict( - vmin = [[self.vmin or d1js['vmin'][0], self.vmin2 or d2js['vmin'][0]]], - vmax = [[self.vmax or d1js['vmax'][0], self.vmax2 or d2js['vmax'][0]]], + vmin = [[d1js['vmin'][0] if self.vmin is None else self.vmin, + d2js['vmin'][0] if self.vmin2 is None else self.vmin2]], + vmax = [[d1js['vmax'][0] if self.vmax is None else self.vmax, + d2js['vmax'][0] if self.vmax2 is None else self.vmax2]], )) if "xfm" in d1js: sdict['xfm'] = [[d1js['xfm'][0], d2js['xfm'][0]]] + if self.alpha is not None: + sdict['alpha'] = [self._alpha_brain.name] + return sdict - def _to_raw(self, data1, data2): + def _to_raw(self, data1, data2, alpha=None): + """Colormap (data1, data2) through the 2D colormap. + + Returns ``(r, g, b, a, nan_mask)`` as uint8 channels. ``a`` is the + colormap alpha multiplied by the (normalized) user ``alpha`` map when + given; it is forced to 0 wherever data1, data2 or alpha is NaN. + """ from matplotlib import pyplot as plt from matplotlib.colors import Normalize cmapdir = options.config.get("webgl", "colormaps") @@ -93,17 +169,16 @@ def _to_raw(self, data1, data2): g.shape = dim1.shape b.shape = dim1.shape a.shape = dim1.shape - # Preserve nan values as alpha = 0 + # NaN in either dimension (or in the alpha map) -> alpha = 0 aidx = np.logical_or(np.isnan(data1), np.isnan(data2)) + if alpha is not None: + alpha = np.asarray(alpha, dtype=float) + aidx = np.logical_or(aidx, np.isnan(alpha)) + user = np.clip(np.nan_to_num(alpha, nan=0.0), 0, 1) + a = np.round(a.astype(float) * user).astype(np.uint8) + aidx = np.broadcast_to(aidx, a.shape) a[aidx] = 0 - # Code from main, to handle alpha input, prob better here but not tested. - # # Possibly move this above setting nans to alpha = 0; - # # Possibly multiply specified alpha by alpha in colormap?? - # if 'alpha' in self.attrs: - # # Over-write alpha from colormap / nans with alpha arg if provided. - # # Question: Might it be important tokeep alpha as an attr? - # a = self.attrs.pop('alpha') - return r, g, b, a + return r, g, b, a, aidx @property def subject(self): @@ -134,13 +209,25 @@ class Volume2D(Dataview2D): Colormap (or colormap name) to use. If not given defaults to the `default_cmap2d` in your pycortex options.cfg file. vmin : float, optional - Minimum value in colormap for dim1. If not given defaults to TODO:WHAT + Minimum value in colormap for dim1. If not given, the ``vmin`` of dim1 + is used (its own ``vmin`` for a Volume/Vertex object, the 1st percentile + of the data for an array). vmax : float, optional - Maximum value in colormap for dim1. If not given defaults to TODO:WHAT + Maximum value in colormap for dim1. If not given, the ``vmax`` of dim1 + is used (its own ``vmax`` for a Volume/Vertex object, the 99th + percentile of the data for an array). vmin2 : float, optional - Minimum value in colormap for dim2. If not given defaults to TODO:WHAT + Minimum value in colormap for dim2. If not given, the ``vmin`` of dim2 + is used (same rule as ``vmin``). vmax2 : float, optional - Maximum value in colormap for dim2. If not given defaults to TODO:WHAT + Maximum value in colormap for dim2. If not given, the ``vmax`` of dim2 + is used (same rule as ``vmax``). + alpha : ndarray or Volume/Vertex, optional + Per-voxel (per-vertex) opacity multiplied into the colormap alpha. + Arrays are taken in [0, 1]; Volume/Vertex objects are normalized by + their own ``vmin``/``vmax``. Honored identically by quickflat and the + WebGL viewer. Wherever dim1, dim2 or alpha is NaN the data is + rendered fully transparent, regardless of this map. **kwargs All additional arguments in kwargs are passed to the VolumeData and Dataview @@ -189,18 +276,29 @@ def raw(self) -> VolumeRGB: if self.dim1.xfmname != self.dim2.xfmname: raise ValueError("Both Volumes must have same xfmname to generate single raw volume") - if ((self.dim1.linear and self.dim2.linear) and - (self.dim1.mask.shape == self.dim2.mask.shape) and - np.all(self.dim1.mask == self.dim2.mask)): - r, g, b, a = self._to_raw(self.dim1.data, self.dim2.data) + def _same_mask(a, b): + return (a.linear and b.linear and a.mask.shape == b.mask.shape + and np.all(a.mask == b.mask)) + + linear = _same_mask(self.dim1, self.dim2) + if linear and self.alpha is not None: + linear = _same_mask(self.dim1, self.alpha) + if linear: + r, g, b, a, nan_mask = self._to_raw( + self.dim1.data, self.dim2.data, self._normalized_alpha()) else: - r, g, b, a = self._to_raw(self.dim1.volume, self.dim2.volume) - # Allow manual override of alpha channel - kws = dict(subject=self.dim1.subject, xfmname=self.dim1.xfmname, - state=self.state, description=self.description, **self.attrs) - if not 'alpha' in self.attrs: - kws['alpha'] = a - return VolumeRGB(r, g, b, **kws) + r, g, b, a, nan_mask = self._to_raw( + self.dim1.volume, self.dim2.volume, + self._normalized_alpha(full_volume=True)) + result = VolumeRGB(r, g, b, subject=self.dim1.subject, + xfmname=self.dim1.xfmname, alpha=a, + state=self.state, description=self.description, + priority=self.priority) + result._nan_mask = nan_mask + return result + + def _wrap_alpha(self, alpha): + return Volume(alpha, self.dim1.subject, self.dim1.xfmname, vmin=0, vmax=1) @property @@ -229,13 +327,25 @@ class Vertex2D(Dataview2D): Colormap (or colormap name) to use. If not given defaults to the `default_cmap2d` in your pycortex options.cfg file. vmin : float, optional - Minimum value in colormap for dim1. If not given defaults to TODO:WHAT + Minimum value in colormap for dim1. If not given, the ``vmin`` of dim1 + is used (its own ``vmin`` for a Volume/Vertex object, the 1st percentile + of the data for an array). vmax : float, optional - Maximum value in colormap for dim1. If not given defaults to TODO:WHAT + Maximum value in colormap for dim1. If not given, the ``vmax`` of dim1 + is used (its own ``vmax`` for a Volume/Vertex object, the 99th + percentile of the data for an array). vmin2 : float, optional - Minimum value in colormap for dim2. If not given defaults to TODO:WHAT + Minimum value in colormap for dim2. If not given, the ``vmin`` of dim2 + is used (same rule as ``vmin``). vmax2 : float, optional - Maximum value in colormap for dim2. If not given defaults to TODO:WHAT + Maximum value in colormap for dim2. If not given, the ``vmax`` of dim2 + is used (same rule as ``vmax``). + alpha : ndarray or Volume/Vertex, optional + Per-voxel (per-vertex) opacity multiplied into the colormap alpha. + Arrays are taken in [0, 1]; Volume/Vertex objects are normalized by + their own ``vmin``/``vmax``. Honored identically by quickflat and the + WebGL viewer. Wherever dim1, dim2 or alpha is NaN the data is + rendered fully transparent, regardless of this map. **kwargs All additional arguments in kwargs are passed to the VolumeData and Dataview @@ -278,12 +388,16 @@ def __repr__(self): def raw(self) -> VertexRGB: """VertexRGB object containing the colormapped data from this object. """ - r, g, b, a = self._to_raw(self.dim1.data, self.dim2.data) - # Allow manual override of alpha channel - kws = dict(subject=self.dim1.subject) - if not 'alpha' in self.attrs: - kws['alpha'] = a - return VertexRGB(r, g, b, **kws) + r, g, b, a, nan_mask = self._to_raw( + self.dim1.data, self.dim2.data, self._normalized_alpha()) + result = VertexRGB(r, g, b, subject=self.dim1.subject, alpha=a, + state=self.state, description=self.description, + priority=self.priority) + result._nan_mask = nan_mask + return result + + def _wrap_alpha(self, alpha): + return Vertex(alpha, self.dim1.subject, vmin=0, vmax=1) @property def vertices(self): diff --git a/cortex/dataset/viewRGB.py b/cortex/dataset/viewRGB.py index 8d2c9d4ea..3c2604190 100644 --- a/cortex/dataset/viewRGB.py +++ b/cortex/dataset/viewRGB.py @@ -73,6 +73,75 @@ def HSV2RGB(color: Color[float] | npt.NDArray) -> Color[int]: return (int(r * 255), int(g * 255), int(b * 255)) +def _warn_alpha_range(alpha): + """Warn when a raw (non-uint8) alpha array lies outside [0, 1].""" + alpha = np.asarray(alpha) + if alpha.dtype == np.uint8 or alpha.size == 0: + return + finite = alpha[np.isfinite(alpha)] + if finite.size and (finite.min() < 0 or finite.max() > 1): + warnings.warn( + "Some alpha values are outside the range of [0, 1]. " + "Consider passing a Volume/Vertex object as alpha with explicit " + "vmin, vmax keyword arguments.", + Warning, + ) + + +def _mask_alpha(alpha, mask): + """Return a copy of ``alpha`` (Volume or Vertex) with ``alpha.vmin`` written + wherever ``mask`` is True. + + ``mask`` may live either in the same space as ``alpha.data`` (masked/linear + or full) or, for volumes, in full ``(z, y, x)`` / ``(t, z, y, x)`` space. + Leading (time) axes are broadcast, so a single-frame alpha combined with a + multi-frame NaN mask yields a multi-frame alpha (#629). + + The masked values are never written into a temporary: ``VolumeData.volume`` + returns a freshly unmasked array for linear volumes, so the previous + ``alpha.volume[mask] = vmin`` silently did nothing for those. + """ + mask = np.asarray(mask, dtype=bool) + if not mask.any(): + return alpha + data = np.asarray(alpha.data) + if data.dtype == np.uint8: + # uint8 alpha bypasses the vmin/vmax normalization later on, and its + # inferred vmin is a percentile of the bytes (255 for a constant map): + # transparent is byte 0. + fill = 0 + else: + fill = 0.0 if alpha.vmin is None else alpha.vmin + try: + shape = np.broadcast_shapes(mask.shape, data.shape) + except ValueError: + shape = None + if shape is not None and shape[-data.ndim:] == data.shape: + new = np.array(np.broadcast_to(data, shape)) # copy, keep dtype + new[np.broadcast_to(mask, shape)] = fill + return alpha.copy(new) + + if hasattr(alpha, "volume"): + vol = np.asarray(alpha.volume) # (t, z, y, x), fresh copy if linear + try: + shape = np.broadcast_shapes(mask.shape, vol.shape) + except ValueError: + shape = None + if shape is not None and shape[-3:] == vol.shape[-3:]: + new = np.array(np.broadcast_to(vol, shape)) # copy, keep dtype + new[np.broadcast_to(mask, shape)] = fill + if new.shape[0] == 1 and not alpha.movie: + new = new[0] + return Volume( + new, alpha.subject, alpha.xfmname, vmin=alpha.vmin, vmax=alpha.vmax + ) + + raise ValueError( + "alpha of shape %s is incompatible with data NaN mask of shape %s" + % (data.shape, mask.shape) + ) + + class DataviewRGB(Dataview): """Abstract base class for RGB data views.""" @@ -119,11 +188,8 @@ def _apply_nan_mask(self, alpha: BrainData): in Dataview.raw and stored as ``_nan_mask``.""" nan_mask = getattr(self, "_nan_mask", None) if nan_mask is None: - return - if nan_mask.shape == alpha.data.shape: - alpha.data[nan_mask] = alpha.vmin - elif hasattr(alpha, "volume") and nan_mask.shape == alpha.volume.shape: - alpha.volume[nan_mask] = alpha.vmin + return alpha + return _mask_alpha(alpha, nan_mask) def _write_hdf(self, h5, name="data", xfmname=None): self._cls._write_hdf(self.red, h5) @@ -345,10 +411,16 @@ def color_voxels( green.flat[i] = this_color[1] blue.flat[i] = this_color[2] - # Now make an alpha volume + # Now make an alpha volume. NaN in any channel forces alpha to its + # minimum. Never write into the caller's array. if alpha is None: alpha = np.ones_like(red, np.uint8) * 255 - alpha[mask] = 0 # TODO: this seems like an actual issue + alpha[mask] = 0 + elif isinstance(alpha, (VolumeData, VertexData)): + alpha = _mask_alpha(alpha, mask) + else: + alpha = np.array(alpha, copy=True) + alpha[mask] = 0 return red, green, blue, alpha @@ -570,21 +642,15 @@ def alpha(self) -> Volume: alpha = np.ones(self.red.volume.shape) alpha = Volume(alpha, self.red.subject, self.red.xfmname, vmin=0, vmax=1) if not isinstance(alpha, Volume): - if alpha.dtype != np.uint8 and (alpha.min() < 0 or alpha.max() > 1): - warnings.warn( - "Some alpha values are outside the range of [0, 1]. " - "Consider passing a Volume object as alpha with explicit vmin, vmax " - "keyword arguments.", - Warning, - ) + _warn_alpha_range(alpha) alpha = Volume(alpha, self.red.subject, self.red.xfmname, vmin=0, vmax=1) + # NaN in any color channel -> alpha at its minimum (transparent) rgb = np.array([self.red.volume, self.green.volume, self.blue.volume]) mask = np.isnan(rgb).any(axis=0) - alpha.volume[mask] = alpha.vmin + alpha = _mask_alpha(alpha, mask) - self._apply_nan_mask(alpha) - return alpha + return self._apply_nan_mask(alpha) @alpha.setter def alpha(self, alpha: Optional[Union[npt.NDArray, Volume]]): @@ -626,6 +692,8 @@ def volume(self) -> np.ndarray[tuple[int, int, int, int, int], np.dtype[np.uint8 else: vol /= dv.vmax - dv.vmin + # NaN (e.g. in a user-supplied alpha map) -> 0, never UB cast + vol = np.nan_to_num(vol, nan=0.0) vol = (np.clip(vol, 0, 1) * 255).astype(np.uint8) else: vol = dv.volume.copy() @@ -846,21 +914,15 @@ def alpha(self) -> Vertex: alpha = np.ones(self.red.vertices.shape[1]) alpha = Vertex(alpha, self.red.subject, vmin=0, vmax=1) if not isinstance(alpha, Vertex): - if alpha.dtype != np.uint8 and (alpha.min() < 0 or alpha.max() > 1): - warnings.warn( - "Some alpha values are outside the range of [0, 1]. " - "Consider passing a Vertex object as alpha with explicit vmin, vmax " - "keyword arguments.", - Warning, - ) + _warn_alpha_range(alpha) alpha = Vertex(alpha, self.red.subject, vmin=0, vmax=1) + # NaN in any color channel -> alpha at its minimum (transparent) rgb = np.array([self.red.data, self.green.data, self.blue.data]) mask = np.isnan(rgb).any(axis=0) - alpha.data[mask] = alpha.vmin + alpha = _mask_alpha(alpha, mask) - self._apply_nan_mask(alpha) - return alpha + return self._apply_nan_mask(alpha) @alpha.setter def alpha(self, alpha: Optional[Union[npt.NDArray, Vertex]]): @@ -887,6 +949,8 @@ def vertices(self) -> npt.NDArray[np.uint8]: else: vert /= dv.vmax - dv.vmin + # NaN (e.g. in a user-supplied alpha map) -> 0, never UB cast + vert = np.nan_to_num(vert, nan=0.0) vert = (np.clip(vert, 0, 1) * 255).astype(np.uint8) else: vert = dv.vertices.copy() diff --git a/cortex/dataset/views.py b/cortex/dataset/views.py index 6402cdea0..8d6c3e5c7 100644 --- a/cortex/dataset/views.py +++ b/cortex/dataset/views.py @@ -132,9 +132,16 @@ def _from_hdf_view( h5, data, xfmname=xfmname, vmin=vmin, vmax=vmax, subject=subject, **kwargs ) - if len(data) == 2: + if len(data) in (2, 3): dim1 = _from_hdf_data(h5, data[0], xfmname=xfmname[0], subject=subject) dim2 = _from_hdf_data(h5, data[1], xfmname=xfmname[1], subject=subject) + # Optional third entry: the alpha map of a 2D view (same xfm as dim1) + alpha = None + if len(data) == 3 and data[2] is not None: + # stored normalized to [0, 1] by Dataview2D._write_hdf + alpha = _from_hdf_data( + h5, data[2], xfmname=xfmname[0], subject=subject, vmin=0, vmax=1 + ) cls = Vertex2D if isinstance(dim1, Vertex) else Volume2D return cls( dim1, @@ -144,6 +151,7 @@ def _from_hdf_view( vmax=vmax[0], vmax2=vmax[1], subject=subject, + alpha=alpha, **kwargs, ) elif len(data) == 4: diff --git a/cortex/quickflat/composite.py b/cortex/quickflat/composite.py index e0e0a14d8..7520579b0 100644 --- a/cortex/quickflat/composite.py +++ b/cortex/quickflat/composite.py @@ -121,7 +121,7 @@ def add_curvature(fig, dataview, extents=None, height=None, threshold=True, cont return cvimg def add_data(fig, braindata, height=1024, thick=32, depth=0.5, pixelwise=True, - sampler='nearest', recache=False, nanmean=False): + sampler='nearest', recache=False, nanmean=True): """Add data to quickflat plot Parameters @@ -143,7 +143,7 @@ def add_data(fig, braindata, height=1024, thick=32, depth=0.5, pixelwise=True, sampler : str Name of sampling function used to sample underlying volume data. Options include 'trilinear','nearest','lanczos'; see functions in cortex.mapper.samplers.py for all options - nanmean : bool, optional (default = False) + nanmean : bool, optional (default = True) If True, NaNs in the data will be ignored when averaging across layers. Returns diff --git a/cortex/quickflat/utils.py b/cortex/quickflat/utils.py index 8bcecafa3..cf0919510 100644 --- a/cortex/quickflat/utils.py +++ b/cortex/quickflat/utils.py @@ -12,7 +12,23 @@ from ..options import config -def make_flatmap_image(braindata, height=1024, recache=False, nanmean=False, **kwargs): +def _rgba_nan_mask(raw, shape): + """Boolean NaN mask (or None) stored by a raw conversion (``Volume.raw``, + ``Volume2D.raw``, ...), in the layout of the uint8 RGBA representation + ``shape`` (time axis included). None for native RGB dataviews: there NaN + has become alpha 0 and is indistinguishable from intentional transparency, + exactly as in the WebGL viewer's RGB textures. + """ + stored = getattr(raw, "_nan_mask", None) + if stored is None: + return None + stored = np.asarray(stored, dtype=bool) + if stored.size != int(np.prod(shape)): + return None + return stored.reshape(shape) + + +def make_flatmap_image(braindata, height=1024, recache=False, nanmean=True, **kwargs): """Generate flatmap image from volumetric brain data This @@ -27,8 +43,14 @@ def make_flatmap_image(braindata, height=1024, recache=False, nanmean=False, **k recache : boolean Whether or not to recache intermediate files. Takes longer to plot this way, potentially resolves some errors. Useful if you've made changes to the alignment. - nanmean : bool, optional (default = False) - If True, NaNs in the data will be ignored when averaging across layers. + nanmean : bool, optional (default = True) + If True, NaN voxels are ignored when averaging across cortical + thickness (mean([1, NaN]) = 1); if False, any NaN voxel contributing + to a pixel makes it NaN / transparent. For 2D dataviews the NaN mask + of the conversion to RGBA is used. For RGB dataviews NaN has already + become alpha 0, so fully transparent voxels count as missing: they + are skipped by ``nanmean=True`` and alpha-weighted otherwise (the + same rule as the WebGL viewer's RGB textures). kwargs : idk idk @@ -49,8 +71,10 @@ def make_flatmap_image(braindata, height=1024, recache=False, nanmean=False, **k **kwargs) if isinstance(braindata, dataset.Vertex2D): - data = braindata.raw.vertices + raw = braindata.raw + data = raw.vertices else: + raw = braindata data = braindata.vertices else: pixmap = get_flatcache(braindata.subject, @@ -59,8 +83,10 @@ def make_flatmap_image(braindata, height=1024, recache=False, nanmean=False, **k recache=recache, **kwargs) if isinstance(braindata, dataset.Volume2D): - data = braindata.raw.volume + raw = braindata.raw + data = raw.volume else: + raw = braindata data = braindata.volume if data.shape[0] > 1: @@ -70,8 +96,43 @@ def make_flatmap_image(braindata, height=1024, recache=False, nanmean=False, **k # Convert data to float to avoid image artifacts data = data.astype(float) if data.dtype == np.uint8: + # RGBA data. Average across cortical thickness in *premultiplied* + # space, as the WebGL viewer does (textures are uploaded with + # premultiplyAlpha=true), then un-premultiply. Averaging straight RGBA + # lets transparent voxels (NaN, alpha=0) bleed their (black) color + # into neighbouring pixels, producing dark halos in quickflat only. + rgba = data.reshape(-1, 4).astype(np.float64) / 255. + alpha = rgba[:, 3:4] + premult = np.concatenate([rgba[:, :3] * alpha, alpha], axis=1) + avg = np.asarray(pixmap.dot(premult)) + + # NaN handling (``nanmean``), as in the float branch below. NaN has + # already become alpha 0 in the RGBA conversion, so the voxel validity + # comes from the NaN mask when the conversion provides one (2D views, + # Volume.raw) and otherwise -- native RGB, like the WebGL RGB + # textures -- from "fully transparent". + nan_mask = _rgba_nan_mask(raw, data.shape[:-1]) + if nan_mask is not None: + valid = ~nan_mask.ravel() + else: + valid = alpha[:, 0] > 0 + w_valid = np.asarray(pixmap.dot(valid.astype(np.float64))).ravel() + if nanmean: + # mean over the valid voxels only: mean([c, NaN]) = c + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + avg = np.where(w_valid[:, None] > 0, avg / w_valid[:, None], 0.) + elif nan_mask is not None: + # any NaN voxel contributing to the pixel hides it + w_nan = np.asarray(pixmap.dot((~valid).astype(np.float64))).ravel() + avg[w_nan > 0] = 0. + + out = np.zeros_like(avg) + opaque = avg[:, 3] > 0 + out[opaque, :3] = avg[opaque, :3] / avg[opaque, 3:4] + out[:, 3] = avg[:, 3] img = np.zeros(mask.shape+(4,), dtype=np.uint8) - img[mask] = pixmap * data.reshape(-1, 4) + img[mask] = np.round(np.clip(out, 0, 1) * 255).astype(np.uint8) img = img.transpose(1,0,2)[::-1] # Make img a c-contiguous array or pil will complain when saving it if not img.flags["C_CONTIGUOUS"]: diff --git a/cortex/quickflat/view.py b/cortex/quickflat/view.py index 1d7102905..0a1c84ae7 100644 --- a/cortex/quickflat/view.py +++ b/cortex/quickflat/view.py @@ -42,7 +42,7 @@ def make_figure(braindata: dataset.Dataview, recache: bool=False, pixelwise: boo labelsize: Optional[str]=None, labelcolor: Optional[ColorType]=None, cutout: Optional[str]=None, curvature_brightness: Optional[float]=None, curvature_contrast: Optional[float]=None, curvature_threshold: Optional[bool]=None, fig: Optional[Union[Figure, Axes]]=None, extra_hatch: Optional[tuple[dataset.Dataview, tuple[float, float, float]]]=None, colorbar_ticks: Optional[npt.ArrayLike]=None, colorbar_location: Union[tuple[float, float, float, float], str]='center', roi_list: Optional[list[str]]=None, sulci_list: Optional[list[str]]=None, - nanmean: bool=False) -> Figure: + nanmean: bool=True) -> Figure: """Show a Volume or Vertex on a flatmap with matplotlib. Parameters @@ -121,7 +121,7 @@ def make_figure(braindata: dataset.Dataview, recache: bool=False, pixelwise: boo vmin, vmax specified in the Volume2D object. fig : figure or ax figure into which to plot flatmap - nanmean : bool, optional (default = False) + nanmean : bool, optional (default = True) If True, NaNs in the data will be ignored when averaging across layers. """ from matplotlib import pyplot as plt @@ -343,8 +343,11 @@ def make_svg(fname, braindata, with_labels=False, with_curvature=True, layers=[' ## Render PNG file & retrieve image data arr, extents = make_flatmap_image(braindata, height=height, **kwargs) # Set nans to alpha = 0. to enable transparency when saving as PNG - mask_nans = np.isnan(arr[..., 3]) - arr[mask_nans, 3] = 0. + if arr.ndim == 3 and np.issubdtype(arr.dtype, np.floating): + # RGBA image: NaN alpha -> fully transparent. (2-D scalar images have + # no alpha channel; matplotlib's bad color handles their NaNs.) + mask_nans = np.isnan(arr[..., 3]) + arr[mask_nans, 3] = 0. if hasattr(braindata, 'cmap'): imsave(fp, arr, cmap=braindata.cmap, vmin=braindata.vmin, vmax=braindata.vmax) diff --git a/cortex/tests/test_nan_alpha.py b/cortex/tests/test_nan_alpha.py new file mode 100644 index 000000000..b0fd404ea --- /dev/null +++ b/cortex/tests/test_nan_alpha.py @@ -0,0 +1,434 @@ +"""NaN / alpha handling must be identical across quickflat, the WebGL data +package and the RGB conversions (browser-free tests). + +Rule: a NaN anywhere at a voxel/vertex -- in either dimension of a 2D view, in +any RGB channel, or in the alpha map itself -- renders as alpha 0. Where there +is no NaN, the alpha (2D ``*_alpha`` colormap, ``alpha=`` kwarg, RGB alpha +channel) is honored. +""" +import json +import os +import tempfile +import warnings + +import numpy as np +import pytest + +import cortex +from cortex import dataset +from cortex.testing_utils import has_installed +from cortex.webgl.data import Package +from cortex.webgl.serve import NPEncode + +subj, xfmname, volshape = "S1", "fullhead", (31, 100, 100) +no_inkscape = not has_installed("inkscape") + + +def _nverts(): + return cortex.db.get_surf(subj, "fiducial", merge=True)[0].shape[0] + + +def _vol_grid(): + zz, yy, xx = np.mgrid[0 : volshape[0], 0 : volshape[1], 0 : volshape[2]] + return zz, yy, xx + + +def _make_2d(kind, d1, d2, alpha=None, **kw): + kw.setdefault("cmap", "RdBu_r_alpha") + kw.update(vmin=-1, vmax=1, vmin2=0, vmax2=1) + if kind == "Volume2D": + return cortex.Volume2D(d1, d2, subj, xfmname, alpha=alpha, **kw) + return cortex.Vertex2D(d1, d2, subj, alpha=alpha, **kw) + + +def _rgba(view): + """uint8 RGBA array of the quickflat/raw representation, time axis dropped.""" + if isinstance(view, (cortex.Volume2D, cortex.VolumeRGB)): + arr = view.volume if isinstance(view, cortex.VolumeRGB) else view.raw.volume + else: + arr = view.vertices if isinstance(view, cortex.VertexRGB) else view.raw.vertices + return arr[0] + + +# --------------------------------------------------------------------------- +# Volume2D / Vertex2D: NaN in either dim or in alpha -> 0; alpha= honored +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("kind", ["Volume2D", "Vertex2D"]) +@pytest.mark.parametrize("nan_in", ["dim1", "dim2", "alpha"]) +def test_2d_nan_anywhere_is_transparent(kind, nan_in): + rng = np.random.default_rng(0) + if kind == "Volume2D": + shape = volshape + region = _vol_grid()[1] < 35 + else: + shape = (_nverts(),) + region = np.arange(shape[0]) < shape[0] // 3 + d1 = rng.uniform(-1, 1, shape) + d2 = np.ones(shape) + alpha = rng.uniform(0.2, 1, shape) + {"dim1": d1, "dim2": d2, "alpha": alpha}[nan_in][region] = np.nan + + rgba = _rgba(_make_2d(kind, d1, d2, alpha=alpha)) + assert rgba[region][..., 3].max() == 0, "NaN must give alpha 0" + assert rgba[~region][..., 3].min() > 0, "non-NaN must keep some alpha" + + +@pytest.mark.parametrize("kind", ["Volume2D", "Vertex2D"]) +def test_2d_alpha_kwarg_multiplies_colormap_alpha(kind): + rng = np.random.default_rng(1) + shape = volshape if kind == "Volume2D" else (_nverts(),) + d1 = rng.uniform(-1, 1, shape) + d2 = rng.uniform(0, 1, shape) + user = rng.uniform(0, 1, shape) + + a_cmap = _rgba(_make_2d(kind, d1, d2)).astype(float)[..., 3] + a_both = _rgba(_make_2d(kind, d1, d2, alpha=user))[..., 3] + expected = np.round(a_cmap * user).astype(np.uint8) + np.testing.assert_array_equal(a_both, expected) + + # alpha=1 everywhere is a no-op + a_one = _rgba(_make_2d(kind, d1, d2, alpha=np.ones(shape)))[..., 3] + np.testing.assert_array_equal(a_one, a_cmap.astype(np.uint8)) + + +def test_2d_alpha_accepts_volume_with_own_range(): + rng = np.random.default_rng(2) + d1 = rng.uniform(-1, 1, volshape) + d2 = np.ones(volshape) + acc = rng.uniform(0, 10, volshape) # e.g. a "confidence" in [0, 10] + v = _make_2d( + "Volume2D", d1, d2, alpha=cortex.Volume(acc, subj, xfmname, vmin=0, vmax=10) + ) + a = _rgba(v)[..., 3].astype(float) + a_cmap = _rgba(_make_2d("Volume2D", d1, d2)).astype(float)[..., 3] + np.testing.assert_array_equal(a, np.round(a_cmap * acc / 10.0)) + + +def test_2d_alpha_not_in_attrs_and_json_serializable(): + """The alpha map used to be stuffed into ``attrs`` as an ndarray, which + crashed the WebGL viewer (500: ndarray is not JSON serializable).""" + rng = np.random.default_rng(3) + v = _make_2d( + "Volume2D", + rng.uniform(-1, 1, volshape), + np.ones(volshape), + alpha=rng.uniform(0, 1, volshape), + ) + assert not any(isinstance(x, np.ndarray) for x in v.attrs.values()) + assert isinstance(v.alpha, cortex.Volume) + js = v.to_json() + json.dumps(js, cls=NPEncode) + json.dumps(js) # plain encoder, as used by the mixer.html handler + assert js["alpha"] == [v.alpha.name] + + v_noalpha = _make_2d("Volume2D", rng.uniform(-1, 1, volshape), np.ones(volshape)) + assert "alpha" not in v_noalpha.to_json() + + +def test_2d_alpha_rejects_other_subject_or_xfm(): + rng = np.random.default_rng(4) + d1 = rng.uniform(-1, 1, volshape) + bad = cortex.Volume(rng.uniform(0, 1, volshape), subj, xfmname) + bad.subject = "not_S1" + with pytest.raises(ValueError): + _make_2d("Volume2D", d1, np.ones(volshape), alpha=bad) + + +def test_package_ships_2d_alpha_as_float_brain(): + rng = np.random.default_rng(5) + alpha = rng.uniform(0, 1, volshape) + alpha[0] = np.nan + v = _make_2d("Volume2D", rng.uniform(-1, 1, volshape), np.ones(volshape), alpha=alpha) + pkg = Package(dataset.Dataset(v=v)) + assert len(pkg.brains) == 3 + assert all(b["raw"] is False for b in pkg.brains.values()) + meta = pkg.metadata() + assert meta["views"][0]["alpha"] == [v.alpha.name] + assert v.alpha.name in meta["images"] + + +def test_2d_to_json_keeps_zero_bounds(): + """``vmin=0``/``vmax=0`` must not fall back to the auto range + (truthiness bug; same class as 5482c8bf).""" + rng = np.random.default_rng(6) + dim1 = cortex.Volume(rng.uniform(-1, 1, volshape), subj, xfmname) # auto range + dim2 = cortex.Volume(rng.uniform(-1, 1, volshape), subj, xfmname) + v = cortex.Volume2D(dim1, dim2, vmin=0, vmax=1, vmin2=-1, vmax2=0) + js = v.to_json() + assert js["vmin"][0] == [0, -1] + assert js["vmax"][0] == [1, 0] + + +def test_2d_alpha_hdf_roundtrip(): + rng = np.random.default_rng(7) + alpha = rng.uniform(0, 1, volshape) + v = _make_2d("Volume2D", rng.uniform(-1, 1, volshape), np.ones(volshape), alpha=alpha) + tf = tempfile.NamedTemporaryFile(suffix=".hdf", delete=False) + tf.close() + os.unlink(tf.name) + try: + dataset.Dataset(twod=v).save(tf.name) + loaded = cortex.load(tf.name) + assert isinstance(loaded.twod.alpha, cortex.Volume) + np.testing.assert_allclose(loaded.twod.alpha.data, alpha, atol=1e-6) + np.testing.assert_array_equal(_rgba(loaded.twod), _rgba(v)) + finally: + if os.path.exists(tf.name): + os.unlink(tf.name) + + +# --------------------------------------------------------------------------- +# VolumeRGB / VertexRGB +# --------------------------------------------------------------------------- + + +def test_volumergb_masked_alpha_nan_channel_is_transparent(): + """A masked (linear) alpha Volume used to be written through a temporary + (``alpha.volume[mask] = vmin``), so NaN voxels stayed opaque.""" + rng = np.random.default_rng(8) + mask = cortex.db.get_mask(subj, xfmname, "thick") + zz, yy, xx = _vol_grid() + r = rng.uniform(0, 1, volshape) + r[yy < 35] = np.nan + alpha_lin = cortex.Volume(np.full(mask.sum(), 0.8), subj, xfmname, vmin=0, vmax=1) + rgb = cortex.VolumeRGB( + cortex.Volume(r, subj, xfmname, vmin=0, vmax=1), + cortex.Volume(rng.uniform(0, 1, volshape), subj, xfmname, vmin=0, vmax=1), + cortex.Volume(rng.uniform(0, 1, volshape), subj, xfmname, vmin=0, vmax=1), + subj, + xfmname, + alpha=alpha_lin, + ) + a = _rgba(rgb)[..., 3] + assert a[np.isnan(r) & mask].max() == 0 + assert a[~np.isnan(r) & mask].min() == 204 # round(0.8 * 255) + # the user's alpha object is untouched + assert np.all(alpha_lin.data == 0.8) + + +@pytest.mark.parametrize("cls", ["VertexRGB", "VolumeRGB"]) +def test_rgb_multiframe_nan_masks_per_frame(cls): + """Regression for #629: multi-frame data + NaN raised IndexError because the + auto alpha was single-frame while the NaN mask was (T, ...).""" + rng = np.random.default_rng(9) + T = 3 + shape = (T, _nverts()) if cls == "VertexRGB" else (T,) + volshape + r, g, b = (rng.uniform(0, 1, shape) for _ in range(3)) + idx = (slice(None),) + (0,) * (len(shape) - 2) + (slice(0, 50),) + r[(1,) + idx[1:]] = np.nan # frame 1 only + if cls == "VertexRGB": + rgb = cortex.VertexRGB(r, g, b, subj) + arr = rgb.vertices + else: + rgb = cortex.VolumeRGB(r, g, b, subj, xfmname) + arr = rgb.volume + assert arr.shape[0] == T + nan_here = np.isnan(r) + assert arr[..., 3][nan_here].max() == 0 + assert arr[..., 3][~nan_here].min() > 0 + # a user-supplied single-frame alpha is broadcast, not rejected + if cls == "VertexRGB": + rgb.alpha = np.full(shape[1:], 0.5) + arr = rgb.vertices + else: + rgb.alpha = np.full(shape[1:], 0.5) + arr = rgb.volume + assert arr[..., 3][nan_here].max() == 0 + assert arr[..., 3][~nan_here].min() == 127 # int(0.5 * 255) + + +@pytest.mark.parametrize("cls", ["VertexRGB", "VolumeRGB"]) +def test_rgb_nan_in_alpha_is_transparent(cls): + rng = np.random.default_rng(10) + shape = (_nverts(),) if cls == "VertexRGB" else volshape + alpha = rng.uniform(0.5, 1, shape) + region = np.arange(shape[0]) < shape[0] // 2 + alpha[region] = np.nan + r, g, b = (rng.uniform(0, 1, shape) for _ in range(3)) + with warnings.catch_warnings(): + warnings.simplefilter("error") # no "invalid value in cast" warnings + if cls == "VertexRGB": + arr = cortex.VertexRGB(r, g, b, subj, alpha=alpha).vertices[0] + else: + arr = cortex.VolumeRGB(r, g, b, subj, xfmname, alpha=alpha).volume[0] + assert arr[region][..., 3].max() == 0 + assert arr[~region][..., 3].min() > 0 + + +def test_color_voxels_does_not_mutate_caller_alpha(): + rng = np.random.default_rng(11) + r = rng.uniform(0, 1, volshape) + r[0] = np.nan + alpha = np.ones(volshape) + rgb = cortex.VolumeRGB( + r, rng.uniform(0, 1, volshape), rng.uniform(0, 1, volshape), subj, xfmname, + vmin=0, vmax=1, alpha=alpha, # vmin/vmax -> color_voxels path + ) + assert np.all(alpha == 1.0) + a = rgb.volume[0][..., 3] + assert a[0].max() == 0 + assert a[1:].min() == 255 + + +# --------------------------------------------------------------------------- +# quickflat +# --------------------------------------------------------------------------- + + +def test_make_flatmap_image_rgb_averages_premultiplied(): + """Transparent (alpha 0 / NaN) voxels must not darken neighbouring pixels. + + In alpha-weighted mode (``nanmean=False`` for RGB data) quickflat averages + RGBA over the voxels of a pixel; the WebGL viewer does this in + premultiplied space. Averaging straight RGBA gave dark halos around + transparent regions in quickflat only. Opaque voxels are pure red, the + others are NaN in the red channel (transparent, and black after the + uint8 conversion), so straight averaging would give R ~ 128. + """ + zz, yy, xx = _vol_grid() + checker = ((xx + yy + zz) % 2).astype(bool) + r = np.where(checker, 1.0, np.nan).astype(np.float32) + red = cortex.VolumeRGB( + cortex.Volume(r, subj, xfmname, vmin=0, vmax=1), + cortex.Volume(np.zeros(volshape, np.float32), subj, xfmname, vmin=0, vmax=1), + cortex.Volume(np.zeros(volshape, np.float32), subj, xfmname, vmin=0, vmax=1), + subj, + xfmname, + ) + img, _ = cortex.quickflat.utils.make_flatmap_image(red, nanmean=False) + a = img[..., 3] + partial = (a > 0) & (a < 255) # pixels mixing transparent and opaque voxels + assert partial.sum() > 100 + # Color must stay pure bright red regardless of partial coverage + assert img[partial][:, 0].min() >= 250 + assert img[partial][:, 1:3].max() <= 5 + # and fully opaque pixels are unchanged + if (a == 255).any(): + assert img[a == 255][:, 0].min() == 255 + + +def test_make_flatmap_image_volume_nan_transparent_when_masked(): + mask = cortex.db.get_mask(subj, xfmname, "thick") + data = np.ones(mask.sum()) + data[: data.size // 2] = np.nan + vol = cortex.Volume(data, subj, xfmname, vmin=0, vmax=1) + img, _ = cortex.quickflat.utils.make_flatmap_image(vol, nanmean=True) + assert np.nanmin(img) == 1 + assert np.isnan(img).any() + + +@pytest.mark.skipif(no_inkscape, reason="Inkscape required") +def test_make_svg_scalar_with_nan(): + """``make_svg`` indexed ``arr[..., 3]`` on a 2-D scalar image.""" + data = np.random.default_rng(12).uniform(0, 1, volshape) + data[:, :35] = np.nan + vol = cortex.Volume(data, subj, xfmname, vmin=0, vmax=1, cmap="viridis") + tf = tempfile.NamedTemporaryFile(suffix=".svg", delete=False) + tf.close() + try: + cortex.quickflat.make_svg(tf.name, vol, with_labels=False) + assert os.path.getsize(tf.name) > 0 + finally: + os.unlink(tf.name) + + +def test_quickflat_nanmean_is_default(): + """quickflat ignores NaN voxels when averaging across thickness by + default, like the WebGL viewer's ``nanmean`` surface toggle.""" + import inspect + + for func in ( + cortex.quickflat.make_figure, + cortex.quickflat.utils.make_flatmap_image, + cortex.quickflat.composite.add_data, + ): + assert inspect.signature(func).parameters["nanmean"].default is True, func + + +def test_rgb_uint8_alpha_nan_channel_is_transparent(): + """A uint8 alpha map bypasses normalization and its inferred vmin is a + percentile of the bytes (255 for a constant map); NaN must still give 0.""" + rng = np.random.default_rng(13) + r = rng.uniform(0, 1, volshape) + r[:, :35] = np.nan + alpha = cortex.Volume(np.full(volshape, 255, np.uint8), subj, xfmname) + rgb = cortex.VolumeRGB( + cortex.Volume(r, subj, xfmname, vmin=0, vmax=1), + cortex.Volume(rng.uniform(0, 1, volshape), subj, xfmname, vmin=0, vmax=1), + cortex.Volume(rng.uniform(0, 1, volshape), subj, xfmname, vmin=0, vmax=1), + subj, xfmname, alpha=alpha, + ) + a = rgb.volume[0][..., 3] + assert a[np.isnan(r)].max() == 0 + assert a[~np.isnan(r)].min() == 255 + + +def _alpha_stats(view, nanmean): + img, _ = cortex.quickflat.utils.make_flatmap_image(view, nanmean=nanmean) + a = img[..., 3] + return a + + +@pytest.mark.parametrize("kind", ["Volume2D", "VolumeRGB"]) +def test_quickflat_nanmean_applies_to_rgba_dataviews(kind): + """quickflat's nanmean must also act on dataviews that reach the RGBA + (uint8) branch. 2D views carry an exact NaN mask (nanmean=False hides + any pixel touched by a NaN voxel); for RGB views NaN has become alpha 0, + so nanmean=False is the alpha-weighted average, like WebGL RGB textures. + """ + zz, yy, xx = _vol_grid() + scattered = (xx + yy + zz) % 3 == 0 # a third of the voxels, everywhere + ones = np.ones(volshape) + d = ones.copy() + d[scattered] = np.nan + kw2d = dict(cmap="RdBu_covar", vmin=-1, vmax=1, vmin2=0, vmax2=1) + if kind == "Volume2D": + clean = cortex.Volume2D(ones, ones, subj, xfmname, **kw2d) + view = cortex.Volume2D(d, ones, subj, xfmname, **kw2d) + else: + zeros = np.zeros(volshape) + clean = cortex.VolumeRGB( + cortex.Volume(ones, subj, xfmname, vmin=0, vmax=1), + cortex.Volume(zeros, subj, xfmname, vmin=0, vmax=1), + cortex.Volume(zeros, subj, xfmname, vmin=0, vmax=1), subj, xfmname) + view = cortex.VolumeRGB( + cortex.Volume(d, subj, xfmname, vmin=0, vmax=1), + cortex.Volume(zeros, subj, xfmname, vmin=0, vmax=1), + cortex.Volume(zeros, subj, xfmname, vmin=0, vmax=1), subj, xfmname) + + a_clean = _alpha_stats(clean, True) + brain = a_clean == 255 # pixels fully covered by (valid) cortex + assert brain.sum() > 10000 + a_mean = _alpha_stats(view, True) + a_strict = _alpha_stats(view, False) + # nanmean: the NaN voxels are ignored -> (almost) every brain pixel is + # still fully opaque, exactly like the NaN-free data + assert (a_mean[brain] == 255).mean() > 0.95 + if kind == "Volume2D": + # any NaN voxel hides the pixel -> most pixels disappear, no middle ground + assert (a_strict[brain] == 0).mean() > 0.5 + assert not ((a_strict[brain] > 0) & (a_strict[brain] < 255)).any() + else: + # alpha-weighted average: roughly a third of the opacity is lost + assert 0.4 < a_strict[brain].mean() / 255. < 0.9 + assert ((a_strict[brain] > 0) & (a_strict[brain] < 255)).mean() > 0.5 + + +def test_quickflat_nanmean_native_rgb_without_nan_is_unchanged(): + """For RGB data with intentional alpha (no NaN), nanmean=False keeps the + alpha-weighted average; nanmean=True skips fully transparent voxels.""" + zz, yy, xx = _vol_grid() + checker = ((xx + yy + zz) % 2).astype(np.float32) + red = cortex.VolumeRGB( + cortex.Volume(np.ones(volshape, np.float32), subj, xfmname, vmin=0, vmax=1), + cortex.Volume(np.zeros(volshape, np.float32), subj, xfmname, vmin=0, vmax=1), + cortex.Volume(np.zeros(volshape, np.float32), subj, xfmname, vmin=0, vmax=1), + subj, xfmname, alpha=cortex.Volume(checker, subj, xfmname, vmin=0, vmax=1), + ) + a_weighted = _alpha_stats(red, False) + a_mean = _alpha_stats(red, True) + partial = (a_weighted > 0) & (a_weighted < 255) + assert partial.sum() > 1000 + assert (a_mean[partial] == 255).mean() > 0.95 diff --git a/cortex/tests/test_webgl_data.py b/cortex/tests/test_webgl_data.py index 35302e1af..65fbf3426 100644 --- a/cortex/tests/test_webgl_data.py +++ b/cortex/tests/test_webgl_data.py @@ -192,3 +192,48 @@ def spy_mosaic(arr, show=False): ) > 5 ), "VolumeRGB Package output looks premultiplied; Three.js will then double-attenuate" + + +def test_package_deduplicates_identical_brains(): + """Two dims/channels with byte-identical data share one content-hash name; + the package must contain that brain once, and ``reorder`` must not choke + on it (it used to re-index the already-serialized bytes).""" + import numpy as np + import cortex + from cortex import utils + from cortex.webgl.data import Package + + subj = "S1" + nverts = cortex.db.get_surf(subj, "fiducial", merge=True)[0].shape[0] + x = np.random.default_rng(0).standard_normal(nverts) + twod = cortex.Vertex2D(x, x, subj) + rgb = cortex.VertexRGB(x, x, x, subj) + pkg = Package(cortex.Dataset(twod=twod, rgb=rgb)) + names = [b.name for b in pkg.uniques] + assert len(names) == len(set(names)) == 2 + + # reorder with the same ctm pack the viewer uses (cortex.webgl.make_static) + ctm = utils.get_ctmpack( + subj, ("inflated",), method="mg2", level=9, recache=False, + external_svg=None, overlays_available=None, + ) + pkg.reorder({subj: ctm}) + for name in names: + assert pkg.images[name][0][1:6] == b"NUMPY" + + +def test_package_rejects_same_bytes_different_metadata(): + """Same data bytes give the same content-hash name; if the metadata + differs the package cannot represent both and must say so.""" + import numpy as np + import pytest + import cortex + from cortex.webgl.data import Package + + nverts = cortex.db.get_surf("S1", "fiducial", merge=True)[0].shape[0] + x = np.zeros(nverts) + a = cortex.Vertex(x, "S1") + b = cortex.Vertex(x, "S1") + b.subject = "S2" # never reaches the database: the check comes first + with pytest.raises(ValueError, match="identical bytes"): + Package(cortex.Dataset(a=a, b=b)) diff --git a/cortex/tests/test_webgl_nan_alpha_parity.py b/cortex/tests/test_webgl_nan_alpha_parity.py new file mode 100644 index 000000000..46428d3ae --- /dev/null +++ b/cortex/tests/test_webgl_nan_alpha_parity.py @@ -0,0 +1,316 @@ +"""quickflat (matplotlib) and the WebGL viewer must agree on NaN and alpha. + +For a set of dataviews covering every NaN / alpha pattern, render the flatmap +with ``cortex.quickshow`` and with the headless WebGL viewer and compare the +fraction of the brain that shows data (red-dominant pixels over all non-white +pixels). Pixel-exact comparison is not possible (different rasterizers, +resolutions, curvature rendering), but the fraction of transparent cortex is, +within a tolerance. + +Skipped if playwright is not installed. +""" +import os +import time + +import numpy as np +import pytest + +import cortex +import cortex.export +from cortex.export.save_views import ( + angle_view_params, + default_view_params, + unfold_view_params, +) +from cortex.tests.testing_utils import has_playwright + +pytestmark = pytest.mark.skipif( + not has_playwright, reason="playwright and chromium are required" +) + +subj, xfmname, volshape = "S1", "fullhead", (31, 100, 100) +FLAT = { + **default_view_params, + **angle_view_params["flatmap"], + **unfold_view_params["flatmap"], +} +TOL = 0.12 # absolute tolerance on the visible-data fraction + + +def _fractions(path): + """(fraction of brain pixels showing data, number of brain pixels).""" + from PIL import Image + + im = Image.open(path).convert("RGBA") + bg = Image.new("RGBA", im.size, (255, 255, 255, 255)) + rgb = np.asarray(Image.alpha_composite(bg, im).convert("RGB")).astype(int) + brain = rgb.min(axis=2) < 235 # anything not (near) white + red = (rgb[..., 0] - np.maximum(rgb[..., 1], rgb[..., 2])) > 40 + n_brain = int(brain.sum()) + return (int((red & brain).sum()) / max(n_brain, 1)), n_brain + + +def _quickshow(view, path): + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + fig = cortex.quickshow( + view, with_curvature=True, with_rois=False, with_labels=False, + with_colorbar=False, with_sulci=False, with_borders=False, height=256, + ) + fig.savefig(path, bbox_inches="tight", pad_inches=0, dpi=80) + plt.close(fig) + return _fractions(path) + + +def _webgl(view, path): + with cortex.export.headless_viewer( + view, viewer_params=dict(labels_visible=[], overlays_visible=[]) + ) as handle: + handle._set_view(**FLAT) + time.sleep(3) + handle.getImage(path, (1024, 768)) + for _ in range(300): + if os.path.exists(path) and os.path.getsize(path) > 0: + break + time.sleep(0.1) + time.sleep(0.3) + errors = [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e] + assert not errors, errors + return _fractions(path) + + +def _cases(): + zz, yy, xx = np.mgrid[0 : volshape[0], 0 : volshape[1], 0 : volshape[2]] + post = yy < 35 # posterior slab + pts = cortex.db.get_surf(subj, "fiducial", merge=True)[0] + nv = pts.shape[0] + vpost = pts[:, 1] < np.percentile(pts[:, 1], 35) + mask = cortex.db.get_mask(subj, xfmname, "thick") + + def V(d, **kw): + return cortex.Volume(d, subj, xfmname, **kw) + + def X(d, **kw): + return cortex.Vertex(d, subj, **kw) + + ones_v, ones_x = np.ones(volshape), np.ones(nv) + nan_v = ones_v.copy(); nan_v[post] = np.nan + nan_x = ones_x.copy(); nan_x[vpost] = np.nan + a0_v = (~post).astype(float) + a0_x = (~vpost).astype(float) + kw2d = dict(cmap="RdBu_r_alpha", vmin=-1, vmax=1, vmin2=0, vmax2=1) + zeros_v, zeros_x = np.zeros(volshape), np.zeros(nv) + + return { + "volume_nan": V(nan_v * 5, cmap="Reds", vmin=0, vmax=1), + "vertex_nan": X(nan_x * 5, cmap="Reds", vmin=0, vmax=1), + "volume_masked_nan": V((nan_v * 5)[mask], cmap="Reds", vmin=0, vmax=1), + "volume2d_nan_dim1": cortex.Volume2D(nan_v, ones_v, subj, xfmname, **kw2d), + "volume2d_nan_dim2": cortex.Volume2D(ones_v, nan_v, subj, xfmname, **kw2d), + "vertex2d_nan_dim2": cortex.Vertex2D(ones_x, nan_x, subj, **kw2d), + "volume2d_alpha0": cortex.Volume2D(ones_v, ones_v, subj, xfmname, alpha=a0_v, **kw2d), + "vertex2d_alpha0": cortex.Vertex2D(ones_x, ones_x, subj, alpha=a0_x, **kw2d), + "vertex2d_alpha_nan": cortex.Vertex2D( + ones_x, ones_x, subj, alpha=np.where(vpost, np.nan, 1.0), **kw2d + ), + "volumergb_nan_channel": cortex.VolumeRGB( + V(nan_v, vmin=0, vmax=1), V(zeros_v, vmin=0, vmax=1), + V(zeros_v, vmin=0, vmax=1), subj, xfmname, + ), + "volumergb_masked_alpha": cortex.VolumeRGB( + V(nan_v, vmin=0, vmax=1), V(zeros_v, vmin=0, vmax=1), + V(zeros_v, vmin=0, vmax=1), subj, xfmname, + alpha=V(np.full(mask.sum(), 1.0), vmin=0, vmax=1), + ), + "volumergb_nan_in_alpha": cortex.VolumeRGB( + V(ones_v, vmin=0, vmax=1), V(zeros_v, vmin=0, vmax=1), + V(zeros_v, vmin=0, vmax=1), subj, xfmname, + alpha=np.where(post, np.nan, 1.0), + ), + "volumergb_color_voxels_nan": cortex.VolumeRGB( + nan_v, zeros_v, zeros_v, subj, xfmname, vmin=0, vmax=1, + ), + "vertexrgb_alpha0": cortex.VertexRGB( + X(ones_x, vmin=0, vmax=1), X(zeros_x, vmin=0, vmax=1), + X(zeros_x, vmin=0, vmax=1), subj, alpha=a0_x, + ), + "vertexrgb_nan_in_alpha": cortex.VertexRGB( + X(ones_x, vmin=0, vmax=1), X(zeros_x, vmin=0, vmax=1), + X(zeros_x, vmin=0, vmax=1), subj, alpha=np.where(vpost, np.nan, 1.0), + ), + } + + +@pytest.fixture(scope="module") +def cases(): + return _cases() + + +@pytest.mark.parametrize("name", [ + "volume_nan", "vertex_nan", "volume_masked_nan", + "volume2d_nan_dim1", "volume2d_nan_dim2", "vertex2d_nan_dim2", + "volume2d_alpha0", "vertex2d_alpha0", "vertex2d_alpha_nan", + "volumergb_nan_channel", "volumergb_masked_alpha", "volumergb_nan_in_alpha", + "volumergb_color_voxels_nan", "vertexrgb_alpha0", "vertexrgb_nan_in_alpha", +]) +def test_visible_fraction_matches(name, cases, tmp_path): + view = cases[name] + f_qs, n_qs = _quickshow(view, str(tmp_path / ("qs_%s.png" % name))) + f_wg, n_wg = _webgl(view, str(tmp_path / ("wg_%s.png" % name))) + assert n_qs > 1000 and n_wg > 1000, "brain not found in one of the renders" + # every case hides roughly the posterior third: neither fully shown nor hidden + assert 0.25 < f_qs < 0.9, "quickshow fraction %.2f" % f_qs + assert 0.25 < f_wg < 0.9, "webgl fraction %.2f" % f_wg + assert abs(f_qs - f_wg) < TOL, ( + "%s: quickshow shows data on %.0f%% of the cortex, WebGL on %.0f%%" + % (name, 100 * f_qs, 100 * f_wg) + ) + + +def test_multilayer_nanmean_toggle(tmp_path): + """With several layers, NaN voxels are left out of the average (like + quickflat's ``nanmean=True``) unless the surface's ``nanmean`` toggle is + off, in which case one NaN at any depth makes the fragment transparent.""" + zz, yy, xx = np.mgrid[0 : volshape[0], 0 : volshape[1], 0 : volshape[2]] + d = np.full(volshape, 5.0) + d[(xx + yy + zz) % 3 == 0] = np.nan # a third of the voxels, scattered + vol = cortex.Volume(d, subj, xfmname, cmap="Reds", vmin=0, vmax=1) + view = { + **default_view_params, + **angle_view_params["lateral_pivot"], + **unfold_view_params["inflated"], + } + + def _red(path): + from PIL import Image + + rgb = np.asarray(Image.open(path).convert("RGB")).astype(int) + return int((rgb[..., 0] - np.maximum(rgb[..., 1], rgb[..., 2]) > 50).sum()) + + counts = {} + with cortex.export.headless_viewer( + vol, viewer_params=dict(labels_visible=[], overlays_visible=[]) + ) as handle: + handle._set_view(**view) + time.sleep(2) + for layers, nanmean in [(1, True), (8, True), (8, False)]: + handle.ui.set("surface.%s.layers" % subj, layers) + handle.ui.set("surface.%s.nanmean" % subj, nanmean) + time.sleep(2.5) + path = str(tmp_path / ("layers%d_nanmean%s.png" % (layers, nanmean))) + handle.getImage(path, (512, 384)) + for _ in range(300): + if os.path.exists(path) and os.path.getsize(path) > 0: + break + time.sleep(0.1) + time.sleep(0.3) + counts[(layers, nanmean)] = _red(path) + errors = [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e] + assert not errors, errors + assert counts[(1, True)] > 5000 + # nanmean: averaging over the valid layers shows at least as much cortex + assert counts[(8, True)] >= 0.9 * counts[(1, True)], counts + # toggle off: any NaN among the 8 layers hides the fragment + assert counts[(8, False)] < 0.6 * counts[(8, True)], counts + + +def test_multilayer_nanmean_toggle_rgb(tmp_path): + """Same as above for RGB data: NaN became alpha 0 in the texture, so with + ``nanmean`` fully transparent layer samples are left out of the average.""" + zz, yy, xx = np.mgrid[0 : volshape[0], 0 : volshape[1], 0 : volshape[2]] + r = np.ones(volshape) + r[(xx + yy + zz) % 3 == 0] = np.nan + zeros = np.zeros(volshape) + vol = cortex.VolumeRGB( + cortex.Volume(r, subj, xfmname, vmin=0, vmax=1), + cortex.Volume(zeros, subj, xfmname, vmin=0, vmax=1), + cortex.Volume(zeros, subj, xfmname, vmin=0, vmax=1), subj, xfmname, + ) + view = { + **default_view_params, + **angle_view_params["lateral_pivot"], + **unfold_view_params["inflated"], + } + + def _red(path): + """Total redness: sum of R - max(G, B) over red-dominant pixels, so + that a partially transparent red (alpha-weighted average) scores + lower than an opaque one covering the same pixels.""" + from PIL import Image + + rgb = np.asarray(Image.open(path).convert("RGB")).astype(int) + redness = rgb[..., 0] - np.maximum(rgb[..., 1], rgb[..., 2]) + return int(redness[redness > 50].sum()) + + counts = {} + with cortex.export.headless_viewer( + vol, viewer_params=dict(labels_visible=[], overlays_visible=[]) + ) as handle: + handle._set_view(**view) + time.sleep(2) + for layers, nanmean in [(1, True), (8, True), (8, False)]: + handle.ui.set("surface.%s.layers" % subj, layers) + handle.ui.set("surface.%s.nanmean" % subj, nanmean) + time.sleep(2.5) + path = str(tmp_path / ("rgb_layers%d_nanmean%s.png" % (layers, nanmean))) + handle.getImage(path, (512, 384)) + for _ in range(300): + if os.path.exists(path) and os.path.getsize(path) > 0: + break + time.sleep(0.1) + time.sleep(0.3) + counts[(layers, nanmean)] = _red(path) + errors = [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e] + assert not errors, errors + assert counts[(1, True)] > 5000 * 50 + assert counts[(8, True)] >= 0.9 * counts[(1, True)], counts + # alpha-weighted average: the NaN layers dilute the red (about a third) + assert counts[(8, False)] < 0.85 * counts[(8, True)], counts + + +def test_vertex_movie_nan_in_next_frame_is_transparent(tmp_path): + """Between two frames the vertex shader mixes frame f and f+1. A vertex + that is NaN in f+1 (replaced by 0 in the GPU buffer) must be masked while + interpolating, not fade towards a fake 0.""" + nverts = cortex.db.get_surf(subj, "fiducial", merge=True)[0].shape[0] + nl = cortex.db.get_surf(subj, "fiducial")[0][0].shape[0] + movie = np.full((2, nverts), 5.0) + movie[1, :nl] = np.nan # left hemisphere undefined in frame 1 only + vtx = cortex.Vertex(movie, subj, cmap="Reds", vmin=0, vmax=1) + view = { + **default_view_params, + **angle_view_params["lateral_pivot"], + **unfold_view_params["inflated"], + } + + def _red(path): + from PIL import Image + + rgb = np.asarray(Image.open(path).convert("RGB")).astype(int) + return int((rgb[..., 0] - np.maximum(rgb[..., 1], rgb[..., 2]) > 50).sum()) + + counts = {} + with cortex.export.headless_viewer( + vtx, viewer_params=dict(labels_visible=[], overlays_visible=[]) + ) as handle: + handle._set_view(**view) + time.sleep(2) + for frame in (0.0, 0.5, 1.0): + handle.setFrame(frame) + time.sleep(2) + path = str(tmp_path / ("frame_%.1f.png" % frame)) + handle.getImage(path, (512, 384)) + for _ in range(300): + if os.path.exists(path) and os.path.getsize(path) > 0: + break + time.sleep(0.1) + time.sleep(0.3) + counts[frame] = _red(path) + errors = [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e] + assert not errors, errors + assert counts[1.0] < 0.8 * counts[0.0], counts # left hemisphere hidden in frame 1 + # while blending towards frame 1 the NaN vertices are already masked + assert abs(counts[0.5] - counts[1.0]) <= 0.1 * counts[1.0], counts diff --git a/cortex/tests/test_webgl_switching.py b/cortex/tests/test_webgl_switching.py new file mode 100644 index 000000000..c9e6b5156 --- /dev/null +++ b/cortex/tests/test_webgl_switching.py @@ -0,0 +1,244 @@ +"""NaN and alpha state must not leak between datasets in the WebGL viewer. + +One headless viewer is loaded with several dataviews that differ only in where +they are NaN / transparent. After every ``setData`` switch the rendered image +must match the image obtained when that dataview is shown alone in a fresh +viewer; otherwise a NaN mask, an alpha map or an RGB alpha channel leaked from +the previously displayed dataset. + +All dataviews render *red* where visible (``Reds`` colormap at a constant high +value, red RGB channels, or the red corner of ``RdBu_r_alpha``), so "visible +data" is simply the number of red-dominant pixels. + +All tests are skipped if playwright is not installed. +""" +import time + +import numpy as np +import pytest + +import cortex +import cortex.export +from cortex.export.save_views import ( + angle_view_params, + default_view_params, + unfold_view_params, +) +from cortex.tests.testing_utils import has_playwright + +pytestmark = pytest.mark.skipif( + not has_playwright, reason="playwright and chromium are required" +) + +subj, xfmname, volshape = "S1", "fullhead", (31, 100, 100) +VIEW = { + **default_view_params, + **angle_view_params["lateral_pivot"], + **unfold_view_params["inflated"], +} +VIEWER_PARAMS = dict(labels_visible=[], overlays_visible=[]) +RTOL = 0.05 # relative tolerance on red-pixel counts + + +def _count_red(path): + from PIL import Image + + rgb = np.asarray(Image.open(path).convert("RGB")).astype(int) + return int((rgb[..., 0] - np.maximum(rgb[..., 1], rgb[..., 2]) > 50).sum()) + + +def _render(handle, path, size=(512, 384)): + import os + + handle.getImage(path, size) + for _ in range(300): + if os.path.exists(path) and os.path.getsize(path) > 0: + break + time.sleep(0.1) + else: + raise RuntimeError("image not written: %s" % path) + time.sleep(0.3) + return _count_red(path) + + +def _pageerrors(handle): + return [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e] + + +def make_views(): + """Dataviews that are red where visible; half of them hide one half. + + Every dataview uses distinct data: two views sharing byte-identical data + get the same content-hash name and the package would serve them twice. + """ + pts_left = cortex.db.get_surf(subj, "fiducial")[0][0] + nl = pts_left.shape[0] + nv = cortex.db.get_surf(subj, "fiducial", merge=True)[0].shape[0] + left = np.arange(nv) < nl # left hemisphere vertices + zz, yy, xx = np.mgrid[0 : volshape[0], 0 : volshape[1], 0 : volshape[2]] + half_vox = xx < volshape[2] // 2 + + def V(data, **kw): + return cortex.Volume(data, subj, xfmname, **kw) + + def X(data, **kw): + return cortex.Vertex(data, subj, **kw) + + views = {} + # scalar vertex / volume (Reds, constant high value -> pure red) + d = np.full(nv, 5.0); d[left] = np.nan + views["vtx_nan"] = X(d, cmap="Reds", vmin=0, vmax=1) + views["vtx_full"] = X(np.full(nv, 4.9), cmap="Reds", vmin=0, vmax=1) + d = np.full(volshape, 5.0); d[half_vox] = np.nan + views["vol_nan"] = V(d, cmap="Reds", vmin=0, vmax=1) + views["vol_full"] = V(np.full(volshape, 4.9), cmap="Reds", vmin=0, vmax=1) + # RGB with alpha channel + views["vtxrgb_a0"] = cortex.VertexRGB( + X(np.ones(nv), vmin=0, vmax=1), X(np.zeros(nv), vmin=0, vmax=1), + X(np.zeros(nv), vmin=0, vmax=1), subj, + alpha=X((~left).astype(float), vmin=0, vmax=1), + ) + views["vtxrgb_full"] = cortex.VertexRGB( + X(np.full(nv, 0.99), vmin=0, vmax=1), X(np.zeros(nv), vmin=0, vmax=1), + X(np.zeros(nv), vmin=0, vmax=1), subj, + ) + views["volrgb_a0"] = cortex.VolumeRGB( + V(np.ones(volshape), vmin=0, vmax=1), V(np.zeros(volshape), vmin=0, vmax=1), + V(np.zeros(volshape), vmin=0, vmax=1), subj, xfmname, + alpha=V((~half_vox).astype(float), vmin=0, vmax=1), + ) + views["volrgb_full"] = cortex.VolumeRGB( + V(np.full(volshape, 0.99), vmin=0, vmax=1), V(np.zeros(volshape), vmin=0, vmax=1), + V(np.zeros(volshape), vmin=0, vmax=1), subj, xfmname, + ) + # 2D views: dim1 = 1 (red corner of RdBu_r_alpha), dim2 = 1 (opaque) + kw2d = dict(cmap="RdBu_r_alpha", vmin=-1, vmax=1, vmin2=0, vmax2=1) + views["vtx2d_alpha"] = cortex.Vertex2D( + np.ones(nv), np.ones(nv), subj, alpha=(~left).astype(float), **kw2d + ) + d2 = np.full(nv, 0.99); d2[left] = np.nan + views["vtx2d_nan_dim2"] = cortex.Vertex2D(np.ones(nv), d2, subj, **kw2d) + views["vtx2d_full"] = cortex.Vertex2D(np.full(nv, 0.98), np.ones(nv), subj, **kw2d) + d2 = np.full(volshape, 0.99); d2[half_vox] = np.nan + views["vol2d_nan_dim2"] = cortex.Volume2D(np.ones(volshape), d2, subj, xfmname, **kw2d) + views["vol2d_alpha"] = cortex.Volume2D( + np.ones(volshape), np.full(volshape, 0.98), subj, xfmname, + alpha=(~half_vox).astype(float), **kw2d + ) + return views + + +SEQUENCES = { + "vertex_nan": ["vtx_nan", "vtx_full", "vtx_nan"], + "vertexrgb_alpha": ["vtxrgb_a0", "vtxrgb_full", "vtxrgb_a0"], + "vertex_scalar_vs_rgb": ["vtx_nan", "vtxrgb_full", "vtx_full", "vtxrgb_a0", "vtx_full"], + "volume_nan": ["vol_nan", "vol_full", "vol_nan"], + "volumergb_alpha": ["volrgb_a0", "volrgb_full", "volrgb_a0"], + "volume_vs_vertex": ["vol_nan", "vtx_full", "vol_full", "vtx_nan", "vol_full"], + "twod_alpha_and_nan": [ + "vtx2d_alpha", "vtx_full", "vtx2d_nan_dim2", "vtx2d_full", "vtx2d_alpha", + "vol2d_alpha", "vol_full", "vol2d_nan_dim2", "vol2d_alpha", + ], +} + + +class TestSwitching: + """One multi-dataset viewer; per-dataset baselines from single viewers.""" + + @pytest.fixture(autouse=True, scope="class") + def _viewer(self, tmp_path_factory): + cls = type(self) + cls.views = make_views() + cls.tmp = tmp_path_factory.mktemp("switching") + cls.baseline = {} + with cortex.export.headless_viewer( + cortex.Dataset(**cls.views), viewer_params=VIEWER_PARAMS + ) as handle: + cls.handle = handle + handle._set_view(**VIEW) + time.sleep(2) + yield + + @classmethod + def _baseline(cls, name): + if name not in cls.baseline: + with cortex.export.headless_viewer( + cls.views[name], viewer_params=VIEWER_PARAMS + ) as handle: + handle._set_view(**VIEW) + time.sleep(2) + cls.baseline[name] = _render( + handle, str(cls.tmp / ("baseline_%s.png" % name)) + ) + return cls.baseline[name] + + def _switch_and_count(self, name, tag): + handle = type(self).handle + handle.setData(name) + time.sleep(2.5) + return _render(handle, str(type(self).tmp / ("%s_%s.png" % (tag, name)))) + + @pytest.mark.parametrize("sequence", sorted(SEQUENCES)) + def test_sequence(self, sequence): + handle = type(self).handle + errors_before = len(_pageerrors(handle)) + for step, name in enumerate(SEQUENCES[sequence]): + count = self._switch_and_count(name, "%s_%d" % (sequence, step)) + expected = self._baseline(name) + assert expected > 500, "baseline for %s renders nothing" % name + assert abs(count - expected) <= RTOL * expected, ( + "%s step %d: %s rendered %d red pixels after %s, expected %d " + "(NaN/alpha state leaked from the previous dataset?)" + % ( + sequence, step, name, count, + SEQUENCES[sequence][step - 1] if step else "load", expected, + ) + ) + assert len(_pageerrors(handle)) == errors_before, _pageerrors(handle) + + def test_hidden_half_is_really_hidden(self): + """Sanity check of the metric: the half-NaN / half-transparent + dataviews show clearly fewer red pixels than their full versions.""" + for hidden, full in [ + ("vtx_nan", "vtx_full"), ("vtxrgb_a0", "vtxrgb_full"), + ("vol_nan", "vol_full"), ("volrgb_a0", "volrgb_full"), + ("vtx2d_alpha", "vtx2d_full"), ("vtx2d_nan_dim2", "vtx2d_full"), + ]: + assert self._baseline(hidden) < 0.8 * self._baseline(full), (hidden, full) + + +def test_addData_does_not_leak_nan_or_alpha(tmp_path): + """Data pushed into a running viewer must not inherit the previous + dataset's NaN mask or alpha.""" + views = make_views() + # Reference for the RGB view shown on its own (shading differs between a + # colormapped and an RGB view, so RGB is only compared with RGB). + with cortex.export.headless_viewer(views["vtxrgb_full"], viewer_params=VIEWER_PARAMS) as handle: + handle._set_view(**VIEW) + time.sleep(2) + n_rgb_full_alone = _render(handle, str(tmp_path / "rgb_full_alone.png")) + + with cortex.export.headless_viewer(views["vtx_nan"], viewer_params=VIEWER_PARAMS) as handle: + handle._set_view(**VIEW) + time.sleep(2) + n_nan = _render(handle, str(tmp_path / "nan.png")) + + handle.addData(full=views["vtx_full"]) + time.sleep(3) + n_full = _render(handle, str(tmp_path / "full.png")) + + handle.addData(rgb_a0=views["vtxrgb_a0"]) + time.sleep(3) + n_a0 = _render(handle, str(tmp_path / "rgb_a0.png")) + + handle.addData(rgb_full=views["vtxrgb_full"]) + time.sleep(3) + n_rgb_full = _render(handle, str(tmp_path / "rgb_full.png")) + + assert not _pageerrors(handle), _pageerrors(handle) + + assert n_full > 1.5 * n_nan, "NaN mask leaked into data added with addData" + assert abs(n_rgb_full - n_rgb_full_alone) <= RTOL * n_rgb_full_alone, ( + "RGB alpha (or a NaN mask) leaked into data added with addData" + ) + assert n_a0 < 0.8 * n_rgb_full, "alpha=0 half is not hidden after addData" diff --git a/cortex/webgl/data.py b/cortex/webgl/data.py index 0dfed7d4c..3fa9d1002 100644 --- a/cortex/webgl/data.py +++ b/cortex/webgl/data.py @@ -23,7 +23,30 @@ class Package(object): def __init__(self, data): self.dataset = dataset.normalize(data) - self.uniques = list(data.uniques(collapse=True)) + # Dataset.uniques() is a set of BrainData objects, but equality is + # identity: two distinct objects with byte-identical data (e.g. + # Vertex2D(x, x) or VertexRGB(r, r, r)) share the same content-hash + # name and would be packaged -- and reordered -- twice. Keep one per + # name. + # The name only hashes the data bytes, so also check that same-name + # brains are interchangeable for the viewer; otherwise fail loudly + # instead of serving one brain's data with another's metadata. + seen = {} + self.uniques = [] + for brain in data.uniques(collapse=True): + sig = _brain_signature(brain) + if brain.name in seen: + if seen[brain.name] != sig: + raise ValueError( + "Two dataviews contain data with identical bytes but " + "different metadata (%s vs %s); pycortex identifies " + "data by its content hash and cannot serve both in one " + "viewer. Make the data differ (e.g. add a tiny offset)." + % (seen[brain.name], sig) + ) + continue + seen[brain.name] = sig + self.uniques.append(brain) self.subjects = set() self.brains = dict() @@ -119,6 +142,24 @@ def image_names(self, fmt="/data/{name}/{frame}/"): return names +def _brain_signature(brain): + """Metadata that must match for two same-name brains to be interchangeable.""" + import hashlib + + mask = getattr(brain, "_mask", None) + if mask is not None and not isinstance(mask, str): + mask = hashlib.sha1(np.ascontiguousarray(mask).tobytes()).hexdigest() + data = getattr(brain, "data", None) # RGB dataviews have channels instead + return ( + type(brain).__name__, + brain.subject, + getattr(brain, "xfmname", None), + None if data is None else tuple(np.shape(data)), + None if data is None else str(np.asarray(data).dtype), + mask, + ) + + def _pack_png(mosaic): from PIL import Image diff --git a/cortex/webgl/resources/js/dataset.js b/cortex/webgl/resources/js/dataset.js index 9a0a71643..069475b1b 100644 --- a/cortex/webgl/resources/js/dataset.js +++ b/cortex/webgl/resources/js/dataset.js @@ -59,6 +59,12 @@ var dataset = (function(module) { this.data.push(module.brains[json.data[i]]); } } + // Optional per-voxel/vertex alpha map (Volume2D/Vertex2D alpha=). + // Shipped as a regular float brain; multiplied into the colormapped + // color by the shaders (#define DATAALPHA). NaN in it -> transparent. + this.alphaData = null; + if (json.alpha !== undefined && json.alpha !== null && json.alpha.length > 0) + this.alphaData = module.brains[json.alpha[0]]; this.name = json.name; this.description = json.desc; this.frames = this.data[0].frames; @@ -105,6 +111,9 @@ var dataset = (function(module) { this.uniforms.mosaic = { type:'v2v', value:[new THREE.Vector2(6, 6), new THREE.Vector2(6, 6)]}; this.uniforms.dshape = { type:'v2v', value:[new THREE.Vector2(100, 100), new THREE.Vector2(100, 100)]}; this.uniforms.volxfm = { type:'m4v', value:[new THREE.Matrix4(), new THREE.Matrix4()] }; + // alpha map textures (frame, next frame); sampled with dim1's + // transform/mosaic since the alpha map shares dim1's xfm. + this.uniforms.dataalpha = { type:'tv', value:[null, null]}; } this._dispatch = this.dispatchEvent.bind(this); @@ -116,8 +125,11 @@ var dataset = (function(module) { // $.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 children = this.data.slice(); + if (this.alphaData !== null) + children.push(this.alphaData); var allready = []; - for (var i = 0; i < this.data.length; i++) { + for (var i = 0; i < children.length; i++) { allready.push(false); } var checkResolve = function() { @@ -131,9 +143,9 @@ var dataset = (function(module) { checkResolve(); } }; - for (var i = 0; i < this.data.length; i++) { + for (var i = 0; i < children.length; i++) { (function(idx) { - this.data[idx].loaded + children[idx].loaded .progress(function(available) { if (available > this.delay) markReady(idx); }.bind(this)) @@ -247,6 +259,8 @@ var dataset = (function(module) { opts.sampler = module.samplers[this.filter]; opts.rgb = this.data[0].raw; opts.twod = this.data.length > 1; + // volumes only: vertex alpha is folded into the nanmask attribute + opts.dataalpha = this.alphaData !== null && !this.vertex; opts.voxline = (viewopts.voxlines==='true'); var shadecode = shaderfunc(opts); var shader = new THREE.ShaderMaterial({ @@ -271,6 +285,8 @@ var dataset = (function(module) { for (var i = 0; i < this.data.length; i++) { this.data[i].init(this.uniforms, i, this.xfm, this.filter); } + if (this.alphaData !== null) + this.alphaData.setFilter(this.filter); this.setFrame(0); }; module.DataView.prototype.setFrame = function(time) { @@ -281,28 +297,64 @@ 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. + // Optional alpha map (volumes): bind the alpha textures for this + // frame and the next. A single-frame alpha is reused for every frame + // of a movie (both slots point at the same texture). + var alpha = this.alphaData; + if (alpha !== null && !this.vertex && alpha.textures.length > 0) { + var at = alpha.textures; + this.uniforms.dataalpha.value[0] = at[fframe.mod(at.length)]; + this.uniforms.dataalpha.value[1] = at[(fframe+1).mod(at.length)]; + } + // Vertex data: build the single shared "nanmask" attribute, which the + // shader multiplies into the color. It is 0 wherever any dim (or the + // alpha map) is NaN, otherwise the frame-mixed alpha value (1 when + // there is no alpha map). Folding alpha into this attribute instead + // of adding attributes matters: a 2D vertex view already uses 15 of + // the 16 guaranteed vertex attributes, and exceeding the limit makes + // the program fail to link silently. if (this.vertex && !this.data[0].raw && this.data[0].nanmasks.length > 0) { - var dim0 = this.data[0].nanmasks[fframe]; + // The shader mixes frame fframe with the next one (framemix), and + // NaN was replaced by 0 in the buffers, so both adjacent frames + // must be valid for every dim (and the alpha map). + var masks = []; + var fmix = frame - fframe; + var pushFrames = function(nanmasks) { + if (nanmasks === undefined || nanmasks.length === 0) return; + var f0 = fframe.mod(nanmasks.length), f1 = (fframe + 1).mod(nanmasks.length); + masks.push(nanmasks[f0]); + // the next frame only matters while actually blending into it + if (fmix > 0 && f1 !== f0) masks.push(nanmasks[f1]); + }; + for (var d = 0; d < this.data.length; d++) + pushFrames(this.data[d].nanmasks); + var a0 = null, a1 = null, amix = 0.0; + if (alpha !== null && alpha.verts.length > 0) { + pushFrames(alpha.nanmasks); + a0 = alpha.verts[fframe.mod(alpha.verts.length)]; + a1 = alpha.verts[(fframe+1).mod(alpha.verts.length)]; + amix = fmix; + } var combined; - if (this.data.length === 1) { - combined = dim0; + if (masks.length === 1 && a0 === null) { + combined = masks[0]; } 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 n = masks[0][side].array.length; + var out = new Float32Array(n); + for (var i = 0; i < n; i++) { + var ok = 1.0; + for (var m = 0; m < masks.length; m++) { + if (masks[m][side].array[i] < 0.5) { ok = 0.0; break; } + } + if (ok > 0 && a0 !== null) + ok = (1.0 - amix) * a0[side].array[i] + amix * a1[side].array[i]; + out[i] = ok; } var attr = new THREE.BufferAttribute(out, 1); attr.needsUpdate = true; return attr; - }.bind(this)); + }); } this._dispatch({type:"attribute", name:"nanmask", value:combined}); } diff --git a/cortex/webgl/resources/js/mriview.js b/cortex/webgl/resources/js/mriview.js index 04a8ca54c..48da5c527 100644 --- a/cortex/webgl/resources/js/mriview.js +++ b/cortex/webgl/resources/js/mriview.js @@ -791,6 +791,11 @@ var mriview = (function(module) { 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) { + // NaN was replaced by 0 in the GPU buffer (dataset.js); report + // NaN from the mask instead of a fake 0. + if (d.nanmasks !== undefined && d.nanmasks.length > 0 && + d.nanmasks[0][hemiIdx].array[vertex] < 0.5) + return NaN return d.verts[0][hemiIdx].array[vertex] }) } @@ -989,6 +994,11 @@ var mriview = (function(module) { 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) { + // NaN was replaced by 0 in the GPU buffer (dataset.js); report + // NaN from the mask instead of a fake 0. + if (d.nanmasks !== undefined && d.nanmasks.length > 0 && + d.nanmasks[0][hemiIdx].array[vertex] < 0.5) + return NaN return d.verts[0][hemiIdx].array[vertex] }) } diff --git a/cortex/webgl/resources/js/mriview_surface.js b/cortex/webgl/resources/js/mriview_surface.js index 6bc8dac8f..8b56c67d8 100644 --- a/cortex/webgl/resources/js/mriview_surface.js +++ b/cortex/webgl/resources/js/mriview_surface.js @@ -30,6 +30,7 @@ var mriview = (function(module) { this.volume = 0; this._layers = 1; this._dither = false; + this._nanmean = true; // average only non-NaN layers (matches quickflat nanmean=True) this._pivot = 0; this._shift = 0; this._specular = parseFloat(viewopts.specularity); @@ -107,6 +108,7 @@ var mriview = (function(module) { layers: {action:[this, "setLayers", {1:1, 4:4, 8:8, 16:16, 32:32}]}, toggleMultipleLayers: {action: this.toggleMultipleLayers.bind(this), key: 'm', hidden: true, help: "Toggle multiple layers"}, dither: {action:[this, "setDither"]}, + nanmean: {action:[this, "setNanmean"]}, sampler: {action:[this, "setSampler", ["nearest", "trilinear"]]}, uniform_illumination: {action:[this, "setUniformIllumination"]}, }); @@ -380,6 +382,7 @@ var mriview = (function(module) { extratex: this.uniforms.extratex.value !== null, halo: false, dither: this._dither, + nanmean: this._nanmean, equivolume: this._equivolume, sampler: this._sampler, }); @@ -665,6 +668,12 @@ var mriview = (function(module) { this._dither = val; this.resetShaders(); } + module.Surface.prototype.setNanmean = function(val) { + if (val === undefined) + return this._nanmean; + this._nanmean = val; + this.resetShaders(); + } module.Surface.prototype.setSampler = function(val) { if (val === undefined) return this._sampler; diff --git a/cortex/webgl/resources/js/shaderlib.js b/cortex/webgl/resources/js/shaderlib.js index f67d9a9f0..65273fd6d 100644 --- a/cortex/webgl/resources/js/shaderlib.js +++ b/cortex/webgl/resources/js/shaderlib.js @@ -213,6 +213,8 @@ var Shaderlib = (function() { header += "#define RGBCOLORS\n"; if (opts.twod) header += "#define TWOD\n"; + if (opts.dataalpha) + header += "#define DATAALPHA\n"; if (!opts.viewspace) header += "#define SAMPLE_WORLD\n"; if (opts.lights !== undefined && !opts.lights) @@ -270,6 +272,9 @@ var Shaderlib = (function() { "uniform vec2 mosaic[2];", "uniform vec2 dshape[2];", "uniform sampler2D data[4];", + "#ifdef DATAALPHA", + "uniform sampler2D dataalpha[2];", + "#endif", "varying vec3 vPos_x;", "varying vec3 vPos_y;", @@ -289,6 +294,9 @@ var Shaderlib = (function() { "#else", "vec4 values = vec4(0.);", "#endif", + "#ifdef DATAALPHA", + "vec2 avals = vec2(0.);", + "#endif", "#ifdef RGBCOLORS", "color[0] += "+sampler+"_x(data[0], vPos_x);", @@ -301,10 +309,21 @@ var Shaderlib = (function() { "values.w += "+sampler+"_y(data[3], vPos_y).r;", "#endif", "#endif", + "#ifdef DATAALPHA", + "avals.x += "+sampler+"_x(dataalpha[0], vPos_x).r;", + "avals.y += "+sampler+"_x(dataalpha[1], vPos_x).r;", + "#endif", "#ifdef RGBCOLORS", "vec4 vColor = mix(color[0], color[1], framemix);", "#else", "vec4 vColor = colorlut(values);", + "#endif", + "#ifdef DATAALPHA", + // alpha map: NaN (neither <=0 nor >0) -> transparent, else + // scale all four (premultiplied) channels by the alpha value. + "bvec2 avalid = notEqual(lessThanEqual(avals, vec2(0.)), lessThan(vec2(0.), avals));", + "float aval = clamp(mix(avals.x, avals.y, framemix), 0., 1.);", + "vColor = all(avalid) ? vColor * aval : vec4(0.);", "#endif", "vColor *= dataAlpha;", @@ -333,6 +352,10 @@ var Shaderlib = (function() { header += "#define RGBCOLORS\n"; if (opts.twod) header += "#define TWOD\n"; + if (opts.dataalpha) + header += "#define DATAALPHA\n"; + if (opts.nanmean === undefined || opts.nanmean) + header += "#define NANMEAN\n"; var sampler = opts.sampler || "nearest"; var morphs = opts.morphs; @@ -478,6 +501,9 @@ var Shaderlib = (function() { "uniform vec2 mosaic[2];", "uniform vec2 dshape[2];", "uniform sampler2D data[4];", + "#ifdef DATAALPHA", + "uniform sampler2D dataalpha[2];", + "#endif", "uniform vec3 slicexn;", // normal vector for the x sliceplane "uniform vec3 sliceyn;", @@ -537,24 +563,69 @@ var Shaderlib = (function() { "vec4 color[2]; color[0] = vec4(0.), color[1] = vec4(0.);", "#else", "vec4 values = vec4(0.);", + "#endif", + "float nvalid = 0.;", // number of layer samples averaged in + "#ifdef DATAALPHA", + "vec2 avals = vec2(0.);", "#endif", "", ].join("\n"); //Create samplers for texture volume sampling var fragMid = ""; - var factor = layers > 1 ? (1/layers).toFixed(6) : "1."; var sampling = [ "#ifdef RGBCOLORS", - "color[0] += "+factor+"*"+sampler+"_x(data[0], coord_x);", - "color[1] += "+factor+"*"+sampler+"_x(data[1], coord_x);", + // RGBA textures: NaN already became alpha 0 on the Python side, + // so with NANMEAN a fully transparent sample counts as missing + // and is left out of the (premultiplied) average. + "{", + "vec4 c0 = "+sampler+"_x(data[0], coord_x);", + "vec4 c1 = "+sampler+"_x(data[1], coord_x);", + "#ifdef NANMEAN", + // current frame only: for single-frame data, data[1] is an + // unbound sampler, which reads as opaque black (alpha 1). + "if (c0.a > 0.) {", + "#else", + "{", + "#endif", + "color[0] += c0;", + "color[1] += c1;", + "nvalid += 1.;", + "}", + "}", "#else", - "values.x += "+factor+"*"+sampler+"_x(data[0], coord_x).r;", - "values.y += "+factor+"*"+sampler+"_x(data[1], coord_x).r;", + // One layer sample. With NANMEAN (default, matching quickflat's + // nanmean=True) a sample containing NaN in any frame/dimension + // (or in the alpha map) is left out of the average and the + // fragment is transparent only if no layer was valid. Without + // it, one NaN at any depth makes the whole fragment transparent. + "{", + "vec4 s = vec4(0.);", + "s.x = "+sampler+"_x(data[0], coord_x).r;", + "s.y = "+sampler+"_x(data[1], coord_x).r;", "#ifdef TWOD", - "values.z += "+factor+"*"+sampler+"_y(data[2], coord_y).r;", - "values.w += "+factor+"*"+sampler+"_y(data[3], coord_y).r;", + "s.z = "+sampler+"_y(data[2], coord_y).r;", + "s.w = "+sampler+"_y(data[3], coord_y).r;", + "#endif", + "#ifdef DATAALPHA", + "vec2 sa = vec2("+sampler+"_x(dataalpha[0], coord_x).r, "+sampler+"_x(dataalpha[1], coord_x).r);", + "#endif", + "#ifdef NANMEAN", + "bool ok = all(notEqual(lessThanEqual(s, vec4(0.)), lessThan(vec4(0.), s)));", + "#ifdef DATAALPHA", + "ok = ok && all(notEqual(lessThanEqual(sa, vec2(0.)), lessThan(vec2(0.), sa)));", + "#endif", + "#else", + "bool ok = true;", "#endif", + "if (ok) {", + "values += s;", + "nvalid += 1.;", + "#ifdef DATAALPHA", + "avals += sa;", + "#endif", + "}", + "}", "#endif", ].join("\n"); @@ -604,6 +675,17 @@ var Shaderlib = (function() { } var fragTail = [ + "if (nvalid > 0.) {", + "#ifdef RGBCOLORS", + "color[0] /= nvalid;", + "color[1] /= nvalid;", + "#else", + "values /= nvalid;", + "#ifdef DATAALPHA", + "avals /= nvalid;", + "#endif", + "#endif", + "}", "#ifdef HALO_RENDER", "if (vMedial < .999) {", "float dweight = gl_FragCoord.w;", @@ -620,7 +702,14 @@ var Shaderlib = (function() { "#ifdef RGBCOLORS", "vec4 vColor = mix(color[0], color[1], framemix);", "#else", - "vec4 vColor = colorlut(values);", + "vec4 vColor = nvalid > 0. ? colorlut(values) : vec4(0.);", + "#endif", + "#ifdef DATAALPHA", + // alpha map: NaN (neither <=0 nor >0) -> transparent, else + // scale all four (premultiplied) channels by the alpha value. + "bvec2 avalid = notEqual(lessThanEqual(avals, vec2(0.)), lessThan(vec2(0.), avals));", + "float aval = clamp(mix(avals.x, avals.y, framemix), 0., 1.);", + "vColor = all(avalid) ? vColor * aval : vec4(0.);", "#endif", "vColor *= dataAlpha;", //"vColor.a = (values.x - vmin[0]) / (vmax[0] - vmin[0]);", @@ -688,6 +777,9 @@ var Shaderlib = (function() { header += "#define RGBCOLORS\n"; if (opts.twod) header += "#define TWOD\n"; + // Note: no DATAALPHA define here. For vertex data the alpha map + // is folded into the nanmask attribute by dataset.js (adding + // attributes would exceed MAX_VERTEX_ATTRIBS for 2D views). var morphs = opts.morphs; var volume = opts.volume || 0; @@ -753,11 +845,12 @@ var Shaderlib = (function() { "cuv.y = (mix(data2, data3, framemix) - vmin[1]) / (vmax[1] - vmin[1]);", "#endif", "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.);", + // Data mask: WebGL drivers sanitize NaN in vertex attributes, + // so dataset.js detects NaN in JavaScript and passes a mask + // (0 = NaN in any dim or in the alpha map, otherwise the + // per-vertex alpha, 1 when there is no alpha map). Scaling + // all four channels keeps the premultiplied convention. + "vColor *= clamp(nanmask, 0., 1.);", "#endif", "#ifdef CORTSHEET", diff --git a/docs/dataset.rst b/docs/dataset.rst index 491eff2c0..cfc1ba172 100644 --- a/docs/dataset.rst +++ b/docs/dataset.rst @@ -95,11 +95,25 @@ If you provided either numpy arrays or :class:`Volume` objects without vmin/vmax 2D dataviews ~~~~~~~~~~~~ -In order to specify 2D data views in webgl, this helper class lets you specify a pair of :class:`Volume` objects to be plotted using a 2D colormap. Currently, quickflat does not yet support 2D colormaps. To declare a 2D dataview:: +This helper class lets you specify a pair of :class:`Volume` (or :class:`Vertex`) objects to be plotted using a 2D colormap, both in webgl and in quickflat. To declare a 2D dataview:: dim1 = cortex.Volume.random(subject, xfmname) dim2 = cortex.Volume.random(subject, xfmname) - twod = cortex.Volume2D(dim1, dim2, subject=None, xfmname=None, vmin2=None, vmax2=None, **kwargs) + twod = cortex.Volume2D(dim1, dim2, subject=None, xfmname=None, vmin2=None, vmax2=None, alpha=None, **kwargs) + +The optional ``alpha`` (an array in [0, 1], or a :class:`Volume`/:class:`Vertex` normalized by its own ``vmin``/``vmax``) is a per-voxel opacity that is multiplied into the alpha channel of the 2D colormap. Many bundled 2D colormaps (``*_alpha``) already encode opacity along their second axis; ``alpha`` is useful when both colormap axes are needed for data. + +NaN and transparency +~~~~~~~~~~~~~~~~~~~~ +A NaN anywhere at a voxel or vertex means that its value is undefined, and pycortex renders it fully transparent so that the curvature shows through. This holds for every dataview type and for both renderers (quickflat and the webgl viewer): + + * :class:`Volume` / :class:`Vertex`: NaN data is transparent. + * :class:`Volume2D` / :class:`Vertex2D`: NaN in *either* dimension, or in the ``alpha`` map, is transparent. + * :class:`VolumeRGB` / :class:`VertexRGB`: NaN in *any* color channel, or in the ``alpha`` channel, is transparent. + +Where there is no NaN, the opacity (2D colormap alpha, ``alpha=`` keyword, RGB alpha channel) is honored as given. In the webgl viewer, NaN masks and opacity belong to each dataview and are not carried over when switching between datasets. + +When several voxels contribute to one pixel (quickflat averages across cortical thickness; the webgl viewer when ``layers`` > 1), NaN voxels are left out of the average by default and the pixel is transparent only if every contributing voxel is NaN. This is ``nanmean=True`` in :func:`cortex.quickflat.make_figure` and the ``nanmean`` toggle in the webgl surface controls; switch either off to make any NaN contribution transparent. For :class:`VolumeRGB` / :class:`VertexRGB` the NaN has already become alpha 0 when the data is converted to RGBA, so there fully transparent voxels count as missing: ``nanmean`` skips them, and switching it off gives the alpha-weighted average instead. Dataset ------- diff --git a/examples/datasets/plot_data_with_alpha.py b/examples/datasets/plot_data_with_alpha.py index 8c349383e..0140a2699 100644 --- a/examples/datasets/plot_data_with_alpha.py +++ b/examples/datasets/plot_data_with_alpha.py @@ -26,6 +26,15 @@ per-voxel/per-vertex array (or a :class:`Volume`/:class:`Vertex`) in ``[0, 1]``. +3. **2D data with an alpha map** -- :class:`Volume2D` / :class:`Vertex2D` + also accept ``alpha=``, which is multiplied into the colormap alpha. This + lets you combine any 2D colormap (e.g. a covariance map) with a separate + opacity map. + +In every case, NaN anywhere at a voxel/vertex (in the data, in either +dimension of a 2D view, in any RGB channel, or in the alpha map itself) +renders fully transparent, so the curvature shows through. + 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 @@ -178,16 +187,49 @@ def _bump(surf, seed, sigma): plt.suptitle("VertexRGB(alpha=accuracy): RGB channels masked by 'accuracy'") plt.show() +# %% +# Pattern 3: 2D data + a separate alpha map via Volume2D(alpha=...) +# ----------------------------------------------------------------- +# +# A 2D colormap can encode two quantities (here a covariance-style map of +# ``data`` against ``accuracy``) while a third map controls opacity. The +# ``alpha=`` map is multiplied into the colormap's own alpha channel. +# Here we fade out the inferior part of the volume. NaNs in the data are +# always transparent, whatever ``alpha`` says: the posterior slab below is +# NaN and renders as bare curvature. +data_vol_nan = data_vol.copy() +data_vol_nan[:, :25, :] = np.nan # posterior slab: undefined +alpha_inferior_fade = np.clip(zz / 30.0, 0, 1) # 0 at the bottom, 1 at the top + +v2d_alpha = cortex.Volume2D( + data_vol_nan, + accuracy_vol, + subject, + xfm, + cmap="RdBu_covar", + vmin=-1, + vmax=1, + vmin2=0, + vmax2=1, + alpha=alpha_inferior_fade, +) +cortex.quickshow(v2d_alpha, with_colorbar=True, with_curvature=True) +plt.suptitle("Volume2D(alpha=...): 2D colormap, separate opacity, NaN slab") +plt.show() + # %% # Notes # ----- # -# * Both patterns produce the same composite formula at the pixel level: +# * All 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). +# based on what the "data" is: scalar (use Pattern 1), RGB (use +# Pattern 2) or two quantities plus an opacity (use Pattern 3). +# * NaN anywhere at a voxel/vertex -- data, either 2D dimension, any RGB +# channel, or the alpha map -- is rendered fully transparent. # * The same objects work in the WebGL viewer: -# ``cortex.webgl.show(v2d)`` etc.; opacity is honored identically. +# ``cortex.webgl.show(v2d)`` etc.; NaN and opacity are honored +# identically, also when switching between datasets. # * 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