Merge pull request #283 from nipreps/enh/mcribs-surface-processing

ENH: Add M-CRIB-S surface processing

Author: mgxd

nibabies/cli/parser.py

@@ -664,6 +664,12 @@ def _slice_time_ref(value, parser):
         action="store_true",
         help="Force traditional FreeSurfer surface reconstruction.",
     )
+    g_baby.add_argument(
+        "--surface-recon-method",
+        choices=("infantfs", "freesurfer", "mcribs"),
+        default="infantfs",
+        help="Method to use for surface reconstruction",
+    )
     return parser
 
 
@@ -674,6 +680,14 @@ def parse_args(args=None, namespace=None):
     parser = _build_parser()
     opts = parser.parse_args(args, namespace)
 
+    # Deprecations
+    if opts.force_reconall:
+        config.loggers.cli.warning(
+            "--force-reconall is deprecated and will be removed in a future release."
+            "To run traditional `recon-all`, use `--surface-recon-method freesurfer` instead."
+        )
+        opts.surface_recon_method = "freesurfer"
+
     if opts.config_file:
         skip = {} if opts.reports_only else {"execution": ("run_uuid",)}
         config.load(opts.config_file, skip=skip)
@@ -745,14 +759,21 @@ def parse_args(args=None, namespace=None):
 
     if config.execution.fs_subjects_dir is None:
         if output_layout == "bids":
-            config.execution.fs_subjects_dir = output_dir / "sourcedata" / "infant-freesurfer"
+            config.execution.fs_subjects_dir = output_dir / "sourcedata" / "freesurfer"
         elif output_layout == "legacy":
-            config.execution.fs_subjects_dir = output_dir / "infant-freesurfer"
+            config.execution.fs_subjects_dir = output_dir / "freesurfer"
     if config.execution.nibabies_dir is None:
         if output_layout == "bids":
             config.execution.nibabies_dir = output_dir
         elif output_layout == "legacy":
             config.execution.nibabies_dir = output_dir / "nibabies"
+    if config.workflow.surface_recon_method == "mcribs":
+        if output_layout == "bids":
+            config.execution.mcribs_dir = output_dir / "sourcedata" / "mcribs"
+        elif output_layout == "legacy":
+            config.execution.mcribs_dir = output_dir / "mcribs"
+        # Ensure the directory is created
+        config.execution.mcribs_dir.mkdir(exist_ok=True, parents=True)
 
     # Wipe out existing work_dir
     if opts.clean_workdir and work_dir.exists():

nibabies/config.py

@@ -394,6 +394,8 @@ class execution(_Config):
     """Output verbosity."""
     low_mem = None
     """Utilize uncompressed NIfTIs and other tricks to minimize memory allocation."""
+    mcribs_dir = None
+    """M-CRIB-S processing and output directory."""
     md_only_boilerplate = False
     """Do not convert boilerplate from MarkDown to LaTex and HTML."""
     nibabies_dir = None
@@ -439,6 +441,7 @@ class execution(_Config):
         "fs_subjects_dir",
         "layout",
         "log_dir",
+        "mcribs_dir",
         "nibabies_dir",
         "output_dir",
         "segmentation_atlases_dir",
@@ -543,8 +546,6 @@ class workflow(_Config):
     """Remove the mean from fieldmaps."""
     force_syn = None
     """Run *fieldmap-less* susceptibility-derived distortions estimation."""
-    force_reconall = False
-    """Force traditional FreeSurfer surface reconstruction instead of infant version."""
     hires = None
     """Run FreeSurfer ``recon-all`` with the ``-hires`` flag."""
     ignore = None
@@ -578,6 +579,8 @@ class workflow(_Config):
     spaces = None
     """Keeps the :py:class:`~niworkflows.utils.spaces.SpatialReferences`
     instance keeping standard and nonstandard spaces."""
+    surface_recon_method = "infantfs"
+    """Method to use for surface reconstruction."""
     topup_max_vols = 5
     """Maximum number of volumes to use with TOPUP, per-series (EPI or BOLD)."""
     use_aroma = None

nibabies/data/boilerplate.bib

@@ -392,3 +392,16 @@ @article{infantfs
 	title = {Infant {FreeSurfer}: An automated segmentation and surface extraction pipeline for T1-weighted neuroimaging data of infants 0{\textendash}2 years},
 	journal = {{NeuroImage}}
 }
