ENH: Replace avscale with non-fsl tools

.zenodo.json

@@ -30,6 +30,11 @@
       "affiliation": "Department of Psychology, University of Oslo",
       "orcid": "0000-0003-2574-1705"
     },
+    {
+      "name": "Humphries, Joseph",
+      "affiliation": "Turing Medical",
+      "orcid": "0000-0002-10257956"
+    },
     {
       "name": "Krause, Michael",
       "affiliation": "Max Planck Institute for Human Development, Berlin, Germany",

qsiprep/interfaces/gradients.py

@@ -16,6 +16,8 @@
 from dipy.reconst.dti import decompose_tensor
 from dipy.core.geometry import normalized_vector
 from sklearn.metrics import r2_score
+from transforms3d.affines import decompose44
+from scipy.spatial.transform import Rotation as R
 from .itk import disassemble_transform
 
 LOGGER = logging.getLogger('nipype.interface')
@@ -618,20 +620,45 @@ def get_fsl_motion_params(itk_file, src_file, ref_file, working_dir):
             src_file=src_file, ref_file=ref_file, itk_file=itk_file,
             fsl_file=tmp_fsl_file)
     os.system(fsl_convert_cmd)
-    proc = Popen(['avscale', '--allparams', tmp_fsl_file, src_file], stdout=PIPE,
-                 stderr=PIPE)
-    stdout, _ = proc.communicate()
 
     def get_measures(line):
         line = line.strip().split()
         return np.array([float(num) for num in line[-3:]])
 
-    lines = stdout.decode("utf-8").split("\n")
+    def get_image_center(src_fname):
+        #returns image center in mm
+        src_img = nb.load(src_fname)
+        src_aff = src_img.affine
+        src_center = (np.array(src_img.shape)-1)/2
+        src_center_mm = nb.affines.apply_affine(src_aff, src_center)
+        src_offsets = src_aff[0:3,3]
+        src_center_mm -= src_offsets
+        return src_center_mm
+
+    def get_trans_from_offset(image_center, rotmat):
+        # offset[0] = trans[0] + center[0] - [rot[0,0]*center[0] +rot[0,1]*center[1] + rot[0,2]*center[2]]
+        trans = np.zeros((3,))
+        offsets = rotmat[0:3,3]
+        for i in range(3):
+            offpart = offsets[i] - image_center[i]
+            rotpart = rotmat[i,0]*image_center[0] + rotmat[i,1]*image_center[1] + rotmat[i,2]*image_center[2]
+            trans[i] = offpart + rotpart
+        return trans
+
+    img_center = get_image_center(src_file)
+    c3d_out_xfm = np.loadtxt(fname=tmp_fsl_file,dtype='float')
+    [T,Rotmat,Z,S] = decompose44(c3d_out_xfm)
+    T = get_trans_from_offset(img_center,c3d_out_xfm)
+
     flip = np.array([1, -1, -1])
-    rotation = get_measures(lines[6]) * flip
-    translation = get_measures(lines[8]) * flip
-    scale = get_measures(lines[10])
-    shear = get_measures(lines[12])
+    negflip = np.array([-1, 1, 1])
+    Rotmat_to_convert = R.from_matrix(Rotmat)
+    Rotvec = Rotmat_to_convert.as_rotvec()
+
+    rotation = Rotvec * negflip
+    translation = T * flip
+    scale = Z
+    shear = S
 
     return np.concatenate([scale, shear, rotation, translation])
 

setup.cfg

@@ -55,6 +55,7 @@ install_requires =
     scikit-image
     SimpleITK
     pyAFQ >= 1.0.1
+    transforms3d
 test_requires =
     coverage
     codecov