Merge pull request #252 from mgxd/enh/t2star-output

ENH: Output T2star maps for multiecho data

Author: mgxd

.circleci/bcp_full_outputs.txt

@@ -54,6 +54,7 @@ sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-confounds_times
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-confounds_timeseries.tsv
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-preproc_bold.json
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-preproc_bold.nii.gz
+sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-scanner_to-boldref_mode-image_xfm.txt
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-scanner_to-T1w_mode-image_xfm.txt
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-T1w_to-scanner_mode-image_xfm.txt
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_space-MNIInfant_cohort-1_boldref.nii.gz

nibabies/interfaces/maths.py

@@ -0,0 +1,84 @@
+"""A module for interfaces """
+import os
+
+import numpy as np
+from nipype.interfaces.base import File, SimpleInterface, TraitedSpec, traits
+from nipype.utils.filemanip import fname_presuffix
+
+
+class ClipInputSpec(TraitedSpec):
+    in_file = File(exists=True, mandatory=True, desc="Input imaging file")
+    out_file = File(desc="Output file name")
+    minimum = traits.Float(
+        -np.inf, usedefault=True, desc="Values under minimum are set to minimum"
+    )
+    maximum = traits.Float(np.inf, usedefault=True, desc="Values over maximum are set to maximum")
+
+
+class ClipOutputSpec(TraitedSpec):
+    out_file = File(desc="Output file name")
+
+
+class Clip(SimpleInterface):
+    """Simple clipping interface that clips values to specified minimum/maximum
+    If no values are outside the bounds, nothing is done and the in_file is passed
+    as the out_file without copying.
+    """
+
+    input_spec = ClipInputSpec
+    output_spec = ClipOutputSpec
+
+    def _run_interface(self, runtime):
+        import nibabel as nb
+
+        img = nb.load(self.inputs.in_file)
+        data = img.get_fdata()
+
+        out_file = self.inputs.out_file
+        if out_file:
+            out_file = os.path.join(runtime.cwd, out_file)
+
+        if np.any((data < self.inputs.minimum) | (data > self.inputs.maximum)):
+            if not out_file:
+                out_file = fname_presuffix(
+                    self.inputs.in_file, suffix="_clipped", newpath=runtime.cwd
+                )
+            np.clip(data, self.inputs.minimum, self.inputs.maximum, out=data)
+            img.__class__(data, img.affine, img.header).to_filename(out_file)
+        elif not out_file:
+            out_file = self.inputs.in_file
+
+        self._results["out_file"] = out_file
+        return runtime
+
+
+class Label2MaskInputSpec(TraitedSpec):
+    in_file = File(exists=True, mandatory=True, desc="Input label file")
+    label_val = traits.Int(mandatory=True, dec="Label value to create mask from")
+
+
+class Label2MaskOutputSpec(TraitedSpec):
+    out_file = File(desc="Output file name")
+
+
+class Label2Mask(SimpleInterface):
+    """Create mask file for a label from a multi-label segmentation"""
+
+    input_spec = Label2MaskInputSpec
+    output_spec = Label2MaskOutputSpec
+
+    def _run_interface(self, runtime):
+        import nibabel as nb
+
+        img = nb.load(self.inputs.in_file)
+
+        mask = np.uint16(img.dataobj) == self.inputs.label_val
+        out_img = img.__class__(mask, img.affine, img.header)
+        out_img.set_data_dtype(np.uint8)
+
+        out_file = fname_presuffix(self.inputs.in_file, suffix="_mask", newpath=runtime.cwd)
+
+        out_img.to_filename(out_file)
+
+        self._results["out_file"] = out_file
+        return runtime

nibabies/interfaces/reports.py

@@ -19,6 +19,7 @@
     isdefined,
     traits,
 )
+from niworkflows.interfaces.reportlets import base as nrb
 from smriprep.interfaces.freesurfer import ReconAll
 
 LOGGER = logging.getLogger("nipype.interface")
@@ -296,6 +297,49 @@ def _generate_segment(self):
         )
 
 