+
+@article{mcribs,
+	doi = {10.1038/s41598-020-61326-2},
+	url = {https://doi.org/10.1038%2Fs41598-020-61326-2},
+	year = 2020,
+	month = {mar},
+	publisher = {Springer Science and Business Media {LLC}},
+	volume = {10},
+	number = {1},
+	author = {Chris L. Adamson and Bonnie Alexander and Gareth Ball and Richard Beare and Jeanie L. Y. Cheong and Alicia J. Spittle and Lex W. Doyle and Peter J. Anderson and Marc L. Seal and Deanne K. Thompson},
+	title = {Parcellation of the neonatal cortex using Surface-based Melbourne Children's Regional Infant Brain atlases (M-{CRIB}-S)},
+	journal = {Scientific Reports}
+}

nibabies/interfaces/mcribs.py

@@ -36,13 +36,14 @@ class MCRIBReconAllInputSpec(CommandLineInputSpec):
     )
     t2w_file = File(
         exists=True,
-        required=True,
         copyfile=True,
         desc='T2w (Isotropic + N4 corrected)',
     )
     segmentation_file = File(
+        exists=True,
         desc='Segmentation file (skips tissue segmentation)',
     )
+    mask_file = File(exists=True, desc='T2w mask')
 
     # MCRIBS options
     conform = traits.Bool(
@@ -54,39 +55,31 @@ class MCRIBReconAllInputSpec(CommandLineInputSpec):
         desc='Perform tissue type segmentation',
     )
     surfrecon = traits.Bool(
-        True,
-        usedefault=True,
         argstr='--surfrecon',
         desc='Reconstruct surfaces',
     )
     surfrecon_method = traits.Enum(
         'Deformable',
         argstr='--surfreconmethod %s',
-        usedefault=True,
+        requires=['surfrecon'],
         desc='Surface reconstruction method',
     )
     join_thresh = traits.Float(
-        1.0,
         argstr='--deformablejointhresh %f',
-        usedefault=True,
+        requires=['surfrecon'],
         desc='Join threshold parameter for Deformable',
     )
     fast_collision = traits.Bool(
-        True,
         argstr='--deformablefastcollision',
-        usedefault=True,
+        requires=['surfrecon'],
         desc='Use Deformable fast collision test',
     )
     autorecon_after_surf = traits.Bool(
-        True,
         argstr='--autoreconaftersurf',
-        usedefault=True,
         desc='Do all steps after surface reconstruction',
     )
     segstats = traits.Bool(
-        True,
         argstr='--segstats',
-        usedefault=True,
         desc='Compute statistics on segmented volumes',
     )
     nthreads = traits.Int(
@@ -97,6 +90,7 @@ class MCRIBReconAllInputSpec(CommandLineInputSpec):
 
 class MCRIBReconAllOutputSpec(TraitedSpec):
     mcribs_dir = Directory(desc='MCRIBS output directory')
+    subjects_dir = Directory(desc='FreeSurfer output directory')
 
 
 class MCRIBReconAll(CommandLine):
@@ -111,10 +105,17 @@ def cmdline(self):
         # Avoid processing if valid
         if self.inputs.outdir:
             sid = self.inputs.subject_id
-            logf = Path(self.inputs.outdir) / sid / 'logs' / f'{sid}.log'
-            if logf.exists():
-                logtxt = logf.read_text().splitlines()[-3:]
-                self._no_run = 'Finished without error' in logtxt
+            # Check MIRTK surface recon deformable
+            if self.inputs.surfrecon:
+                surfrecon_dir = Path(self.inputs.outdir) / sid / 'SurfReconDeformable' / sid
+                if self._verify_surfrecon_outputs(surfrecon_dir, error=False):
+                    self._no_run = True
+            # Check FS directory population
+            elif self.inputs.autorecon_after_surf:
+                fs_dir = Path(self.inputs.outdir) / sid / 'freesurfer' / sid
+                if self._verify_autorecon_outputs(fs_dir, error=False):
+                    self._no_run = True
+
             if self._no_run:
                 return "echo MCRIBSReconAll: nothing to do"
         return cmd
@@ -147,21 +148,22 @@ def _setup_directory_structure(self, mcribs_dir: Path) -> None:
         root.mkdir(**mkdir_kw)
 
         # T2w operations
-        t2w = root / 'RawT2' / f'{sid}.nii.gz'
-        t2w.parent.mkdir(**mkdir_kw)
-        if not t2w.exists():
-            shutil.copy(self.inputs.t2w_file, str(t2w))
-
-        if not self.inputs.conform:
-            t2wiso = root / 'RawT2RadiologicalIsotropic' / f'{sid}.nii.gz'
-            t2wiso.parent.mkdir(**mkdir_kw)
-            if not t2wiso.exists():
-                t2wiso.symlink_to(f'../RawT2/{sid}.nii.gz')
-
-            n4 = root / 'TissueSegDrawEM' / sid / 'N4' / f'{sid}.nii.gz'
-            n4.parent.mkdir(**mkdir_kw)
-            if not n4.exists():
-                n4.symlink_to(f'../../../RawT2/{sid}.nii.gz')
+        if self.inputs.t2w_file:
+            t2w = root / 'RawT2' / f'{sid}.nii.gz'
+            t2w.parent.mkdir(**mkdir_kw)
+            if not t2w.exists():
+                shutil.copy(self.inputs.t2w_file, str(t2w))
+
+            if not self.inputs.conform:
+                t2wiso = root / 'RawT2RadiologicalIsotropic' / f'{sid}.nii.gz'
+                t2wiso.parent.mkdir(**mkdir_kw)
+                if not t2wiso.exists():
+                    t2wiso.symlink_to(f'../RawT2/{sid}.nii.gz')
+
+                n4 = root / 'TissueSegDrawEM' / sid / 'N4' / f'{sid}.nii.gz'
+                n4.parent.mkdir(**mkdir_kw)
+                if not n4.exists():
+                    n4.symlink_to(f'../../../RawT2/{sid}.nii.gz')
 
         # Segmentation
         if self.inputs.segmentation_file:
@@ -184,6 +186,11 @@ def _setup_directory_structure(self, mcribs_dir: Path) -> None:
                 if not surfrec.exists():
                     surfrec.symlink_to(f'../../../RawT2/{sid}.nii.gz')
 
+                if self.inputs.mask_file:
+                    surfrec_mask = surfrec.parent / 'brain-mask.nii.gz'
+                    if not surfrec_mask.exists():
+                        shutil.copy(self.inputs.mask_file, str(surfrec_mask))
+
         if self.inputs.surfrecon:
             # Create FreeSurfer layout to safeguard against cd-ing into missing directories
             for d in ('surf', 'mri', 'label', 'scripts', 'stats'):
@@ -196,21 +203,82 @@ def _run_interface(self, runtime):
         # if users wish to preserve their runs
         mcribs_dir = self.inputs.outdir or Path(runtime.cwd) / 'mcribs'
         self._mcribs_dir = Path(mcribs_dir)
-        self._setup_directory_structure(self._mcribs_dir)
+        if self.inputs.surfrecon:
+            assert self.inputs.t2w_file, "Missing T2w input"
+            self._setup_directory_structure(self._mcribs_dir)
         # overwrite CWD to be in MCRIB subject's directory
         runtime.cwd = str(self._mcribs_dir / self.inputs.subject_id)
         return super()._run_interface(runtime)
 
     def _list_outputs(self):
         outputs = self._outputs().get()
-        outputs['mcribs_dir'] = str(self._mcribs_dir)
-
-        # Copy freesurfer directory into FS subjects dir
         sid = self.inputs.subject_id
-        mcribs_fs = self._mcribs_dir / sid / 'freesurfer' / sid
-        if mcribs_fs.exists() and self.inputs.subjects_dir:
+        if self.inputs.surfrecon:
+            # verify surface reconstruction was successful
+            surfrecon_dir = self._mcribs_dir / sid / 'SurfReconDeformable' / sid
+            self._verify_surfrecon_outputs(surfrecon_dir, error=True)
+
+        outputs['mcribs_dir'] = str(self._mcribs_dir)
+        if self.inputs.autorecon_after_surf and self.inputs.subjects_dir:
+            mcribs_fs = self._mcribs_dir / sid / 'freesurfer' / sid
+            self._verify_autorecon_outputs(mcribs_fs, error=True)
             dst = Path(self.inputs.subjects_dir) / self.inputs.subject_id
             if not dst.exists():
                 shutil.copytree(mcribs_fs, dst)
+            outputs['subjects_dir'] = self.inputs.subjects_dir
 
         return outputs
+
+    @staticmethod
+    def _verify_surfrecon_outputs(surfrecon_dir: Path, error: bool) -> bool:
+        """
+        Sanity check to ensure the surface reconstruction was successful.
+
+        MCRIBReconAll does not return a failing exit code if a step failed, which leads
+        this interface to be marked as completed without error in such cases.
+        """
+        # fmt:off
+        surfrecon_files = {
+            'meshes': (
+                'pial-lh-reordered.vtp',
+                'pial-rh-reordered.vtp',
+                'white-rh.vtp',
+                'white-lh.vtp',
+            )
+        }
+        # fmt:on
+        for d, fls in surfrecon_files.items():
+            for fl in fls:
+                if not (surfrecon_dir / d / fl).exists():
+                    if error:
+                        raise FileNotFoundError(f"SurfReconDeformable missing: {fl}")
+                    return False
+        return True
+
+    @staticmethod
+    def _verify_autorecon_outputs(fs_dir: Path, error: bool) -> bool:
+        """
+        Sanity check to ensure the necessary FreeSurfer files have been created.
+
+        MCRIBReconAll does not return a failing exit code if a step failed, which leads
+        this interface to be marked as completed without error in such cases.
+        """
+        # fmt:off
+        fs_files = {
+            'mri': ('T2.mgz', 'aseg.presurf.mgz', 'ribbon.mgz', 'brain.mgz'),
+            'label': ('lh.cortex.label', 'rh.cortex.label'),
+            'stats': ('aseg.stats', 'brainvol.stats', 'lh.aparc.stats', 'rh.curv.stats'),
+            'surf': (
+                'lh.pial', 'rh.pial',
+                'lh.white', 'rh.white',
+                'lh.curv', 'rh.curv',
+                'lh.thickness', 'rh.thickness'),
+        }
+        # fmt:on
+        for d, fls in fs_files.items():
+            for fl in fls:
+                if not (fs_dir / d / fl).exists():
+                    if error:
+                        raise FileNotFoundError(f"FreeSurfer directory missing: {fl}")
+                    return False
+        return True