Merge branch 'dev' of github.com:freesurfer/freesurfer into dev

CMakeLists.txt

@@ -659,7 +659,7 @@ if(NOT MINIMAL)
     qdec_glmfit
     resurf
     spline3
-    subfields
+    subregions
     stat_normalize
     stim_polar
     mri_synthseg

include/GradUnwarp.h

@@ -3,7 +3,7 @@
 
 #include "gcamorph.h"
 #include "mri.h"
-#include "vol_geom.h"
+//#include "vol_geom.h"
 
 typedef struct
 {

include/transform.h

@@ -24,11 +24,8 @@
 #include "const.h"
 #include "float.h"
 
-#include "vol_geom.h"
-
 typedef enum { MINC, TKREG, GENERIC, UNKNOWN=-1 } TransformType;
 
-#if 0
 typedef struct
 {
   int           valid;   /* whether this is a
@@ -46,7 +43,6 @@ typedef struct
   char          fname[STRLEN];  // volume filename
 }
 VOL_GEOM, VG;
-#endif
 
 typedef struct
 {

include/vol_geom.h

@@ -1,7 +1,6 @@
 #ifndef VOL_GEOM_H
 #define VOL_GEOM_H
 
-//#include "gca.h"
 #include "mri.h"
 
 struct VOL_GEOM;
@@ -68,4 +67,5 @@ struct VOL_GEOM
 };
 
 typedef VOL_GEOM VG;
+
 #endif

mris_gradient_unwarp/mris_gradient_unwarp.cpp

@@ -159,18 +159,22 @@ int main(int argc, char *argv[])
     vox2ras_orig = extract_i_to_r(origvol);  //MRIxfmCRS2XYZ(origvol, 0);
     inv_vox2ras_orig = MatrixInverse(vox2ras_orig, NULL);  //extract_r_to_i(origvol);
 
-    vg.copyFromMRI(origvol);
+    //vg.copyFromMRI(origvol);  // vol_geom.cpp
+    getVolGeom(origvol, &vg);
 
     printVolInfo(origvol, vox2ras_orig, inv_vox2ras_orig);
   }
   else if (insurf != NULL)
   {
     origsurf = MRISread(insurf);
 
-    vox2ras_orig     = origsurf->vg.getVox2RAS();
-    inv_vox2ras_orig = origsurf->vg.getRAS2Vox();
+    //vox2ras_orig     = origsurf->vg.getVox2RAS();
+    //inv_vox2ras_orig = origsurf->vg.getRAS2Vox();
+    //vg.copy(&origsurf->vg);  // vol_geom.cpp
 
-    vg.copy(&origsurf->vg);
+    vox2ras_orig     = vg_i_to_r(&origsurf->vg);
+    inv_vox2ras_orig = vg_r_to_i(&origsurf->vg);
+    copyVolGeom(&origsurf->vg, &vg);
   }
 
   GradUnwarp *gradUnwarp = new GradUnwarp(nthreads);

python/freesurfer/subregions/__init__.py

@@ -0,0 +1,4 @@
+from freesurfer.subregions.process import structure_names
+from freesurfer.subregions.process import run_cross_sectional
+from freesurfer.subregions.process import run_longitudinal
+from freesurfer.samseg.gemsbindings import setGlobalDefaultNumberOfThreads as set_thread_count

python/freesurfer/subfields/brainstem.py → python/freesurfer/subregions/brainstem.py

@@ -5,8 +5,8 @@
 import freesurfer as fs
 
 from freesurfer import samseg
-from freesurfer.subfields import utils
-from freesurfer.subfields.core import MeshModel
+from freesurfer.subregions import utils
+from freesurfer.subregions.core import MeshModel
 
 
 class BrainstemSubstructures(MeshModel):

python/freesurfer/subfields/core.py → python/freesurfer/subregions/core.py

@@ -6,7 +6,7 @@
 import freesurfer as fs
 
 from freesurfer import samseg
-from freesurfer.subfields import utils
+from freesurfer.subregions import utils
 
 
 class MeshModel:

python/freesurfer/subfields/hippocampus.py → python/freesurfer/subregions/hippocampus.py

@@ -6,8 +6,8 @@
 import freesurfer as fs
 
 from freesurfer import samseg
-from freesurfer.subfields import utils
-from freesurfer.subfields.core import MeshModel
+from freesurfer.subregions import utils
+from freesurfer.subregions.core import MeshModel
 
 
 class HippoAmygdalaSubfields(MeshModel):

python/freesurfer/subfields/process.py → python/freesurfer/subregions/process.py

@@ -3,10 +3,10 @@
 import freesurfer as fs
 
 from freesurfer import samseg
-from freesurfer.subfields import utils
-from freesurfer.subfields.thalamus import ThalamicNuclei
-from freesurfer.subfields.brainstem import BrainstemSubstructures
-from freesurfer.subfields.hippocampus import HippoAmygdalaSubfields
+from freesurfer.subregions import utils
+from freesurfer.subregions.thalamus import ThalamicNuclei
+from freesurfer.subregions.brainstem import BrainstemSubstructures
+from freesurfer.subregions.hippocampus import HippoAmygdalaSubfields
 
 
 model_lookup = {

python/freesurfer/subfields/thalamus.py → python/freesurfer/subregions/thalamus.py

@@ -5,8 +5,8 @@
 import freesurfer as fs
 
 from freesurfer import samseg
-from freesurfer.subfields import utils
-from freesurfer.subfields.core import MeshModel
+from freesurfer.subregions import utils
+from freesurfer.subregions.core import MeshModel
 
 
 class ThalamicNuclei(MeshModel):

samseg/gems_train_mesh

@@ -0,0 +1,165 @@
+#!/usr/bin/env python3
+
+###########################
+#
+# Takes inputs required for `kvlBuildAtlasMesh` and a training schedule
+# and makes repeated calls to `kvlBuildAtlasMesh` according to the training
+# schedule, using the output of one call to `kvlBuildAtlasMesh` as input to the
+# next (via `explicitStartCollection`)
+#
+# The file `training-schedule-example.txt` contains an example training schedule
+# every line should contain 3 numbers
+#   - The number of interations to run for this epoch (unsigned int)
+#   - The 'stiffness' factor to pass to kvlBuildAtlasMesh (float)
+#   - The 'edgeCollapseEncouragementFactor' to pass to kvlBuildAtlasMesh (float)
+#
+# Example:
+#  /autofs/cluster/gerenuk/pwighton/fs/freesurfer/samseg/gems_train_mesh \
+#    --num-upsamples 1 \
+#    --mesh-size 3 3 3 \
+#    --schedule-file /autofs/cluster/gerenuk/pwighton/fs/freesurfer/samseg/training-schedule-example.txt \
+#    --work-dir /autofs/cluster/gerenuk/pwighton/samseg/test-gems-train-mesh-02 \
+#    --binary /autofs/cluster/gerenuk/pwighton/samseg/install/gems/bin/kvlBuildAtlasMesh \
+#    --label-files \
+#        seg_1.mgz \
+#        seg_2.mgz \
+#        seg_3.mgz
+###########################
+
+import os
+import sys
+import argparse
+import tempfile
+import errno
+import shutil
+
+def parse_args(args):
+  parser = argparse.ArgumentParser()
+  parser.add_argument('-n','--num-upsamples', required=True, \
+                      help='The number of upsapmling steps for `kvlBuildAtlasMesh` to perform.')
+  parser.add_argument('-m','--mesh-size', nargs=3, required=True, \
+                      help='The mesh size (x, y, z) to pass to `kvlBuildAtlasMesh`.')
+  parser.add_argument('-s','--schedule-file', required=False, default=None, \
+                      help='Filename containing the training schedule for sucessive calls to `kvlBuildAtlasMesh`')
+  parser.add_argument('-w', '--work-dir', required=False, default=None, \
+                      help='Directory under which to keep output and intermediary files (will be created)')
+  parser.add_argument('-b', '--binary', required=False, default=None, \
+                      help='Location of `kvlBuildAtlasMesh` binary')
+  parser.add_argument('-l', '--label-files', nargs='+', required=True, \
+                      help='The ground truth lables to train against')
+  #parser.add_argument('-t', '--lookup-table', required=False, default=None, \
+  #                    help='The FreeSurfer-like lookup table')
+  args = parser.parse_args()
+
+  if args.schedule_file is None:
+    args.schedule = default_training_schedule()
+  else:
+    args.schedule = read_schedule_file(args.schedule_file)
+
+  # Make sure we can find the kvlBuildAtlasMesh binary    
+  if args.binary is None:
+    if os.environ.get('FREESURFER_HOME'):
+      args.binary = os.path.join(os.environ.get('FREESURFER_HOME'),'gems/bin/kvlBuildAtlasMesh')
+  if args.binary is None or not os.path.exists(args.binary):
+    print("gems_train_mesh:  ERROR: Can't find kvlBuildAtlasMesh, either set the FREESURFER_HOME env var or use -b")
+    raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), args.binary)
+
+  # Make sure we can find the FreeSurfer lookup table
+  #if args.lookup_table is None:
+  #  if os.environ.get('FREESURFER_HOME'):
+  #    args.lookup_table = os.path.join(os.environ.get('FREESURFER_HOME'),'FreeSurferColorLUT.txt')
+  #if args.lookup_table is None or not os.path.exists(args.lookup_table)
+  #  print("gems_train_mesh:  ERROR: Can't find FreeSurfer-like lookup table, either set the FREESURFER_HOME env var or use -t")
+  #  raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), args.binary)
+
+  # Make sure all input label files exist
+  for file in args.label_files:
+    if not os.path.exists(file):
+      print("gems_train_mesh:  ERROR: Can't find label file "+file)
+      raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), file)
+
+  # Make the working dir if it doesn't exist
+  if args.work_dir is None:
+    args.work_dir = tempfile.mkdtemp()
+    print("gems_train_mesh:  No workdir specified so using "+args.work_dir)
+  else:
+    os.makedirs(args.work_dir, exist_ok=True)
+
+  return args
+
+# The training shedule format is a list of lists, eg:
+# [ [num_itr, stiffness, edgeCollaseFactor],
+#   [num_itr, stiffness, edgeCollaseFactor] ]
+def default_training_schedule():
+  return [ [1, 1.0, 1.0],
+           [10, 0.1, 1.25] ]
+
+def read_schedule_file(filename):
+  with open(filename) as f:
+    schedule = [[float(x) for x in line.split()] for line in f]
+  # Every line of the training schedule should have exactly three values:
+  #   - num itr (uint)
+  #   - stiffness (float)
+  #   - edge collapse factor (float)
+  for epoch in schedule:
+    assert(len(epoch) == 3)  
+  return schedule
+
+# Runs kvl_BuildAtlasMesh and returns a filename string with the last meshCollection created   
+def run_kvlBuildAtlasMesh(binary, work_dir, epoch_num, num_upsamples, mesh_size, label_files, num_itr=5, stiffness=1.0, edgeCollapse=1.0, startCollection=None):
+  launch_dir = os.path.join(work_dir, 'epoch_'+'{:04d}'.format(epoch_num)+'_launch_dir')
+  out_dir = os.path.join(work_dir, 'epoch_'+'{:04d}'.format(epoch_num)+'_out_dir')
+
+  os.makedirs(launch_dir, exist_ok=True)
+  os.makedirs(out_dir, exist_ok=True)
+
+  print(f"gems_train_mesh: epoch: {epoch_num}\n  launch dir: {launch_dir}\n  out dir: {out_dir}")
+
+  # Copy explicitStartCollection.gz to launch dir if specified
+  if startCollection is not None:
+    startCollection_dest = os.path.join(launch_dir, 'explicitStartCollection.gz')
+    shutil.copy(startCollection, startCollection_dest)
+
+  # launch  
+  orig_dir = os.getcwd()
+  os.chdir(launch_dir)
+  launch_cmd = f"{binary} {num_upsamples} {mesh_size[0]} {mesh_size[1]} {mesh_size[2]} {stiffness} {num_itr} {edgeCollapse} {out_dir} "+' '.join(label_files)
+  print("  launch command: "+launch_cmd)
+  os.system(launch_cmd)
+  os.chdir(orig_dir)
+
+  # kvlBuildAtlasMesh will create meshCollections in out_dir up to and including num_itr.
+  # If we can't find the last meshCollection, then something went wrong with kvlBuildAtlasMesh 
+  lastMeshCollection_file = os.path.join(out_dir, f"CurrentMeshCollection{num_itr}.gz")
+  if not os.path.exists(lastMeshCollection_file):
+    print(f"gems_train_mesh:  ERROR: Can't find the file {lastMeshCollection_file} that we were expecting kvlBuildAtlasMesh to create")
+    raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), lastMeshCollection_file)
+  return lastMeshCollection_file
+
+def main(argv):
+  args = parse_args(argv)
+  orig_cwd = os.getcwd()
+  startCollection = None
+  try:
+    for epoch_num, epoch in enumerate(args.schedule):
+      kvl_num_itr = int(epoch[0])
+      kvl_stiffness = float(epoch[1])
+      kvl_edge_collapse = float(epoch[2])
+      startCollection = run_kvlBuildAtlasMesh(\
+                              binary=args.binary, \
+                              work_dir=args.work_dir, \
+                              epoch_num=epoch_num, \
+                              num_upsamples=args.num_upsamples, \
+                              mesh_size=args.mesh_size, \
+                              label_files=args.label_files, \
+                              num_itr=kvl_num_itr, \
+                              stiffness=kvl_stiffness, \
+                              edgeCollapse=kvl_edge_collapse, \
+                              startCollection=startCollection)
+  except Exception as e:
+    os.chdir(orig_cwd)
+    raise e
+
+if __name__ == "__main__":
+  sys.exit(main(sys.argv))
+

samseg/training-schedule-example.txt

@@ -0,0 +1,4 @@
+1  1.0   1.0
+2  0.5   1.05
+4  0.25  1.10
+10 0.1   1.25

subregions/CMakeLists.txt

@@ -0,0 +1,3 @@
+project(subregions)
+
+install_pyscript(segment_subregions)

subfields/segment_subfields → subregions/segment_subregions

@@ -2,28 +2,28 @@
 
 import os
 import freesurfer as fs
-from freesurfer import subfields
+from freesurfer import subregions
 
 
 description = f'''
 Cross-sectional and longitudinal segmentation for the following
-structures: {", ".join(subfields.structure_names)}. To segment
+structures: {", ".join(subregions.structure_names)}. To segment
 the thalamic nuclei, for example, in a cross-sectional analysis:
 
-    segment_subfields thalamus --cross subj
+    segment_subregions thalamus --cross subj
 
 Similarly, for a longitudinal analysis, run:
 
-    segment_subfields thalamus --long-base base \\ 
-                      --long-tps m00.long.base m24.long.base
+    segment_subregions thalamus --long-base base \\ 
+                        --long-tps m00.long.base m24.long.base
 
 where each argument to --long-tps represents a subject timepoint.
 Output segmentations and computed structure volumes will be
 saved to the subject's `mri` subdirectory.
 '''
 
 parser = fs.utils.ArgumentParser(description=description)
-parser.add_argument('structure', help=f'Structure to segment. Options are: {", ".join(subfields.structure_names)}.')
+parser.add_argument('structure', help=f'Structure to segment. Options are: {", ".join(subregions.structure_names)}.')
 parser.add_argument('--cross', help='Subject to segment in cross-sectional analysis.')
 parser.add_argument('--long-tps', nargs='+', help='Subject timepoints to segment in longitudinal analysis.')
 parser.add_argument('--long-base', help='Base subject for longitudinal analysis.')
@@ -40,7 +40,7 @@ if not fs.fshome():
     fs.fatal('FREESURFER_HOME must be set!')
 
 # Specify the maximum number of threads the GEMS code will use
-subfields.set_thread_count(args.threads)
+subregions.set_thread_count(args.threads)
 
 # Sanity check on process type
 if args.cross and (args.long_base or args.long_tps):
@@ -92,16 +92,16 @@ for side in sides:
             'tempDir': args.temp_dir,
         })
         if args.structure == 'hippo-amygdala':
-            # Provide some extra data for hippocampal subfields
+            # Provide some extra data for hippocampal subregions
             subjParameters['wmParcFileName'] = os.path.join(subjdir, 'mri', 'wmparc.mgz')
         return subjParameters
 
     if args.cross:
         # Cross-sectional analysis
         parameters = config_subj_params(args.cross, out_dir=args.out_dir)
-        subfields.run_cross_sectional(args.structure, parameters)
+        subregions.run_cross_sectional(args.structure, parameters)
     else:
         # Longitudinal analysis
         baseParameters = config_subj_params(args.long_base, temp_subdir=f'base')
         tpParameters = [config_subj_params(subj, temp_subdir=f'tp{t:02d}') for t, subj in enumerate(args.long_tps)]
-        subfields.run_longitudinal(args.structure, baseParameters, tpParameters)
+        subregions.run_longitudinal(args.structure, baseParameters, tpParameters)

utils/CMakeLists.txt

@@ -145,7 +145,7 @@ add_library(utils STATIC
   MRISrigidBodyAlignGlobal.cpp
   mris_sphshapepvf.cpp
   GradUnwarp.cpp
-  vol_geom.cpp
+  #vol_geom.cpp
   mrisurf.cpp
   mrisurf_base.cpp
   mrisurf_compute_dxyz.cpp