[RTM] ENH: Conformation transforms

.circle/data/ds005_outputs.txt

@@ -24,6 +24,7 @@ ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-MNI152NLin2009cAsym_class-WM_pro
 ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-MNI152NLin2009cAsym_dtissue.nii.gz
 ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-MNI152NLin2009cAsym_preproc.nii.gz
 ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-MNI152NLin2009cAsym_target-T1w_warp.h5
+ds005/out/fmriprep/sub-01/anat/sub-01_T1w_space-orig_target-T1w_affine.txt
 ds005/out/fmriprep/sub-01/anat/sub-01_T1w_target-fsnative_affine.txt
 ds005/out/fmriprep/sub-01/anat/sub-01_T1w_target-MNI152NLin2009cAsym_warp.h5
 ds005/out/fmriprep/sub-01/func

.circle/data/ds054_outputs.txt

@@ -16,6 +16,7 @@ ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-MNI152NLin2009cAsym_clas
 ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-MNI152NLin2009cAsym_dtissue.nii.gz
 ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-MNI152NLin2009cAsym_preproc.nii.gz
 ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-MNI152NLin2009cAsym_target-T1w_warp.h5
+ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_space-orig_target-T1w_affine.txt
 ds054/out/fmriprep/sub-100185/anat/sub-100185_T1w_target-MNI152NLin2009cAsym_warp.h5
 ds054/out/fmriprep/sub-100185/func
 ds054/out/fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_bold_confounds.tsv

fmriprep/interfaces/__init__.py

@@ -13,7 +13,7 @@
 )
 from .surf import NormalizeSurf, GiftiNameSource, GiftiSetAnatomicalStructure
 from .reports import SubjectSummary, FunctionalSummary, AboutSummary
-from .utils import TPM2ROI, AddTPMs, AddTSVHeader
+from .utils import TPM2ROI, AddTPMs, AddTSVHeader, ConcatAffines
 from .fmap import FieldEnhance
 from .confounds import GatherConfounds, ICAConfounds
 from .itk import MCFLIRT2ITK, MultiApplyTransforms

fmriprep/interfaces/freesurfer.py

