Merge branch 'bf/coefficient_field' into 'master'

Bf/coefficient field

See merge request fsl/fslpy!252

CHANGELOG.rst

@@ -2,8 +2,8 @@ This document contains the ``fslpy`` release history in reverse chronological
 order.
 
 
-3.3.0 (Under development)
--------------------------
+3.3.0 (Tuesday 22nd September 2020)
+-----------------------------------
 
 
 Added
@@ -13,6 +13,7 @@ Added
   ``imtest``, ``fsl_abspath``, ``remove_ext``, ``Text2Vest``, and
   ``Vest2Text``.
 * New :func:`.gps` function, wrapping the FSL ``gps`` command.
+* New :func:`.vest.loadVestFile` and :func:`.vest.generateVest` functions.
 
 
 Changed
@@ -22,11 +23,13 @@ Changed
 * Updates to the BIDS filetree specification.
 
 
-Changed
-^^^^^^^
+Fixed
+^^^^^
 
 
-* Adjustments to tests and documentation
+* The :class:`.CoefficientField` class now works with alternate reference
+  images (i.e. a reference image with different dimensions to that which
+  was originally used when the non-linear transformation was calculated).
 
 
 3.2.2 (Thursday 9th July 2020)

fsl/transform/fnirt.py

@@ -225,6 +225,27 @@ def _readFnirtCoefficientField(fname, img, src, ref):
     # image voxel coordinates.
     fieldToRefMat = affine.scaleOffsetXform(knotSpacing, 0)
 
+    # But if the provided reference has
+    # different resolution to the
+    # reference that was originally
+    # used to generate the warp field,
+    # we need to adjust the field
+    # accordingly. We assume that the
+    # references are aligned in the FSL
+    # coordinate system, so simply apply
+    # a scaling factor calculated by
+    # dividing the original reference
+    # pixdims by the provided reference
+    # pixdims.
+    refPixdims = np.array([img.header['intent_p1'],
+                           img.header['intent_p2'],
+                           img.header['intent_p3']])
+
+    if not np.all(np.isclose(ref.pixdim[:3], refPixdims)):
+        fieldToRefMat = affine.concat(
+            affine.scaleOffsetXform(refPixdims / ref.pixdim[:3], 0),
+            fieldToRefMat)
+
     return nonlinear.CoefficientField(fname,
                                       src,
                                       ref,
@@ -236,7 +257,7 @@ def _readFnirtCoefficientField(fname, img, src, ref):
                                       fieldToRefMat=fieldToRefMat)
 
 
-def readFnirt(fname, src, ref, defType=None):
+def readFnirt(fname, src, ref, defType=None, intent=None):
     """Reads a non-linear FNIRT transformation image, returning
     a :class:`.DeformationField` or :class:`.CoefficientField` depending
     on the file type.
@@ -247,28 +268,38 @@ def readFnirt(fname, src, ref, defType=None):
     :arg defType:  Deformation type - either ``'absolute'`` or ``'relative'``.
                    Only used if the file is a deformation field. If not
                    provided, is automatically inferred from the data.
+    :arg intent:   NIFTI intent code of ``fname``. e.g.
+                   :attr:`.constants.FSL_FNIRT_DISPLACEMENT_FIELD`. If not
+                   provided, the intent is read from the image header.
     """
 
-    # Figure out whether the file
-    # is a deformation field or
-    # a coefficient field
-
-    img   = fslimage.Image(fname, loadData=False)
-    disps = (constants.FSL_FNIRT_DISPLACEMENT_FIELD,
-             constants.FSL_TOPUP_FIELD)
-    coefs = (constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
-             constants.FSL_DCT_COEFFICIENTS,
-             constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS,
-             constants.FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS,
-             constants.FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS)
+    if defType not in (None, 'absolute', 'relative'):
+        raise ValueError('defType must be None, "absolute" or "relative" '
+                         '(passed in as {})'.format(defType))
+
+    # Figure out whether the file is a
+    # deformation field or a coefficient
+    # field by checking the intent code.
+    # If the intent is provided, assume
+    # that the caller knows the type of
+    # the field.
+    img    = fslimage.Image(fname, loadData=False)
+    intent = intent or img.intent
+    disps  = (constants.FSL_FNIRT_DISPLACEMENT_FIELD,
+              constants.FSL_TOPUP_FIELD)
+    coefs  = (constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
+              constants.FSL_DCT_COEFFICIENTS,
+              constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS,
+              constants.FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS,
+              constants.FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS)
 
-    if img.intent in disps:
+    if intent in disps:
         return _readFnirtDeformationField(fname, img, src, ref, defType)
-    elif img.intent in coefs:
+    elif intent in coefs:
         return _readFnirtCoefficientField(fname, img, src, ref)
     else:
-        raise ValueError('Cannot determine type of nonlinear '
-                         'file {}'.format(fname))
+        raise ValueError('Cannot determine type of nonlinear warp field '
+                         '{} (intent code: {})'.format(fname, intent))
 
 
 def toFnirt(field):

fsl/transform/nonlinear.py

@@ -496,13 +496,11 @@ class CoefficientField(NonLinearTransform):
 
         fdata      = self.data
         nx, ny, nz = self.shape[:3]
