Bumpy flatmap redesign

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)
 

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.

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]

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.

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,