diff --git a/cortex/brainctm.py b/cortex/brainctm.py index 8c76c986f..b67c88b57 100644 --- a/cortex/brainctm.py +++ b/cortex/brainctm.py @@ -29,6 +29,7 @@ class BrainCTM: def __init__(self, subject, decimate=False): self.subject = subject self.types = [] + self.has_volume = False left, right = db.get_surf(subject, "fiducial") try: @@ -42,6 +43,7 @@ def __init__(self, subject, decimate=False): self.left = DecimatedHemi(left[0], left[1], fleft[1], pia=pleft[0]) self.right = DecimatedHemi(right[0], right[1], fright[1], pia=pright[0]) self.addSurf("wm", addtype=False, renorm=False) + self.has_volume = True except IOError: self.left = DecimatedHemi(left[0], left[1], fleft[1]) self.right = DecimatedHemi(right[0], right[1], fright[1]) @@ -52,6 +54,7 @@ def __init__(self, subject, decimate=False): self.left = Hemi(pleft[0], left[1]) self.right = Hemi(pright[0], right[1]) self.addSurf("wm", addtype=False, renorm=False) + self.has_volume = True except IOError: self.left = Hemi(left[0], left[1]) self.right = Hemi(right[0], right[1]) @@ -112,6 +115,43 @@ def addCurvature(self, **kwargs): self.left.aux[:,1] = npz.left self.right.aux[:,1] = npz.right + def addEquivolumeAreas(self, **kwargs): + """Load the vertex areas the viewer's equivolume depth sampling needs. + + They ride in the two spare components of ``auxdat`` -- ``x`` is the + medial wall mask and ``y`` the curvature -- rather than taking vertex + attribute slots of their own, which the surface shaders are already + right up against. + """ + if not self.has_volume: + return + npz = db.get_surfinfo(self.subject, type='equivolume_areas', **kwargs) + try: + self.left.aux[:,2] = npz['wm_left'][self.left.mask] + self.left.aux[:,3] = npz['pia_left'][self.left.mask] + self.right.aux[:,2] = npz['wm_right'][self.right.mask] + self.right.aux[:,3] = npz['pia_right'][self.right.mask] + except AttributeError: + self.left.aux[:,2] = npz['wm_left'] + self.left.aux[:,3] = npz['pia_left'] + self.right.aux[:,2] = npz['wm_right'] + self.right.aux[:,3] = npz['pia_right'] + npz.close() + + def addBumpyFlat(self, **kwargs): + """Load the relief that gives the flatmap its bumps. + + Each vertex gets a (0, 0, height) offset from its position on the flat + white matter surface, in flatmap units. + """ + if self.flatlims is None or not self.has_volume: + return + + npz = db.get_surfinfo(self.subject, type='bumpy_flatmap', **kwargs) + self.left.setBump(npz['bump_left']) + self.right.setBump(npz['bump_right']) + npz.close() + def save(self, path, method='mg2', external_svg=None, overlays_available=None, **kwargs): """Save CTM file for static html display. @@ -212,6 +252,7 @@ def __init__(self, pts, polys, norms=None): self.pts = pts self.polys = polys self.flat = None + self.bump = None self.surfs = {} self.aux = np.zeros((len(self.ctm), 4)) @@ -232,8 +273,24 @@ def setFlat(self, pts): self.ctm.addUV(pts[:,:2].astype(float), 'uv') self.flat = pts[:,:2] + def setBump(self, offsets): + '''Bumpy flatmap offsets, padded to the four components a ctm attribute + map always has.''' + self.bump = np.hstack([offsets, np.zeros((len(offsets), 1))]) + def save(self, **kwargs): self.ctm.addAttrib(self.aux, 'auxdat') + if self.bump is not None: + self.ctm.addAttrib(self.bump, 'flatoffset') + + # OpenCTM only has eight attribute map slots. Two go to auxdat and the + # bumpy flatmap offsets and one to the white matter surface, so a viewer + # can carry at most five extra surfaces. + if len(self.ctm.attribs) > 8: + raise ValueError( + "too many surfaces for one ctm file: %d attribute maps, and " + "OpenCTM allows 8. Pass fewer entries in `types`. (%s)" + % (len(self.ctm.attribs), ", ".join(self.ctm.attribs))) self.ctm.save(**kwargs) ctm = CTMfile(self.tfName) mesh = ctm.getMesh() @@ -276,6 +333,9 @@ def setFlat(self, pts): def addSurf(self, pts, **kwargs): super().addSurf(pts[self.mask], **kwargs) + def setBump(self, offsets): + super().setBump(offsets[self.mask]) + def make_pack(outfile, subj, types=("inflated",), method='raw', level=0, decimate=False, disp_layers=['rois'], external_svg=None, overlays_available=None,): @@ -288,6 +348,8 @@ def make_pack(outfile, subj, types=("inflated",), method='raw', level=0, ctm = BrainCTM(subj, decimate=decimate) ctm.addCurvature() + ctm.addEquivolumeAreas() + ctm.addBumpyFlat() for name in types: ctm.addSurf(name) diff --git a/cortex/database.py b/cortex/database.py index cbcbc60d1..ca18e63a3 100644 --- a/cortex/database.py +++ b/cortex/database.py @@ -251,7 +251,8 @@ def get_anat(self, subject: str, type: Literal['raw', 'brainmask', 'whitematter' from . import volume return volume.anat2epispace(anatnib.get_fdata().T.astype(float), subject, xfmname, order=order) - def get_surfinfo(self, subject: str, type: str="curvature", recache: bool=False, **kwargs) -> Vertex: + def get_surfinfo(self, subject: str, type: str="curvature", + recache: bool=False, **kwargs) -> Optional[Vertex]: """Return auxiliary surface information from the filestore. Surface info is defined as anatomical information specific to a subject in surface space. A Vertex class will be returned as necessary. Info not found in the filestore will be automatically generated. diff --git a/cortex/defaults.cfg b/cortex/defaults.cfg index 3bd21f41b..a753e4276 100644 --- a/cortex/defaults.cfg +++ b/cortex/defaults.cfg @@ -151,6 +151,15 @@ specularity = 1.0 overlayscale = 1 anim_speed = 2 bumpy_flatmap = false +# How much to exaggerate the relief of a bumpy flatmap. 1.0 is the computed +# geometry at its true scale; larger values make the gyri and sulci easier to +# read at the cost of no longer being to scale. The default exaggerates a +# little, because a 2-5 mm slab is subtle next to a whole flatmap and the +# relief is there to be read. Purely a display setting -- it does not change +# the cached geometry, so it can be adjusted freely. This is where the viewer's +# bumpy_flatmap_scale slider starts; the slider runs to five times true scale, +# or to twice this value if that is higher. +bumpy_flatmap_scale = 1.55 allow_tilt = false [curvature] diff --git a/cortex/freesurfer.py b/cortex/freesurfer.py index f520b9b7b..ac3e74705 100644 --- a/cortex/freesurfer.py +++ b/cortex/freesurfer.py @@ -391,6 +391,31 @@ def import_flat(fs_subject, patch, hemis=['lh', 'rh'], cx_subject=None, os.unlink(overlays_file) # Regenerate it? + # clear_cache only empties the cache/ directory, but surface-info/ holds + # derived data too, and some of it is computed from the flatmap that has + # just been replaced. Those files would otherwise stay stale forever. + _clear_flat_surfinfo(cx_subject) + + +def _clear_flat_surfinfo(cx_subject): + """Delete the surface info files that are derived from the flatmap. + + Called when a flatmap is (re)imported. Anything computed from the flat + surface -- distortion, the flatmap border, the bumpy flatmap -- + describes the *old* flatmap once a new one is imported. + """ + surfiform = database.db.get_paths(cx_subject)['surfinfo'] + directory = os.path.dirname(surfiform) + if not os.path.isdir(directory): + return + flat_derived = ("distortion", "flat_border", "bumpy_flatmap") + for fname in os.listdir(directory): + # Options are appended to the type as "type[opt=val].npz", so match on + # the leading type name rather than the whole filename. + stem = fname.split("[")[0].rsplit(".", 1)[0] + if stem in flat_derived: + os.unlink(os.path.join(directory, fname)) + def _remove_disconnected_polys(polys): """Remove polygons that are not in the main connected component. diff --git a/cortex/polyutils/__init__.py b/cortex/polyutils/__init__.py index 3806c6517..8444ba9f0 100644 --- a/cortex/polyutils/__init__.py +++ b/cortex/polyutils/__init__.py @@ -1,4 +1,9 @@ +from .bumpy import ( + FlatSlab, + folding_height, + naive_prism_height, +) from .distortion import Distortion from .misc import ( _memo, diff --git a/cortex/polyutils/bumpy.py b/cortex/polyutils/bumpy.py new file mode 100644 index 000000000..60f2122bf --- /dev/null +++ b/cortex/polyutils/bumpy.py @@ -0,0 +1,251 @@ +"""Bumpy flatmaps: giving a flatmap the relief of the cortical slab. + +Cortex is 2 to 5 mm thick, so if you peeled the slab off the white matter and +laid it down the white side would end up flat and the pial side would sit some +distance above it. That relief is a folding cue which survives even when the +flatmap is covered in data. + +The height is ``V_frustum / A_wm``, the folded volume of each column over the +*folded* white matter area beneath it, which with ``r = sqrt(A_pia / A_wm)`` is +``thickness * (1 + r + r**2) / 3``. `r` is what carries the folding: the pia has +more area than the white matter under a gyral crown and less in a fundus. + +Dividing by the *flattened* area instead is the more obvious model of a slab +laid flat, and it does not work. A flatmap's area distortion measures +essentially uncorrelated with both mean and Gaussian curvature, so as a +denominator it contributes no folding and injects the flattening algorithm's +artifacts in its place; what is left is close to a map of cortical thickness, +which is blobby and reads as round knobs. `naive_prism_height` computes that +version, for comparison. +""" + + +import numpy as np +from scipy import sparse + +try: + from scipy.sparse.linalg import factorized as _factorized +except ImportError: + from scipy.sparse.linalg.dsolve import factorized as _factorized + +from .misc import _memo, face_area, face_volume +from .surface import Surface + +__all__ = ["FlatSlab", "folding_height", "naive_prism_height"] + + +def _flat_plane(flat): + """Flatmap coordinates as a 3D point set lying in the z = 0 plane. + + Only the first two columns of a pycortex flat surface are the flatmap + coordinates; the third is left over from the flattening and is ignored + everywhere else too (see `cortex.brainctm`, which stores `pts[:, :2]`). + """ + plane = np.zeros((len(flat), 3)) + plane[:, :2] = np.asarray(flat)[:, :2] + return plane + + +def _lumped(values, polys, nverts): + """A per-face quantity gathered onto vertices, a third to each corner.""" + return np.bincount(np.asarray(polys).ravel(), + weights=np.repeat(values, 3), minlength=nverts) / 3.0 + + +def _regularise_log_height(flat, polys, values, correlation_length): + """Smooth a positive field in the log, on the flat mesh. + + `Surface.smooth` solves ``(M + t L) y = M x`` with `L` the cotangent + stiffness and `M` the lumped mass, which is the screened-Poisson system + wanted here with ``t = lc**2``. Working in the log keeps a twenty-fold + compression as an offset of three rather than a twenty-fold spike, and + makes this a geometric rather than an arithmetic mean -- the right + averaging for a ratio. + """ + good = values > 0 + target = np.zeros(len(flat)) + target[good] = np.log(values[good]) + if good.any(): + target[~good] = np.median(target[good]) + + surf = Surface(flat, polys) + smoothed = surf.smooth(target, correlation_length ** 2) + # `Surface.smooth` returns zero for vertices in no triangle, which would + # come back as a height of one; leave those at the median instead. + isolated = np.asarray(surf.connected.sum(1)).ravel() == 0 + smoothed[isolated] = target[isolated] + return np.exp(smoothed) + + +def folding_height(flat, wm, pia, polys, correlation_length=0.5): + """Slab height driven by the pial flare, not by the flattening. + + ``V_frustum / A_wm``, which with ``r = sqrt(A_pia / A_wm)`` is + ``thickness * (1 + r + r**2) / 3``. Over a gyral crown the pia has more area + than the white matter beneath it, so `r` rises; in a fundus it falls. See + the module docstring for why the denominator is the folded area. + """ + nverts = len(wm) + polys = np.asarray(polys) + awm = _lumped(face_area(wm[polys]), polys, nverts) + apia = _lumped(face_area(pia[polys]), polys, nverts) + thickness = np.sqrt(((pia - wm) ** 2).sum(1)) + + good = awm > 0 + r = np.zeros(nverts) + r[good] = np.sqrt(np.maximum(apia[good], 0.0) / awm[good]) + height = thickness * (1.0 + r + r ** 2) / 3.0 + + # Same log-space regularisation as the prism height -- this is a ratio too, + # and a vertex whose white matter area nearly vanishes would otherwise + # spike. Smoothed on the flat mesh, not the folded one, because that is + # where the relief is going to be looked at and where `correlation_length` + # is measured for every other field here. + return _regularise_log_height(flat, polys, height, correlation_length) + + +def naive_prism_height(flat, wm, pia, polys): + """Height of a volume-preserving vertical prism over the flatmap. + + The obvious reading of the bumpy flatmap idea: each column keeps its folded + volume over a base of the flattened triangle. Kept for comparison -- it + spikes wherever flattening compressed a triangle hard, and smoothing a field + of ratios afterwards does not remove spikes that dominate the mean. See the + module docstring for the deeper problem with the denominator. + + Parameters + ---------- + flat : 2D ndarray, shape (total_verts, 3) + Location of each vertex in flatmap space. + wm : 2D ndarray, shape (total_verts, 3) + Location of each vertex on the white matter surface. + pia : 2D ndarray, shape (total_verts, 3) + Location of each vertex on the pial surface. + polys : 2D ndarray, shape (total_polys, 3) + Triangle vertex indices, shared by all three surfaces. + + Returns + ------- + height : 1D ndarray, shape (total_verts,) + Height of the pial surface above the flatmap at each vertex, in the + units of the input surfaces. + """ + nverts = len(wm) + polys = np.asarray(polys) + vol = face_volume(wm, pia, polys) + area = face_area(_flat_plane(flat)[polys]) + + # One height per triangle, then averaged onto the vertices. Note this is a + # mean of ratios, and that is the whole problem: a triangle the flattening + # crushed has a tiny area in the denominator, so it contributes an enormous + # height that dominates every vertex it touches. + faceheight = np.zeros(len(polys)) + good = area > 0 + faceheight[good] = vol[good] / area[good] + + counts = np.bincount(polys.ravel(), minlength=nverts) + height = np.bincount(polys.ravel(), weights=np.repeat(faceheight, 3), + minlength=nverts) + return np.where(counts > 0, height / np.maximum(counts, 1), 0.0) + + +class FlatSlab(object): + """The cortical slab's relief, as an offset from the flat white surface. + + `folding_height`, smoothed and with the whole-map swell taken out. + + Parameters + ---------- + flat : 2D ndarray, shape (total_verts, 3) + Location of each vertex in flatmap space. Only the first two columns + are used; the flatmap is taken to lie in the z = 0 plane. + wm, pia : 2D ndarray, shape (total_verts, 3) + The white matter and pial surfaces. + polys : 2D ndarray, shape (total_polys, 3) + Triangle vertex indices of the *flat* surface. Vertices in no triangle + -- the medial wall, which is cut away from the flatmap -- get a zero + offset. + correlation_length : float, optional + Length scale over which the height ratio is regularised -- enough to + stop a vertex whose white matter area nearly vanishes from spiking. + polish : float, optional + Smoothing of the finished relief, as a diffusion time. Beware the + scale: `t` halves a wavelength of ``2 * pi * sqrt(t)``, so the default + cuts at 11 mm and not, as it looks, at 3. Much more takes the gyri too. + detrend : float, optional + Wavelength above which the relief is flattened. Cortex is regionally + thicker in some lobes than others, and on S1 that swell is a third of + the relief's variance in one map-spanning mode -- real, but not folding, + and the viewer's scale setting multiplies it along with everything else. + Pass 0 to keep it. + + Attributes + ---------- + info : dict + Height statistics from the last computed relief. + """ + def __init__(self, flat, wm, pia, polys, correlation_length=0.5, + polish=3.0, detrend=64.0): + self.flat = np.asarray(flat, dtype=np.double) + self.wm = np.asarray(wm, dtype=np.double) + self.pia = np.asarray(pia, dtype=np.double) + self.polys = np.asarray(polys) + self.correlation_length = correlation_length + self.polish = polish + self.detrend = detrend + self.info = {} + self._cache = {} + + if not (len(self.flat) == len(self.wm) == len(self.pia)): + raise ValueError("flat, wm and pia must have the same number of " + "vertices, got %d, %d and %d" + % (len(self.flat), len(self.wm), len(self.pia))) + + @property + @_memo + def _submesh(self): + """Restrict to the vertices that are actually on the flatmap. + + Returns ``(mask, subflat, subwm, subpia, subpolys)``, with `subflat` in + the z = 0 plane and the triangles reindexed to match. + """ + mask = np.zeros(len(self.wm), dtype=bool) + mask[self.polys.ravel()] = True + + vmap = np.zeros(len(self.wm), dtype=np.int64) + vmap[mask] = np.arange(mask.sum()) + + return (mask, _flat_plane(self.flat)[mask], self.wm[mask], + self.pia[mask], vmap[self.polys]) + + @property + @_memo + def relaxed(self): + """The pial surface's height above the flat white matter surface. + + Returns + ------- + offsets : 2D ndarray, shape (total_verts, 3) + The first two columns are zero -- the relief is purely vertical -- + and the third is the height. Vertices off the flatmap are zero. + """ + mask, flat, wm, pia, polys = self._submesh + surf = Surface(flat, polys) + + relief = folding_height(flat, wm, pia, polys, self.correlation_length) + if self.polish: + relief = surf.smooth(relief, self.polish) + if self.detrend: + swell = surf.smooth(relief, (self.detrend / (2 * np.pi)) ** 2) + # about its own mean, so the sheet keeps its average thickness + relief = relief - (swell - swell.mean()) + + self.info = dict(n_verts=int(mask.sum()), + height_mean=float(relief.mean()), + height_std=float(relief.std()), + height_min=float(relief.min()), + height_max=float(relief.max())) + + offsets = np.zeros((len(self.wm), 3)) + offsets[mask, 2] = relief + return offsets diff --git a/cortex/polyutils/misc.py b/cortex/polyutils/misc.py index 20369b467..4660014c8 100644 --- a/cortex/polyutils/misc.py +++ b/cortex/polyutils/misc.py @@ -44,12 +44,20 @@ def face_area(pts): return 0.5 * np.sqrt((np.cross(pts[:,1]-pts[:,0], pts[:,2]-pts[:,0])**2).sum(1)) def face_volume(pts1, pts2, polys): - '''Volume of each face in a polyhedron sheet''' - vols = np.zeros((len(polys),)) - for i, face in enumerate(polys): - vols[i] = brick_vol(np.append(pts1[face], pts2[face], axis=0)) - if i % 1000 == 0: - print(i) + '''Volume of each face in a polyhedron sheet + + The same three-tetrahedron decomposition as `brick_vol`, over every face at + once. Nodes 0-2 are the `pts1` triangle and 3-5 the `pts2` one. + ''' + polys = np.asarray(polys) + corners = np.concatenate([np.asarray(pts1)[polys], + np.asarray(pts2)[polys]], axis=1) + vols = np.zeros(len(polys)) + for a, b, c, d in ((0, 1, 2, 4), (0, 2, 3, 4), (2, 3, 4, 5)): + edges = np.stack([corners[:, b] - corners[:, a], + corners[:, c] - corners[:, a], + corners[:, d] - corners[:, a]], axis=1) + vols += np.abs(np.linalg.det(edges)) / 6 return vols def decimate(pts, polys): diff --git a/cortex/surfinfo.py b/cortex/surfinfo.py index f59525648..465fd512f 100644 --- a/cortex/surfinfo.py +++ b/cortex/surfinfo.py @@ -219,3 +219,90 @@ def make_surface_graph(tris): lines.append(pts[pbnd,:2]) np.savez(outfile, lines=lines, ismwalls=ismwalls) + +def bumpy_flatmap(outfile, subject, **kwargs): + """Give the flatmap the relief of the cortical slab, and cache it. + + Hemispheres with no flat surface get an array of zeros rather than an error. + + Parameters + ---------- + outfile : str + Path where the offsets will be saved as an npz file. + subject : str + Subject in the pycortex database whose flatmap gets the relief. + **kwargs + Passed to `cortex.polyutils.FlatSlab`. + + Notes + ----- + Stored under ``bump_left`` and ``bump_right`` rather than ``left`` and + ``right`` on purpose: `get_surfinfo` turns a file with the latter into a + `Vertex` by concatenating them, which assumes one value per vertex and + would mangle these three-component offsets. + """ + offsets = [] + for hemi in ("lh", "rh"): + wm, _ = db.get_surf(subject, "wm", hemi) + pia, _ = db.get_surf(subject, "pia", hemi) + try: + flat, flatpolys = db.get_surf(subject, "flat", hemi) + except IOError: + offsets.append(np.zeros_like(wm)) + continue + offsets.append(polyutils.FlatSlab(flat, wm, pia, flatpolys, + **kwargs).relaxed) + + # Compressed, and single precision: these are millimetre offsets that end + # up in a float32 vertex attribute anyway, so the extra digits are only + # taking up room in the filestore. + np.savez_compressed(outfile, + bump_left=offsets[0].astype(np.float32), + bump_right=offsets[1].astype(np.float32)) + + +def equivolume_areas(outfile, subject, smooth=1.0): + """ + Compute smoothed vertex areas on the white matter and pial surfaces. + + These are what the webgl viewer's equivolume depth sampling needs in order to + turn a requested volume fraction through the cortical sheet into a position + between the two surfaces. They used to be recomputed in javascript on every + viewer load with a uniform "umbrella" smoothing; computing them here instead + means they are cached, and lets them be smoothed with the cotangent-weighted + operator in `cortex.polyutils.Surface.smooth`, which respects the varying + size of the triangles rather than treating every neighbour equally. + + Parameters + ---------- + outfile : str + Path where the areas will be saved as an npz file. + subject : str + Subject in the pycortex database for whom the areas will be computed. + smooth : float, optional + Amount of smoothing to apply. Default 1.0. Pass 0 for the raw + barycentric vertex areas. + + Notes + ----- + Stored under ``wm_left`` / ``wm_right`` / ``pia_left`` / ``pia_right``; see + the note in `bumpy_flatmap` for why ``left`` and ``right`` are avoided. + + Changing the smoothing operator does move the depths the viewer samples at, + slightly: against the five umbrella iterations the javascript used, the + depth equivolume sampling picks for a requested fraction of 0.5 shifts by + about 0.01 of the cortical thickness at the median on S1. The tail is larger, + but it is concentrated where the two areas are nearly equal and the depth is + poorly determined either way. + """ + areas = dict() + for hemi, side in zip(["lh", "rh"], ["left", "right"]): + for name in ["wm", "pia"]: + pts, polys = db.get_surf(subject, name, hemi) + surf = polyutils.Surface(pts, polys) + # The lumped mass matrix of the Laplace-Beltrami operator is exactly + # the barycentric vertex area, i.e. a third of each incident face. + _, vertex_area, _, _ = surf.laplace_operator + areas["%s_%s" % (name, side)] = surf.smooth(vertex_area, smooth) + + np.savez(outfile, **areas) diff --git a/cortex/tests/test_bumpy.py b/cortex/tests/test_bumpy.py new file mode 100644 index 000000000..f85e4b20b --- /dev/null +++ b/cortex/tests/test_bumpy.py @@ -0,0 +1,332 @@ +"""Tests for the bumpy flatmap in `cortex.polyutils.bumpy`. + +These check the relief against properties that can be stated in advance -- +folding produces relief at constant thickness, the height ignores flatmap +distortion, a ridge stays a ridge -- rather than against a stored result, so +they say whether the answer is right and not just whether it changed. +""" + +import numpy as np +from scipy.spatial import Delaunay + +from cortex import polyutils +from cortex.polyutils import bumpy + + +def slab_grid(n=13, spacing=1.0, thickness=2.5, stretch=1.0): + """A flat rectangular slab, optionally flattened with a uniform stretch. + + Returns ``(flat, wm, pia, polys, index)`` where `index` maps a grid position + to its vertex number, so a test can pick out the middle of the patch and + stay clear of the free edges. + """ + xs = np.arange(n) * spacing + X, Y = np.meshgrid(xs, xs, indexing='ij') + wm = np.stack([X.ravel(), Y.ravel(), np.zeros(X.size)], axis=1) + pia = wm + np.array([0., 0., thickness]) + + flat = wm.copy() + flat[:, :2] *= stretch + + index = np.arange(n * n).reshape(n, n) + polys = [] + for i in range(n - 1): + for j in range(n - 1): + a, b = index[i, j], index[i + 1, j] + c, d = index[i + 1, j + 1], index[i, j + 1] + polys += [[a, b, c], [a, c, d]] + return flat, wm, pia, np.array(polys), index + + +def test_prism_volume_matches_brick_vol(): + """`face_volume` is `brick_vol` over every face at once.""" + rng = np.random.default_rng(0) + wm = rng.normal(size=(30, 3)) + pia = wm + rng.normal(size=(30, 3)) * 0.3 + np.array([0., 0., 2.]) + polys = np.array([[0, 1, 2], [1, 2, 3], [4, 5, 6], [7, 8, 9]]) + + expected = np.array([polyutils.brick_vol(np.append(wm[f], pia[f], axis=0)) + for f in polys]) + got = polyutils.face_volume(wm, pia, polys) + assert np.allclose(got, expected) + + +def test_vertices_off_the_flatmap_are_untouched(): + """Vertices in no flat triangle -- the medial wall -- get a zero offset.""" + flat, wm, pia, polys, index = slab_grid(n=9, stretch=1.1) + # Cut a corner out of the flatmap, the way the medial wall is cut away. + corner = {index[0, 0], index[0, 1], index[1, 0]} + keep = np.array([not (set(f) & corner) for f in polys]) + + offsets = bumpy.FlatSlab(flat, wm, pia, polys[keep]).relaxed + + off_map = np.ones(len(wm), bool) + off_map[polys[keep].ravel()] = False + assert off_map.sum() > 0 + assert np.abs(offsets[off_map]).max() == 0.0 + assert np.abs(offsets[~off_map, 2]).min() > 0.0 + + +def _s1_patch(radius=32): + """A patch of S1's flatmap, with its white, pial and curvature data.""" + from cortex import db + wm, _ = db.get_surf("S1", "wm", "lh") + pia, _ = db.get_surf("S1", "pia", "lh") + flat, flatpolys = db.get_surf("S1", "flat", "lh") + curv = db.get_surfinfo("S1", type="curvature").data[:len(wm)] + + surf = polyutils.Surface(bumpy._flat_plane(flat), flatpolys) + seed = flatpolys[len(flatpolys) // 2, 0] + inside = np.zeros(len(flat), bool) + inside[surf.get_euclidean_patch(seed, radius)["vertex_mask"]] = True + keep = inside[flatpolys].all(1) + + sub = np.zeros(len(flat), bool) + sub[flatpolys[keep].ravel()] = True + remap = np.zeros(len(flat), np.int64) + remap[sub] = np.arange(sub.sum()) + return (flat[sub], wm[sub], pia[sub], remap[flatpolys[keep]], curv[sub]) + + +def _silk(pts, polys): + """RMS angle, in degrees, between the normals of neighbouring triangles. + + This is what "smooth" means to a shader: the shading normal is the + derivative of the surface, so it is the angle between neighbours and not + the height itself that decides whether the relief reads as silk or as + crumpled foil. + """ + tri = pts[polys] + n = np.cross(tri[:, 1] - tri[:, 0], tri[:, 2] - tri[:, 0]) + n /= np.maximum(np.linalg.norm(n, axis=1), 1e-12)[:, None] + + edges = np.sort(np.vstack([polys[:, [0, 1]], polys[:, [1, 2]], + polys[:, [2, 0]]]), axis=1) + face = np.tile(np.arange(len(polys)), 3) + order = np.lexsort((edges[:, 1], edges[:, 0])) + e, f = edges[order], face[order] + shared = np.all(e[:-1] == e[1:], axis=1) + a, b = f[:-1][shared], f[1:][shared] + + dot = np.clip((n[a] * n[b]).sum(1), -1, 1) + return np.degrees(np.sqrt((np.arccos(dot) ** 2).mean())) + + +def test_the_relief_is_smooth_enough_to_shade(): + """The bumped flatmap has to be smooth at the scale of a triangle. + + A shading normal is the derivative of the height field, so noise that is + invisible in the heights is glaring in the lighting. + """ + flat, wm, pia, polys, _ = _s1_patch() + plane = bumpy._flat_plane(flat) + + slab = bumpy.FlatSlab(flat, wm, pia, polys) + relief = slab.relaxed + # Two thresholds, because a bare angle is not scale-free -- the angle + # between neighbouring faces grows with height for the same shape. The + # absolute number is what the eye sees; the ratio is smoothness per unit + # relief, and it is the one that says the surface itself is not rough. + assert _silk(plane + relief, polys) < 4.5 + assert _silk(plane + relief, polys) / relief[:, 2].std() < 5.0 + + rough = bumpy.FlatSlab(flat, wm, pia, polys, polish=0) + assert _silk(plane + rough.relaxed, polys) > _silk(plane + slab.relaxed, + polys) + + +def _band(surf, x, wavelength): + """Low-pass at a given wavelength. + + `Surface.smooth`'s `factor` is a diffusion time, and one backward-Euler step + has transfer function 1/(1 + k^2 t), so `factor = t` is half power at a + wavelength of 2*pi*sqrt(t) -- not sqrt(t), and not the sqrt(2t) the + Gaussian-equivalent sigma suggests. It is easy to be wrong here by 4.4x. + """ + return surf.smooth(x.copy(), (wavelength / (2 * np.pi)) ** 2) + + +def test_the_relief_carries_gyral_scale_signal(): + """The relief has to have energy at the scale of gyri, not just be smooth. + + Every other measure here is local -- `_silk` compares neighbouring faces, + and a correlation against raw curvature is dominated by whichever band has + most variance. None of them can see a relief that is beautifully smooth and + anatomically empty, which is what over-smoothing produces. + """ + flat, wm, pia, polys, curv = _s1_patch() + surf = polyutils.Surface(bumpy._flat_plane(flat), polys) + height = bumpy.FlatSlab(flat, wm, pia, polys).relaxed[:, 2] + + def gyral(x): + """The 8-16 mm band, where the folding lives.""" + return _band(surf, x, 16.0) - _band(surf, x, 8.0) + + band, cband = gyral(height), gyral(curv) + assert np.corrcoef(band, cband)[0, 1] > 0.35, ( + "the relief has no folding signal left at gyral scale") + + # and it must not have been smoothed into a featureless sheet: the gyral + # band has to hold real amplitude, not a rounding error on the mean + assert band.std() / height.std() > 0.1, ( + "the gyral band has been smoothed away relative to the whole relief") + + +def _elongation(surf, field, scale=8.0): + """Structure-tensor coherence: 1 for a ridge, 0 for a round bump. + + Measured from the field itself, which is what keeps it an independent check + even when the smoothing has been oriented by curvature. + """ + g = surf.surface_gradient(np.asarray(field, float), at_verts=False) + idx = np.array([(0, 0), (0, 1), (1, 1)]) + packed = (g[:, :, None] * g[:, None, :])[:, idx[:, 0], idx[:, 1]] + + area = surf.face_areas + w = np.maximum(np.asarray(surf.connected.dot(area)).ravel(), 1e-20) + v = surf.connected.dot(area[:, None] * packed) / w[:, None] + v = np.column_stack([surf.smooth(v[:, k].copy(), + (scale / (2 * np.pi)) ** 2) + for k in range(v.shape[1])]) + packed = v[surf.polys].mean(1) + + sxx, sxy, syy = packed[:, 0], packed[:, 1], packed[:, 2] + tr = sxx + syy + disc = np.sqrt(np.maximum((sxx - syy) ** 2 + 4 * sxy ** 2, 0.0)) + ok = tr > 1e-20 + return float((disc[ok]).sum() / tr[ok].sum()) + + +def test_the_elongation_measure_tells_ridges_from_knobs(): + """Pin the instrument before trusting what it says about cortex.""" + g = np.linspace(0, 48, 120) + X, Y = np.meshgrid(g, g) + xy = np.column_stack([X.ravel(), Y.ravel()]) + surf = polyutils.Surface(np.column_stack([xy, np.zeros(len(xy))]), + Delaunay(xy).simplices) + x, y = xy[:, 0], xy[:, 1] + + ridges = np.sin(2 * np.pi * x / 8.0) + knobs = ridges * np.sin(2 * np.pi * y / 8.0) + assert _elongation(surf, ridges) > 0.9 + assert _elongation(surf, knobs) < 0.4 + + +def test_a_ridge_stays_a_ridge(): + """The whole point: an elongated thickness ridge must not come out beaded. + + Every other synthetic case in this file is a square grid stretched equally + in both directions, so none of them can tell an elongated relief from a + chain of round bumps -- which is exactly the defect this guards. The slab is + given a thickness ridge running along y, and the relief has to come back at + least as elongated as an honest measurement of that ridge. + """ + n, spacing = 41, 1.0 + flat, wm, pia, polys, index = slab_grid(n=n, spacing=spacing, thickness=2.0) + + # a ridge along y: thicker in a 6 mm-wide band, uniform down its length + x = wm[:, 0] - wm[:, 0].mean() + bump = 1.2 * np.exp(-(x / 3.0) ** 2) + pia = pia + np.column_stack([np.zeros(len(x)), np.zeros(len(x)), bump]) + + surf = polyutils.Surface(bumpy._flat_plane(flat), polys) + relief = bumpy.FlatSlab(flat, wm, pia, polys).relaxed[:, 2] + + assert _elongation(surf, relief) > 0.85, ( + "the relief lost the ridge's direction; it is beaded rather than " + "elongated") + # and it has to still be a ridge in the right place, not just elongated + assert np.corrcoef(relief, bump)[0, 1] > 0.8 + + +def _curved_ridge_slab(n=41, thickness=2.0, height=2.0, width=5.0): + """A slab of *exactly* constant thickness, folded into a ridge along y. + + The pia is offset along the surface normal, so on the convex crown it has + more area than the white matter beneath it and in the flanks less. Thickness + is uniform to machine precision, so any relief this produces comes from the + folding and from nothing else. + """ + g = np.arange(n) * 1.0 + X, Y = np.meshgrid(g, g, indexing='ij') + z = height * np.exp(-((X - g[-1] / 2) / width) ** 2) + wm = np.stack([X.ravel(), Y.ravel(), z.ravel()], axis=1) + + gx, gy = np.gradient(z, 1.0, 1.0, axis=(0, 1)) + normal = np.stack([-gx.ravel(), -gy.ravel(), np.ones(n * n)], axis=1) + normal /= np.linalg.norm(normal, axis=1)[:, None] + pia = wm + thickness * normal + + index = np.arange(n * n).reshape(n, n) + polys = [] + for i in range(n - 1): + for j in range(n - 1): + a, b = index[i, j], index[i + 1, j] + c, d = index[i + 1, j + 1], index[i, j + 1] + polys += [[a, b, c], [a, c, d]] + flat = np.column_stack([X.ravel(), Y.ravel(), np.zeros(n * n)]) + return flat, wm, pia, np.array(polys), z.ravel() + + +def test_folding_alone_produces_relief(): + """Constant thickness, folded: there must still be a bump on the crown. + + This is the property the relief was missing. Cortical thickness is a fairly + blobby field and on its own it gives a relief of round knobs; what makes a + flatmap look like gyri is the pial flare, the pia carrying more area than + the white matter beneath a crown. Here thickness is uniform to machine + precision, so the flare is the *only* signal available and a quantity that + ignores it would return a flat sheet. + """ + flat, wm, pia, polys, crown = _curved_ridge_slab() + + thickness = np.linalg.norm(pia - wm, axis=1) + assert thickness.std() < 1e-12, "the fixture is supposed to be uniform" + + height = bumpy.folding_height(flat, wm, pia, polys) + assert height.std() > 0.02 * thickness.mean() + assert np.corrcoef(height, crown)[0, 1] > 0.6 + + +def test_folding_height_ignores_the_flatmap_distortion(): + """The point of the new denominator, stated as a property. + + `folding_height` divides by the *folded* white matter area, so distorting + the flatmap must not change it -- the flatmap enters only as the mesh the + regularisation is solved on. `naive_prism_height` divides by the + *flattened* area, so the same distortion moves it a great deal. That + difference is the whole reason for the change: a flatmap's area distortion + measures essentially uncorrelated with curvature, so as a denominator it + contributes no folding and injects the flattening algorithm's artifacts + instead. + """ + flat, wm, pia, polys, _ = _curved_ridge_slab(n=31) + + # a smooth, folding-unrelated area distortion of the flatmap + warped = flat.copy() + warped[:, 0] *= 1.0 + 0.3 * np.sin(2 * np.pi * flat[:, 1] / 30.0) + + def shift(fn): + a, b = fn(flat), fn(warped) + return np.abs(b - a).mean() / a.mean() + + folding = shift(lambda f: bumpy.folding_height(f, wm, pia, polys)) + prism = shift(lambda f: bumpy.naive_prism_height(f, wm, pia, polys)) + assert folding < 0.02 + assert prism > 5 * folding + + +def test_folding_height_is_the_frustum_over_the_white_area(): + """Pin the algebra: the height really is V_frustum / A_wm.""" + flat, wm, pia, polys, _ = _curved_ridge_slab(n=21) + + awm = bumpy._lumped(polyutils.face_area(wm[polys]), polys, len(wm)) + apia = bumpy._lumped(polyutils.face_area(pia[polys]), polys, len(wm)) + r = np.sqrt(apia / awm) + thickness = np.linalg.norm(pia - wm, axis=1) + expected = thickness * (1 + r + r ** 2) / 3.0 + + # unregularised comparison, so use a correlation length short enough that + # the smoothing is not what is being tested + got = bumpy.folding_height(flat, wm, pia, polys, correlation_length=1e-4) + np.testing.assert_allclose(got, expected, rtol=2e-3) diff --git a/cortex/tests/test_webgl_headless.py b/cortex/tests/test_webgl_headless.py index f97340ebc..d4753ed4a 100644 --- a/cortex/tests/test_webgl_headless.py +++ b/cortex/tests/test_webgl_headless.py @@ -72,6 +72,24 @@ def _wait_for_file(path, timeout=30): raise RuntimeError(f"File {path!r} not written within {timeout}s") +def _assert_not_blank(path): + """Fail if the render came out as a single flat color. + + A shader that fails to compile or link (or geometry that never made it to + the GPU) leaves the canvas showing nothing but the background, which is + otherwise indistinguishable from a successful render: the png is written, + and no javascript exception is raised. + """ + from PIL import Image + + rgb = np.asarray(Image.open(path).convert("RGB")).reshape(-1, 3) + ncolors = len(np.unique(rgb, axis=0)) + assert ncolors > 10, ( + f"{path} has only {ncolors} distinct color(s); the brain was probably " + "never drawn." + ) + + # --------------------------------------------------------------------------- # Group 1: Data type smoke tests # --------------------------------------------------------------------------- @@ -90,6 +108,7 @@ def test_datatype_renders(dtype_name, tmp_path): _wait_for_file(outfile) assert os.path.isfile(outfile) assert os.path.getsize(outfile) > 0 + _assert_not_blank(outfile) # No uncaught JS errors pageerrors = [e for e in handle._pw_thread.browser_errors if "[pageerror]" in e] assert len(pageerrors) == 0, f"JS errors: {pageerrors}" @@ -1054,3 +1073,88 @@ def _bump(surf, seed, sigma): print(f"\nVisual comparison saved to:\n {out_path}\n") assert out_path.exists() assert out_path.stat().st_size > 0 + + +# --------------------------------------------------------------------------- +# Group 6: Bumpy flatmap +# --------------------------------------------------------------------------- + + +def _ensure_bumpy_flatmap(): + """Put a bumpy flatmap in the database for the test subject if there is + none. It is a few seconds per hemisphere, so just ask for it.""" + cortex.db.get_surfinfo(subj, type='bumpy_flatmap').close() + + +@pytest.mark.timeout(900) +def test_bumpy_flatmap_changes_the_render(tmp_path): + """The relief reaches the shader and visibly changes the flatmap. + + This is the end-to-end check on the whole path: the height computed in + `cortex.polyutils.FlatSlab`, the cached surface info, the ``flatoffset`` + attribute in the ctm, the height javascript packs into ``flatbump.w``, and + the displacement the vertex shader applies. Any break in that chain shows + up here as two identical images. + """ + from PIL import Image + + _ensure_bumpy_flatmap() + # The pack may predate the offsets, in which case it has no flatoffset map. + cortex.utils.get_ctmpack(subj, recache=True) + + vol = cortex.Volume(np.random.randn(*volshape), subj, xfmname) + original = cortex.options.config.get("webgl_viewopts", "bumpy_flatmap") + images = {} + + def capture(handle, name): + outfile = str(tmp_path / ("%s.png" % name)) + handle.getImage(outfile, (512, 384)) + _wait_for_file(outfile) + _assert_not_blank(outfile) + pageerrors = [e for e in handle._pw_thread.browser_errors + if "[pageerror]" in e] + assert len(pageerrors) == 0, f"JS errors: {pageerrors}" + images[name] = np.asarray(Image.open(outfile).convert("RGB")) + + try: + for bumpy in (False, True): + cortex.options.config.set("webgl_viewopts", "bumpy_flatmap", + "true" if bumpy else "false") + with cortex.export.headless_viewer(vol, viewer_params={}) as handle: + handle._set_view(**{**default_view_params, + **unfold_view_params["inflated"]}) + capture(handle, "inflated_%s" % bumpy) + handle._set_view(**{**default_view_params, + **unfold_view_params["flatmap"]}) + capture(handle, "bumpy_%s" % bumpy) + if bumpy: + # And the same viewer with the relief exaggerated, which is + # what the bumpy_flatmap_scale slider drives. + # Only the rendered image is checked, not a read-back of + # the value: reading any surface menu property through the + # javascript proxy returns an empty dict, for unfold and + # depth just as much as for this one. + handle.ui.set("surface.%s.bumpy_flatmap_scale" % subj, 4.0) + time.sleep(0.3) + capture(handle, "bumpy_scaled") + finally: + cortex.options.config.set("webgl_viewopts", "bumpy_flatmap", original) + + assert not np.array_equal(images["bumpy_True"], images["bumpy_False"]), ( + "the bumpy flatmap rendered identically to the flat one; the relief " + "never reached the shader" + ) + assert not np.array_equal(images["bumpy_scaled"], images["bumpy_True"]), ( + "exaggerating the relief changed nothing; the bumpy_flatmap_scale " + "slider is not reaching the shader" + ) + # The offsets are in flatmap coordinates, so they mean nothing until the + # surface is flat: the displacement ramps in over inflated-to-flat, and at + # the inflated state it must not have started. This once regressed the + # other way -- the displacement ramped over anatomical-to-inflated while + # the shading normal ramped over inflated-to-flat, so geometry and + # lighting disagreed across the whole first half of the unfold. + assert np.array_equal(images["inflated_True"], images["inflated_False"]), ( + "the bumpy flatmap changed the inflated surface; a flatmap-space " + "offset is leaking into the folded surfaces" + ) diff --git a/cortex/tests/test_webgl_shaders.py b/cortex/tests/test_webgl_shaders.py new file mode 100644 index 000000000..9d4d8df25 --- /dev/null +++ b/cortex/tests/test_webgl_shaders.py @@ -0,0 +1,195 @@ +"""Tests that every webgl shader variant compiles and links. + +WebGL only guarantees 16 vertex attribute slots (``MAX_VERTEX_ATTRIBS``), and +the surface shaders use nearly all of them. A shader that asks for one slot too +many still compiles: it fails at *link* time, which three.js only reports on +the browser console and which shows up in the viewer as an unexplained black +screen. That is how ``Vertex2D`` data broke (gh-714), so every combination of +options the viewer generates shaders with is linked here. + +These tests only need Chromium; no subject database or viewer is involved. +""" + +import json +import os + +import pytest + +import cortex.webgl +from cortex.tests.testing_utils import has_playwright + +pytestmark = pytest.mark.skipif( + not has_playwright, reason="playwright and chromium are required" +) + +JS_PATH = os.path.join(os.path.dirname(cortex.webgl.__file__), "resources", "js") + +# The declarations THREE.WebGLProgram prepends to every shader it builds +# (three.js r69, resources/js/three.js). Only the ones the surface shaders +# actually rely on are listed; a missing one shows up as a compile error rather +# than as a silently passing test. +VERTEX_PREFIX = """ +precision highp float; +precision highp int; +#define MAX_DIR_LIGHTS 3 +#define MAX_POINT_LIGHTS 0 +#define MAX_SPOT_LIGHTS 0 +#define MAX_HEMI_LIGHTS 0 +#define MAX_SHADOWS 0 +uniform mat4 modelMatrix; +uniform mat4 modelViewMatrix; +uniform mat4 projectionMatrix; +uniform mat4 viewMatrix; +uniform mat3 normalMatrix; +uniform vec3 cameraPosition; +attribute vec3 position; +attribute vec3 normal; +attribute vec2 uv; +attribute vec2 uv2; +""" + +FRAGMENT_PREFIX = """ +precision highp float; +precision highp int; +#define MAX_DIR_LIGHTS 3 +#define MAX_POINT_LIGHTS 0 +#define MAX_SPOT_LIGHTS 0 +#define MAX_HEMI_LIGHTS 0 +#define MAX_SHADOWS 0 +uniform mat4 viewMatrix; +uniform vec3 cameraPosition; +""" + +# Loads the viewer's shader library into a page and exposes a hook that builds +# one shader variant and links it, the way THREE.WebGLProgram does. +PAGE = """ + + + + +""" + +# The options the viewer generates surface shaders with. ``morphs`` is the +# number of surfaces to mix between (anatomical, inflated and flat), ``volume`` +# says the subject has a white matter surface; the rest come from the dataview +# and from the surface menu. +SURFACE_OPTS = dict(morphs=3, volume=1, layers=1, rois=True, extratex=False, + halo=False, dither=False, voxline=False, sampler="nearest") + + +def _surface_variants(): + """Every (shader, opts) pair the viewer can ask for a surface shader.""" + for shader in ("surface_vertex", "surface_pixel"): + for rgb in (False, True): + for twod in (False, True): + if rgb and twod: + continue # RGB data has no second dimension + for hasflat in (False, True): + for equivolume in (False, True): + opts = dict(SURFACE_OPTS, rgb=rgb, twod=twod, + hasflat=hasflat, equivolume=equivolume) + name = "%s-%s%s%s%s" % ( + shader, + "rgb" if rgb else "cmap", + "-2d" if twod else "", + "-flat" if hasflat else "", + "-equivolume" if equivolume else "", + ) + yield pytest.param(shader, opts, id=name) + + +def _variants(): + yield from _surface_variants() + # The shaders the picker renders with; they morph the same geometry but + # carry no data. + yield pytest.param("pick", dict(morphs=3, volume=1), id="pick") + yield pytest.param("depth", dict(morphs=3, volume=1), id="depth") + + +@pytest.fixture(scope="module") +def link_shader(tmp_path_factory): + """Return a function linking one shader variant in a real GL context.""" + from playwright.sync_api import sync_playwright + + page_path = tmp_path_factory.mktemp("shaders") / "shaders.html" + page_path.write_text( + PAGE.replace("__JSDIR__", JS_PATH) + .replace("__VERTEX_PREFIX__", json.dumps(VERTEX_PREFIX)) + .replace("__FRAGMENT_PREFIX__", json.dumps(FRAGMENT_PREFIX)) + ) + + with sync_playwright() as playwright: + browser = playwright.chromium.launch( + headless=True, + args=["--enable-webgl", "--use-gl=swiftshader", "--no-sandbox", + "--disable-dev-shm-usage"], + ) + page = browser.new_page() + page.goto("file://%s" % page_path, wait_until="load", timeout=60000) + if not page.evaluate("() => !!window.linkShader"): + browser.close() + pytest.skip("no WebGL context available in this browser") + yield lambda shader, opts: page.evaluate( + "args => window.linkShader(args[0], args[1])", [shader, opts] + ) + browser.close() + + +@pytest.mark.parametrize("shader,opts", list(_variants())) +def test_shader_links(shader, opts, link_shader): + """Each shader variant has to compile *and* link. + + A variant that uses more vertex attributes than the driver has slots for + compiles fine and fails to link, which leaves the viewer showing nothing at + all. + """ + result = link_shader(shader, opts) + assert result["compiled"], "%s did not compile:\n%s" % (shader, result["log"]) + assert result["linked"], ( + "%s compiled but did not link, using %d of the %d available vertex " + "attributes:\n%s" % (shader, len(result["attributes"]), + result["max_attributes"], result["log"]) + ) + assert len(result["attributes"]) <= result["max_attributes"] diff --git a/cortex/utils.py b/cortex/utils.py index 8d231d3d8..dcfe5ccb5 100644 --- a/cortex/utils.py +++ b/cortex/utils.py @@ -113,7 +113,10 @@ def get_ctmpack(subject, types=("inflated",), method="raw", level=0, recache=Fal """ lvlstr = ("%dd" if decimate else "%d")%level # Generates different cache files for each combination of disp_layers - ctmcache = "%s_[{types}]_{method}_{level}_v3.json"%subject + # v4: the ctm now also carries the equivolume vertex areas and the bumpy + # flatmap offsets, so packs cached by an older version are missing data the + # shaders expect. + ctmcache = "%s_[{types}]_{method}_{level}_v4.json"%subject ctmcache = ctmcache.format(types=','.join(types), method=method, level=lvlstr) diff --git a/cortex/webgl/resources/js/mriview_surface.js b/cortex/webgl/resources/js/mriview_surface.js index fdbea89b2..e2d4a20bc 100644 --- a/cortex/webgl/resources/js/mriview_surface.js +++ b/cortex/webgl/resources/js/mriview_surface.js @@ -64,6 +64,13 @@ var mriview = (function(module) { this.object = new THREE.Group(); this.object.name = 'Surface'; + //Exaggeration of the bumpy flatmap's relief. A missing or unparseable + //setting falls back to 1, the true scale; zero is a legitimate value + //and has to survive, so this cannot just be `|| 1`. + var bumpscale = parseFloat(viewopts.bumpy_flatmap_scale); + if (isNaN(bumpscale)) + bumpscale = 1.0; + this.uniforms = THREE.UniformsUtils.merge( [ THREE.UniformsLib[ "lights" ], { @@ -76,6 +83,7 @@ var mriview = (function(module) { thickmix: { type:'f', value:0.5}, surfmix: { type:'f', value:0}, bumpyflat: { type:'i', value:viewopts.bumpy_flatmap == 'true'}, + bumpyflat_scale: { type:'f', value:bumpscale}, allowtilt: { type:'i', value:viewopts.allow_tilt == 'true'}, // equivolume: { type:'i', value:viewopts.equivolume == 'true'}, @@ -108,6 +116,12 @@ var mriview = (function(module) { "fiducial surface": {action: this.to_fiducial_surface.bind(this), key: 'u', help: "Fiducial surface"}, "WM surface": {action: this.to_white_matter_surface.bind(this), key: 'y', help: "White matter surface"}, bumpy_flatmap: {action:[this, "setBumpyFlat"]}, + //The slider starts wherever the config file put it. Its range goes + //to 5x true scale, or to twice the configured value if that is + //already higher, so a deliberately large setting is not clamped + //away the first time the menu is drawn. + bumpy_flatmap_scale: {action:[this.uniforms.bumpyflat_scale, "value", + 0, Math.max(5, 2 * bumpscale)]}, allow_tilt: {action:[this, "setAllowTilt"]}, equivolume: {action:[this, "setEquivolume"]}, changeDepth: {action: this.changeDepth.bind(this), wheel: true, modKeys: ['altKey'], hidden: true, help:'Change depth'}, @@ -168,17 +182,6 @@ var mriview = (function(module) { right:{positions:[], normals:[]}, }; - // Smoothing parameters for the bumpy flatmap. These must be declared - // outside the per-hemisphere loop below: they used to live inside it, - // and `var` hoisting meant they were still undefined on the first - // iteration, so the left hemisphere was left unsmoothed while the - // right one picked up the values assigned during the left pass. - var areasmoothfactor = 0.1; - var areasmoothiter = 5; - - var distsmoothfactor = 0.1; - var distsmoothiter = 20; - for (var name in names) { var hemi = geometries[names[name]]; posdata[name].map = hemi.indexMap; @@ -201,23 +204,15 @@ var mriview = (function(module) { posdata[name].normals.push(hemi.attributes['mixNorms'+i]) delete hemi.attributes[json.names[i]]; } - //Setup flatmap mix - var wmareas = module.computeAreas(hemi.attributes.wm, hemi.attributes.index, hemi.offsets); - wmareas = module.iterativelySmoothVertexData(hemi.attributes.wm, hemi.attributes.index, hemi.offsets, wmareas, areasmoothfactor, areasmoothiter); - hemi.wmareas = wmareas; - - var pialareas = module.computeAreas(hemi.attributes.position, hemi.attributes.index, hemi.offsets); - pialareas = module.iterativelySmoothVertexData(hemi.attributes.position, hemi.attributes.index, hemi.offsets, pialareas, areasmoothfactor, areasmoothiter); - hemi.pialareas = pialareas; - - var pialarea_attr = new THREE.BufferAttribute(pialareas, 1); - pialarea_attr.needsUpdate = true; - var wmarea_attr = new THREE.BufferAttribute(wmareas, 1); - wmarea_attr.needsUpdate = true; - - hemi.addAttribute('pialarea', pialarea_attr); - hemi.addAttribute('wmarea', wmarea_attr); + //The white matter and pial vertex areas that equivolume depth + //sampling needs already arrived in auxdat.zw -- x is the medial + //wall mask and y the curvature -- computed and smoothed in + //python by cortex.surfinfo.equivolume_areas. They ride in the + //spare components of an attribute that is already here because + //WebGL only guarantees 16 vertex attribute slots and these + //shaders use every one of them. + //Setup flatmap mix if (this.flatlims !== undefined) { var flats = this._makeFlat(hemi.attributes.uv.array, json.flatlims, names[name]); hemi.addAttribute('mixSurfs'+json.names.length, new THREE.BufferAttribute(flats.pos, 4)); @@ -227,51 +222,55 @@ var mriview = (function(module) { posdata[name].positions.push(hemi.attributes['mixSurfs'+json.names.length]); posdata[name].normals.push(hemi.attributes['mixNorms'+json.names.length]); - // var flatareas = module.computeAreas(hemi.attributes.mixSurfs1, hemi.culled.index, hemi.culled.offsets); - // var flatareascale = flatscale ** 2; - // flatareas = flatareas.map(function (a) { return a / flatareascale;}); - // flatareas = module.iterativelySmoothVertexData(hemi.attributes.position, hemi.attributes.index, hemi.offsets, flatareas, smoothfactor, smoothiter); - // hemi.flatareas = flatareas; - - var dists = module.computeDist(hemi.attributes.position, hemi.attributes.wm); - dists = module.iterativelySmoothVertexData(hemi.attributes.position, hemi.attributes.index, hemi.offsets, dists, distsmoothfactor, distsmoothiter); - - var vertexvolumes = module.computeVertexPrismVolume(wmareas, pialareas, dists); - hemi.vertexvolumes = vertexvolumes; - var flatheights = module.computeFlatVolumeHeight(wmareas, vertexvolumes); - flatheights.array = flatheights.array.map(function (h) {return h * flatscale;}); - // flatheights.array = flatheights.array.map(Math.sqrt); - - var flat_offset_verts; - if ( name == "left" ) { - flat_offset_verts = module.offsetVerts(hemi.attributes.mixSurfs1, flatheights, 0, -1); - } else { - flat_offset_verts = module.offsetVerts(hemi.attributes.mixSurfs1, flatheights, 0, 1); + //The bumpy flatmap: the height of the pial surface above + //the flat white matter surface, computed in python (see + //cortex.polyutils.FlatSlab). The relief is purely vertical + //and the flatmap plane is (y, z) in viewer space, so the + //height goes along x. + var flatsurf = hemi.attributes['mixSurfs'+json.names.length]; + var nverts = flatsurf.array.length / 4; + //Absent for a subject with no white matter surface, and + //for a ctm cached before this moved into python; either way + //the flatmap just stays flat. + var offsets = (hemi.attributes.wm !== undefined) + ? hemi.attributes.flatoffset : undefined; + var mirror = (name == "left") ? -1 : 1; + + var displaced = new Float32Array(flatsurf.array); + var height = new Float32Array(nverts); + for (var v = 0; v < nverts; v++) { + if (offsets !== undefined) + height[v] = mirror * flatscale * offsets.array[v*4+2]; + displaced[v*4] += height[v]; + } + + //The medial wall is in no triangle and would get a zero + //normal, which the shader's normalize turns into a NaN. + //Fall back to the sheet's own normal, read off a real + //triangle so the sign follows the winding. + var flatnorm = module.flatSheetNormal( + flatsurf.array, 4, hemi.attributes.index, hemi.offsets); + var bumpnorms = module.computeNormal( + new THREE.BufferAttribute(displaced, 4), + hemi.attributes.index, hemi.offsets, flatnorm); + //flatbump.xyz is the shading normal and its w the height, + //so the whole relief rides in one attribute -- these + //shaders use all 16 slots WebGL guarantees. + var flatbump = new Float32Array(nverts * 4); + for (var v = 0; v < nverts; v++) { + flatbump[v*4] = bumpnorms.array[v*3]; + flatbump[v*4+1] = bumpnorms.array[v*3+1]; + flatbump[v*4+2] = bumpnorms.array[v*3+2]; + flatbump[v*4+3] = height[v]; } - // // hemi.addAttribute('offsetflat', flat_offset_verts); - // var flatoff_geom = new THREE.BufferGeometry(); - // flatoff_geom.addAttribute('position', flat_offset_verts); - // flatoff_geom.addAttribute('index', hemi.attributes.index); - // flatoff_geom.computeVertexNormals(); - - // console.log(flatoff_geom); - // this.flatoff = flatoff_geom; - - // hemi.addAttribute('flatBumpNorms', flatoff_geom.attributes.normal); - hemi.addAttribute('flatheight', flatheights); - hemi.addAttribute('flatBumpNorms', module.computeNormal(flat_offset_verts, hemi.attributes.index, hemi.offsets) ); - } else { - // Fill these attributes so the shader doesn't choke, even though - // there's no flatmap - // just set flatheight to 1 everywhere. - // var flatheight_arr = new Float32Array(hemi.attributes.position.position / hemi.attributes.position.itemSize); - // flatheight_arr = flatheight_arr.map(function (x) {return 1.0;}); - // var flatheight = new THREE.BufferAttribute(flatheight_arr, 1); - // hemi.addAttribute('flatheight', flatheight); - - // // and set the flatBumpNorms to the - // var flatBumpNorms = module.computeNormal() + + var flatbump_attr = new THREE.BufferAttribute(flatbump, 4); + flatbump_attr.needsUpdate = true; + hemi.addAttribute('flatbump', flatbump_attr); } + //With no flatmap there is nothing to add: the shader only + //declares flatbump under #ifdef HASFLAT, so it never looks for + //an attribute that is not here. //Generate an index list that has culled non-flatmap vertices var culled = module._cull_flatmap_vertices(hemi.attributes.index.array, hemi.attributes.auxdat.array, hemi.offsets); diff --git a/cortex/webgl/resources/js/mriview_utils.js b/cortex/webgl/resources/js/mriview_utils.js index 04f6b2a41..fd5f5ae58 100644 --- a/cortex/webgl/resources/js/mriview_utils.js +++ b/cortex/webgl/resources/js/mriview_utils.js @@ -210,7 +210,39 @@ var mriview = (function(module) { return {pos:pos, norm:norm, base:base}; } - module.computeNormal = function(vertices, index, offsets) { + //The unit normal of a planar sheet, taken from its first triangle with any + //area. Every triangle of a flatmap gives the same answer up to the sign, + //and reading it off the mesh gets that sign from the winding. + module.flatSheetNormal = function(positions, stride, index, offsets) { + var indices = index.array; + var pA = new THREE.Vector3(), pB = new THREE.Vector3(), + pC = new THREE.Vector3(), cb = new THREE.Vector3(), + ab = new THREE.Vector3(); + + for (var j = 0; j < offsets.length; j++) { + var start = offsets[j].start, count = offsets[j].count, + base = offsets[j].index; + for (var i = start; i < start + count; i += 3) { + var vA = base + indices[i], vB = base + indices[i+1], + vC = base + indices[i+2]; + pA.set(positions[vA*stride], positions[vA*stride+1], positions[vA*stride+2]); + pB.set(positions[vB*stride], positions[vB*stride+1], positions[vB*stride+2]); + pC.set(positions[vC*stride], positions[vC*stride+1], positions[vC*stride+2]); + cb.subVectors(pC, pB); + ab.subVectors(pA, pB); + cb.cross(ab); + if (cb.length() > 1e-8) + return cb.normalize().toArray(); + } + } + return [0, 0, 0]; + } + + //`fallback`, if given, is a 3-element array used for vertices whose + //accumulated face normals vanish -- a vertex in no triangle at all, or one + //whose incident triangles have zero area. Without it those vertices get a + //zero normal, which the fragment shader then normalizes into a NaN. + module.computeNormal = function(vertices, index, offsets, fallback) { var i, il; var j, jl; @@ -290,9 +322,9 @@ var mriview = (function(module) { normals[ i + 2 ] *= n; if (isNaN(normals[i])) { - normals[i] = 0; - normals[i+1] = 0; - normals[i+2] = 0; + normals[i] = fallback === undefined ? 0 : fallback[0]; + normals[i+1] = fallback === undefined ? 0 : fallback[1]; + normals[i+2] = fallback === undefined ? 0 : fallback[2]; } } @@ -301,205 +333,6 @@ var mriview = (function(module) { return attr; } - module.computeAreas = function(vertices, index, offsets) { - var i, il; - var j, jl; - - var indices = index.array; - var positions = vertices.array; - var stride = vertices.itemSize; - - var areas = new Float32Array( vertices.array.length / vertices.itemSize ); - - var vA, vB, vC, x, y, z, area, - pA = new THREE.Vector3(), - pB = new THREE.Vector3(), - pC = new THREE.Vector3(), - - cb = new THREE.Vector3(), - ab = new THREE.Vector3(), - - tri = new THREE.Triangle(); - - for ( j = 0, jl = offsets.length; j < jl; ++ j ) { - - var start = offsets[ j ].start; - var count = offsets[ j ].count; - var index = offsets[ j ].index; - - for ( i = start, il = start + count; i < il; i += 3 ) { - - vA = index + indices[ i ]; - vB = index + indices[ i + 1 ]; - vC = index + indices[ i + 2 ]; - - x = positions[ vA * stride ]; - y = positions[ vA * stride + 1 ]; - z = positions[ vA * stride + 2 ]; - pA.set( x, y, z ); - - x = positions[ vB * stride ]; - y = positions[ vB * stride + 1 ]; - z = positions[ vB * stride + 2 ]; - pB.set( x, y, z ); - - x = positions[ vC * stride ]; - y = positions[ vC * stride + 1 ]; - z = positions[ vC * stride + 2 ]; - pC.set( x, y, z ); - - tri.set( pA, pB, pC ); - varea = tri.area() / 3; - - - areas[vA] += varea; - areas[vB] += varea; - areas[vC] += varea; - - } - } - - return areas; - } - - module.smoothVertexData = function(vertices, index, offsets, data, factor) { - var i, il; - var j, jl; - - var indices = index.array; - var positions = vertices.array; - var stride = vertices.itemSize; - - var smoothdata = new Float32Array( data.length ), - counts = new Uint16Array( data.length ); - - var vA, vB, vC, x, y, z, area, - pA = new THREE.Vector3(), - pB = new THREE.Vector3(), - pC = new THREE.Vector3(), - - cb = new THREE.Vector3(), - ab = new THREE.Vector3(), - - tri = new THREE.Triangle(); - - for ( j = 0, jl = offsets.length; j < jl; ++ j ) { - - var start = offsets[ j ].start; - var count = offsets[ j ].count; - var index = offsets[ j ].index; - - for ( i = start, il = start + count; i < il; i += 3 ) { - - vA = index + indices[ i ]; - vB = index + indices[ i + 1 ]; - vC = index + indices[ i + 2 ]; - - // x = positions[ vA * stride ]; - // y = positions[ vA * stride + 1 ]; - // z = positions[ vA * stride + 2 ]; - // pA.set( x, y, z ); - - // x = positions[ vB * stride ]; - // y = positions[ vB * stride + 1 ]; - // z = positions[ vB * stride + 2 ]; - // pB.set( x, y, z ); - - // x = positions[ vC * stride ]; - // y = positions[ vC * stride + 1 ]; - // z = positions[ vC * stride + 2 ]; - // pC.set( x, y, z ); - - // tri.set( pA, pB, pC ); - // varea = tri.area() / 3; - - smoothdata[vA] += data[vB] + data[vC]; - smoothdata[vB] += data[vA] + data[vC]; - smoothdata[vC] += data[vA] + data[vB]; - - counts[vA] += 2; - counts[vB] += 2; - counts[vC] += 2; - - } - } - - for ( var k = 0; k < data.length; k ++ ) { - smoothdata[k] = smoothdata[k] / counts[k] * factor + (data[k] * (1.0 - factor)); - } - - return smoothdata; - } - - module.iterativelySmoothVertexData = function(vertices, index, offsets, data, factor, iters) { - var smoothed = data; - for ( var i = 0; i < iters; i ++ ) { - smoothed = module.smoothVertexData(vertices, index, offsets, smoothed, factor); - } - return smoothed; - } - - module.computeVertexPrismVolume = function(areas1, areas2, dists) { - var num = areas1.length; - - var volumes = new Float32Array(num); - - for ( var i = 0; i < num; i ++ ) { - volumes[i] = dists[i] / 3 * ( areas1[i] + areas2[i] + Math.sqrt(areas1[i] * areas2[i])); - } - - return volumes; - } - - module.computeFlatVolumeHeight = function(areas, volumes) { - var num = areas.length; - var flatheights = new Float32Array(num); - - for ( var i = 0; i < num; i ++ ) { - flatheights[i] = volumes[i] / areas[i]; - } - - var attr = new THREE.BufferAttribute(flatheights, 1); - attr.needsUpdate = true; - return attr; - } - - module.computeDist = function(verts1, verts2) { - var stride1 = verts1.itemSize; - var stride2 = verts2.itemSize; - var num = verts1.length / stride1; - var dists = new Float32Array(num); - - var tv1 = new THREE.Vector3(), tv2 = new THREE.Vector3(); - - for ( var i = 0; i < num; i ++ ) { - tv1.set(verts1.array[i * stride1], - verts1.array[i * stride1 + 1], - verts1.array[i * stride1 + 2]); - tv2.set(verts2.array[i * stride2], - verts2.array[i * stride2 + 1], - verts2.array[i * stride2 + 2]); - dists[i] = tv1.distanceTo(tv2); - } - return dists; - // var attr = new THREE.BufferAttribute(dists, 1); - // attr.needsUpdate = true; - // return attr; - } - - module.offsetVerts = function(verts, offsets, axis, factor) { - var num = verts.length / verts.itemSize; - var newpos = new Float32Array(verts.array); - - for ( var i = 0; i < num; i ++ ) { - newpos[i*verts.itemSize + axis] += factor * offsets.array[i]; - } - - var attr = new THREE.BufferAttribute(newpos, verts.itemSize); - attr.needsUpdate = true; - return attr; - } - //Generates a hatch texture in canvas module.makeHatch = function(size, linewidth, spacing) { //Creates a cross-hatching pattern in canvas for the dropout shading diff --git a/cortex/webgl/resources/js/shaderlib.js b/cortex/webgl/resources/js/shaderlib.js index dc4deb82e..cea96af4a 100644 --- a/cortex/webgl/resources/js/shaderlib.js +++ b/cortex/webgl/resources/js/shaderlib.js @@ -174,19 +174,24 @@ var Shaderlib = (function() { return glsl; }, - // thickmixer: header code that loads the uniforms and attributes needed to - // do equivolume sampling, for vertex shaders + // thickmixer: header code that loads the uniforms needed to do + // equivolume sampling, for vertex shaders. The white matter and pial + // vertex areas it needs ride along in auxdat.zw, which mriview_surface + // fills in when the surfaces load, rather than in attributes of their + // own: WebGL only guarantees 16 vertex attribute slots and these + // shaders are right up against that limit, so whatever fits in the + // spare components of an attribute that is already there goes there. thickmixer: [ "uniform float thickmix;", "uniform int equivolume;", - "attribute float wmarea;", - "attribute float pialarea;", ].join("\n"), // thickmixer_main: translates a desired volume fraction into linear mixing - // parameter. + // parameter. Requires auxdat to be declared by the including shader. thickmixer_main: [ "#ifdef EQUIVOLUME", + "float wmarea = auxdat.z;", + "float pialarea = auxdat.w;", "float use_thickmix = 1. - (1. / (pialarea - wmarea) * (-1. * wmarea + sqrt((1. - thickmix) * pialarea * pialarea + thickmix * wmarea * wmarea)));", "#else", "float use_thickmix = thickmix;", @@ -361,6 +366,7 @@ var Shaderlib = (function() { // "uniform float thickmix;", utils.thickmixer, "uniform int bumpyflat;", + "uniform float bumpyflat_scale;", "float f_bumpyflat = float(bumpyflat);", "attribute vec4 wm;", @@ -368,8 +374,14 @@ var Shaderlib = (function() { "attribute vec4 auxdat;", "#ifdef HASFLAT", - "attribute vec3 flatBumpNorms;", - "attribute float flatheight;", + //xyz: normal of the bump-displaced flatmap, w: the x component + //of the bump offset. The other two components ride in the + //unused w of `wm` and of the flat morph target -- see + //`bumpvector` below. The surface shaders use all 16 of the + //vertex attribute slots webgl guarantees, so a bumpy flatmap + //has to fit in the spare components of attributes that are + //already there rather than adding one of its own. + "attribute vec4 flatbump;", "#endif", // "attribute float dropout;", @@ -425,20 +437,28 @@ var Shaderlib = (function() { "vec3 pos, norm;", "mixfunc(mpos, mnorm, pos, norm);", - // "norm = mix(flatBumpNorms, normalize(onorm), thickmix);", - // "norm = normalize(flatBumpNorms);", - "#ifdef CORTSHEET", // "#ifdef HASFLAT", - "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * mix(1., 0., use_thickmix) * flatheight * f_bumpyflat;", + //The relief is purely vertical and the flatmap's + //out-of-plane axis is x, so the scale setting is vertical + //exaggeration, as on a topographic map. Javascript bakes + //the per-hemisphere mirroring and the flatmap scale in. + "vec3 bumpvector = vec3(flatbump.w * bumpyflat_scale, 0., 0.);", + //Only over inflated-to-flat: a height in flatmap + //coordinates means nothing on a folded surface. + "pos += clamp(surfmix*"+(morphs-1)+". - "+(morphs-2)+"., 0., 1.) * mix(1., 0., use_thickmix) * f_bumpyflat * bumpvector;", "#else", "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * .62 * distance(position, wm.xyz) * mix(1., 0., use_thickmix);", "#endif", "#endif", "#ifdef HASFLAT", - "vNormal = normalMatrix * mix(norm, flatBumpNorms, (1.0 - use_thickmix) * clamp(surfmix*"+(morphs-1)+". - 1.0, 0., 1.) * f_bumpyflat);", + //Scaling a height field by s scales the normal's in-plane + //components by s and leaves the out-of-plane one alone -- + //exact, and at scale 0 it gives back the flat normal. + "vec3 bumpnorm = normalize(vec3(flatbump.x, bumpyflat_scale * flatbump.y, bumpyflat_scale * flatbump.z));", + "vNormal = normalMatrix * mix(norm, bumpnorm, (1.0 - use_thickmix) * clamp(surfmix*"+(morphs-1)+". - "+(morphs-2)+"., 0., 1.) * f_bumpyflat);", "#else", "vNormal = normalMatrix * norm;", "#endif", @@ -665,14 +685,9 @@ var Shaderlib = (function() { wm: { type: 'v4', value:null }, wmnorm: { type: 'v3', value:null }, auxdat: { type: 'v4', value:null }, - wmarea: { type: 'f', value:null }, - pialarea: { type: 'f', value:null }, - // flatBumpNorms: { type: 'v3', value:null }, - // flatheight: { type: 'f', value:null }, }; if (opts.hasflat) { - attributes.flatBumpNorms = { type: 'v3', value:null }; - attributes.flatheight = { type: 'f', value:null }; + attributes.flatbump = { type: 'v4', value:null }; } for (var i = 0; i < morphs-1; i++) { attributes['mixSurfs'+i] = { type:'v4', value:null}; @@ -711,6 +726,7 @@ var Shaderlib = (function() { // "uniform float thickmix;", utils.thickmixer, "uniform int bumpyflat;", + "uniform float bumpyflat_scale;", "float f_bumpyflat = float(bumpyflat);", "varying vec4 vColor;", @@ -730,8 +746,14 @@ var Shaderlib = (function() { "attribute vec4 auxdat;", "#ifdef HASFLAT", - "attribute vec3 flatBumpNorms;", - "attribute float flatheight;", + //xyz: normal of the bump-displaced flatmap, w: the x component + //of the bump offset. The other two components ride in the + //unused w of `wm` and of the flat morph target -- see + //`bumpvector` below. The surface shaders use all 16 of the + //vertex attribute slots webgl guarantees, so a bumpy flatmap + //has to fit in the spare components of attributes that are + //already there rather than adding one of its own. + "attribute vec4 flatbump;", "#endif", // "attribute float dropout;", @@ -788,14 +810,25 @@ var Shaderlib = (function() { "#ifdef CORTSHEET", "#ifdef HASFLAT", - "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * mix(1., 0., use_thickmix) * flatheight * f_bumpyflat;", + //The relief is purely vertical and the flatmap's + //out-of-plane axis is x, so the scale setting is vertical + //exaggeration, as on a topographic map. Javascript bakes + //the per-hemisphere mirroring and the flatmap scale in. + "vec3 bumpvector = vec3(flatbump.w * bumpyflat_scale, 0., 0.);", + //Only over inflated-to-flat: a height in flatmap + //coordinates means nothing on a folded surface. + "pos += clamp(surfmix*"+(morphs-1)+". - "+(morphs-2)+"., 0., 1.) * mix(1., 0., use_thickmix) * f_bumpyflat * bumpvector;", "#else", "pos += clamp(surfmix*"+(morphs-1)+"., 0., 1.) * normalize(norm) * .62 * distance(position, wm.xyz) * mix(1., 0., use_thickmix);", "#endif", "#endif", "#ifdef HASFLAT", - "vNormal = normalMatrix * mix(norm, flatBumpNorms, (1.0 - use_thickmix) * clamp(surfmix*"+(morphs-1)+". - 1.0, 0., 1.) * f_bumpyflat);", + //Scaling a height field by s scales the normal's in-plane + //components by s and leaves the out-of-plane one alone -- + //exact, and at scale 0 it gives back the flat normal. + "vec3 bumpnorm = normalize(vec3(flatbump.x, bumpyflat_scale * flatbump.y, bumpyflat_scale * flatbump.z));", + "vNormal = normalMatrix * mix(norm, bumpnorm, (1.0 - use_thickmix) * clamp(surfmix*"+(morphs-1)+". - "+(morphs-2)+"., 0., 1.) * f_bumpyflat);", "#else", "vNormal = normalMatrix * norm;", "#endif", @@ -888,13 +921,10 @@ var Shaderlib = (function() { wm: { type: 'v4', value:null }, wmnorm: { type: 'v3', value:null }, auxdat: { type: 'v4', value:null }, - wmarea: { type: 'f', value:null }, - pialarea: { type: 'f', value:null }, }; if (opts.hasflat) { - attributes.flatBumpNorms = { type: 'v3', value:null }; - attributes.flatheight = { type: 'f', value:null }; + attributes.flatbump = { type: 'v4', value:null }; } for (var i = 0; i < 4; i++) diff --git a/docs/api_reference_flat.rst b/docs/api_reference_flat.rst index ff356b9cb..284f7413c 100644 --- a/docs/api_reference_flat.rst +++ b/docs/api_reference_flat.rst @@ -175,6 +175,15 @@ polyutils Surface Distortion + FlatSlab + +.. autosummary:: + :toctree:generated/ + + coarsen_flat_mesh + lame_parameters + naive_prism_height + prolongation_matrix segment @@ -204,6 +213,8 @@ surfinfo thickness tissots_indicatrix flat_border + bumpy_flatmap + equivolume_areas utils diff --git a/docs/database.rst b/docs/database.rst index f47209093..859edb361 100644 --- a/docs/database.rst +++ b/docs/database.rst @@ -266,6 +266,10 @@ Surface info The filestore also manages several important quantifications about the surfaces. These include Tissot's Indicatrix and the flatmap surface distortion. There are stored in the ``/surface-info`` directory. This is also where the per-vertex curvature, sulcal depth and thickness imported from Freesurfer_ are stored (see :ref:`database-freesurfer-import`). Each file is an ``.npz`` file holding one array per hemisphere, under the keys ``left`` and ``right``, and can be loaded with ``cortex.db.get_surfinfo``. +Two of them cannot be returned as a ``Vertex`` object, because that conversion concatenates arrays stored under the keys ``left`` and ``right`` and assumes one value per vertex. ``bumpy_flatmap.npz`` holds the bumpy flatmap's relief (see :class:`cortex.polyutils.FlatSlab`) as a three-component offset per vertex, of which only the third is nonzero, under the keys ``bump_left`` and ``bump_right``, and ``equivolume_areas.npz`` holds four scalar maps rather than two -- the white matter and pial vertex areas the webgl viewer's equivolume depth sampling needs -- under ``wm_left``, ``wm_right``, ``pia_left`` and ``pia_right``. For these, ``cortex.db.get_surfinfo`` hands back the ``.npz`` file itself rather than a ``Vertex``; remember to close it. + +Anything computed from the flatmap -- the distortion maps, the flatmap border and the bumpy flatmap -- is deleted and regenerated when a new flatmap is imported with ``cortex.freesurfer.import_flat``, since it describes a flatmap that no longer exists. + Views ----- @@ -311,9 +315,11 @@ Here is an example entry into the filestore... ├── overlays.svg ├── rois.svg ├── surface-info + │ ├── bumpy_flatmap.npz │ ├── curvature.npz │   ├── distortion[dist_type=areal].npz │ ├── distortion[dist_type=metric].npz + │ ├── equivolume_areas.npz │ ├── sulcaldepth.npz │ └── thickness.npz ├── surfaces diff --git a/docs/userguide/webgl.rst b/docs/userguide/webgl.rst index 442e66539..82eabd2be 100644 --- a/docs/userguide/webgl.rst +++ b/docs/userguide/webgl.rst @@ -105,16 +105,18 @@ flatten flatten cortical surface Surface Controls **************** -======== ============================ -name description -======== ============================ -unfold level of unfolding -pivot angle between hemispheres -shift distance between hemispheres -depth cortical depth -left toggle left hemisphere -right toggle right hemisphere -======== ============================ +=================== ================================== +name description +=================== ================================== +unfold level of unfolding +pivot angle between hemispheres +shift distance between hemispheres +depth cortical depth +bumpy_flatmap give the flatmap relief, see below +bumpy_flatmap_scale how much to exaggerate that relief +left toggle left hemisphere +right toggle right hemisphere +=================== ================================== Lighting Controls @@ -139,6 +141,34 @@ that reads as shaded relief. Like ``pivot``, the sliders move to show the values in effect and can still be dragged afterwards; the next change to ``unfold`` or ``bumpy_flatmap`` drives them back from the configured defaults. +The relief is the shape a slab of cortex would take if it were peeled off the +white matter and laid flat: thicker where the pia carries more area than the +white matter beneath it, which is over a gyral crown, and thinner in a sulcal +fundus. See :class:`cortex.polyutils.FlatSlab` for the quantity and +:mod:`cortex.polyutils.bumpy` for why it is that one. It costs a few seconds +per hemisphere and is cached in the subject's database entry like any other +surface info. + +``bumpy_flatmap_scale`` exaggerates the relief, and is a slider in the surface +controls as well as a setting. At 1.0 the bumps are at their true scale, which +for a 2-5 mm slab is subtle next to a whole flatmap, so the default exaggerates +somewhat; larger values make the folding easier still to read. The relief is +purely vertical, so this is vertical exaggeration as on a topographic map and +turning it up does not move anything sideways relative to the data drawn +underneath. The shading follows, since a height +field scaled by *s* has a normal whose in-plane components scale by *s* and +whose out-of-plane component does not. The slider starts wherever the +configuration file put it, and runs from 0 to five times true scale -- or to +twice the configured value if that is already higher. It does nothing while +``bumpy_flatmap`` is off, since there is no relief to scale. Being a display +setting, it does not invalidate the cached geometry, so it can be dragged +around freely. + +The relief appears over the second half of the unfold, between the inflated +surface and the flatmap, and the anatomical and inflated surfaces are untouched +by it. The offsets are in the flatmap's own coordinates, which have no meaning +on a folded surface. + Overlay Controls **************** diff --git a/examples/surface_analyses/plot_bumpy_flatmap.py b/examples/surface_analyses/plot_bumpy_flatmap.py new file mode 100644 index 000000000..1464bc857 --- /dev/null +++ b/examples/surface_analyses/plot_bumpy_flatmap.py @@ -0,0 +1,113 @@ +""" +============== +Bumpy Flatmaps +============== + +A flatmap throws away the folding of the cortical surface, which is a shame: +where the gyri and sulci were is genuinely useful context, and once a flatmap is +covered in data there is nothing left to show it. The usual fix is to shade the +flatmap by curvature underneath the data, but that only works where the data is +transparent. + +The bumpy flatmap puts the folding back as relief instead of as color. Cortex is +2 to 5 mm thick, so imagine peeling the cortical slab off the white matter, +making the flattening cuts, and laying it down. The white matter side ends up +flat, and the pial side sits some distance above it. + +How far above turns out to depend a great deal on which area you divide the +column's volume by. The *flattened* area is the more obvious choice and it does +not work: a flatmap's area distortion measures essentially uncorrelated with +curvature, so that denominator contributes no folding and injects the +flattening algorithm's artifacts instead, besides giving enormous heights +wherever flattening crushed a triangle. What is left is close to a map of +cortical thickness, which is blobby and reads as round knobs. + +`cortex.polyutils.folding_height` divides by the *folded* white matter area +instead, which with ``r = sqrt(A_pia / A_wm)`` is ``thickness * (1 + r + r**2) +/ 3``. `r` is what carries the folding: the pia has more area than the white +matter beneath a gyral crown and less in a fundus. + +The naive volume-over-flattened-area field is plotted below for comparison. +""" + +import numpy as np +import matplotlib.pyplot as plt + +import cortex +from cortex.polyutils import naive_prism_height + +subject = "S1" + +# Cached in the database, so this is a lookup; computing it takes a few seconds +# per hemisphere. See `cortex.polyutils.FlatSlab` for the parameters. +npz = cortex.db.get_surfinfo(subject, type="bumpy_flatmap") +offsets = np.vstack([npz["bump_left"], npz["bump_right"]]) +npz.close() + +# The two things it is being compared against, computed per hemisphere. +naive, thickness, onmap = [], [], [] +for hemi in ["lh", "rh"]: + wm, polys = cortex.db.get_surf(subject, "wm", hemi) + pia, _ = cortex.db.get_surf(subject, "pia", hemi) + flat, flatpolys = cortex.db.get_surf(subject, "flat", hemi) + + naive.append(naive_prism_height(flat, wm, pia, flatpolys)) + thickness.append(np.linalg.norm(pia - wm, axis=1)) + + # Vertices cut away from the flatmap -- the medial wall -- have no relief. + mask = np.zeros(len(wm), bool) + mask[flatpolys.ravel()] = True + onmap.append(mask) + +naive = np.concatenate(naive) +thickness = np.concatenate(thickness) +onmap = np.concatenate(onmap) +relief = offsets[:, 2] + +############################################################################### +# The relief itself, shown alongside the naive volume-preserving height on the +# same color scale, which is what makes the difference obvious: the naive map is +# mostly flat with a scatter of very bright spikes, because its range is set by +# a handful of crushed triangles. + +vmax = float(np.percentile(relief[onmap], 99)) +for name, field in [("folding height", relief), + ("naive V/A_flat", naive)]: + vertex = cortex.Vertex(field, subject, vmin=0, vmax=vmax, cmap="viridis") + cortex.quickshow(vertex, with_rois=False, with_labels=False, + with_curvature=False) + plt.title("bumpy flatmap height, %s (mm)" % name) + +############################################################################### +# Where the spikes are. Plotted as distributions, on a log axis because the +# naive height's problem is entirely in its tail. Cortex is 2-5 mm thick, so a +# sensible bumpy flatmap should sit in roughly that range; the naive height runs +# well past it. + +fig, ax = plt.subplots(figsize=(7, 4)) +bins = np.logspace(-1, 2, 120) +for name, field in [("cortical thickness", thickness), + ("naive V/A_flat", naive), ("folding height", relief)]: + ax.hist(np.clip(field[onmap], 1e-2, None), bins=bins, histtype="step", + label=name, linewidth=1.5) +ax.set_xscale("log") +ax.set_yscale("log") +ax.set_xlabel("height above the flatmap (mm)") +ax.set_ylabel("vertices") +ax.legend() +ax.set_title("the naive height's problem is its tail") + +############################################################################### +# Finally, the check that the relief means what it should: height against mean +# curvature. Gyral crowns, where the pia carries more area than the white matter +# beneath it, end up thicker; sulcal fundi end up thinner. + +curv = cortex.db.get_surfinfo(subject, type="curvature") +fig, ax = plt.subplots(figsize=(6, 4)) +ax.hexbin(curv.data[onmap], relief[onmap], gridsize=60, bins="log", + cmap="Blues") +ax.set_xlabel("mean curvature (sulci < 0 < gyri)") +ax.set_ylabel("relief height (mm)") +ax.set_title("relief follows the folding") + +plt.show() diff --git a/filestore/db/S1/surface-info/bumpy_flatmap.npz b/filestore/db/S1/surface-info/bumpy_flatmap.npz new file mode 100644 index 000000000..72fe0fb69 Binary files /dev/null and b/filestore/db/S1/surface-info/bumpy_flatmap.npz differ