@@ -48,19 +48,33 @@ class StructuralReference(fs.RobustTemplate):
         .../sub-01_ses-test_T1w.nii.gz --template mri_robust_template_out.mgz'
 
     """
-    @property
-    def cmdline(self):
-        cmd = super(StructuralReference, self).cmdline
-        if len(self.inputs.in_files) > 1:
-            return cmd
+
+    def _num_vols(self):
+        n_files = len(self.inputs.in_files)
+        if n_files != 1:
+            return n_files
 
         img = nb.load(self.inputs.in_files[0])
-        if len(img.shape) > 3 and img.shape[3] > 1:
-            return cmd
+        if len(img.shape) == 3:
+            return 1
+
+        return img.shape[3]
 
-        out_file = self._list_outputs()['out_file']
-        copyfile(self.inputs.in_files[0], out_file)
-        return "echo Only one time point!"
+    @property
+    def cmdline(self):
+        if self._num_vols() == 1:
+            return "echo Only one time point!"
+        return super(StructuralReference, self).cmdline
+
+    def _list_outputs(self):
+        outputs = super(StructuralReference, self)._list_outputs()
+        if self._num_vols() == 1:
+            in_file = self.inputs.in_files[0]
+            transform_file = outputs['transform_outputs'][0]
+            outputs['out_file'] = in_file
+            fs.utils.LTAConvert(in_lta='identity.nofile', source_file=in_file,
+                                target_file=in_file, out_lta=transform_file).run()
+        return outputs
 
 
 class MakeMidthicknessInputSpec(fs.utils.MRIsExpandInputSpec):

fmriprep/interfaces/images.py

@@ -198,15 +198,16 @@ def _run_interface(self, runtime):
 
 
 class ConformInputSpec(BaseInterfaceInputSpec):
-    in_file = File(exists=True, mandatory=True, desc='Input T1w image')
+    in_file = File(exists=True, mandatory=True, desc='Input image')
     target_zooms = traits.Tuple(traits.Float, traits.Float, traits.Float,
                                 desc='Target zoom information')
     target_shape = traits.Tuple(traits.Int, traits.Int, traits.Int,
                                 desc='Target shape information')
 
 
 class ConformOutputSpec(TraitedSpec):
-    out_file = File(exists=True, desc='Conformed T1w image')
+    out_file = File(exists=True, desc='Conformed image')
+    transform = File(exists=True, desc='Conformation transform')
 
 
 class Conform(SimpleInterface):
@@ -234,6 +235,13 @@ def _run_interface(self, runtime):
         zooms = np.array(reoriented.header.get_zooms()[:3])
         shape = np.array(reoriented.shape[:3])
 
+        # Reconstruct transform from orig to reoriented image
+        ornt_xfm = nb.orientations.inv_ornt_aff(
+            nb.io_orientation(orig_img.affine), orig_img.shape)
+        # Identity unless proven otherwise
+        target_affine = reoriented.affine.copy()
+        conform_xfm = np.eye(4)
+
         xyz_unit = reoriented.header.get_xyzt_units()[0]
         if xyz_unit == 'unknown':
             # Common assumption; if we're wrong, unlikely to be the only thing that breaks
@@ -247,12 +255,9 @@ def _run_interface(self, runtime):
         rescale = not np.allclose(zooms, target_zooms, atol=atol)
         resize = not np.all(shape == target_shape)
         if rescale or resize:
-            target_affine = np.eye(4, dtype=reoriented.affine.dtype)
             if rescale:
                 scale_factor = target_zooms / zooms
                 target_affine[:3, :3] = reoriented.affine[:3, :3].dot(np.diag(scale_factor))
-            else:
-                target_affine[:3, :3] = reoriented.affine[:3, :3]
 
             if resize:
                 # The shift is applied after scaling.
@@ -261,10 +266,9 @@ def _run_interface(self, runtime):
                 # Use integer shifts to avoid unnecessary interpolation
                 offset = (reoriented.affine[:3, 3] * size_factor - reoriented.affine[:3, 3])
                 target_affine[:3, 3] = reoriented.affine[:3, 3] + offset.astype(int)
-            else:
-                target_affine[:3, 3] = reoriented.affine[:3, 3]
 
             data = nli.resample_img(reoriented, target_affine, target_shape).get_data()
+            conform_xfm = np.linalg.inv(reoriented.affine).dot(target_affine)
             reoriented = reoriented.__class__(data, target_affine, reoriented.header)
 
         # Image may be reoriented, rescaled, and/or resized
@@ -274,7 +278,14 @@ def _run_interface(self, runtime):
         else:
             out_name = fname
 
+        transform = ornt_xfm.dot(conform_xfm)
+        assert np.allclose(orig_img.affine.dot(transform), target_affine)
+
+        mat_name = fname_presuffix(fname, suffix='.mat', newpath=runtime.cwd, use_ext=False)
+        np.savetxt(mat_name, transform, fmt='%.08f')
+
         self._results['out_file'] = out_name
+        self._results['transform'] = mat_name
 
         return runtime
 
@@ -286,6 +297,7 @@ class ReorientInputSpec(BaseInterfaceInputSpec):
 
 class ReorientOutputSpec(TraitedSpec):
     out_file = File(exists=True, desc='Reoriented T1w image')
+    transform = File(exists=True, desc='Reorientation transform')
 
 
 class Reorient(SimpleInterface):
@@ -299,14 +311,22 @@ def _run_interface(self, runtime):
         orig_img = nb.load(fname)
         reoriented = nb.as_closest_canonical(orig_img)
 
+        # Reconstruct transform from orig to reoriented image
+        ornt_xfm = nb.orientations.inv_ornt_aff(
+            nb.io_orientation(orig_img.affine), orig_img.shape)
+
         # Image may be reoriented
         if reoriented is not orig_img:
             out_name = fname_presuffix(fname, suffix='_ras', newpath=runtime.cwd)
             reoriented.to_filename(out_name)
         else:
             out_name = fname
 
+        mat_name = fname_presuffix(fname, suffix='.mat', newpath=runtime.cwd, use_ext=False)
+        np.savetxt(mat_name, ornt_xfm, fmt='%.08f')
+
         self._results['out_file'] = out_name
+        self._results['transform'] = mat_name
 
         return runtime
 

fmriprep/interfaces/utils.py

@@ -8,16 +8,18 @@
 
 """
 
