BUG: fix gs_qc to work with N(streamlines)/mm^3

gyral_structure/qc/surface.py

@@ -24,13 +24,14 @@ def faces_to_vertices(arr, surface: CorticalMesh):
     return np.array(mapping.dot(arr)).flatten() / nneigh
 
 
-def plot_all(surface: CorticalMesh, mask, basis: BasisFunc,
+def plot_all(surface: CorticalMesh, mask, nstream, basis: BasisFunc,
              params=(), directory='.', skip_nsplit=False):
     """
     Plots all the QC plots on the surface
 
     :param surface: cortical mesh to be shown
     :param mask: which vertices to include
+    :param nstream: number of streamlines expected
     :param basis: basis function
     :param params: parameters
     :param directory: directory to store the plots
@@ -44,22 +45,28 @@ def plot_all(surface: CorticalMesh, mask, basis: BasisFunc,
     density = np.sqrt((field_faces ** 2).sum(-1))
     radiality = ncrossing / density
 
+    expected = nstream / masked_surface.size_faces()
+
     if not skip_nsplit:
         import matplotlib.pyplot as plt
         fig, axes = plt.subplots(1, 1)
         nsplit.plot_offset(request, basis, params, axes=axes)
         fig.savefig(op.join(directory, 'nsplit.png'))
 
-    scl = cifti2.ScalarAxis(['density (crossing)', 'density (total)', 'radiality'])
+    scl = cifti2.ScalarAxis([
+        'normed density (crossing)', 'normed density (total)',
+        'density (crossing)', 'density (total)', 'radiality'
+    ])
     bm = cifti2.BrainModelAxis.from_mask(mask, name=surface.anatomy.cifti)
     cifti2.Cifti2Image(
-        np.stack([faces_to_vertices(arr, masked_surface) for arr in [ncrossing, density, radiality]], axis=0),
+        np.stack([faces_to_vertices(arr, masked_surface) for arr in [ncrossing / expected, density / expected,
+                                                                     ncrossing, density, radiality]], axis=0),
         header=cifti2.Cifti2Header.from_axes((scl, bm))
     ).to_filename(op.join(directory, 'field.dscalar.nii'))
 
     for toplot, name, label in [
-        (ncrossing, 'density_crossing', r'$\rho$(streamlines crossing)'),
-        (density, 'density_total',  r'$\rho$(streamlines)'),
+        (ncrossing / expected, 'density_crossing', r'$\rho$(streamlines crossing)'),
+        (density / expected, 'density_total',  r'$\rho$(streamlines)'),
         (radiality, 'radiality',  'radial alignment'),
     ]:
         plot_array(surface, mask, toplot, op.join(directory, name), label)

gyral_structure/scripts/gs_qc

@@ -4,10 +4,14 @@ parser = argparse.ArgumentParser("""Produces QC plots for a fitted dataset""")
 
 parser.add_argument("fit", help='HDF5 file with the fit (from `gs_fit`)')
 parser.add_argument("out", help='output directory name')
-parser.add_argument("-s", '--surface', nargs=3, action='append', help=(
-    'expects three arguments: (<name of surface> <filename of surface> <filename of surface mask>). ' +
-    'Can be repeated'
-))
+parser.add_argument("-w", "--white", required=True,
+                    help="GIFTI file with WM/GM boundary")
+parser.add_argument("-m", "--midthickness",
+                    help="GIFTI file with mid-thickness surface")
+parser.add_argument("-p", "--pial", required=True,
+                    help="GIFTI file with pial surface")
+parser.add_argument("-s", "--surface-mask",
+                    help="GIFTI file with surface mask of which vertices to consider")
 parser.add_argument("-v", '--volume', nargs=2, action='append', help=(
         'expects two arguments: (<name of volume> <filename of volume>). ' +
         'Can be repeated'
@@ -18,24 +22,38 @@ parser.add_argument("-sn", "--skip_nsplit", action='store_true',
 args = parser.parse_args()
 
 from gyral_structure.basis import read
-from mcutils.surface import CorticalMesh
+from mcutils.surface import CorticalMesh, Cortex
 import os
 from gyral_structure import qc
 import nibabel as nib
+import numpy as np
 
 bf, params = read(args.fit)
 
-if args.surface is not None:
-    for surf_name, surf_fn, mask_surf_fn in args.surface:
-        out_directory = os.path.join(args.out, surf_name)
-        os.makedirs(out_directory, exist_ok=True)
-        qc.surface.plot_all(
-            CorticalMesh.read(surf_fn),
-            nib.load(mask_surf_fn).darrays[0].data > 0,
-            bf, params,
-            out_directory,
-            skip_nsplit=args.skip_nsplit,
-        )
+if args.surface_mask is None:
+    mask = np.ones(CorticalMesh.read(args.white).nvertices, dtype='bool')
+else:
+    mask = nib.load(args.surface_mask).darrays[0].data > 0
+
+cortex = Cortex([
+    CorticalMesh.read(args.white)[mask],
+    CorticalMesh.read(args.pial)[mask],
+])
+target_streamlines = cortex.wedge_volume()
+
+for surf_name in ('white', 'midthickness', 'pial'):
+    surf_fn = getattr(args, surf_name)
+    if surf_fn is None:
+        continue
+    out_directory = os.path.join(args.out, surf_name)
+    os.makedirs(out_directory, exist_ok=True)
+    qc.surface.plot_all(
+        CorticalMesh.read(surf_fn),
+        mask, target_streamlines,
+        bf, params,
+        out_directory,
+        skip_nsplit=args.skip_nsplit,
+    )
 
 if args.volume is not None:
     for vol_name, mask_vol_fn in args.volume: