[WIP] Generating of TVB compatible datasets

bin/aseg2srf

@@ -0,0 +1,119 @@
+#!/bin/bash
+
+# Print usage if no argument is given
+if [ -z "$1" ]; then
+cat <<EOU
+Generate surfaces for the subcortical structures
+segmented with FreeSurfer.
+
+Usage:
+sub2srf -s <list of subjects> [-l <list of labels>] [-d]
+
+Options:
+-s <list> : Specify a list of subjects between quotes,
+            e.g. -s "john bill mary mark" or a text file
+            containing one subject per line.
+-l <list> : Specify a list of labels between quotes,
+            e.g. "10 11 14 52 253", or a text file
+            containing one label per line, or ignore
+            this option to convert all labels.
+-d          Debug mode. Leave all temporary files.
+
+Requirements:
+FreeSurfer must have been configured and the variables
+FREESURFER_HOME and SUBJECTS_DIR must have been correctly set.
+
+_____________________________________
+Anderson M. Winkler
+Institute of Living / Yale University
+Jul/2009
+http://brainder.org
+EOU
+exit
+fi
+
+# List of labels to be converted if no list is specified
+LABLIST="4 5 7 8 10 11 12 13 14 15 16 17 18 26 28 43 44 46 47 49 50 51 52 53 54 58 60 251 252 253 254 255"
+
+# Check and accept arguments
+SBJLIST=""
+DEBUG=N
+while getopts 's:l:d' OPTION
+do
+  case ${OPTION} in
+    s) SBJLIST=$( [[ -f ${OPTARG} ]] && cat ${OPTARG} || echo "${OPTARG}" ) ;;
+    l) LABLIST=$( [[ -f ${OPTARG} ]] && cat ${OPTARG} || echo "${OPTARG}" ) ;;
+    d) DEBUG=Y ;;
+  esac
+done
+
+# Prepare a random string to save temporary files
+#if hash md5 2> /dev/null ; then
+#  RND0=$(head -n 1 /dev/random | md5)
+#elif hash md5sum 2> /dev/null ; then
+#  RND0=$(head -n 1 /dev/random | md5sum)
+#fi
+#RNDSTR=${RND0:0:12}
+RNDSTR=${RANDOM:0:12}
+
+# Define a function for Ctrl+C as soon as the RNDSTR is defined
+trap bashtrap INT
+bashtrap()
+{
+  break ; break
+  [[ "${s}" != "" ]] && rm -rf ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR} # try to delete temp files
+  exit 1
+}
+
+# For each subject
+for s in ${SBJLIST} ; do
+
+  # Create directories for temp files and results
+  mkdir -p ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}
+  mkdir -p ${SUBJECTS_DIR}/${s}/ascii
+
+  # For each label
+  for lab in ${LABLIST} ; do
+
+    # Label string
+    lab0=$(printf %03d ${lab})
+
+    # Pre-tessellate
+    echo "==> Pre-tessellating: ${s}, ${lab0}"
+    ${FREESURFER_HOME}/bin/mri_pretess \
+           ${SUBJECTS_DIR}/${s}/mri/aseg.mgz ${lab} \
+           ${SUBJECTS_DIR}/${s}/mri/norm.mgz \
+           ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}_filled.mgz
+
+    # Tessellate
+    echo "==> Tessellating: ${s}, ${lab0}"
+    ${FREESURFER_HOME}/bin/mri_tessellate \
+           ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}_filled.mgz \
+           ${lab} ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}_notsmooth
+
+    # Smooth
+    echo "==> Smoothing: ${s}, ${lab0}"
+    ${FREESURFER_HOME}/bin/mris_smooth -nw \
+           ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}_notsmooth \
+           ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}
+
+    # Convert to ASCII
+    echo "==> Converting to ASCII: ${s}, ${lab0}"
+    ${FREESURFER_HOME}/bin/mris_convert \
+           ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0} \
+           ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}.asc
+    mv ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}.asc \
+       ${SUBJECTS_DIR}/${s}/ascii/aseg_${lab0}.srf
+  done
+
+  # Get rid of temp files
+  if [ "${DEBUG}" == "Y" ] ; then
+    echo "==> Temporary files for ${s} saved at:"
+    echo "${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}"
+  else
+    echo "==> Removing temporary files for ${s}"
+    rm -rf ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}
+  fi
+  echo "==> Done: ${s}"
+done
+exit 0