+class LabeledHistogramInputSpec(nrb._SVGReportCapableInputSpec):
+    in_file = traits.File(exists=True, mandatory=True, desc="Image containing values to plot")
+    label_file = traits.File(
+        exists=True,
+        desc="Mask or label image where non-zero values will be used to extract data from in_file",
+    )
+    mapping = traits.Dict(desc="Map integer label values onto names of voxels")
+    xlabel = traits.Str("voxels", usedefault=True, desc="Description of values plotted")
+
+
+class LabeledHistogram(nrb.ReportingInterface):
+    input_spec = LabeledHistogramInputSpec
+
+    def _generate_report(self):
+        import nibabel as nb
+        import numpy as np
+        import seaborn as sns
+        from matplotlib import pyplot as plt
+        from nilearn.image import resample_to_img
+
+        report_file = self._out_report
+        img = nb.load(self.inputs.in_file)
+        data = img.get_fdata(dtype=np.float32)
+
+        if self.inputs.label_file:
+            label_img = nb.load(self.inputs.label_file)
+            if label_img.shape != img.shape[:3] or not np.allclose(label_img.affine, img.affine):
+                label_img = resample_to_img(label_img, img, interpolation="nearest")
+            labels = np.uint16(label_img.dataobj)
+        else:
+            labels = np.uint8(data > 0)
+
+        uniq_labels = np.unique(labels[labels > 0])
+        label_map = self.inputs.mapping or {label: label for label in uniq_labels}
+
+        rois = {label_map.get(label, label): data[labels == label] for label in label_map}
+        with sns.axes_style('whitegrid'):
+            fig = sns.histplot(rois, bins=50)
+            fig.set_xlabel(self.inputs.xlabel)
+        plt.savefig(report_file)
+        plt.close()
+
+
 def get_world_pedir(ornt, pe_direction):
     """Return world direction of phase encoding"""
     # TODO: Move to niworkflows

nibabies/workflows/bold/base.py

@@ -69,7 +69,7 @@
     init_bold_surf_wf,
 )
 from .stc import init_bold_stc_wf
-from .t2s import init_bold_t2s_wf
+from .t2s import init_bold_t2s_wf, init_t2s_reporting_wf
 
 
 def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=None):
@@ -346,10 +346,9 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
     outputnode = pe.Node(
         niu.IdentityInterface(
             fields=[
+                "anat2bold_xfm",
                 "bold_anat",
                 "bold_anat_ref",
-                "bold2anat_xfm",
-                "anat2bold_xfm",
                 "bold_mask_anat",
                 "bold_aseg_anat",
                 "bold_aparc_anat",
@@ -363,6 +362,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
                 "bold_mask_native",
                 "bold_echos_native",
                 "bold_cifti",
+                "bold2anat_xfm",
                 "cifti_variant",
                 "cifti_metadata",
                 "cifti_density",
@@ -372,6 +372,10 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
                 "melodic_mix",
                 "nonaggr_denoised_file",
                 "confounds_metadata",
+                "t2star_bold",
+                "t2star_t1",
+                "t2star_std",
+                "hmc_xforms",
             ]
         ),
         name="outputnode",
@@ -433,6 +437,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
             ('bold_anat_ref', 'inputnode.bold_t1_ref'),
             ('bold2anat_xfm', 'inputnode.bold2anat_xfm'),
             ('anat2bold_xfm', 'inputnode.anat2bold_xfm'),
+            ('hmc_xforms', 'inputnode.hmc_xforms'),
             ('bold_aseg_anat', 'inputnode.bold_aseg_t1'),
             ('bold_aparc_anat', 'inputnode.bold_aparc_t1'),
             ('bold_mask_anat', 'inputnode.bold_mask_t1'),
@@ -449,6 +454,9 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
             ('cifti_variant', 'inputnode.cifti_variant'),
             ('cifti_metadata', 'inputnode.cifti_metadata'),
             ('cifti_density', 'inputnode.cifti_density'),
+            ('t2star_bold', 'inputnode.t2star_bold'),
+            ('t2star_t1', 'inputnode.t2star_t1'),
+            ('t2star_std', 'inputnode.t2star_std'),
             ('confounds_metadata', 'inputnode.confounds_metadata'),
             ('acompcor_masks', 'inputnode.acompcor_masks'),
             ('tcompcor_mask', 'inputnode.tcompcor_mask'),
@@ -547,8 +555,31 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
             name="bold_t2smap_wf",
         )
 
+        t2s_reporting_wf = init_t2s_reporting_wf()
+
+        ds_report_t2scomp = pe.Node(
+            DerivativesDataSink(
+                desc="t2scomp",
+                datatype="figures",
+                dismiss_entities=("echo",),
+            ),
+            name="ds_report_t2scomp",
+            run_without_submitting=True,
+        )
+
+        ds_report_t2star_hist = pe.Node(
+            DerivativesDataSink(
+                desc="t2starhist",
+                datatype="figures",
+                dismiss_entities=("echo",),
+            ),
+            name="ds_report_t2star_hist",
+            run_without_submitting=True,
+        )
+
     bold_final = pe.Node(
-        niu.IdentityInterface(fields=["bold", "boldref", "mask", "bold_echos"]), name="bold_final"
+        niu.IdentityInterface(fields=["bold", "boldref", "mask", "bold_echos", "t2star"]),
+        name="bold_final",
     )
 
     # Mask input BOLD reference image
@@ -580,6 +611,9 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
             ('bold_ref', 'inputnode.raw_ref_image')]),
         (validate_bolds, bold_hmc_wf, [
             (("out_file", pop_file), 'inputnode.bold_file')]),
