diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index aafbf6355b4af1fad83324291892b86829b642b1..2630bdf90fb8990aa2cb351c377482cbc11923ae 100644
--- a/CHANGELOG.rst
+++ b/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)
diff --git a/fsl/transform/fnirt.py b/fsl/transform/fnirt.py
index 4fa4051c86cd63524a4cca57120c5c1adf5266f9..36dbda8c66f270d03148dd1f5d9a3d521115a46b 100644
--- a/fsl/transform/fnirt.py
+++ b/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):
diff --git a/fsl/transform/nonlinear.py b/fsl/transform/nonlinear.py
index 77bc212b87d0fc71dfd3f4155fd5ed36a0d12f95..74788a3dc0ebf8b2bcd98a053837858975b027f9 100644
--- a/fsl/transform/nonlinear.py
+++ b/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):
diff --git a/tests/test_transform/test_fnirt.py b/tests/test_transform/test_fnirt.py
index eefc68b64e471dc052ff4fb92b637f56a599a8fa..27ceb40703e5e2506dffb762c59df562153feaf1 100644
--- a/tests/test_transform/test_fnirt.py
+++ b/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):
diff --git a/tests/test_transform/test_nonlinear.py b/tests/test_transform/test_nonlinear.py
index 323b6aff94f3e741221a000c96d051fb367476b9..8bfb0c3aba0d7d0746cf48d56b19631dfc1c2536 100644
--- a/tests/test_transform/test_nonlinear.py
+++ b/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')