bin/main.sh

@@ -13,6 +13,9 @@ sh $CODE/reconALL.sh
 source $CODE/postreconALL.sh
 #<-It returns all surfaces and volumes...
 
+# Generate surfaces for subcortical structures
+$CODE/aseg2srf -s ${SUBJECT}
+
 
 #Now, for Steps 2 we split to three parallel streams:
 

bin/run_surfaces_to_struct_dataset.py

@@ -0,0 +1,33 @@
+#!/usr/bin/env python3
+
+import logging
+import os
+import sys
+
+from tvb.recon.flow.surfaces_to_structural_datasets import SurfacesToStructuralDataset
+from tvb.recon.cli.runner import SimpleRunner
+
+def main():
+    subjects_dir, subjid, source_lut, target_lut, weights_file, tract_lengths_file, out_file, out_surf = sys.argv[1:]
+
+    logging.basicConfig(level=logging.INFO)
+    runner = SimpleRunner()
+
+    flow = SurfacesToStructuralDataset(
+        os.path.join(subjects_dir, subjid, "surf"),
+        os.path.join(subjects_dir, subjid, "label"),
+        os.path.join(subjects_dir, subjid, "ascii"),
+        source_lut,
+        target_lut,
+        weights_file,
+        tract_lengths_file,
+        out_file,
+        out_surf
+    )
+
+    flow.run(runner)
+
+
+
+if __name__ == "__main__":
+    main()

tvb/recon/algo/structural_dataset.py

@@ -0,0 +1,70 @@
+
+import numpy as np
+import os
+import tempfile
+from typing import List
+from zipfile import ZipFile
+
+
+class StructuralDataset:
+    def __init__(self,
+                 orientations: np.ndarray,
+                 areas: np.ndarray,
+                 centers: np.ndarray,
+                 cortical: np.ndarray,
+                 weights_ut: np.ndarray,
+                 tract_lengths_ut: np.ndarray,
+                 names: List[str]):
+
+        nregions = len(names)
+
+        assert orientations.shape == (nregions, 3)
+        assert areas.shape == (nregions,)
+        assert centers.shape == (nregions, 3)
+        assert cortical.shape == (nregions,)
+        assert weights_ut.shape == (nregions, nregions)
+        assert tract_lengths_ut.shape == (nregions, nregions)
+
+        # Upper triangular -> symmetric matrices
+        assert np.sum(np.tril(weights_ut, -1)) == 0
+        assert np.sum(np.tril(tract_lengths_ut, -1)) == 0
+        self.weights = weights_ut + weights_ut.transpose() - np.diag(np.diag(weights_ut))
+        self.tract_lengths = tract_lengths_ut + tract_lengths_ut.transpose() - np.diag(np.diag(tract_lengths_ut))
+
+        self.orientations = orientations
+        self.areas = areas
+        self.centers = centers
+        self.cortical = cortical
+        self.names = names
+
+    def save_to_txt_zip(self, filename: os.PathLike):
+
+        tmpdir = tempfile.TemporaryDirectory()
+
+        file_areas = os.path.join(tmpdir.name, 'areas.txt')
+        file_orientations = os.path.join(tmpdir.name, 'average_orientations.txt')
+        file_centres = os.path.join(tmpdir.name, 'centres.txt')
+        file_cortical = os.path.join(tmpdir.name, 'cortical.txt')
+        file_weights = os.path.join(tmpdir.name, 'weights.txt')
+        file_tract_lengths = os.path.join(tmpdir.name, 'tract_lengths.txt')
+
+        np.savetxt(file_areas, self.areas, fmt='%.2f')
+        np.savetxt(file_orientations, self.orientations, fmt='%.2f %.2f %.2f')
+        np.savetxt(file_cortical, self.cortical, fmt='%d')
+        np.savetxt(file_weights, self.weights, fmt='%d')
+        np.savetxt(file_tract_lengths, self.tract_lengths, fmt='%.3f')
+
+        with open(file_centres, 'w') as f:
+            for i, name in enumerate(self.names):
+                f.write('%s %.4f %.4f %.4f\n' % (name, self.centers[i, 0], self.centers[i, 1], self.centers[i, 2]))
+
+        with ZipFile(filename, 'w') as zip_file:
+            zip_file.write(file_areas, os.path.basename(file_areas))
+            zip_file.write(file_orientations, os.path.basename(file_orientations))
+            zip_file.write(file_centres, os.path.basename(file_centres))
+            zip_file.write(file_cortical, os.path.basename(file_cortical))
+            zip_file.write(file_weights, os.path.basename(file_weights))
+            zip_file.write(file_tract_lengths, os.path.basename(file_tract_lengths))
+
+    def save_to_h5(self, filename):
+        pass

