From 670811c691853bc39d099211665018b26de7a98d Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauld.mccarthy@gmail.com>
Date: Mon, 16 Jan 2017 16:34:55 +0000
Subject: [PATCH] Image has new voxelsToScaledVoxels method, which returns
 scaled voxel xform. Transform module has new flirtMatrixToSform function
 which turns a FLIRT matrix into a src -> struc sform transformation.

---
 fsl/data/image.py      | 23 +++++++++++++++++++++++
 fsl/utils/transform.py | 35 +++++++++++++++++++++++++++++++++++
 2 files changed, 58 insertions(+)

diff --git a/fsl/data/image.py b/fsl/data/image.py
index 0fc1f6ff0..ec8b9fd39 100644
--- a/fsl/data/image.py
+++ b/fsl/data/image.py
@@ -422,6 +422,29 @@ class Nifti(notifier.Notifier):
         return npla.det(self.__voxToWorldMat) > 0
 
 
+    @memoize.Instanceify(memoize.memoize)
+    def voxelsToScaledVoxels(self):
+        """Returns a transformation matrix which transforms from voxel
+        coordinates into scaled voxel coordinates, with a left-right flip
+        if the image appears to be stored in neurological order.
+
+        See http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT/FAQ#What_is_the\
+        _format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to\
+        _the_transformation_parameters.3F
+        """
+
+        shape          = list(self.shape[ :3])
+        pixdim         = list(self.pixdim[:3])
+        voxToPixdimMat = np.diag(pixdim + [1.0])
+        
+        if self.isNeurological():
+            x              = (shape[0] - 1) * pixdim[0]
+            flip           = transform.scaleOffsetXform([-1, 1, 1], [x, 0, 0])
+            voxToPixdimMat = transform.concat(voxToPixdimMat, flip)
+
+        return voxToPixdimMat
+
+
     def sameSpace(self, other):
         """Returns ``True`` if the ``other`` image (assumed to be a
         :class:`Nifti` instance) has the same dimensions and is in the
diff --git a/fsl/utils/transform.py b/fsl/utils/transform.py
index df82ddc55..f7bd2f21f 100644
--- a/fsl/utils/transform.py
+++ b/fsl/utils/transform.py
@@ -19,6 +19,7 @@ spaces. The following functions are provided:
    rotMatToAxisAngles
    axisAnglesToRotMat
    axisBounds
+   flirtMatrixToSform
 """
 
 import numpy        as np
@@ -378,3 +379,37 @@ def _fillPoints(p, axes):
         newp[:, ax] = p[:, i]
 
     return newp
+
+
+def flirtMatrixToSform(flirtMat, srcImage, refImage):
+    """Converts the given ``FLIRT`` transformation matrix into a
+    transformation from the source image voxel coordinate system to
+    the reference image world coordinate system.
+
+    FLIRT transformation matrices transform from the source image scaled voxel
+    coordinate system into the reference image scaled voxel coordinate system
+    (voxels scaled by pixdims, with a left-right flip if the image sform has a
+    positive determinant).
+
+    So to construct a transformation from source image voxel coordinates
+    into reference image world coordinates, we need to combine the following:
+
+      1. Source voxels -> Source scaled voxels
+      2. Source scaled voxels -> Reference scaled voxels (the FLIRT matrix)
+      3. Reference scaled voxels -> Reference voxels
+      4. Reference voxels -> Reference world (the reference image sform)
+
+    :arg flirtMat: A ``(4, 4)`` transformation matrix
+    :arg srcImage: Source :class:`.Image`
+    :arg refImage: Reference :class:`.Image`
+    """
+    
+    srcScaledVoxelMat    = srcImage.voxelsToScaledVoxels()
+    refScaledVoxelMat    = refImage.voxelsToScaledVoxels()
+    refVoxToWorldMat     = refImage.voxToWorldMat
+    refInvScaledVoxelMat = invert(refScaledVoxelMat)
+
+    return concat(refVoxToWorldMat,
+                  refInvScaledVoxelMat,
+                  flirtMat,
+                  srcScaledVoxelMat)
-- 
GitLab