+import os
 import numpy as np
 import nibabel as nb
 import scipy.ndimage as nd
 
 from niworkflows.nipype import logging
 from niworkflows.nipype.utils.filemanip import fname_presuffix
 from niworkflows.nipype.interfaces.base import (
-    traits, isdefined, TraitedSpec, BaseInterfaceInputSpec, File, InputMultiPath,
-    SimpleInterface
+    traits, isdefined, File, InputMultiPath,
+    TraitedSpec, DynamicTraitedSpec, BaseInterfaceInputSpec, SimpleInterface
 )
+from niworkflows.nipype.interfaces.io import add_traits
 
 IFLOGGER = logging.getLogger('interfaces')
 
@@ -187,6 +189,44 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class ConcatAffinesInputSpec(DynamicTraitedSpec, BaseInterfaceInputSpec):
+    invert = traits.Bool(False, usedefault=True, desc='Invert output transform')
+
+
+class ConcatAffinesOutputSpec(TraitedSpec):
+    out_mat = File(exists=True, desc='Output transform')
+
+
+class ConcatAffines(SimpleInterface):
+    input_spec = ConcatAffinesInputSpec
+    output_spec = ConcatAffinesOutputSpec
+
+    def __init__(self, num_affines=0, *args, **kwargs):
+        super(ConcatAffines, self).__init__(*args, **kwargs)
+        self._num_affines = num_affines
+        trait_type = File(exists=True)
+        if num_affines == 0:
+            add_traits(self.inputs, ['mat_list'], trait_type)
+        elif num_affines < 26:
+            add_traits(self.inputs, self._get_names(num_affines), trait_type)
+
+    @staticmethod
+    def _get_names(num_affines):
+        A = ord('A') - 1
+        return ['mat_{}to{}'.format(chr(X), chr(X + 1))
+                for X in range(A + num_affines, A, -1)]
+
+    def _run_interface(self, runtime):
+        out_mat = os.path.join(runtime.cwd, 'concat.mat')
+        in_list = [self.inputs.get()[name] for name in self._get_names(self._num_affines)]
+
+        out_xfm = _concat_xfms(in_list, invert=self.inputs.invert)
+        np.savetxt(out_mat, out_xfm, fmt=str('%.12g'))
+
+        self._results['out_mat'] = out_mat
+        return runtime
+
+
 def _tpm2roi(in_tpm, in_mask, mask_erosion_mm=None, erosion_mm=None,
              mask_erosion_prop=None, erosion_prop=None, pthres=0.95):
     """
@@ -237,3 +277,15 @@ def _tpm2roi(in_tpm, in_mask, mask_erosion_mm=None, erosion_mm=None,
     roi_img.set_data_dtype(np.uint8)
     roi_img.to_filename(roi_fname)
     return roi_fname, eroded_mask_file or in_mask
+
+
+def _concat_xfms(in_list, invert):
+    transforms = [np.loadtxt(in_mat) for in_mat in in_list]
+    out_xfm = transforms.pop(0)
+    for xfm in transforms:
+        out_xfm = out_xfm.dot(xfm)
+
+    if invert:
+        out_xfm = np.linalg.inv(out_xfm)
+
+    return out_xfm