tvb/recon/cli/core.py

@@ -31,8 +31,7 @@ def __new__(cls: type, name: str, bases: Tuple[type],
                 attrs: Dict[str, object]) -> type:
         for key, value in list(attrs.items()):
             if not key.startswith('_'):
-                value:
-                    str
+                value: str
                 attrs[key] = Flag(key, value)
         return super().__new__(cls, name, bases, attrs)
 
@@ -55,8 +54,7 @@ def __new__(cls: type, name: str, bases: Tuple[type],
                 attrs: Dict[str, object]) -> type:
         for key, value in list(attrs.items()):
             if not key.startswith('_'):
-                value:
-                    str
+                value: str
                 attrs[key] = EnvVar(key, value)
         return super().__new__(cls, name, bases, attrs)
 

tvb/recon/cli/fs.py

@@ -204,6 +204,14 @@ class mris_decimate(BaseFsCLI):
     exe = 'mris_decimate'
 
 
+class mris_info(BaseFsCLI):
+    """
+    The mris_info command from the fs package
+
+    """
+    exe = 'mris_info'
+
+
 class mris_smooth(BaseFsCLI):
     """
     The mris_smooth command from the fs package.
@@ -276,17 +284,17 @@ def convert(in_, out,
 
 # TODO in/out typed Union[Volume, VolFile] 
 def binarize(in_file: str, out_file: str,
-             min_val: Union[int, float, NoneType]=None,
+             min_val: Union[int, float, None]=None,
              erode: int=0,
              dilate: int=0,
-             mask_file: Union[str, NoneType]=None,
+             mask_file: Union[str, None]=None,
             ):
     "Binarize image, possibly with dilation or erosion."
     cmd = mri_binarize
     args = [
         cmd.exe,
         cmd.Flags.in_file, in_file,
-        cmd.Flags.out_file, out_file,
+        cmd.Flags.out_file, out_file]
 
     if min_val is not None:
         args += [cmd.Flags.min_val, str(min_val)]
@@ -298,3 +306,10 @@ def binarize(in_file: str, out_file: str,
         args += [cmd.Flags.mask_file, mask_file]
     return args
 
+
+def mris_convert_run(in_: os.PathLike, out: os.PathLike):
+    return [mris_convert.exe, in_, out]
+
+
+def mris_info_run(in_: os.PathLike):
+    return [mris_info.exe, in_]

tvb/recon/cli/runner.py

@@ -45,7 +45,7 @@ def tmp_fname(self, fname: str): pass
     def fname(self, fname: str): pass
 
     @abc.abstractmethod
-    def run(self, args): pass
+    def run(self, args, capture_output: bool): pass
 
 
 class SimpleRunner(Runner):
@@ -110,10 +110,18 @@ def stringify_args(self, args) -> Sequence[str]:
             str_args.append(val)
         return str_args
 
-    def run(self, args, **kwargs):
+    def run(self, args, capture_output=False, **kwargs):
         str_args = self.stringify_args(args)
         for i, arg in enumerate(str_args):
             self.logger.debug('args[%d]: %r', i, arg)
         tic = time.time()
-        subprocess.check_call(str_args, **kwargs)
+
+        if not capture_output:
+            subprocess.check_call(str_args, **kwargs)
+            output = None
+        else:
+            output = subprocess.check_output(str_args, stderr=subprocess.STDOUT, universal_newlines=True, **kwargs)
+
         self.logger.debug('elapsed %.2fs' % (time.time() - tic, ))
+        return output
+