-        ix, iy, iz = self.ref.shape[:3]
 
         # Convert the given voxel coordinates
         # into the corresponding coefficient
         # field voxel coordinates
-        x, y, z    = coords.T
-        i, j, k    = affine.transform(coords, self.refToFieldMat).T
+        i, j, k = affine.transform(coords, self.refToFieldMat).T
 
         # i, j, k: coefficient field indices
         # u, v, w: position of the ref voxel
@@ -582,6 +580,10 @@ def convertDeformationType(field, defType=None):
         if field.deformationType == 'absolute': defType = 'relative'
         else:                                   defType = 'absolute'
 
+    if defType not in ('absolute', 'relative'):
+        raise ValueError('defType must be "absolute" or "relative" '
+                         '("{}" passed)'.format(defType))
+
     # Regardless of the conversion direction,
     # we need the coordinates of every voxel
     # in the reference coordinate system.
@@ -601,8 +603,8 @@ def convertDeformationType(field, defType=None):
     # assumed to be) the relative shift. Or,
     # to convert from absolute to relative,
     # we subtract the reference image voxels.
-    if   defType == 'absolute': return field.data + coords
-    elif defType == 'relative': return field.data - coords
+    if defType == 'absolute': return field.data + coords
+    else:                     return field.data - coords
 
 
 def convertDeformationSpace(field, from_, to):

tests/test_transform/test_fnirt.py

@@ -9,10 +9,13 @@ import os.path as op
 import itertools as it
 
 import numpy as np
+import nibabel as nib
 
 import pytest
 
 import fsl.data.image          as fslimage
+import fsl.utils.tempdir       as tempdir
+import fsl.data.constants      as constants
 import fsl.transform.affine    as affine
 import fsl.transform.nonlinear as nonlinear
 import fsl.transform.fnirt     as fnirt
@@ -51,6 +54,39 @@ def test_readFnirt():
     assert disp.refSpace == 'fsl'
 
 
+def test_readFnirt_defType_intent():
+    src  = op.join(datadir, 'src.nii.gz')
+    ref  = op.join(datadir, 'ref.nii.gz')
+    coef = op.join(datadir, 'coefficientfield.nii.gz')
+    disp = op.join(datadir, 'displacementfield.nii.gz')
+
+    src  = fslimage.Image(src)
+    ref  = fslimage.Image(ref)
+
+    field = fnirt.readFnirt(disp, src, ref, defType='absolute')
+    assert field.deformationType == 'absolute'
+    field = fnirt.readFnirt(disp, src, ref, defType='relative')
+    assert field.deformationType == 'relative'
+
+    img = nib.load(coef)
+    img.header['intent_code'] = 0
+    with tempdir.tempdir():
+        img.to_filename('field.nii.gz')
+
+        with pytest.raises(ValueError):
+            fnirt.readFnirt('field', src, ref)
+
+        field = fnirt.readFnirt(
+            'field', src, ref,
+            intent=constants.FSL_CUBIC_SPLINE_COEFFICIENTS)
+        assert isinstance(field, nonlinear.CoefficientField)
+
+        field = fnirt.readFnirt(
+            'field', src, ref,
+            intent=constants.FSL_FNIRT_DISPLACEMENT_FIELD)
+        assert isinstance(field, nonlinear.DeformationField)
+
+
 def test_toFnirt():
 
     def check(got, exp):

tests/test_transform/test_nonlinear.py

@@ -311,6 +311,38 @@ def test_CoefficientField_transform():
         assert np.all(np.isclose(gotnp, srccoordsnp[srcspace], **tol))
 
 
+def test_coefficientField_transform_altref():
+
+    # test coordinates (manually determined).
+    # original ref image is 2mm isotropic,
+    # resampled is 1mm. Each tuple contains:
+    #
+    # (src, ref2mm, ref1mm)
+    coords = [
+        ((18.414, 26.579, 25.599), (11, 19, 11), (22, 38, 22)),
+        ((14.727, 22.480, 20.340), ( 8, 17,  8), (16, 34, 16)),
+        ((19.932, 75.616, 27.747), (11, 45,  5), (22, 90, 10))
+    ]
+
+    nldir  = op.join(datadir, 'nonlinear')
+    src    = op.join(nldir, 'src.nii.gz')
+    ref    = op.join(nldir, 'ref.nii.gz')
+    cf     = op.join(nldir, 'coefficientfield.nii.gz')
+
+    src      = fslimage.Image(src)
+    ref2mm   = fslimage.Image(ref)
+    ref1mm   = ref2mm.adjust((1, 1, 1))
+    cfref2mm = fnirt.readFnirt(cf, src, ref2mm)
+    cfref1mm = fnirt.readFnirt(cf, src, ref1mm)
+
+    for srcc, ref2mmc, ref1mmc in coords:
+        ref2mmc = cfref2mm.transform(ref2mmc, 'voxel', 'voxel')
+        ref1mmc = cfref1mm.transform(ref1mmc, 'voxel', 'voxel')
+
+        assert np.all(np.isclose(ref2mmc, srcc, 1e-4))
+        assert np.all(np.isclose(ref1mmc, srcc, 1e-4))
+
+
 def test_coefficientFieldToDeformationField():
 
     nldir = op.join(datadir, 'nonlinear')