+        (bold_hmc_wf, outputnode, [
+            ("outputnode.xforms", "hmc_xforms"),
+        ]),
         # Native-space BOLD files
         (final_boldref_wf, final_boldref_mask, [('outputnode.epi_ref_file', 'in_file')]),
         (final_boldref_wf, bold_final, [('outputnode.epi_ref_file', 'boldref')]),
@@ -589,6 +623,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
             ('boldref', 'bold_native_ref'),
             ('mask', 'bold_mask_native'),
             ('bold_echos', 'bold_echos_native'),
+            ('t2star', 't2star_bold'),
         ]),
         # EPI-T1 registration workflow
         (inputnode, bold_reg_wf, [
@@ -668,15 +703,22 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
         workflow.connect([
             # update name source for optimal combination
             (inputnode, func_derivatives_wf, [
-                (('bold_file', combine_meepi_source), 'inputnode.source_file')]),
-            (join_echos, bold_t2s_wf, [
-                ('bold_files', 'inputnode.bold_file')]),
-            (bold_t2s_wf, split_opt_comb, [
-                ('outputnode.bold', 'in_file')]),
-            (split_opt_comb, bold_t1_trans_wf, [
-                ("out_files", "inputnode.bold_split")]),
+                (("bold_file", combine_meepi_source), "inputnode.source_file"),
+            ]),
+            (join_echos, bold_t2s_wf, [("bold_files", "inputnode.bold_file")]),
             (join_echos, bold_final, [("bold_files", "bold_echos")]),
-            (bold_t2s_wf, bold_final, [("outputnode.bold", "bold")]),
+            (bold_t2s_wf, split_opt_comb, [("outputnode.bold", "in_file")]),
+            (split_opt_comb, bold_t1_trans_wf, [("out_files", "inputnode.bold_split")]),
+            (bold_t2s_wf, bold_final, [("outputnode.bold", "bold"),
+                                       ("outputnode.t2star_map", "t2star")]),
+            (inputnode, t2s_reporting_wf, [("anat_dseg", "inputnode.label_file")]),
+            (bold_reg_wf, t2s_reporting_wf, [
+                ("outputnode.itk_t1_to_bold", "inputnode.label_bold_xform")
+            ]),
+            (bold_final, t2s_reporting_wf, [("t2star", "inputnode.t2star_file"),
+                                            ("boldref", "inputnode.boldref")]),
+            (t2s_reporting_wf, ds_report_t2scomp, [('outputnode.t2s_comp_report', 'in_file')]),
+            (t2s_reporting_wf, ds_report_t2star_hist, [("outputnode.t2star_hist", "in_file")]),
         ])
         # fmt:on
 
@@ -704,6 +746,23 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
         ])
         # fmt:on
 
+        if multiecho:
+            t2star_to_t1w = pe.Node(
+                ApplyTransforms(interpolation="LanczosWindowedSinc", float=True),
+                name="t2star_to_t1w",
+                mem_gb=0.1,
+            )
+            # fmt:off
+            workflow.connect([
+                (bold_reg_wf, t2star_to_t1w, [("outputnode.itk_bold_to_t1", "transforms")]),
+                (bold_t1_trans_wf, t2star_to_t1w, [
+                    ("outputnode.bold_mask_t1", "reference_image")
+                ]),
+                (bold_final, t2star_to_t1w, [("t2star", "input_image")]),
+                (t2star_to_t1w, outputnode, [("output_image", "t2star_t1")]),
+            ])
+            # fmt:on
+
     if spaces.get_spaces(nonstandard=False, dim=(3,)):
         # Apply transforms in 1 shot
         # Only use uncompressed output if AROMA is to be run
@@ -712,6 +771,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
             mem_gb=mem_gb["resampled"],
             omp_nthreads=omp_nthreads,
             spaces=spaces,
+            multiecho=multiecho,
             name="bold_std_trans_wf",
             use_compression=not config.execution.low_mem,
         )
@@ -727,6 +787,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
                 ('anat_aparc', 'inputnode.bold_aparc')]),
             (bold_final, bold_std_trans_wf, [
                 ('mask', 'inputnode.bold_mask'),
+                ('t2star', 'inputnode.t2star'),
             ]),
             (bold_reg_wf, bold_std_trans_wf, [
                 ('outputnode.itk_bold_to_t1', 'inputnode.itk_bold_to_t1')]),
@@ -761,8 +822,8 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non
         else:
             # fmt:off
             workflow.connect([
-                (split_opt_comb, bold_std_trans_wf, [
-                    ('out_files', 'inputnode.bold_split')])
+                (split_opt_comb, bold_std_trans_wf, [('out_files', 'inputnode.bold_split')]),
+                (bold_std_trans_wf, outputnode, [("outputnode.t2star_std", "t2star_std")]),
             ])
             # fmt:on