diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index ebb88428d9544339a06c4ba56797bbcfeb2bb763..54fcdefa352f5ea5e00a1fde425a01c3ed368d65 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -12,6 +12,18 @@ Added
* New :mod:`.image.roi` module, for exracting an ROI of an image, or expanding
its field-of-view.
+* New :meth:`.Image.getAffine` method, for retrieving an affine between any of
+ the voxel, FSL, or world coordinate systems.
+* New :mod:`fsl.transforms` package, which contains classes and functions for
+ working with linear and non-linear FLIRT and FNIRT transformations.
+* New static methods :meth:`.Nifti.determineShape`,
+ :meth:`.Nifti.determineAffine`, :meth:`.Nifti.generateAffines`, and
+ :meth:`.Nifti.identifyAffine`.
+* New prototype :mod:`fsl.transforms.x5` module, for reading/writing linear
+ and non-linear X5 files (*preliminary release, subject to change*).
+* New prototype :mod:`.fsl_convert_x5` :mod:`.fsl_apply_x5` programs, for
+ working with X5 transformations (*preliminary release, subject to change*).
+
Changed
@@ -20,6 +32,16 @@ Changed
* The :mod:`.resample_image` script has been updated to support resampling of
images with more than 3 dimensions.
+* `h5py <https://www.h5py.org/>`_ has been added to the ``fslpy`` dependencies.
+
+
+Deprecated
+^^^^^^^^^^
+
+
+* The :mod:`fsl.utils.transform` module has been deprecated its functions can
+ now be found in the :mod:`fsl.transforms.affine` and
+ :mod:`fsl.transform.flirt` modules.
2.3.1 (Friday July 5th 2019)
@@ -169,7 +191,7 @@ Added
* Simple built-in :mod:`.deprecated` decorator.
* New :mod:`fsl.data.utils` module, which currently contains one function
:func:`.guessType`, which guesses the data type of a file/directory path.
-* New `.commonBase` function for finding the common prefix of a set of
+* New :func:`.commonBase` function for finding the common prefix of a set of
file/directory paths.
@@ -209,7 +231,7 @@ Fixed
^^^^^
-* Fixed an issue with the `.dicom.loadSeries` using memory-mapping for
+* Fixed an issue with the :func:`.dicom.loadSeries` using memory-mapping for
image files that would subsequently be deleted.
* Fixed an issue in the :class:`.GiftiMesh` class, where
``numpy``/``nibabel`` was returning read-only index arrays.
diff --git a/README.rst b/README.rst
index ac2f6b895b93c2ab54cd6db2ec719754e7200bfe..6b0f069bb8212acae7a37cab9ba5dc18b0ba4ca7 100644
--- a/README.rst
+++ b/README.rst
@@ -91,6 +91,10 @@ The ``rtree`` library assumes that ``libspatialindex`` is installed on
your system.
+The :mod:`fsl.transform.x5` module uses `h5py <https://www.h5py.org/>`_,
+which requires ``libhdf5``.
+
+
Documentation
-------------
diff --git a/doc/_static/theme_overrides.css b/doc/_static/theme_overrides.css
new file mode 100644
index 0000000000000000000000000000000000000000..7c1a5202231086af91e2bd4b388d8de4b064f27f
--- /dev/null
+++ b/doc/_static/theme_overrides.css
@@ -0,0 +1,10 @@
+/* override table width restrictions */
+.wy-table-responsive table td, .wy-table-responsive table th {
+ white-space: normal;
+}
+
+.wy-table-responsive {
+ margin-bottom: 24px;
+ max-width: 100%;
+ overflow: visible;
+}
diff --git a/doc/conf.py b/doc/conf.py
index a7fb2b4d202b8f8256f827ccad8cf34c42860d32..90254c22f1cee57d5fc65a2cdfe395e0a27f3a93 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -33,7 +33,8 @@ date = datetime.date.today()
# ones.
extensions = [
'sphinx.ext.autodoc',
- 'sphinx.ext.autosummary',
+ 'sphinx.ext.viewcode',
+ 'sphinx.ext.autosummary',
'sphinx.ext.mathjax',
'sphinx.ext.graphviz',
'sphinx.ext.todo',
@@ -121,6 +122,7 @@ pygments_style = 'sphinx'
# a list of builtin themes.
html_theme = 'sphinx_rtd_theme'
+
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
# documentation.
@@ -148,7 +150,13 @@ html_theme = 'sphinx_rtd_theme'
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = []
+html_static_path = ['_static']
+
+html_context = {
+ 'css_files': [
+ '_static/theme_overrides.css', # overrides for wide tables in RTD theme
+ ],
+ }
# Add any extra paths that contain custom files (such as robots.txt or
# .htaccess) here, relative to this directory. These files are copied
@@ -354,7 +362,7 @@ epub_exclude_files = ['search.html']
autoclass_content = 'class'
# Document private members and special members (e.g. __init__)
-autodocsourc_default_flags = ['private-members', 'special-members']
+autodoc_default_flags = ['private-members', 'special-members']
# Documentation for python modules is in the same order
# as the source code.
@@ -368,7 +376,7 @@ def autodoc_skip_member(app, what, name, obj, skip, options):
if what == 'class':
attName = name.split('.')[-1]
return skip or attName.startswith('_sync_')
-
+
return skip or False
diff --git a/doc/contributing.rst b/doc/contributing.rst
index 4e76856b7505fe98a204090b96afa0dcb6cd825a..ef710683b31c8094be90fe39a44d85d2e3e20be4 100644
--- a/doc/contributing.rst
+++ b/doc/contributing.rst
@@ -34,7 +34,7 @@ the following labels (this convention has been inherited from `nibabel
* *BF* : bug fix
* *RF* : refactoring
- * *NF* : new feature
+ * *ENH*: enhancement/new feature
* *BW* : addresses backward-compatibility
* *OPT* : optimization
* *BK* : breaks something and/or tests fail
@@ -102,7 +102,7 @@ Unit and integration tests are currently run with ``py.test`` and
``coverage``.
- Aim for 100% code coverage.
-- Tests must pass on python 2.7, 3.4, 3.5, and 3.6
+- Tests must pass on python 3.5, 3.6, and 3.7.
Coding conventions
diff --git a/doc/fsl.rst b/doc/fsl.rst
index fb702bca8580962ba9fd1da2a06c2e895efa767a..c92138a86ddb95c109a747fa3a4624d6d74613d5 100644
--- a/doc/fsl.rst
+++ b/doc/fsl.rst
@@ -7,6 +7,7 @@
fsl.data
fsl.scripts
fsl.utils
+ fsl.transform
fsl.version
fsl.wrappers
diff --git a/doc/fsl.scripts.fsl_convert_x5.rst b/doc/fsl.scripts.fsl_convert_x5.rst
new file mode 100644
index 0000000000000000000000000000000000000000..827d483b40e1ef773866243d763abd32df530c27
--- /dev/null
+++ b/doc/fsl.scripts.fsl_convert_x5.rst
@@ -0,0 +1,7 @@
+``fsl.scripts.fsl_convert_x5``
+==============================
+
+.. automodule:: fsl.scripts.fsl_convert_x5
+ :members:
+ :undoc-members:
+ :show-inheritance:
diff --git a/doc/fsl.scripts.rst b/doc/fsl.scripts.rst
index ece0e2b6e0731daa11612a583244a955bfcf0f99..743a3c12db30828d94c0ec0559a03b1248ae5744 100644
--- a/doc/fsl.scripts.rst
+++ b/doc/fsl.scripts.rst
@@ -5,6 +5,7 @@
:hidden:
fsl.scripts.atlasq
+ fsl.scripts.fsl_convert_x5
fsl.scripts.fsl_ents
fsl.scripts.imcp
fsl.scripts.imglob
diff --git a/doc/fsl.transform.affine.rst b/doc/fsl.transform.affine.rst
new file mode 100644
index 0000000000000000000000000000000000000000..93792be10f033d6b36801012d2a223dd05d75409
--- /dev/null
+++ b/doc/fsl.transform.affine.rst
@@ -0,0 +1,7 @@
+``fsl.transform.affine``
+========================
+
+.. automodule:: fsl.transform.affine
+ :members:
+ :undoc-members:
+ :show-inheritance:
diff --git a/doc/fsl.utils.transform.rst b/doc/fsl.transform.flirt.rst
similarity index 57%
rename from doc/fsl.utils.transform.rst
rename to doc/fsl.transform.flirt.rst
index b21d270e979bf113d60e5966879b1586d4292d65..836a82f49332eecf65e5ce70ca1778619a1f20d0 100644
--- a/doc/fsl.utils.transform.rst
+++ b/doc/fsl.transform.flirt.rst
@@ -1,7 +1,7 @@
-``fsl.utils.transform``
+``fsl.transform.flirt``
=======================
-.. automodule:: fsl.utils.transform
+.. automodule:: fsl.transform.flirt
:members:
:undoc-members:
:show-inheritance:
diff --git a/doc/fsl.transform.fnirt.rst b/doc/fsl.transform.fnirt.rst
new file mode 100644
index 0000000000000000000000000000000000000000..4a1efe50fc5ed17c2844bf49a8838bd467dc245e
--- /dev/null
+++ b/doc/fsl.transform.fnirt.rst
@@ -0,0 +1,7 @@
+``fsl.transform.fnirt``
+=======================
+
+.. automodule:: fsl.transform.fnirt
+ :members:
+ :undoc-members:
+ :show-inheritance:
diff --git a/doc/fsl.transform.nonlinear.rst b/doc/fsl.transform.nonlinear.rst
new file mode 100644
index 0000000000000000000000000000000000000000..d137a81f30c54cb5e436047307feca854353cea6
--- /dev/null
+++ b/doc/fsl.transform.nonlinear.rst
@@ -0,0 +1,7 @@
+``fsl.transform.nonlinear``
+===========================
+
+.. automodule:: fsl.transform.nonlinear
+ :members:
+ :undoc-members:
+ :show-inheritance:
diff --git a/doc/fsl.transform.rst b/doc/fsl.transform.rst
new file mode 100644
index 0000000000000000000000000000000000000000..fdbd859adaa5828574327c49463d3084d0cd9769
--- /dev/null
+++ b/doc/fsl.transform.rst
@@ -0,0 +1,16 @@
+``fsl.transform``
+=================
+
+.. toctree::
+ :hidden:
+
+ fsl.transform.affine
+ fsl.transform.flirt
+ fsl.transform.fnirt
+ fsl.transform.nonlinear
+ fsl.transform.x5
+
+.. automodule:: fsl.transform
+ :members:
+ :undoc-members:
+ :show-inheritance:
diff --git a/doc/fsl.transform.x5.rst b/doc/fsl.transform.x5.rst
new file mode 100644
index 0000000000000000000000000000000000000000..9b6978cf742e613b850e3157aa78cc489661066e
--- /dev/null
+++ b/doc/fsl.transform.x5.rst
@@ -0,0 +1,7 @@
+``fsl.transform.x5``
+====================
+
+.. automodule:: fsl.transform.x5
+ :members:
+ :undoc-members:
+ :show-inheritance:
diff --git a/doc/fsl.utils.rst b/doc/fsl.utils.rst
index 5722c7eb4f726f3b17c1ce32f7edee2f85779e4a..f253bac7e11914f3d5ccebcd152a576cd326ce66 100644
--- a/doc/fsl.utils.rst
+++ b/doc/fsl.utils.rst
@@ -23,7 +23,6 @@
fsl.utils.run
fsl.utils.settings
fsl.utils.tempdir
- fsl.utils.transform
fsl.utils.weakfuncref
.. automodule:: fsl.utils
diff --git a/doc/images/fnirt_coefficient_field.png b/doc/images/fnirt_coefficient_field.png
new file mode 100644
index 0000000000000000000000000000000000000000..b4d5f548dfae121874b4e24b7a058c8e16fc4938
Binary files /dev/null and b/doc/images/fnirt_coefficient_field.png differ
diff --git a/doc/images/fnirt_coefficient_field.xcf b/doc/images/fnirt_coefficient_field.xcf
new file mode 100644
index 0000000000000000000000000000000000000000..0470fe8bd4bd418e9c2ffcd720ef80bcc3198a14
Binary files /dev/null and b/doc/images/fnirt_coefficient_field.xcf differ
diff --git a/doc/images/fnirt_warp_field.png b/doc/images/fnirt_warp_field.png
new file mode 100644
index 0000000000000000000000000000000000000000..734f30b5aa97e0c09f1512119daa5c8620b999f9
Binary files /dev/null and b/doc/images/fnirt_warp_field.png differ
diff --git a/doc/images/fnirt_warp_field.xcf b/doc/images/fnirt_warp_field.xcf
new file mode 100644
index 0000000000000000000000000000000000000000..2d1a9cca1b475652801c2344b54ab225c999ab61
Binary files /dev/null and b/doc/images/fnirt_warp_field.xcf differ
diff --git a/doc/images/nonlinear_registration_process.png b/doc/images/nonlinear_registration_process.png
new file mode 100644
index 0000000000000000000000000000000000000000..3232f02f0805a34c5a49f8bf353308f6b2fc13b8
Binary files /dev/null and b/doc/images/nonlinear_registration_process.png differ
diff --git a/doc/images/nonlinear_registration_process.xcf b/doc/images/nonlinear_registration_process.xcf
new file mode 100644
index 0000000000000000000000000000000000000000..6bd91d952323fe9135105878c442879adc87eed3
Binary files /dev/null and b/doc/images/nonlinear_registration_process.xcf differ
diff --git a/doc/mock_modules.txt b/doc/mock_modules.txt
index 8bb6c230a3289007c59d160f6b6ec1cb7572b809..a3d52e045b97b85ef44d551318e5bb84a5b2b6ff 100644
--- a/doc/mock_modules.txt
+++ b/doc/mock_modules.txt
@@ -1,4 +1,5 @@
deprecation
+h5py
nibabel
nibabel.fileslice
nibabel.freesurfer
@@ -6,3 +7,4 @@ numpy
numpy.linalg
scipy
scipy.ndimage
+scipy.ndimage.interpolation
diff --git a/fsl/__init__.py b/fsl/__init__.py
index d0a5a517ec659fafefaa8b1d582ff5369c2e55fc..1a3917dc763ec2e28bb5f21fe00ebe71154590fe 100644
--- a/fsl/__init__.py
+++ b/fsl/__init__.py
@@ -13,6 +13,7 @@ following sub-packages:
fsl.data
fsl.utils
fsl.scripts
+ fsl.transform
fsl.version
fsl.wrappers
@@ -22,4 +23,4 @@ following sub-packages:
for details.
"""
-__path__ = __import__('pkgutil').extend_path(__path__, __name__)
+__path__ = __import__('pkgutil').extend_path(__path__, __name__) # noqa
diff --git a/fsl/data/atlases.py b/fsl/data/atlases.py
index b28eb952c98ccbda29f2097ca8eda9ae660cba4c..8d5dad9bd679b4bd8786de4a5a8e41c21db3d310 100644
--- a/fsl/data/atlases.py
+++ b/fsl/data/atlases.py
@@ -54,7 +54,7 @@ import fsl.data.image as fslimage
import fsl.data.constants as constants
from fsl.utils.platform import platform as platform
import fsl.utils.image.resample as resample
-import fsl.utils.transform as transform
+import fsl.transform.affine as affine
import fsl.utils.notifier as notifier
import fsl.utils.settings as fslsettings
@@ -572,7 +572,7 @@ class AtlasDescription(object):
# Load the appropriate transformation matrix
# and transform all those voxel coordinates
# into world coordinates
- coords = transform.transform(coords, self.xforms[0])
+ coords = affine.transform(coords, self.xforms[0])
# Update the coordinates
# in our label objects
@@ -810,7 +810,7 @@ class LabelAtlas(Atlas):
"""
if not voxel:
- loc = transform.transform([loc], self.worldToVoxMat)[0]
+ loc = affine.transform([loc], self.worldToVoxMat)[0]
loc = [int(v) for v in loc.round()]
if loc[0] < 0 or \
@@ -983,7 +983,7 @@ class ProbabilisticAtlas(Atlas):
"""
if not voxel:
- loc = transform.transform([loc], self.worldToVoxMat)[0]
+ loc = affine.transform([loc], self.worldToVoxMat)[0]
loc = [int(v) for v in loc.round()]
if loc[0] < 0 or \
diff --git a/fsl/data/featanalysis.py b/fsl/data/featanalysis.py
index 7c38eb7773465fa79cc952e4038c38451e7e3b97..718ef679b0502337ecf5aba68a481b336fedd0e6 100644
--- a/fsl/data/featanalysis.py
+++ b/fsl/data/featanalysis.py
@@ -42,16 +42,15 @@ The following functions return the names of various files of interest:
"""
-import collections
-import logging
-import os.path as op
-import numpy as np
+import collections
+import logging
+import os.path as op
+import numpy as np
-import fsl.utils.path as fslpath
-import fsl.utils.transform as transform
-
-from . import image as fslimage
-from . import featdesign
+import fsl.utils.path as fslpath
+import fsl.transform.affine as affine
+from . import image as fslimage
+from . import featdesign
log = logging.getLogger(__name__)
@@ -467,9 +466,9 @@ def loadClusterResults(featdir, settings, contrast):
zcog = [c.zcogx, c.zcogy, c.zcogz]
copemax = [c.copemaxx, c.copemaxy, c.copemaxz]
- zmax = transform.transform([zmax], coordXform)[0].round()
- zcog = transform.transform([zcog], coordXform)[0].round()
- copemax = transform.transform([copemax], coordXform)[0].round()
+ zmax = affine.transform([zmax], coordXform)[0].round()
+ zcog = affine.transform([zcog], coordXform)[0].round()
+ copemax = affine.transform([copemax], coordXform)[0].round()
c.zmaxx, c.zmaxy, c.zmaxz = zmax
c.zcogx, c.zcogy, c.zcogz = zcog
diff --git a/fsl/data/image.py b/fsl/data/image.py
index 1bf901c287d5c2ba3ea033e11f08fdb79cc55e3d..dd423d4b5f8fc5341ea104378e5cf1d58e3aed1f 100644
--- a/fsl/data/image.py
+++ b/fsl/data/image.py
@@ -34,6 +34,7 @@ and file names:
import os
import os.path as op
+import itertools as it
import string
import logging
import tempfile
@@ -46,7 +47,7 @@ import nibabel as nib
import nibabel.fileslice as fileslice
import fsl.utils.meta as meta
-import fsl.utils.transform as transform
+import fsl.transform.affine as affine
import fsl.utils.notifier as notifier
import fsl.utils.memoize as memoize
import fsl.utils.path as fslpath
@@ -135,13 +136,32 @@ class Nifti(notifier.Notifier, meta.Meta):
methods.
- **The affine transformation**
+ **Affine transformations**
- The :meth:`voxToWorldMat` and :meth:`worldToVoxMat` attributes contain
- transformation matrices for transforming between voxel and world
- coordinates. The ``Nifti`` class follows the same process as ``nibabel``
- in selecting the affine (see
+ The ``Nifti`` class is aware of three coordinate systems:
+
+ - The ``voxel`` coordinate system, used to access image data
+
+ - The ``world`` coordinate system, where voxel coordinates are
+ transformed into a millimetre coordinate system, defined by the
+ ``sform`` and/or ``qform`` elements of the NIFTI header.
+
+ - The ``fsl`` coordinate system, where voxel coordinates are scaled by
+ the ``pixdim`` values in the NIFTI header, and the X axis is inverted
+ if the voxel-to-world affine has a positive determinant.
+
+
+ The :meth:`getAffine` method is a simple means of acquiring an affine
+ which will transform between any of these coordinate systems.
+
+
+ See `here <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>`_
+ for more details on the ``fsl`` coordinate system..
+
+
+ The ``Nifti`` class follows the same process as ``nibabel`` in determining
+ the ``voxel`` to ``world`` affine (see
http://nipy.org/nibabel/nifti_images.html#the-nifti-affines):
@@ -204,9 +224,11 @@ class Nifti(notifier.Notifier, meta.Meta):
=============== ========================================================
``'transform'`` The affine transformation matrix has changed. This topic
- will occur when the ``voxToWorldMat`` is changed.
+ will occur when the :meth:`voxToWorldMat` is changed.
+ ``'header'`` A header field has changed. This will occur when the
+ :meth:`intent` is changed.
=============== ========================================================
- """
+ """ # noqa
def __init__(self, header):
@@ -225,22 +247,59 @@ class Nifti(notifier.Notifier, meta.Meta):
raise ValueError('Unrecognised header: {}'.format(header))
header = header
- origShape, shape, pixdim = self.__determineShape(header)
+ origShape, shape, pixdim = Nifti.determineShape(header)
+ voxToWorldMat = Nifti.determineAffine(header)
+ affines, isneuro = Nifti.generateAffines(voxToWorldMat,
+ shape,
+ pixdim)
+
+ self.header = header
+ self.__shape = shape
+ self.__origShape = origShape
+ self.__pixdim = pixdim
+ self.__affines = affines
+ self.__isNeurological = isneuro
+
+
+ @staticmethod
+ def determineShape(header):
+ """This method is called by :meth:`__init__`. It figures out the actual
+ shape of the image data, and the zooms/pixdims for each data axis. Any
+ empty trailing dimensions are squeezed, but the returned shape is
+ guaranteed to be at least 3 dimensions. Returns:
+
+ - A sequence/tuple containing the image shape, as reported in the
+ header.
+ - A sequence/tuple containing the effective image shape.
+ - A sequence/tuple containing the zooms/pixdims.
+ """
- voxToWorldMat = self.__determineTransform(header)
- worldToVoxMat = transform.invert(voxToWorldMat)
+ # The canonicalShape function figures out
+ # the data shape that we should use.
+ origShape = list(header.get_data_shape())
+ shape = canonicalShape(origShape)
+ pixdims = list(header.get_zooms())
- self.header = header
- self.__shape = shape
- self.__intent = header.get('intent_code',
- constants.NIFTI_INTENT_NONE)
- self.__origShape = origShape
- self.__pixdim = pixdim
- self.__voxToWorldMat = voxToWorldMat
- self.__worldToVoxMat = worldToVoxMat
+ # if get_zooms() doesn't return at
+ # least len(shape) values, we'll
+ # fall back to the pixdim field in
+ # the header.
+ if len(pixdims) < len(shape):
+ pixdims = header['pixdim'][1:]
+ pixdims = pixdims[:len(shape)]
+
+ # should never happen, but if we only
+ # have zoom values for the original
+ # (< 3D) shape, pad them with 1s.
+ if len(pixdims) < len(shape):
+ pixdims = pixdims + [1] * (len(shape) - len(pixdims))
- def __determineTransform(self, header):
+ return origShape, shape, pixdims
+
+
+ @staticmethod
+ def determineAffine(header):
"""Called by :meth:`__init__`. Figures out the voxel-to-world
coordinate transformation matrix that is associated with this
``Nifti`` instance.
@@ -250,33 +309,40 @@ class Nifti(notifier.Notifier, meta.Meta):
# specially, as FNIRT clobbers the
# sform section of the NIFTI header
# to store other data.
- #
- # TODO Nibabel <= 2.1.0 has a bug in header.get
- # for fields with a value of 0. When this
- # bug gets fixed, we can replace the if-else
- # block below with this:
- #
- # intent = header.get('intent_code', -1)
- # qform = header.get('qform_code', -1)
- # sform = header.get('sform_code', -1)
- #
- # TODO Change this in fslpy 2.0.0
- #
- if isinstance(header, nib.nifti1.Nifti1Header):
- intent = header['intent_code']
- qform = header['qform_code']
- sform = header['sform_code']
- else:
- intent = -1
- qform = -1
- sform = -1
-
- if intent in (constants.FSL_FNIRT_DISPLACEMENT_FIELD,
+ intent = header.get('intent_code', -1)
+ qform = header.get('qform_code', -1)
+ sform = header.get('sform_code', -1)
+
+ # FNIRT non-linear coefficient files
+ # clobber the sform/qform/intent
+ # and pixdims of the nifti header,
+ # so we can't correctly place it in
+ # the world coordinate system. See
+ # $FSLDIR/src/fnirt/fnirt_file_writer.cpp
+ # and fsl.transform.nonlinear for more
+ # details.
+ if intent in (constants.FSL_DCT_COEFFICIENTS,
constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
- constants.FSL_DCT_COEFFICIENTS,
- constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS):
- log.debug('FNIRT output image detected - using qform matrix')
- voxToWorldMat = np.array(header.get_qform())
+ constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS,
+ constants.FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS,
+ constants.FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS):
+
+ log.debug('FNIRT coefficient field detected - generating affine')
+
+ # Knot spacing is stored in the pixdims
+ # (specified in terms of reference image
+ # voxels), and reference image pixdims
+ # are stored as intent code parameters.
+ # If we combine the two, we can at least
+ # get the shape/size of the coefficient
+ # field about right
+ knotpix = header.get_zooms()[:3]
+ refpix = (header.get('intent_p1', 1) or 1,
+ header.get('intent_p2', 1) or 1,
+ header.get('intent_p3', 1) or 1)
+ voxToWorldMat = affine.concat(
+ affine.scaleOffsetXform(refpix, 0),
+ affine.scaleOffsetXform(knotpix, 0))
# If the qform or sform codes are unknown,
# then we can't assume that the transform
@@ -290,7 +356,7 @@ class Nifti(notifier.Notifier, meta.Meta):
# should just be a straight scaling matrix.
elif qform == 0 and sform == 0:
pixdims = header.get_zooms()
- voxToWorldMat = transform.scaleOffsetXform(pixdims, 0)
+ voxToWorldMat = affine.scaleOffsetXform(pixdims, 0)
# Otherwise we let nibabel decide
# which transform to use.
@@ -300,34 +366,102 @@ class Nifti(notifier.Notifier, meta.Meta):
return voxToWorldMat
- def __determineShape(self, header):
- """This method is called by :meth:`__init__`. It figures out the actual
- shape of the image data, and the zooms/pixdims for each data axis. Any
- empty trailing dimensions are squeezed, but the returned shape is
- guaranteed to be at least 3 dimensions. Returns:
+ @staticmethod
+ def generateAffines(voxToWorldMat, shape, pixdim):
+ """Called by :meth:`__init__`, and the :meth:`voxToWorldMat` setter.
+ Generates and returns a dictionary containing affine transformations
+ between the ``voxel``, ``fsl``, and ``world`` coordinate
+ systems. These affines are accessible via the :meth:`getAffine`
+ method.
+
+ :arg voxToWorldMat: The voxel-to-world affine transformation
+ :arg shape: Image shape (number of voxels along each dimension
+ :arg pixdim: Image pixdims (size of one voxel along each
+ dimension)
+ :returns: A tuple containing:
+
+ - a dictionary of affine transformations between
+ each pair of coordinate systems
+ - ``True`` if the image is to be considered
+ "neurological", ``False`` otherwise - see the
+ :meth:`isNeurological` method.
+ """
- - A sequence/tuple containing the image shape, as reported in the
- header.
- - A sequence/tuple containing the effective image shape.
- - A sequence/tuple containing the zooms/pixdims.
+ import numpy.linalg as npla
+
+ affines = {}
+ shape = list(shape[ :3])
+ pixdim = list(pixdim[:3])
+ voxToScaledVoxMat = np.diag(pixdim + [1.0])
+ isneuro = npla.det(voxToWorldMat) > 0
+
+ if isneuro:
+ x = (shape[0] - 1) * pixdim[0]
+ flip = affine.scaleOffsetXform([-1, 1, 1],
+ [ x, 0, 0])
+ voxToScaledVoxMat = affine.concat(flip, voxToScaledVoxMat)
+
+ affines['fsl', 'fsl'] = np.eye(4)
+ affines['voxel', 'voxel'] = np.eye(4)
+ affines['world', 'world'] = np.eye(4)
+ affines['voxel', 'world'] = voxToWorldMat
+ affines['world', 'voxel'] = affine.invert(voxToWorldMat)
+ affines['voxel', 'fsl'] = voxToScaledVoxMat
+ affines['fsl', 'voxel'] = affine.invert(voxToScaledVoxMat)
+ affines['fsl', 'world'] = affine.concat(affines['voxel', 'world'],
+ affines['fsl', 'voxel'])
+ affines['world', 'fsl'] = affine.concat(affines['voxel', 'fsl'],
+ affines['world', 'voxel'])
+
+ return affines, isneuro
+
+
+ @staticmethod
+ def identifyAffine(image, xform, from_=None, to=None):
+ """Attempt to identify the source or destination space for the given
+ affine.
+
+ ``xform`` is assumed to be an affine transformation which can be used
+ to transform coordinates between two coordinate systems associated with
+ ``image``.
+
+ If one of ``from_`` or ``to`` is provided, the other will be derived.
+ If neither are provided, both will be derived. See the
+ :meth:`.Nifti.getAffine` method for details on the valild values that
+ ``from_`` and ``to`` may take.
+
+ :arg image: :class:`.Nifti` instance associated with the affine.
+
+ :arg xform: ``(4, 4)`` ``numpy`` array encoding an affine
+ transformation
+
+ :arg from_: Label specifying the coordinate system which ``xform``
+ takes as input
+
+ :arg to: Label specifying the coordinate system which ``xform``
+ produces as output
+
+ :returns: A tuple containing:
+ - A label for the ``from_`` coordinate system
+ - A label for the ``to`` coordinate system
"""
- # The canonicalShape function figures out
- # the data shape that we should use.
- origShape = list(header.get_data_shape())
- shape = canonicalShape(origShape)
- pixdims = list(header.get_zooms())
+ if (from_ is not None) and (to is not None):
+ return from_, to
- # if get_zooms() doesn't return at
- # least len(shape) values, we'll
- # fall back to the pixdim field in
- # the header.
- if len(pixdims) < len(shape):
- pixdims = header['pixdim'][1:]
+ if from_ is not None: froms = [from_]
+ else: froms = ['voxel', 'fsl', 'world']
+ if to is not None: tos = [to]
+ else: tos = ['voxel', 'fsl', 'world']
- pixdims = pixdims[:len(shape)]
+ for from_, to in it.product(froms, tos):
- return origShape, shape, pixdims
+ candidate = image.getAffine(from_, to)
+
+ if np.all(np.isclose(candidate, xform)):
+ return from_, to
+
+ raise ValueError('Could not identify affine')
def strval(self, key):
@@ -382,6 +516,16 @@ class Nifti(notifier.Notifier, meta.Meta):
return tuple(self.__shape)
+ @property
+ def ndim(self):
+ """Returns the number of dimensions in this image. This number may not
+ match the number of dimensions specified in the NIFTI header, as
+ trailing dimensions of length 1 are ignored. But it is guaranteed to be
+ at least 3.
+ """
+ return len(self.__shape)
+
+
@property
def pixdim(self):
"""Returns a tuple containing the image pixdims (voxel sizes)."""
@@ -390,9 +534,17 @@ class Nifti(notifier.Notifier, meta.Meta):
@property
def intent(self):
- """Returns the NIFTI intent code of this image.
- """
- return self.__intent
+ """Returns the NIFTI intent code of this image. """
+ return self.header.get('intent_code', constants.NIFTI_INTENT_NONE)
+
+
+ @intent.setter
+ def intent(self, val):
+ """Sets the NIFTI intent code of this image. """
+ # analyze has no intent
+ if (self.niftiVersion > 0) and (val != self.intent):
+ self.header.set_intent(val, allow_unknown=True)
+ self.notify(topic='header')
@property
@@ -426,12 +578,39 @@ class Nifti(notifier.Notifier, meta.Meta):
return units
+ def getAffine(self, from_, to):
+ """Return an affine transformation which can be used to transform
+ coordinates from ``from_`` to ``to``.
+
+ Valid values for the ``from_`` and ``to`` arguments are:
+
+ - ``'voxel'``: The voxel coordinate system
+ - ``'world'``: The world coordinate system, as defined by the image
+ sform/qform
+ - ``'fsl'``: The FSL coordinate system (scaled voxels, with a
+ left-right flip if the sform/qform has a positive determinant)
+
+ :arg from_: Source coordinate system
+ :arg to: Destination coordinate system
+ :returns: A ``numpy`` array of shape ``(4, 4)``
+ """
+ from_ = from_.lower()
+ to = to .lower()
+
+ if from_ not in ('voxel', 'fsl', 'world') or \
+ to not in ('voxel', 'fsl', 'world'):
+ raise ValueError('Invalid source/reference spaces: '
+ '{} -> {}'.format(from_, to))
+
+ return np.copy(self.__affines[from_, to])
+
+
@property
def worldToVoxMat(self):
"""Returns a ``numpy`` array of shape ``(4, 4)`` containing an
affine transformation from world coordinates to voxel coordinates.
"""
- return np.array(self.__worldToVoxMat)
+ return self.getAffine('world', 'voxel')
@property
@@ -439,7 +618,7 @@ class Nifti(notifier.Notifier, meta.Meta):
"""Returns a ``numpy`` array of shape ``(4, 4)`` containing an
affine transformation from voxel coordinates to world coordinates.
"""
- return np.array(self.__voxToWorldMat)
+ return self.getAffine('voxel', 'world')
@voxToWorldMat.setter
@@ -463,8 +642,12 @@ class Nifti(notifier.Notifier, meta.Meta):
header.set_sform(xform, code=sformCode)
- self.__voxToWorldMat = self.__determineTransform(header)
- self.__worldToVoxMat = transform.invert(self.__voxToWorldMat)
+ affines, isneuro = Nifti.generateAffines(xform,
+ self.shape,
+ self.pixdim)
+
+ self.__affines = affines
+ self.__isNeurological = isneuro
log.debug('Affine changed:\npixdims: '
'%s\nsform: %s\nqform: %s',
@@ -475,6 +658,27 @@ class Nifti(notifier.Notifier, meta.Meta):
self.notify(topic='transform')
+ @property
+ def voxToScaledVoxMat(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
+ """
+ return self.getAffine('voxel', 'fsl')
+
+
+ @property
+ def scaledVoxToVoxMat(self):
+ """Returns a transformation matrix which transforms from scaled voxels
+ into voxels, the inverse of the :meth:`voxToScaledVoxMat` transform.
+ """
+ return self.getAffine('fsl', 'voxel')
+
+
def mapIndices(self, sliceobj):
"""Adjusts the given slice object so that it may be used to index the
underlying ``nibabel`` NIFTI image object.
@@ -490,16 +694,6 @@ class Nifti(notifier.Notifier, meta.Meta):
return fileslice.canonical_slicers(sliceobj, self.__origShape)
- @property
- def ndim(self):
- """Returns the number of dimensions in this image. This number may not
- match the number of dimensions specified in the NIFTI header, as
- trailing dimensions of length 1 are ignored. But it is guaranteed to be
- at least 3.
- """
- return len(self.__shape)
-
-
def getXFormCode(self, code=None):
"""This method returns the code contained in the NIFTI header,
indicating the space to which the (transformed) image is oriented.
@@ -546,6 +740,7 @@ class Nifti(notifier.Notifier, meta.Meta):
return int(code)
+
# TODO Check what has worse performance - hashing
# a 4x4 array (via memoizeMD5), or the call
# to aff2axcodes. I'm guessing the latter,
@@ -563,7 +758,6 @@ class Nifti(notifier.Notifier, meta.Meta):
return nib.orientations.aff2axcodes(xform, inaxes)
- @memoize.Instanceify(memoize.memoize)
def isNeurological(self):
"""Returns ``True`` if it looks like this ``Nifti`` object has a
neurological voxel orientation, ``False`` otherwise. This test is
@@ -582,51 +776,7 @@ class Nifti(notifier.Notifier, meta.Meta):
_format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to\
_the_transformation_parameters.3F
"""
- import numpy.linalg as npla
- return npla.det(self.__voxToWorldMat) > 0
-
-
- @property
- def voxToScaledVoxMat(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
- """
- return self.__voxToScaledVoxMat()
-
-
- @memoize.Instanceify(memoize.memoize)
- def __voxToScaledVoxMat(self):
- """See :meth:`voxToScaledVoxMat`. """
-
- 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(flip, voxToPixdimMat)
-
- return voxToPixdimMat
-
-
- @property
- def scaledVoxToVoxMat(self):
- """Returns a transformation matrix which transforms from scaled voxels
- into voxels, the inverse of the :meth:`voxToScaledVoxMat` transform.
- """
- return self.__scaledVoxToVoxMat()
-
-
- @memoize.Instanceify(memoize.memoize)
- def __scaledVoxToVoxMat(self):
- """See :meth:`scaledVoxToVoxMat`. """
- return transform.invert(self.voxToScaledVoxMat)
+ return self.__isNeurological
def sameSpace(self, other):
@@ -879,7 +1029,13 @@ class Image(Nifti):
nibImage = ctr(image, xform, header=header)
- # otherwise, we assume that it is a nibabel image
+ # If it's an Image object, we
+ # just take the nibabel image
+ elif isinstance(image, Image):
+ nibImage = image.nibImage
+
+ # otherwise, we assume that
+ # it is a nibabel image
else:
nibImage = image
@@ -914,9 +1070,10 @@ class Image(Nifti):
threaded=threaded)
# Listen to ourself for changes
- # to the voxToWorldMat, so we
+ # to header attributse so we
# can update the saveState.
- self.register(self.name, self.__transformChanged, topic='transform')
+ self.register(self.name, self.__headerChanged, topic='transform')
+ self.register(self.name, self.__headerChanged, topic='header')
if calcRange:
self.calcRange()
@@ -1055,8 +1212,8 @@ class Image(Nifti):
self.__nibImage.set_sform(xform, code)
- def __transformChanged(self, *args, **kwargs):
- """Called when the ``voxToWorldMat`` of this :class:`Nifti` instance
+ def __headerChanged(self, *args, **kwargs):
+ """Called when header properties of this :class:`Nifti` instance
changes. Updates the :attr:`saveState` accordinbgly.
"""
if self.__saveState:
@@ -1137,7 +1294,9 @@ class Image(Nifti):
# We save the image out to a temp file,
# then close the old image, move the
# temp file to the real destination,
- # then re-open the file.
+ # then re-open the file. This is done
+ # to ensure that all references to the
+ # old file are destroyed.
tmphd, tmpfname = tempfile.mkstemp(suffix=getExt(filename))
os.close(tmphd)
diff --git a/fsl/data/mesh.py b/fsl/data/mesh.py
index f3d1fb93bfdeaff3a0fb3afac2cd2dd9a0bb38cf..c15e49ced5f972e39e01af4c520053ecdb6b809a 100644
--- a/fsl/data/mesh.py
+++ b/fsl/data/mesh.py
@@ -33,9 +33,9 @@ import collections
import os.path as op
import numpy as np
-import fsl.utils.meta as meta
-import fsl.utils.notifier as notifier
-import fsl.utils.transform as transform
+import fsl.utils.meta as meta
+import fsl.utils.notifier as notifier
+import fsl.transform.affine as affine
log = logging.getLogger(__name__)
@@ -693,7 +693,7 @@ def calcFaceNormals(vertices, indices):
v2 = vertices[indices[:, 2]]
fnormals = np.cross((v1 - v0), (v2 - v0))
- fnormals = transform.normalise(fnormals)
+ fnormals = affine.normalise(fnormals)
return fnormals
@@ -724,7 +724,7 @@ def calcVertexNormals(vertices, indices, fnormals):
vnormals[v2, :] += fnormals[i]
# normalise to unit length
- return transform.normalise(vnormals)
+ return affine.normalise(vnormals)
def needsFixing(vertices, indices, fnormals, loBounds, hiBounds):
@@ -769,4 +769,4 @@ def needsFixing(vertices, indices, fnormals, loBounds, hiBounds):
# vertex to the camera is positive
# If it isn't, we need to flip the
# triangle winding order.
- return np.dot(n, transform.normalise(camera - vert)) < 0
+ return np.dot(n, affine.normalise(camera - vert)) < 0
diff --git a/fsl/data/mghimage.py b/fsl/data/mghimage.py
index 6c36340d923c8df61ddc8a5ee0ca2d308bde54dc..a143e06ce1e3ce9a2367444c7e67408459e3cd36 100644
--- a/fsl/data/mghimage.py
+++ b/fsl/data/mghimage.py
@@ -14,9 +14,9 @@ import os.path as op
import six
import nibabel as nib
-import fsl.utils.path as fslpath
-import fsl.utils.transform as transform
-import fsl.data.image as fslimage
+import fsl.utils.path as fslpath
+import fsl.transform.affine as affine
+import fsl.data.image as fslimage
ALLOWED_EXTENSIONS = ['.mgz', '.mgh']
@@ -55,12 +55,12 @@ class MGHImage(fslimage.Image):
filename = None
data = image.get_data()
- affine = image.affine
+ xform = image.affine
vox2surf = image.header.get_vox2ras_tkr()
fslimage.Image.__init__(self,
data,
- xform=affine,
+ xform=xform,
name=name,
dataSource=filename)
@@ -68,9 +68,9 @@ class MGHImage(fslimage.Image):
self.setMeta('mghImageFile', filename)
self.__voxToSurfMat = vox2surf
- self.__surfToVoxMat = transform.invert(vox2surf)
- self.__surfToWorldMat = transform.concat(affine, self.__surfToVoxMat)
- self.__worldToSurfMat = transform.invert(self.__surfToWorldMat)
+ self.__surfToVoxMat = affine.invert(vox2surf)
+ self.__surfToWorldMat = affine.concat(xform, self.__surfToVoxMat)
+ self.__worldToSurfMat = affine.invert(self.__surfToWorldMat)
def save(self, filename=None):
diff --git a/fsl/scripts/fsl_apply_x5.py b/fsl/scripts/fsl_apply_x5.py
new file mode 100644
index 0000000000000000000000000000000000000000..0981d960863ca18673501ebe4b7422fb1e68ded1
--- /dev/null
+++ b/fsl/scripts/fsl_apply_x5.py
@@ -0,0 +1,148 @@
+#!/usr/bin/env python
+#
+# fsl_apply_x5.py - Apply an X5 transformation to an image.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""The ``fsl_apply_x5`` script can be used to apply an X5 transformation file
+to resample an image.
+"""
+
+
+import functools as ft
+import sys
+import argparse
+
+import fsl.transform.x5 as x5
+import fsl.transform.nonlinear as nonlinear
+import fsl.utils.parse_data as parse_data
+import fsl.utils.image.resample as resample
+import fsl.data.image as fslimage
+
+
+def parseArgs(args):
+ """Parses command-line arguments.
+
+ :arg args: Sequence of command-line arguments. If ``None``, ``sys.argv``
+ is used
+ :returns: An ``argparse.Namespace`` object containing parsed arguments
+ """
+
+ parser = argparse.ArgumentParser('fsl_apply_x5')
+ flags = {
+ 'input' : ('input',),
+ 'xform' : ('xform',),
+ 'output' : ('output',),
+ 'interp' : ('-i', '--interp'),
+ 'ref' : ('-r', '--ref'),
+ }
+
+ helps = {
+ 'input' : 'Input image',
+ 'xform' : 'X5 transformation file',
+ 'output' : 'Output image',
+ 'interp' : 'Interpolation (default: linear)',
+ 'ref' : 'Alternate reference image (default: '
+ 'reference specified in X5 file)',
+ }
+ opts = {
+ 'input' : dict(help=helps['input'],
+ type=parse_data.Image),
+ 'xform' : dict(help=helps['xform']),
+ 'output' : dict(help=helps['output'],
+ type=parse_data.ImageOut),
+ 'interp' : dict(help=helps['interp'],
+ choices=('nearest', 'linear', 'cubic'),
+ default='linear'),
+ 'ref' : dict(help=helps['ref'],
+ type=ft.partial(parse_data.Image, loadData=False)),
+ }
+
+ parser.add_argument(*flags['input'], **opts['input'])
+ parser.add_argument(*flags['xform'], **opts['xform'])
+ parser.add_argument(*flags['output'], **opts['output'])
+ parser.add_argument(*flags['interp'], **opts['interp'])
+ parser.add_argument(*flags['ref'], **opts['ref'])
+
+ if len(args) == 0:
+ parser.print_help()
+ sys.exit(0)
+
+ args = parser.parse_args(args)
+
+ if args.interp == 'nearest': args.interp = 0
+ elif args.interp == 'linear': args.interp = 1
+ elif args.interp == 'cubic': args.interp = 3
+
+ return args
+
+
+def applyLinear(args):
+ """Applies a linear X5 transformation file to the input.
+
+ :arg args: ``argparse.Namespace`` object
+ :returns: The transformed input as an :class:`.Image` object
+ """
+
+ input = args.input
+ xform, src, ref = x5.readLinearX5(args.xform)
+
+ if args.ref is not None:
+ ref = args.ref
+
+ res, xform = resample.resampleToReference(input,
+ ref,
+ matrix=xform,
+ order=args.interp)
+
+ return fslimage.Image(res, xform=xform, header=ref.header)
+
+
+def applyNonlinear(args):
+ """Applies a non-linear X5 transformation file to the input.
+
+ :arg args: ``argparse.Namespace`` object
+ :returns: The transformed input as an :class:`.Image` object
+ """
+
+ field = x5.readNonLinearX5(args.xform)
+
+ if args.ref is None: ref = field.ref
+ else: ref = args.ref
+
+ result = nonlinear.applyDeformation(args.input,
+ field,
+ ref=ref,
+ order=args.interp,
+ mode='constant')
+
+ return fslimage.Image(result, header=ref.header)
+
+
+def main(args=None):
+ """Entry point. Parse command-line arguments, then calls
+ :func:`applyLinear` or :func:`applyNonlinear` depending on the x5 file
+ type.
+ """
+
+ if args is None:
+ args = sys.argv[1:]
+
+ print()
+ print('Warning: this version of fsl_apply_x5 is a development release.\n'
+ 'Interface, behaviour, and input/output formats of future versions\n'
+ 'may differ substantially from this version.')
+ print()
+
+ args = parseArgs(args)
+
+ if x5.inferType(args.xform) == 'linear':
+ result = applyLinear(args)
+ else:
+ result = applyNonlinear(args)
+
+ result.save(args.output)
+
+
+if __name__ == '__main__':
+ sys.exit(main())
diff --git a/fsl/scripts/fsl_convert_x5.py b/fsl/scripts/fsl_convert_x5.py
new file mode 100644
index 0000000000000000000000000000000000000000..d70b78d670091ba54787e009e8ea4f9810cbb5b4
--- /dev/null
+++ b/fsl/scripts/fsl_convert_x5.py
@@ -0,0 +1,214 @@
+#!/usr/bin/env python
+#
+# fsl_convert_x5.py - Convert between FSL and X5 transformation files.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This script can be used to convert between FSL and X5 linear and non-linear
+transformation file formats.
+"""
+
+
+import os.path as op
+import functools as ft
+import textwrap as tw
+import sys
+import shutil
+import logging
+import argparse
+
+import fsl.data.image as fslimage
+import fsl.utils.parse_data as parse_data
+import fsl.transform.flirt as flirt
+import fsl.transform.fnirt as fnirt
+import fsl.transform.x5 as x5
+
+
+log = logging.getLogger(__name__)
+
+
+def parseArgs(args):
+ """Create an argument parser and parse the given ``args``.
+
+ :arg args: Sequence of commane line arguments.
+ :return: An ``argparse.Namespace`` object containing parsed arguments.
+ """
+
+ helps = {
+ 'input' : 'Input file',
+ 'output' : 'Output file',
+ 'source' : 'Source image',
+ 'reference' : 'Reference image',
+ 'input_format' : 'Input format - if not provided, the input format '
+ 'is automatically inferred from the input file '
+ 'name.',
+ 'output_format' : 'Output format - if not provided, the output format '
+ 'is automatically inferred from the output file '
+ 'name.',
+ }
+
+ epilog = tw.dedent("""
+ If the input and output file formats are the same, the input file
+ is simply copied to the output file.
+ """).strip()
+
+ parser = argparse.ArgumentParser('fsl_convert_x5')
+ subparsers = parser.add_subparsers(dest='ctype')
+ flirt = subparsers.add_parser('flirt', epilog=epilog)
+ fnirt = subparsers.add_parser('fnirt', epilog=epilog)
+ imgtype = ft.partial(parse_data.Image, loadData=False)
+
+ flirt.add_argument('input', help=helps['input'])
+ flirt.add_argument('output', help=helps['output'])
+ flirt.add_argument('-s', '--source', help=helps['source'],
+ type=imgtype)
+ flirt.add_argument('-r', '--reference', help=helps['reference'],
+ type=imgtype)
+ flirt.add_argument('-if', '--input_format', help=helps['input_format'],
+ choices=('x5', 'mat'))
+ flirt.add_argument('-of', '--output_format', help=helps['output_format'],
+ choices=('x5', 'mat'))
+
+ fnirt .add_argument('input', help=helps['input'])
+ fnirt .add_argument('output', help=helps['output'])
+ fnirt .add_argument('-s', '--source', help=helps['source'],
+ type=imgtype)
+ fnirt .add_argument('-r', '--reference', help=helps['reference'],
+ type=imgtype)
+ fnirt .add_argument('-if', '--input_format', help=helps['input_format'],
+ choices=('x5', 'nii'))
+ fnirt .add_argument('-of', '--output_format', help=helps['output_format'],
+ choices=('x5', 'nii'))
+
+ if len(args) == 0:
+ parser.print_help()
+ sys.exit(0)
+
+ if len(args) == 1:
+ if args[0] == 'flirt':
+ flirt.print_help()
+ sys.exit(0)
+ elif args[0] == 'fnirt':
+ fnirt.print_help()
+ sys.exit(0)
+
+ args = parser.parse_args(args)
+
+ # If input/output formats were not
+ # specified, infer them from the
+ # file names
+ def getfmt(arg, fname, exist):
+ ext = op.splitext(fname)[1]
+ if ext in ('.mat', '.x5'):
+ return ext[1:]
+
+ fname = fslimage.addExt(fname, mustExist=exist)
+
+ if fslimage.looksLikeImage(fname):
+ return 'nii'
+ parser.error('Could not infer format from '
+ 'filename: {}'.format(args.input))
+
+ if args.input_format is None:
+ args.input_format = getfmt('input', args.input, True)
+ if args.output_format is None:
+ args.output_format = getfmt('output', args.output, False)
+
+ # The source and reference arguments are
+ # required if the input is a FLIRT matrix
+ # or a FNIRT displacement/coefficient field.
+ if args.input_format in ('mat', 'nii') and \
+ (args.source is None or args.reference is None):
+ parser.error('You must specify a source and reference '
+ 'when the input is not an X5 file!')
+
+ return args
+
+
+def flirtToX5(args):
+ """Convert a linear FLIRT transformation matrix to an X5 transformation
+ file.
+ """
+ src = args.source
+ ref = args.reference
+ xform = flirt.readFlirt(args.input)
+ xform = flirt.fromFlirt(xform, src, ref, 'world', 'world')
+ x5.writeLinearX5(args.output, xform, src, ref)
+
+
+def X5ToFlirt(args):
+ """Convert a linear X5 transformation file to a FLIRT matrix. """
+ xform, src, ref = x5.readLinearX5(args.input)
+ xform = flirt.toFlirt(xform, src, ref, 'world', 'world')
+ flirt.writeFlirt(xform, args.output)
+
+
+def fnirtToX5(args):
+ """Convert a non-linear FNIRT transformation into an X5 transformation
+ file.
+ """
+ src = args.source
+ ref = args.reference
+ field = fnirt.readFnirt(args.input, src=src, ref=ref)
+ field = fnirt.fromFnirt(field, 'world', 'world')
+ x5.writeNonLinearX5(args.output, field)
+
+
+def X5ToFnirt(args):
+ """Convert a non-linear X5 transformation file to a FNIRT-style
+ transformation file.
+ """
+ field = x5.readNonLinearX5(args.input)
+ field = fnirt.toFnirt(field)
+ field.save(args.output)
+
+
+def doFlirt(args):
+ """Converts a linear FIRT transformation file to or from the X5 format. """
+ infmt = args.input_format
+ outfmt = args.output_format
+
+ if (infmt, outfmt) == ('x5', 'mat'): X5ToFlirt(args)
+ elif (infmt, outfmt) == ('mat', 'x5'): flirtToX5(args)
+ else: shutil.copy(args.input, args.output)
+
+
+def doFnirt(args):
+ """Converts a non-linear FNIRT transformation file to or from the X5
+ format.
+ """
+ infmt = args.input_format
+ outfmt = args.output_format
+
+ if (infmt, outfmt) == ('x5', 'nii'): X5ToFnirt(args)
+ elif (infmt, outfmt) == ('nii', 'x5'): fnirtToX5(args)
+ else: shutil.copy(args.input, args.output)
+
+
+def main(args=None):
+ """Entry point. Calls :func:`doFlirt` or :func:`doFnirt` depending on
+ the sub-command specified in the arguments.
+
+ :arg args: Sequence of command-line arguments. If not provided,
+ ``sys.argv`` is used.
+ """
+
+ if args is None:
+ args = sys.argv[1:]
+
+ print()
+ print('Warning: this version of fsl_convert_x5 is a development release.\n'
+ 'Interface, behaviour, and input/output formats of future versions\n'
+ 'may differ substantially from this version.')
+ print()
+
+ args = parseArgs(args)
+
+ if args.ctype == 'flirt': doFlirt(args)
+ elif args.ctype == 'fnirt': doFnirt(args)
+
+ return 0
+
+
+if __name__ == '__main__':
+ sys.exit(main())
diff --git a/fsl/transform/__init__.py b/fsl/transform/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..31bad1d68dcd5676ffa319f42a906ce92170fca4
--- /dev/null
+++ b/fsl/transform/__init__.py
@@ -0,0 +1,19 @@
+#!/usr/bin/env python
+#
+# __init__.py - Functions for working with linear and non-linear FSL
+# transformations.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains functions for working with linear and non-linear FSL
+transformations. Functionality is split across the following modules:
+
+.. autosummary::
+ :nosignatures:
+
+ ~fsl.transform.affine
+ ~fsl.transform.flirt
+ ~fsl.transform.fnirt
+ ~fsl.transform.nonlinear
+ ~fsl.transform.x5
+"""
diff --git a/fsl/transform/affine.py b/fsl/transform/affine.py
new file mode 100644
index 0000000000000000000000000000000000000000..e0364bf94030c76d6df99a5222ce6439bed8df38
--- /dev/null
+++ b/fsl/transform/affine.py
@@ -0,0 +1,584 @@
+#!/usr/bin/env python
+#
+# affine.py - Utility functions for working with affine transformations.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains utility functions for working with affine
+transformations. The following functions are available:
+
+.. autosummary::
+ :nosignatures:
+
+ transform
+ scaleOffsetXform
+ invert
+ concat
+ compose
+ decompose
+ rotMatToAffine
+ rotMatToAxisAngles
+ axisAnglesToRotMat
+ axisBounds
+ rmsdev
+
+And a few more functions are provided for working with vectors:
+
+.. autosummary::
+ :nosignatures:
+
+ veclength
+ normalise
+ transformNormal
+"""
+
+
+import collections.abc as abc
+import numpy as np
+import numpy.linalg as linalg
+
+
+def invert(x):
+ """Inverts the given matrix using ``numpy.linalg.inv``. """
+ return linalg.inv(x)
+
+
+def concat(*xforms):
+ """Combines the given matrices (returns the dot product)."""
+
+ result = xforms[0]
+
+ for i in range(1, len(xforms)):
+ result = np.dot(result, xforms[i])
+
+ return result
+
+
+def veclength(vec):
+ """Returns the length of the given vector(s).
+
+ Multiple vectors may be passed in, with a shape of ``(n, 3)``.
+ """
+ vec = np.array(vec, copy=False).reshape(-1, 3)
+ return np.sqrt(np.einsum('ij,ij->i', vec, vec))
+
+
+def normalise(vec):
+ """Normalises the given vector(s) to unit length.
+
+ Multiple vectors may be passed in, with a shape of ``(n, 3)``.
+ """
+ vec = np.array(vec, copy=False).reshape(-1, 3)
+ n = (vec.T / veclength(vec)).T
+
+ if n.size == 3:
+ n = n[0]
+
+ return n
+
+
+def scaleOffsetXform(scales, offsets):
+ """Creates and returns an affine transformation matrix which encodes
+ the specified scale(s) and offset(s).
+
+
+ :arg scales: A tuple of up to three values specifying the scale factors
+ for each dimension. If less than length 3, is padded with
+ ``1.0``.
+
+ :arg offsets: A tuple of up to three values specifying the offsets for
+ each dimension. If less than length 3, is padded with
+ ``0.0``.
+
+ :returns: A ``numpy.float32`` array of size :math:`4 \\times 4`.
+ """
+
+ oktypes = (abc.Sequence, np.ndarray)
+
+ if not isinstance(scales, oktypes): scales = [scales]
+ if not isinstance(offsets, oktypes): offsets = [offsets]
+ if not isinstance(scales, list): scales = list(scales)
+ if not isinstance(offsets, list): offsets = list(offsets)
+
+ lens = len(scales)
+ leno = len(offsets)
+
+ if lens < 3: scales = scales + [1.0] * (3 - lens)
+ if leno < 3: offsets = offsets + [0.0] * (3 - leno)
+
+ xform = np.eye(4, dtype=np.float64)
+
+ xform[0, 0] = scales[0]
+ xform[1, 1] = scales[1]
+ xform[2, 2] = scales[2]
+
+ xform[0, 3] = offsets[0]
+ xform[1, 3] = offsets[1]
+ xform[2, 3] = offsets[2]
+
+ return xform
+
+
+def compose(scales, offsets, rotations, origin=None):
+ """Compose a transformation matrix out of the given scales, offsets
+ and axis rotations.
+
+ :arg scales: Sequence of three scale values.
+
+ :arg offsets: Sequence of three offset values.
+
+ :arg rotations: Sequence of three rotation values, in radians, or
+ a rotation matrix of shape ``(3, 3)``.
+
+ :arg origin: Origin of rotation - must be scaled by the ``scales``.
+ If not provided, the rotation origin is ``(0, 0, 0)``.
+ """
+
+ preRotate = np.eye(4)
+ postRotate = np.eye(4)
+
+ rotations = np.array(rotations)
+
+ if rotations.shape == (3,):
+ rotations = axisAnglesToRotMat(*rotations)
+
+ if origin is not None:
+ preRotate[ 0, 3] = -origin[0]
+ preRotate[ 1, 3] = -origin[1]
+ preRotate[ 2, 3] = -origin[2]
+ postRotate[0, 3] = origin[0]
+ postRotate[1, 3] = origin[1]
+ postRotate[2, 3] = origin[2]
+
+ scale = np.eye(4, dtype=np.float64)
+ offset = np.eye(4, dtype=np.float64)
+ rotate = np.eye(4, dtype=np.float64)
+
+ scale[ 0, 0] = scales[ 0]
+ scale[ 1, 1] = scales[ 1]
+ scale[ 2, 2] = scales[ 2]
+ offset[ 0, 3] = offsets[0]
+ offset[ 1, 3] = offsets[1]
+ offset[ 2, 3] = offsets[2]
+
+ rotate[:3, :3] = rotations
+
+ return concat(offset, postRotate, rotate, preRotate, scale)
+
+
+def decompose(xform, angles=True):
+ """Decomposes the given transformation matrix into separate offsets,
+ scales, and rotations, according to the algorithm described in:
+
+ Spencer W. Thomas, Decomposing a matrix into simple transformations, pp
+ 320-323 in *Graphics Gems II*, James Arvo (editor), Academic Press, 1991,
+ ISBN: 0120644819.
+
+ It is assumed that the given transform has no perspective components. Any
+ shears in the affine are discarded.
+
+ :arg xform: A ``(3, 3)`` or ``(4, 4)`` affine transformation matrix.
+
+ :arg angles: If ``True`` (the default), the rotations are returned
+ as axis-angles, in radians. Otherwise, the rotation matrix
+ is returned.
+
+ :returns: The following:
+
+ - A sequence of three scales
+ - A sequence of three translations (all ``0`` if ``xform``
+ was a ``(3, 3)`` matrix)
+ - A sequence of three rotations, in radians. Or, if
+ ``angles is False``, a rotation matrix.
+ """
+
+ # The inline comments in the code below are taken verbatim from
+ # the referenced article, [except for notes in square brackets].
+
+ # The next step is to extract the translations. This is trivial;
+ # we find t_x = M_{4,1}, t_y = M_{4,2}, and t_z = M_{4,3}. At this
+ # point we are left with a 3*3 matrix M' = M_{1..3,1..3}.
+ xform = xform.T
+
+ if xform.shape == (4, 4):
+ translations = xform[ 3, :3]
+ xform = xform[:3, :3]
+ else:
+ translations = np.array([0, 0, 0])
+
+ M1 = xform[0]
+ M2 = xform[1]
+ M3 = xform[2]
+
+ # The process of finding the scaling factors and shear parameters
+ # is interleaved. First, find s_x = |M'_1|.
+ sx = np.sqrt(np.dot(M1, M1))
+
+ # Then, compute an initial value for the xy shear factor,
+ # s_xy = M'_1 * M'_2. (this is too large by the y scaling factor).
+ sxy = np.dot(M1, M2)
+
+ # The second row of the matrix is made orthogonal to the first by
+ # setting M'_2 = M'_2 - s_xy * M'_1.
+ M2 = M2 - sxy * M1
+
+ # Then the y scaling factor, s_y, is the length of the modified
+ # second row.
+ sy = np.sqrt(np.dot(M2, M2))
+
+ # The second row is normalized, and s_xy is divided by s_y to
+ # get its final value.
+ M2 = M2 / sy
+ sxy = sxy / sy
+
+ # The xz and yz shear factors are computed as in the preceding,
+ sxz = np.dot(M1, M3)
+ syz = np.dot(M2, M3)
+
+ # the third row is made orthogonal to the first two rows,
+ M3 = M3 - sxz * M1 - syz * M2
+
+ # the z scaling factor is computed,
+ sz = np.sqrt(np.dot(M3, M3))
+
+ # the third row is normalized, and the xz and yz shear factors are
+ # rescaled.
+ M3 = M3 / sz
+ sxz = sxz / sz
+ syz = sxz / sz
+
+ # The resulting matrix now is a pure rotation matrix, except that it
+ # might still include a scale factor of -1. If the determinant of the
+ # matrix is -1, negate the matrix and all three scaling factors. Call
+ # the resulting matrix R.
+ #
+ # [We do things different here - if the rotation matrix has negative
+ # determinant, the flip is encoded in the x scaling factor.]
+ R = np.array([M1, M2, M3])
+ if linalg.det(R) < 0:
+ R[0] = -R[0]
+ sx = -sx
+
+ # Finally, we need to decompose the rotation matrix into a sequence
+ # of rotations about the x, y, and z axes. [This is done in the
+ # rotMatToAxisAngles function]
+ if angles: rotations = rotMatToAxisAngles(R.T)
+ else: rotations = R.T
+
+ return np.array([sx, sy, sz]), translations, rotations
+
+
+def rotMatToAffine(rotmat, origin=None):
+ """Convenience function which encodes the given ``(3, 3)`` rotation
+ matrix into a ``(4, 4)`` affine.
+ """
+ return compose([1, 1, 1], [0, 0, 0], rotmat, origin)
+
+
+def rotMatToAxisAngles(rotmat):
+ """Given a ``(3, 3)`` rotation matrix, decomposes the rotations into
+ an angle in radians about each axis.
+ """
+
+ yrot = np.sqrt(rotmat[0, 0] ** 2 + rotmat[1, 0] ** 2)
+
+ if np.isclose(yrot, 0):
+ xrot = np.arctan2(-rotmat[1, 2], rotmat[1, 1])
+ yrot = np.arctan2(-rotmat[2, 0], yrot)
+ zrot = 0
+ else:
+ xrot = np.arctan2( rotmat[2, 1], rotmat[2, 2])
+ yrot = np.arctan2(-rotmat[2, 0], yrot)
+ zrot = np.arctan2( rotmat[1, 0], rotmat[0, 0])
+
+ return [xrot, yrot, zrot]
+
+
+def axisAnglesToRotMat(xrot, yrot, zrot):
+ """Constructs a ``(3, 3)`` rotation matrix from the given angles, which
+ must be specified in radians.
+ """
+
+ xmat = np.eye(3)
+ ymat = np.eye(3)
+ zmat = np.eye(3)
+
+ xmat[1, 1] = np.cos(xrot)
+ xmat[1, 2] = -np.sin(xrot)
+ xmat[2, 1] = np.sin(xrot)
+ xmat[2, 2] = np.cos(xrot)
+
+ ymat[0, 0] = np.cos(yrot)
+ ymat[0, 2] = np.sin(yrot)
+ ymat[2, 0] = -np.sin(yrot)
+ ymat[2, 2] = np.cos(yrot)
+
+ zmat[0, 0] = np.cos(zrot)
+ zmat[0, 1] = -np.sin(zrot)
+ zmat[1, 0] = np.sin(zrot)
+ zmat[1, 1] = np.cos(zrot)
+
+ return concat(zmat, ymat, xmat)
+
+
+def axisBounds(shape,
+ xform,
+ axes=None,
+ origin='centre',
+ boundary='high',
+ offset=1e-4):
+ """Returns the ``(lo, hi)`` bounds of the specified axis/axes in the
+ world coordinate system defined by ``xform``.
+
+ If the ``origin`` parameter is set to ``centre`` (the default),
+ this function assumes that voxel indices correspond to the voxel
+ centre. For example, the voxel at ``(4, 5, 6)`` covers the space:
+
+ ``[3.5 - 4.5, 4.5 - 5.5, 5.5 - 6.5]``
+
+ So the bounds of the specified shape extends from the corner at
+
+ ``(-0.5, -0.5, -0.5)``
+
+ to the corner at
+
+ ``(shape[0] - 0.5, shape[1] - 0.5, shape[1] - 0.5)``
+
+ If the ``origin`` parameter is set to ``corner``, this function
+ assumes that voxel indices correspond to the voxel corner. In this
+ case, a voxel at ``(4, 5, 6)`` covers the space:
+
+ ``[4 - 5, 5 - 6, 6 - 7]``
+
+ So the bounds of the specified shape extends from the corner at
+
+ ``(0, 0, 0)``
+
+ to the corner at
+
+ ``(shape[0], shape[1], shape[1])``.
+
+
+ If the ``boundary`` parameter is set to ``high``, the high voxel bounds
+ are reduced by a small amount (specified by the ``offset`` parameter)
+ before they are transformed to the world coordinate system. If
+ ``boundary`` is set to ``low``, the low bounds are increased by a small
+ amount. The ``boundary`` parameter can also be set to ``'both'``, or
+ ``None``. This option is provided so that you can ensure that the
+ resulting bounds will always be contained within the image space.
+
+ :arg shape: The ``(x, y, z)`` shape of the data.
+
+ :arg xform: Transformation matrix which transforms voxel coordinates
+ to the world coordinate system.
+
+ :arg axes: The world coordinate system axis bounds to calculate.
+
+ :arg origin: Either ``'centre'`` (the default) or ``'corner'``.
+
+ :arg boundary: Either ``'high'`` (the default), ``'low'``, ''`both'``,
+ or ``None``.
+
+ :arg offset: Amount by which the boundary voxel coordinates should be
+ offset. Defaults to ``1e-4``.
+
+ :returns: A tuple containing the ``(low, high)`` bounds for each
+ requested world coordinate system axis.
+ """
+
+ origin = origin.lower()
+
+ # lousy US spelling
+ if origin == 'center':
+ origin = 'centre'
+
+ if origin not in ('centre', 'corner'):
+ raise ValueError('Invalid origin value: {}'.format(origin))
+ if boundary not in ('low', 'high', 'both', None):
+ raise ValueError('Invalid boundary value: {}'.format(boundary))
+
+ scalar = False
+
+ if axes is None:
+ axes = [0, 1, 2]
+
+ elif not isinstance(axes, abc.Iterable):
+ scalar = True
+ axes = [axes]
+
+ x, y, z = shape[:3]
+
+ points = np.zeros((8, 3), dtype=np.float32)
+
+ if origin == 'centre':
+ x0 = -0.5
+ y0 = -0.5
+ z0 = -0.5
+ x -= 0.5
+ y -= 0.5
+ z -= 0.5
+ else:
+ x0 = 0
+ y0 = 0
+ z0 = 0
+
+ if boundary in ('low', 'both'):
+ x0 += offset
+ y0 += offset
+ z0 += offset
+
+ if boundary in ('high', 'both'):
+ x -= offset
+ y -= offset
+ z -= offset
+
+ points[0, :] = [x0, y0, z0]
+ points[1, :] = [x0, y0, z]
+ points[2, :] = [x0, y, z0]
+ points[3, :] = [x0, y, z]
+ points[4, :] = [x, y0, z0]
+ points[5, :] = [x, y0, z]
+ points[6, :] = [x, y, z0]
+ points[7, :] = [x, y, z]
+
+ tx = transform(points, xform)
+
+ lo = tx[:, axes].min(axis=0)
+ hi = tx[:, axes].max(axis=0)
+
+ if scalar: return (lo[0], hi[0])
+ else: return (lo, hi)
+
+
+def transform(p, xform, axes=None, vector=False):
+ """Transforms the given set of points ``p`` according to the given affine
+ transformation ``xform``.
+
+
+ :arg p: A sequence or array of points of shape :math:`N \\times 3`.
+
+ :arg xform: A ``(4, 4)`` affine transformation matrix with which to
+ transform the points in ``p``.
+
+ :arg axes: If you are only interested in one or two axes, and the source
+ axes are orthogonal to the target axes (see the note below),
+ you may pass in a 1D, ``N*1``, or ``N*2`` array as ``p``, and
+ use this argument to specify which axis/axes that the data in
+ ``p`` correspond to.
+
+ :arg vector: Defaults to ``False``. If ``True``, the points are treated
+ as vectors - the translation component of the transformation
+ is not applied. If you set this flag, you pass in a ``(3, 3)``
+ transformation matrix.
+
+ :returns: The points in ``p``, transformed by ``xform``, as a ``numpy``
+ array with the same data type as the input.
+
+
+ .. note:: The ``axes`` argument should only be used if the source
+ coordinate system (the points in ``p``) axes are orthogonal
+ to the target coordinate system (defined by the ``xform``).
+
+ In other words, you can only use the ``axes`` argument if
+ the ``xform`` matrix consists solely of translations and
+ scalings.
+ """
+
+ p = _fillPoints(p, axes)
+ t = np.dot(xform[:3, :3], p.T).T
+
+ if not vector:
+ t = t + xform[:3, 3]
+
+ if axes is not None:
+ t = t[:, axes]
+
+ if t.size == 1: return t[0]
+ else: return t
+
+
+def transformNormal(p, xform, axes=None):
+ """Transforms the given point(s), under the assumption that they
+ are normal vectors. In this case, the points are transformed by
+ ``invert(xform[:3, :3]).T``.
+ """
+ return transform(p, invert(xform[:3, :3]).T, axes, vector=True)
+
+
+def _fillPoints(p, axes):
+ """Used by the :func:`transform` function. Turns the given array p into
+ a ``N*3`` array of ``x,y,z`` coordinates. The array p may be a 1D array,
+ or an ``N*2`` or ``N*3`` array.
+ """
+
+ if not isinstance(p, abc.Iterable): p = [p]
+
+ p = np.array(p)
+
+ if axes is None: return p
+
+ if not isinstance(axes, abc.Iterable): axes = [axes]
+
+ if p.ndim == 1:
+ p = p.reshape((len(p), 1))
+
+ if p.ndim != 2:
+ raise ValueError('Points array must be either one or two '
+ 'dimensions')
+
+ if len(axes) != p.shape[1]:
+ raise ValueError('Points array shape does not match specified '
+ 'number of axes')
+
+ newp = np.zeros((len(p), 3), dtype=p.dtype)
+
+ for i, ax in enumerate(axes):
+ newp[:, ax] = p[:, i]
+
+ return newp
+
+
+def rmsdev(T1, T2, R=None, xc=None):
+ """Calculates the RMS deviation of the given affine transforms ``T1`` and
+ ``T2``. This can be used as a measure of the 'distance' between two
+ affines.
+
+ The ``T1`` and ``T2`` arguments may be either full ``(4, 4)`` affines, or
+ ``(3, 3)`` rotation matrices.
+
+ See FMRIB technical report TR99MJ1, available at:
+
+ https://www.fmrib.ox.ac.uk/datasets/techrep/
+
+ :arg T1: First affine
+ :arg T2: Second affine
+ :arg R: Sphere radius
+ :arg xc: Sphere centre
+ :returns: The RMS deviation between ``T1`` and ``T2``.
+ """
+
+ if R is None:
+ R = 1
+
+ if xc is None:
+ xc = np.zeros(3)
+
+ # rotations only
+ if T1.shape == (3, 3):
+ M = np.dot(T2, invert(T1)) - np.eye(3)
+ A = M[:3, :3]
+ t = np.zeros(3)
+
+ # full affine
+ else:
+ M = np.dot(T2, invert(T1)) - np.eye(4)
+ A = M[:3, :3]
+ t = M[:3, 3]
+
+ Axc = np.dot(A, xc)
+
+ erms = np.dot((t + Axc).T, t + Axc)
+ erms = 0.2 * R ** 2 * np.dot(A.T, A).trace() + erms
+ erms = np.sqrt(erms)
+
+ return erms
diff --git a/fsl/transform/flirt.py b/fsl/transform/flirt.py
new file mode 100644
index 0000000000000000000000000000000000000000..e1c84b7617e23eda40aac4692fe586d4942de372
--- /dev/null
+++ b/fsl/transform/flirt.py
@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+#
+# flirt.py - Functions for working with FLIRT matrices.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains functions for working with FLIRT affine transformation
+matrices. The following functions are available:
+
+.. autosummary::
+ :nosignatures:
+
+ readFlirt
+ writeFlirt
+ fromFlirt
+ toFlirt
+ flirtMatrixToSform
+ sformToFlirtMatrix
+
+
+
+FLIRT transformation matrices are affine matrices of shape ``(4, 4)`` which
+encode a linear transformation from a source image to a reference image. FLIRT
+matrices are defined in terms of *FSL coordinates*, which is a coordinate
+system where voxels are scaled by pixdims, and with a left-right flip if the
+image ``sform`` has a positive determinant.
+
+
+FLIRT matrices thus encode a transformation from source image FSL coordinates
+to reference image FSL coordinates.
+"""
+
+
+import numpy as np
+
+from .affine import concat
+
+
+def readFlirt(fname):
+ """Reads a FLIRT matrix from a file. """
+ return np.loadtxt(fname)
+
+
+def writeFlirt(xform, fname):
+ """Writes the given FLIRT matrix to a file. """
+ np.savetxt(fname, xform, fmt='%1.15g')
+
+
+def fromFlirt(xform, src, ref, from_='voxel', to='world'):
+ """Convert a FLIRT affine matrix into another affine.
+
+ Given a FLIRT matrix, and the source and reference images with which the
+ FLIRT matrix is associated, converts the matrix into an affine matrix
+ which can transform coordinates from the source image ``from_`` coordinate
+ system to the reference image ``to`` coordinate system.
+
+ The ``from_`` and ``to`` arguments specify the desired spaces that the
+ returned affine should transform between. The default values
+ (``from_='voxel'`` and ``to='world'`` will generate an affine which
+ transforms from voxels in the source image to world-coordinates in the
+ reference image.
+
+ Valid values for the ``from_`` and ``to`` arguments are:
+
+ - ``voxel``: The voxel coordinate system
+ - ``fsl``: The FSL coordiante system (voxels scaled by pixdims, with the
+ X axis inverted if the image sform/qform has a positive determinant)
+ - ``world`` The world coordinate system
+
+ See the :class:`.Nifti` class documentation and the
+ :meth:`.Nifti.getAffine` method for more details.
+
+ :arg xform: ``numpy`` array of shape ``(4, 4)`` containing a FLIRT
+ transformation matrix.
+ :arg src: :class:`.Nifti` object, the ``xform`` source image
+ :arg ref: :class:`.Nifti` object, the ``xform`` reference image
+ :arg from_: Desired source coordinate system
+ :arg to: Desired reference coordinate system
+ :returns: ``numpy`` array of shape ``(4, 4)`` containing a matrix
+ encoding a transformation from the source ``from_`` to
+ the reference ``to`` coordinate systems.
+ """
+ premat = src.getAffine(from_, 'fsl')
+ postmat = ref.getAffine('fsl', to)
+ return concat(postmat, xform, premat)
+
+
+def toFlirt(xform, src, ref, from_='voxel', to='world'):
+ """Convert an affine matrix into a FLIRT matrix.
+
+ :returns: ``numpy`` array of shape ``(4, 4)`` containing a matrix
+ encoding a transformation from the source ``from_`` to
+ the reference ``to`` coordinate systems.
+ :arg src: :class:`.Nifti` object, the ``xform`` source image
+ :arg ref: :class:`.Nifti` object, the ``xform`` reference image
+ :arg from_: ``xform`` source coordinate system
+ :arg to: ``xform`` target coordinate system
+ :returns: A ``numpy`` array of shape ``(4, 4)`` containing a FLIRT
+ transformation matrix from ``src`` to ``ref``.
+ """
+ premat = src.getAffine('fsl', from_)
+ postmat = ref.getAffine(to, 'fsl')
+ return concat(postmat, xform, premat)
+
+
+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.
+
+ 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`
+ """
+ return fromFlirt(flirtMat, srcImage, refImage, 'voxel', 'world')
+
+
+def sformToFlirtMatrix(srcImage, refImage, srcXform=None):
+ """Under the assumption that the given ``srcImage`` and ``refImage`` share a
+ common world coordinate system (defined by their
+ :attr:`.Nifti.voxToWorldMat` attributes), this function will calculate and
+ return a transformation matrix from the ``srcImage`` FSL coordinate system
+ to the ``refImage`` FSL coordinate system, that can be saved to disk and
+ used with FLIRT, to resample the source image to the space of the
+ reference image.
+
+ :arg srcImage: Source :class:`.Image`
+ :arg refImage: Reference :class:`.Image`
+ :arg srcXform: Optionally used in place of the ``srcImage``
+ :attr:`.Nifti.voxToWorldMat`
+ """
+
+ srcScaledVoxToVoxMat = srcImage.scaledVoxToVoxMat
+ srcVoxToWorldMat = srcImage.voxToWorldMat
+ refWorldToVoxMat = refImage.worldToVoxMat
+ refVoxToScaledVoxMat = refImage.voxToScaledVoxMat
+
+ if srcXform is not None:
+ srcVoxToWorldMat = srcXform
+
+ return concat(refVoxToScaledVoxMat,
+ refWorldToVoxMat,
+ srcVoxToWorldMat,
+ srcScaledVoxToVoxMat)
diff --git a/fsl/transform/fnirt.py b/fsl/transform/fnirt.py
new file mode 100644
index 0000000000000000000000000000000000000000..2a2351a0ed8f7cd1aa5ff2d48b62626361845964
--- /dev/null
+++ b/fsl/transform/fnirt.py
@@ -0,0 +1,388 @@
+#!/usr/bin/env python
+#
+# fnirt.py - Functions for working with FNIRT non-linear transformations.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains functions for working with FNIRT non-linear
+transformations. The following functions are available:
+
+.. autosummary::
+ :nosignatures:
+
+ readFnirt
+ toFnirt
+ fromFnirt
+
+
+Non-linear registration using FNIRT
+-----------------------------------
+
+
+FNIRT is used to calculate a non-linear registration from a source image to a
+reference image. FNIRT outputs the resulting non-linear transformation as
+either:
+
+ - A deformation/warp field which contains displacements or coordinates.
+ - A coefficient field which can be used to generate a warp field.
+
+
+Non-linear registration using FNIRT generally follows the process depicted
+here:
+
+
+.. image:: images/nonlinear_registration_process.png
+ :width: 80%
+ :align: center
+
+
+First, an initial linear registration is performed from the source image to
+the reference image using FLIRT; this provides an initial global alignment
+which can be used as the starting point for the non-linear registration. Next,
+FNIRT is used to non-linearly register the aligned source image to the
+reference image. Importantly, both of these steps are performed using FSL
+coordinates.
+
+
+Note that we have three spaces, and three sets of coordinate systems, to
+consider:
+
+ 1. Source image space - the source image, before initial linear registration
+ to the reference image
+
+ 2. Aligned-source image space - the source image, after it has been linearly
+ transformed to the reference image space
+
+ 3. Reference image space
+
+
+The initial affine registration calculates a transformation between spaces 1
+and 2, and the non-linear registration calculates a transformation between
+spaces 2 and 3. Note that the fields-of-view for spaces 2 and 3 are
+equivalent.
+
+
+The non-linear transformation file generated by FNIRT will contain the initial
+linear registration, with it either being encoded directly into the warps (for
+a warp field), or being stored in the NIfTI header (for a coefficient field).
+
+
+FNIRT warp fields
+^^^^^^^^^^^^^^^^^
+
+
+A FNIRT deformation field (a.k.a. warp field) is defined in the same space as
+the reference image, and may contain:
+
+ - *relative displacements*, where each voxel in the warp field contains an
+ offset which can be added to the reference image coordinates for that
+ voxel, in order to calculate the corresponding source image coordinates.
+
+ - *absolute coordinates*, where each voxel in the warp field simply contains
+ the source image coordinates which correspond to those reference image
+ coordinates.
+
+
+If an initial linear registration was used as the starting point for FNIRT,
+this is encoded into the displacements/coordinates themselves, so they can be
+used to transform from the reference image to the *original* source image
+space.
+
+
+.. image:: images/fnirt_warp_field.png
+ :width: 80%
+ :align: center
+
+
+FNIRT coefficient fields
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+
+A FNIRT coefficient field contains the coefficients of a set of quadratic or
+cubic B-spline functions defined on a regular 3D grid overlaid on the
+reference image voxel coordinate system. Each coefficient in this grid may be
+referred to as a *control point* or a *knot*.
+
+
+Evaluating the spline functions at a particular location in the grid will
+result in a relative displacement which can be added to that location's
+reference image coordinates, in order to determine the corresponding source
+image coordinates.
+
+
+If an initial linear registration was used as the starting point for FNIRT,
+the generated displacement field will encode a transformation to *aligned*
+source image coordinates, and the initial affine will be stored in the NIfTI
+header of the coefficient field file.
+
+
+.. image:: images/fnirt_coefficient_field.png
+ :width: 80%
+ :align: center
+"""
+
+
+import logging
+
+import nibabel as nib
+import numpy as np
+
+import fsl.data.constants as constants
+import fsl.data.image as fslimage
+from . import affine
+from . import nonlinear
+
+
+log = logging.getLogger(__name__)
+
+
+def _readFnirtDeformationField(fname, img, src, ref, defType=None):
+ """Loads ``fname``, assumed to be a FNIRT deformation field.
+
+ :arg fname: File name of FNIRT deformation field
+
+ :arg img: ``fname`` loaded as an :class:`.Image`
+
+ :arg src: Source image
+
+ :arg ref: Reference image
+
+ :arg defType: Deformation type - either ``'absolute'`` or ``'relative'``.
+ If not provided, is automatically inferred from the data.
+
+ :return: A :class:`.DeformationField` object
+ """
+ return nonlinear.DeformationField(fname,
+ src,
+ ref,
+ srcSpace='fsl',
+ refSpace='fsl',
+ defType=defType)
+
+
+def _readFnirtCoefficientField(fname, img, src, ref):
+ """Loads ``fname``, assumed to be a FNIRT coefficient field.
+
+ :arg fname: File name of FNIRT coefficient field
+ :arg img: ``fname`` loaded as an :class:`.Image`
+ :arg src: Source image
+ :arg ref: Reference image
+ :return: A :class:`.CoefficientField`
+ """
+
+ # FNIRT uses NIFTI header fields in
+ # non-standard ways to store some
+ # additional information about the
+ # coefficient field. See
+ # $FSLDIR/src/fnirt/fnirt_file_writer.cpp
+ # for more details.
+
+ # The field type (quadratic, cubic,
+ # or discrete-cosine-transform) is
+ # inferred from the intent. There is
+ # no support in this implementation
+ # for DCT fields
+ cubics = (constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
+ constants.FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS)
+ quads = (constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS,
+ constants.FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS)
+
+ if img.intent in cubics:
+ fieldType = 'cubic'
+ elif img.intent in quads:
+ fieldType = 'quadratic'
+ else:
+ fieldType = 'cubic'
+ log.warning('Unrecognised/unsupported coefficient '
+ 'field type (assuming cubic b-spline): '
+ '{}'.format(img.intent))
+
+ # Knot spacing (in voxels) is
+ # stored in the pixdims
+ knotSpacing = img.pixdim[:3]
+
+ # The sform contains an initial
+ # global src-to-ref affine
+ # (the starting point for the
+ # non-linear registration). This
+ # is encoded as a flirt matrix,
+ # i.e. it transforms from
+ # source-scaled-voxels to
+ # ref-scaled-voxels
+ srcToRefMat = img.header.get_sform()
+
+ # The fieldToRefMat affine tells
+ # the CoefficientField class how
+ # to transform coefficient field
+ # voxel coordinates into
+ # displacement field/reference
+ # image voxel coordinates.
+ fieldToRefMat = affine.scaleOffsetXform(knotSpacing, 0)
+
+ return nonlinear.CoefficientField(fname,
+ src,
+ ref,
+ srcSpace='fsl',
+ refSpace='fsl',
+ fieldType=fieldType,
+ knotSpacing=knotSpacing,
+ srcToRefMat=srcToRefMat,
+ fieldToRefMat=fieldToRefMat)
+
+
+def readFnirt(fname, src, ref, defType=None):
+ """Reads a non-linear FNIRT transformation image, returning
+ a :class:`.DeformationField` or :class:`.CoefficientField` depending
+ on the file type.
+
+ :arg fname: File name of FNIRT transformation
+ :arg src: Source image
+ :arg ref: Reference image
+ :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.
+ """
+
+ # 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 img.intent in disps:
+ return _readFnirtDeformationField(fname, img, src, ref, defType)
+ elif img.intent in coefs:
+ return _readFnirtCoefficientField(fname, img, src, ref)
+ else:
+ raise ValueError('Cannot determine type of nonlinear '
+ 'file {}'.format(fname))
+
+
+def toFnirt(field):
+ """Convert a :class:`.NonLinearTransform` to a FNIRT-compatible
+ :class:`.DeformationField` or :class:`.CoefficientField`.
+
+ :arg field: :class:`.NonLinearTransform` to convert
+ :return: A FNIRT-compatible :class:`.DeformationField` or
+ :class:`.CoefficientField`.
+ """
+
+ # If we have a coefficient field
+ # which transforms between fsl
+ # space, we can just create a copy.
+ if isinstance(field, nonlinear.CoefficientField) and \
+ (field.srcSpace == 'fsl' and field.refSpace == 'fsl'):
+
+ # We start with a nibabel image,
+ # as we need to mess with the header
+ # fields to store all of the FNIRT
+ # coefficient field information
+ fieldBase = nib.nifti1.Nifti1Image(field.data, None)
+
+ # Set the intent
+ if field.fieldType == 'cubic':
+ intent = constants.FSL_CUBIC_SPLINE_COEFFICIENTS
+ elif field.fieldType == 'quadratic':
+ intent = constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS
+ fieldBase.header['intent_code'] = intent
+
+ # Reference image pixdims are
+ # stored in the intent code
+ # parameters.
+ fieldBase.header['intent_p1'] = field.ref.pixdim[0]
+ fieldBase.header['intent_p2'] = field.ref.pixdim[1]
+ fieldBase.header['intent_p3'] = field.ref.pixdim[2]
+
+ # Pixdims are used to
+ # store the knot spacing,
+ pixdims = list(field.knotSpacing) + [1]
+ qform = np.diag(pixdims)
+
+ # The sform is used to store the
+ # initial src-to-ref affine
+ if field.srcToRefMat is not None: sform = field.srcToRefMat
+ else: sform = np.eye(4)
+
+ # The qform offsets are
+ # used to store the
+ # reference image shape
+ qform[:3, 3] = field.ref.shape[:3]
+
+ fieldBase.header.set_zooms(pixdims)
+ fieldBase.set_sform(sform, 1)
+ fieldBase.set_qform(qform, 1)
+ fieldBase.update_header()
+
+ field = nonlinear.CoefficientField(
+ fieldBase,
+ src=field.src,
+ ref=field.ref,
+ srcSpace='fsl',
+ refSpace='fsl',
+ fieldType=field.fieldType,
+ knotSpacing=field.knotSpacing,
+ fieldToRefMat=field.fieldToRefMat,
+ srcToRefMat=field.srcToRefMat)
+
+ # Otherwise we have a non-FSL coefficient
+ # field, or a deformation field.
+ else:
+
+ # We can't convert a CoefficientField
+ # which doesn't transform in FSL
+ # coordinates, because the coefficients
+ # will have been calculated between some
+ # other source/reference coordinate
+ # systems, and we can't adjust the
+ # coefficients to encode an FSL->FSL
+ # deformation.
+ if isinstance(field, nonlinear.CoefficientField):
+ field = nonlinear.coefficientFieldToDeformationField(field)
+
+ # Again, if we have a displacement
+ # field which transforms between
+ # fsl spaces, we can just take a copy
+ if field.srcSpace == 'fsl' and field.refSpace == 'fsl':
+ field = nonlinear.DeformationField(
+ field.data,
+ header=field.header,
+ src=field.src,
+ ref=field.ref,
+ defType=field.deformationType)
+
+ # Otherwise we have to adjust the
+ # displacements so they transform
+ # between fsl coordinates.
+ field = nonlinear.convertDeformationSpace(
+ field, from_='fsl', to='fsl')
+
+ field.header['intent_code'] = constants.FSL_FNIRT_DISPLACEMENT_FIELD
+
+ return field
+
+
+def fromFnirt(field, from_='world', to='world'):
+ """Convert a FNIRT-style :class:`.NonLinearTransform` to a generic
+ :class:`.DeformationField`.
+
+ :arg field: A FNIRT-style :class:`.CoefficientField` or
+ :class:`.DeformationField`
+
+ :arg from_: Desired reference image coordinate system
+
+ :arg to: Desired source image coordinate system
+
+ :return: A :class:`.DeformationField` which contains displacements
+ from the reference image ``from_`` cordinate system to the
+ source image ``to`` coordinate syste.
+ """
+ if isinstance(field, nonlinear.CoefficientField):
+ field = nonlinear.coefficientFieldToDeformationField(field)
+ return nonlinear.convertDeformationSpace(field, from_=from_, to=to)
diff --git a/fsl/transform/nonlinear.py b/fsl/transform/nonlinear.py
new file mode 100644
index 0000000000000000000000000000000000000000..3010cda86f246e1e8697a4cfaee412a000fc91b5
--- /dev/null
+++ b/fsl/transform/nonlinear.py
@@ -0,0 +1,886 @@
+#!/usr/bin/env python
+#
+# nonlinear.py - Functions/classes for non-linear transformations.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains data structures and functions for working with
+FNIRT-style nonlinear transformations.
+
+
+The :class:`DeformationField` and :class:`CoefficientField` can be used to
+load and interact with FNIRT-style transformation images. The following
+utility functions are also available:
+
+
+.. autosummary::
+ :nosignatures:
+
+ detectDeformationType
+ convertDeformationType
+ convertDeformationSpace
+ applyDeformation
+ coefficientFieldToDeformationField
+"""
+
+
+import logging
+import itertools as it
+
+import numpy as np
+import scipy.ndimage.interpolation as ndinterp
+
+import fsl.utils.memoize as memoize
+import fsl.data.image as fslimage
+import fsl.utils.image.resample as resample
+from . import affine
+
+
+log = logging.getLogger(__name__)
+
+
+class NonLinearTransform(fslimage.Image):
+ """Class which represents a nonlinear transformation. This is just a base
+ class for the :class:`DeformationField` and :class:`CoefficientField`
+ classes.
+
+
+ A nonlinear transformation is an :class:`.Image` which contains
+ some mapping between a source image coordinate system and a reference image
+ coordinate system.
+
+
+ In FSL, non-linear transformations are defined in terms of the reference
+ image coordinate system. At a given location in the reference image
+ space, the non-linear mapping at that location can be used to calculate
+ the corresponding location in the source image space. Therefore, these
+ non-linear transformation effectively encode a transformation *from* the
+ reference image *to* the (unwarped) source image.
+ """
+
+
+ def __init__(self,
+ image,
+ src,
+ ref,
+ srcSpace=None,
+ refSpace=None,
+ **kwargs):
+ """Create a ``NonLinearTransform``. See the :meth:`.Nifti.getAffine`
+ method for an overview of the values that ``srcSpace`` and ``refSpace``
+ may take.
+
+ :arg image: A string containing the name of an image file to
+ load, or a :mod:`numpy` array, or a :mod:`nibabel`
+ image object.
+
+ :arg src: :class:`.Nifti` representing the source image.
+
+ :arg ref: :class:`.Nifti` representing the reference image.
+
+ :arg srcSpace: Coordinate system in the source image that this
+ ``NonLinearTransform`` maps to. Defaults to
+ ``'fsl'``.
+
+ :arg refSpace: Coordinate system in the reference image that this
+ ``NonLinearTransform`` maps from. Defaults to
+ ``'fsl'``.
+
+ All other arguments are passed through to :meth:`.Image.__init__`.
+ """
+
+ # TODO Could make more general by replacing
+ # srcSpace and refSpace with src/ref affines,
+ # which transform tofrom (e.g.) source/ref
+ # voxels to/from the space required by the
+ # deformation field
+
+ if srcSpace is None: srcSpace = 'fsl'
+ if refSpace is None: refSpace = 'fsl'
+
+ if srcSpace not in ('fsl', 'voxel', 'world') or \
+ refSpace not in ('fsl', 'voxel', 'world'):
+ raise ValueError('Invalid source/reference space: {} -> {}'.format(
+ srcSpace, refSpace))
+
+ fslimage.Image.__init__(self, image, **kwargs)
+
+ self.__src = fslimage.Nifti(src.header.copy())
+ self.__ref = fslimage.Nifti(ref.header.copy())
+ self.__srcSpace = srcSpace
+ self.__refSpace = refSpace
+
+
+ @property
+ def src(self):
+ """Return a reference to the :class:`.Nifti` instance representing
+ the source image.
+ """
+ return self.__src
+
+
+ @property
+ def ref(self):
+ """Return a reference to the :class:`.Nifti` instance representing
+ the reference image.
+ """
+ return self.__ref
+
+
+ @property
+ def srcSpace(self):
+ """Return the source image coordinate system this
+ ``NonLinearTransform`` maps from - see :meth:`.Nifti.getAffine`.
+ """
+ return self.__srcSpace
+
+
+ @property
+ def refSpace(self):
+ """Return the reference image coordinate system this
+ ``NonLinearTransform`` maps to - see :meth:`.Nifti.getAffine`.
+ """
+ return self.__refSpace
+
+
+ def transform(self, coords, from_=None, to=None):
+ """Transform coordinates from the reference image space to the source
+ image space. Implemented by sub-classes.
+
+ See the :meth:`.Nifti.getAffine` method for an overview of the values
+ that ``from_`` and ``to`` may take.
+
+ :arg coords: A sequence of XYZ coordinates, or ``numpy`` array of shape
+ ``(n, 3)`` containing ``n`` sets of coordinates in the
+ reference space.
+
+ :arg from_: Reference image space that ``coords`` are defined in
+
+ :arg to: Source image space to transform ``coords`` into
+
+ :returns: The corresponding coordinates in the source image space.
+ """
+ raise NotImplementedError()
+
+
+class DeformationField(NonLinearTransform):
+ """Class which represents a deformation (a.k.a. warp) field which, at each
+ voxel, contains either:
+
+ - a relative displacement from the reference image space to the source
+ image space, or
+ - absolute coordinates in the source space
+
+
+ It is assumed that the a ``DeformationField`` is aligned with the
+ reference image in their world coordinate systems (i.e. their ``sform``
+ affines project the reference image and the deformation field into
+ alignment).
+ """
+
+
+ def __init__(self, image, src, ref=None, **kwargs):
+ """Create a ``DisplacementField``.
+
+ :arg ref: Optional. If not provided, it is assumed that the
+ reference is defined in the same space as ``image``.
+
+ :arg defType: Either ``'absolute'`` or ``'relative'``, indicating
+ the type of this displacement field. If not provided,
+ will be inferred via the :func:`detectDeformationType`
+ function.
+
+ All other arguments are passed through to
+ :meth:`NonLinearTransform.__init__`.
+ """
+
+ if ref is None:
+ ref = self
+
+ defType = kwargs.pop('defType', None)
+
+ if defType not in (None, 'relative', 'absolute'):
+ raise ValueError('Invalid value for defType: {}'.format(defType))
+
+ NonLinearTransform.__init__(self, image, src, ref, **kwargs)
+
+ self.__defType = defType
+
+
+ @property
+ def deformationType(self):
+ """The type of this ``DeformationField`` - ``'absolute'`` or
+ ``'relative'``.
+ """
+ if self.__defType is None:
+ self.__defType = detectDeformationType(self)
+ return self.__defType
+
+
+ @property
+ def absolute(self):
+ """``True`` if this ``DeformationField`` contains absolute
+ coordinates.
+ """
+ return self.deformationType == 'absolute'
+
+
+ @property
+ def relative(self):
+ """``True`` if this ``DeformationField`` contains relative
+ displacements.
+ """
+ return self.deformationType == 'relative'
+
+
+ def transform(self, coords, from_=None, to=None):
+ """Transform the given XYZ coordinates from the reference image space
+ to the source image space.
+
+ :arg coords: A sequence of XYZ coordinates, or ``numpy`` array of shape
+ ``(n, 3)`` containing ``n`` sets of coordinates in the
+ reference space.
+
+ :arg from_: Reference image space that ``coords`` are defined in
+
+ :arg to: Source image space to transform ``coords`` into
+
+ :returns: ``coords``, transformed into the source image space
+ """
+
+ if from_ is None: from_ = self.refSpace
+ if to is None: to = self.srcSpace
+
+ coords = np.asanyarray(coords)
+
+ # We may need to pre-transform the
+ # coordinates so they are in the
+ # same reference image space as the
+ # displacements
+ if from_ != self.refSpace:
+ xform = self.ref.getAffine(from_, self.refSpace)
+ coords = affine.transform(coords, xform)
+
+ # We also need to get the coordinates
+ # in field voxels, so we can look up
+ # the displacements/coordinates. We
+ # can get this through the assumption
+ # that field and ref are aligned in
+ # the world coordinate system
+ xform = affine.concat(self .getAffine('world', 'voxel'),
+ self.ref.getAffine(self.refSpace, 'world'))
+
+ if np.all(np.isclose(xform, np.eye(4))):
+ voxels = coords
+ else:
+ voxels = affine.transform(coords, xform)
+
+ # Mask out the coordinates
+ # that are out of bounds of
+ # the deformation field
+ voxels = np.round(voxels).astype(np.int)
+ voxmask = (voxels >= [0, 0, 0]) & (voxels < self.shape[:3])
+ voxmask = voxmask.all(axis=1)
+ voxels = voxels[voxmask]
+
+ xs, ys, zs = voxels.T
+
+ if self.absolute: disps = self.data[xs, ys, zs, :]
+ else: disps = self.data[xs, ys, zs, :] + coords[voxmask]
+
+ # Make sure the coordinates are in
+ # the requested source image space
+ if to != self.srcSpace:
+ xform = self.src.getAffine(self.srcSpace, to)
+ disps = affine.transform(disps, xform)
+
+ # Nans for input coordinates which
+ # were outside of the field
+ outcoords = np.full(coords.shape, np.nan)
+ outcoords[voxmask] = disps
+
+ return outcoords
+
+
+class CoefficientField(NonLinearTransform):
+ """Class which represents a B-spline coefficient field generated by FNIRT.
+
+ The :meth:`displacements` method can be used to calculate relative
+ displacements for a set of reference space voxel coordinates.
+
+
+ A FNIRT nonlinear transformation often contains a *premat*, a global
+ affine transformation from the source space to the reference space, which
+ was calculated with FLIRT, and used as the starting point for the
+ non-linear optimisation performed by FNIRT.
+
+
+ This affine may be provided when creating a ``CoefficientField`` as the
+ ``srcToRefMat`` argument to :meth:`__init__`, and is subsequently accessed
+ via the :meth:`srcToRefMat` attribute.
+ """
+
+
+ def __init__(self,
+ image,
+ src,
+ ref,
+ srcSpace=None,
+ refSpace=None,
+ fieldType='cubic',
+ knotSpacing=None,
+ fieldToRefMat=None,
+ srcToRefMat=None,
+ **kwargs):
+ """Create a ``CoefficientField``.
+
+ :arg fieldType: Must currently be ``'cubic'``
+
+ :arg knotSpacing: A tuple containing the spline knot spacings along
+ each axis.
+
+ :arg fieldToRefMat: Affine transformation which can transform reference
+ image voxel coordinates into coefficient field
+ voxel coordinates.
+
+ :arg srcToRefMat: Optional initial global affine transformation from
+ the source image to the reference image. This is
+ assumed to be a FLIRT-style matrix, i.e. it
+ transforms from source image ``srcSpace``
+ coordinates into reference image ``refSpace``
+ coordinates (typically ``'fsl'`` coordinates, i.e.
+ scaled voxels potentially with a left-right flip).
+
+ See the :class:`NonLinearTransform` class for details on the other
+ arguments.
+ """
+
+ if fieldType not in ('cubic',):
+ raise ValueError('Unsupported field type: {}'.format(fieldType))
+
+ if srcToRefMat is not None: srcToRefMat = np.copy(srcToRefMat)
+ if fieldToRefMat is None: fieldToRefMat = np.eye(4)
+
+ NonLinearTransform.__init__(self,
+ image,
+ src,
+ ref,
+ srcSpace,
+ refSpace,
+ **kwargs)
+
+ self.__fieldType = fieldType
+ self.__knotSpacing = tuple(knotSpacing)
+ self.__refToSrcMat = None
+ self.__srcToRefMat = srcToRefMat
+ self.__fieldToRefMat = np.copy(fieldToRefMat)
+ self.__refToFieldMat = affine.invert(self.__fieldToRefMat)
+
+ if srcToRefMat is not None:
+ self.__refToSrcMat = affine.invert(srcToRefMat)
+
+
+ @property
+ def fieldType(self):
+ """Return the type of the coefficient field, currently always
+ ``'cubic'``.
+ """
+ return self.__fieldType
+
+
+ @property
+ def knotSpacing(self):
+ """Return a tuple containing spline knot spacings along the x, y, and
+ z axes.
+ """
+ return self.__knotSpacing
+
+
+ @property
+ def fieldToRefMat(self):
+ """Return an affine transformation which can transform coefficient
+ field voxel coordinates into reference image voxel coordinates.
+ """
+ return np.copy(self.__fieldToRefMat)
+
+
+ @property
+ def refToFieldMat(self):
+ """Return an affine transformation which can transform reference
+ image voxel coordinates into coefficient field voxel coordinates.
+ """
+ return np.copy(self.__refToFieldMat)
+
+
+ @property
+ def srcToRefMat(self):
+ """Return the initial source-to-reference affine, or ``None`` if
+ there isn't one.
+ """
+ return self.__srcToRefMat
+
+
+ @property
+ def refToSrcMat(self):
+ """Return the inverse of the initial source-to-reference affine, or
+ ``None`` if there isn't one.
+ """
+ return self.__refToSrcMat
+
+
+ @memoize.Instanceify(memoize.memoize)
+ def asDeformationField(self, defType='relative', premat=True):
+ """Convert this ``CoefficientField`` to a :class:`DeformationField`.
+
+ This method is a wrapper around
+ :func:`coefficientFieldToDeformationField`
+ """
+ return coefficientFieldToDeformationField(self, defType, premat)
+
+
+ def transform(self, coords, from_=None, to=None, premat=True):
+ """Overrides :meth:`NonLinearTransform.transform`. Transforms the
+ given ``coords`` from the reference image space into the source image
+ space.
+
+ :arg coords: A sequence of XYZ coordinates, or ``numpy`` array of shape
+ ``(n, 3)`` containing ``n`` sets of coordinates in the
+ reference space.
+
+ :arg from_: Reference image space that ``coords`` are defined in
+
+ :arg to: Source image space to transform ``coords`` into
+
+ :arg premat: If ``True``, the inverse :meth:`srcToRefMat` is applied
+ to the coordinates after the displacements have been
+ added.
+
+ :returns: ``coords``, transformed into the source image space
+ """
+ df = self.asDeformationField(premat=premat)
+ return df.transform(coords, from_, to)
+
+
+ def displacements(self, coords):
+ """Calculate the relative displacements for the given coordinates.
+
+ :arg coords: ``(N, 3)`` array of reference image voxel coordinates.
+ :return: A ``(N, 3)`` array of relative displacements to the
+ source image for ``coords``
+ """
+
+ if self.fieldType != 'cubic':
+ raise NotImplementedError()
+
+ # See
+ # https://www.cs.jhu.edu/~cis/cista/746/papers/\
+ # RueckertFreeFormBreastMRI.pdf
+ # https://www.fmrib.ox.ac.uk/datasets/techrep/tr07ja2/tr07ja2.pdf
+
+ # Cubic b-spline basis functions
+ def b0(u):
+ return ((1 - u) ** 3) / 6
+
+ def b1(u):
+ return (3 * (u ** 3) - 6 * (u ** 2) + 4) / 6
+
+ def b2(u):
+ return (-3 * (u ** 3) + 3 * (u ** 2) + 3 * u + 1) / 6
+
+ def b3(u):
+ return (u ** 3) / 6
+
+ b = [b0, b1, b2, b3]
+
+ 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: coefficient field indices
+ # u, v, w: position of the ref voxel
+ # on the current spline
+ u = np.remainder(i, 1)
+ v = np.remainder(j, 1)
+ w = np.remainder(k, 1)
+ i = np.floor(i).astype(np.int)
+ j = np.floor(j).astype(np.int)
+ k = np.floor(k).astype(np.int)
+
+ disps = np.zeros(coords.shape)
+
+ for l, m, n in it.product(range(4), range(4), range(4)):
+ il = i + l
+ jm = j + m
+ kn = k + n
+ mask = (il >= 0) & \
+ (il < nx) & \
+ (jm >= 0) & \
+ (jm < ny) & \
+ (kn >= 0) & \
+ (kn < nz)
+
+ il = il[mask]
+ jm = jm[mask]
+ kn = kn[mask]
+ uu = u[ mask]
+ vv = v[ mask]
+ ww = w[ mask]
+
+ cx, cy, cz = fdata[il, jm, kn, :].T
+ c = b[l](uu) * b[m](vv) * b[n](ww)
+
+ disps[mask, 0] += c * cx
+ disps[mask, 1] += c * cy
+ disps[mask, 2] += c * cz
+
+ return disps
+
+
+def detectDeformationType(field):
+ """Attempt to automatically determine whether a deformation field is
+ specified in absolute or relative coordinates.
+
+ :arg field: A :class:`DeformationField`
+
+ :returns: ``'absolute'`` if it looks like ``field`` contains absolute
+ coordinates, ``'relative'`` otherwise.
+ """
+
+ # This test is based on the assumption
+ # that a deformation field containing
+ # absolute coordinates will have a
+ # greater standard deviation than one
+ # which contains relative coordinates.
+ absdata = field[:]
+ reldata = convertDeformationType(field, 'relative')
+ stdabs = absdata.std(axis=(0, 1, 2)).sum()
+ stdrel = reldata.std(axis=(0, 1, 2)).sum()
+
+ if stdabs > stdrel: return 'absolute'
+ else: return 'relative'
+
+
+def convertDeformationType(field, defType=None):
+ """Convert a deformation field between storing absolute coordinates or
+ relative displacements.
+
+ :arg field: A :class:`DeformationField` instance
+ :arg defType: Either ``'absolute'`` or ``'relative'``. If not provided,
+ the opposite type to ``field.deformationType`` is used.
+ :returns: A ``numpy.array`` containing the adjusted deformation field.
+ """
+
+ if defType is None:
+ if field.deformationType == 'absolute': defType = 'relative'
+ else: defType = 'absolute'
+
+ # Regardless of the conversion direction,
+ # we need the coordinates of every voxel
+ # in the reference coordinate system.
+ dx, dy, dz = field.shape[:3]
+ xform = affine.concat(field.ref.getAffine('world', field.refSpace),
+ field .getAffine('voxel', 'world'))
+
+ coords = np.meshgrid(np.arange(dx),
+ np.arange(dy),
+ np.arange(dz), indexing='ij')
+ coords = np.array(coords).transpose((1, 2, 3, 0))
+ coords = affine.transform(coords.reshape((-1, 3)), xform)
+ coords = coords.reshape((dx, dy, dz, 3))
+
+ # If converting from relative to absolute,
+ # we just add the coordinates to (what is
+ # 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
+
+
+def convertDeformationSpace(field, from_, to):
+ """Adjust the source and/or reference spaces of the given deformation
+ field. See the :meth:`.Nifti.getAffine` method for the valid values for
+ the ``from_`` and ``to`` arguments.
+
+ :arg field: A :class:`DeformationField` instance
+ :arg from_: New reference image coordinate system
+ :arg to: New source image coordinate system
+
+ :returns: A new :class:`DeformationField` which transforms between
+ the reference ``from_`` coordinate system and the source ``to``
+ coordinate system.
+ """
+
+ if field.srcSpace == to and field.refSpace == from_:
+ return field
+
+ # Get the field in absolute coordinates
+ # if necessary - these are our source
+ # coordinates in the original "to" space.
+ fieldcoords = field.data
+ if field.relative: srccoords = convertDeformationType(field)
+ else: srccoords = fieldcoords
+
+ srccoords = srccoords.reshape((-1, 3))
+
+ # Now transform those source coordinates
+ # from the original source space to the
+ # source space specified by "to"
+ if to != field.srcSpace:
+
+ srcmat = field.src.getAffine(field.srcSpace, to)
+ srccoords = affine.transform(srccoords, srcmat)
+
+ # If we have been asked to return
+ # absolute coordinates, the
+ # reference "from_" coordinate
+ # system is irrelevant - we're done.
+ if field.absolute:
+ fieldcoords = srccoords
+
+ # Otherwise our deformation field
+ # will contain relative displacements
+ # between the reference image "from_"
+ # coordinate system and the source
+ # image "to" coordinate system. We
+ # need to re-calculate the relative
+ # displacements between the new
+ # reference "from_" space and source
+ # "to" space.
+ else:
+ refcoords = np.meshgrid(np.arange(field.shape[0]),
+ np.arange(field.shape[1]),
+ np.arange(field.shape[2]), indexing='ij')
+ refcoords = np.array(refcoords)
+ refcoords = refcoords.transpose((1, 2, 3, 0)).reshape((-1, 3))
+
+ xform = affine.concat(
+ field.ref.getAffine('world', from_),
+ field .getAffine('voxel', 'world'))
+
+ if not np.all(np.isclose(xform, np.eye(4))):
+ refcoords = affine.transform(refcoords, xform)
+
+ fieldcoords = srccoords - refcoords
+
+ return DeformationField(
+ fieldcoords.reshape(field.shape),
+ header=field.header,
+ src=field.src,
+ ref=field.ref,
+ srcSpace=to,
+ refSpace=from_,
+ defType=field.deformationType)
+
+
+def applyDeformation(image, field, ref=None, order=1, mode=None, cval=None):
+ """Applies a :class:`DeformationField` to an :class:`.Image`.
+
+ The image is transformed into the space of the field's reference image
+ space. See the ``scipy.ndimage.interpolation.map_coordinates`` function
+ for details on the ``order``, ``mode`` and ``cval`` options.
+
+ If an alternate reference image is provided via the ``ref`` argument,
+ the deformation field is resampled into its space, and then applied to
+ the input image. It is therefore assumed that an alternate ``ref`` is
+ aligned in world coordinates with the field's actual reference image.
+
+ :arg image: :class:`.Image` to be transformed
+
+ :arg field: :class:`DeformationField` to use
+
+ :arg ref: Alternate reference image - if not provided, ``field.ref``
+ is used
+
+ :arg order: Spline interpolation order, passed through to the
+ ``scipy.ndimage.affine_transform`` function - ``0``
+ corresponds to nearest neighbour interpolation, ``1``
+ (the default) to linear interpolation, and ``3`` to
+ cubic interpolation.
+
+ :arg mode: How to handle regions which are outside of the image FOV.
+ Defaults to `''nearest'``.
+
+ :arg cval: Constant value to use when ``mode='constant'``.
+
+ :return: ``numpy.array`` containing the transformed image data.
+ """
+
+ if order is None: order = 1
+ if mode is None: mode = 'nearest'
+ if cval is None: cval = 0
+ if ref is None: ref = field.ref
+
+ # We need the field to contain
+ # absolute source image voxel
+ # coordinates
+ field = convertDeformationSpace(field, 'voxel', 'voxel')
+ if field.deformationType != 'absolute':
+ field = DeformationField(convertDeformationType(field, 'absolute'),
+ header=field.header,
+ src=field.src,
+ ref=field.ref,
+ srcSpace='voxel',
+ refSpace='voxel',
+ defType='absolute')
+
+ # If the field is not voxel-aligned
+ # to the reference, we need to
+ # resample the field itself into the
+ # reference image space (assumed to
+ # be world-aligned). If field and ref
+ # are not not world aligned, regions
+ # of the field outside of the
+ # reference image space will contain
+ # -1s, so will be detected as out of
+ # bounds by map_coordinates below.
+ #
+ # This will potentially result in
+ # truncation at the field boundaries,
+ # but there's nothing we can do about
+ # that.
+ src = field.src
+
+ if not field.sameSpace(ref):
+ field = resample.resampleToReference(field,
+ ref,
+ order=order,
+ mode='constant',
+ cval=-1)[0]
+ else:
+ field = field.data
+
+ # If the input image is in a
+ # different space to the field
+ # source space, we need to
+ # adjust the resampling matrix.
+ # We assume world-world alignment
+ # between the original source
+ # and the image to be resampled
+ if not image.sameSpace(src):
+ post = affine.concat(image.getAffine('world', 'voxel'),
+ src .getAffine('voxel', 'world'))
+ shape = field.shape
+ field = field.reshape((-1, 3))
+ field = affine.transform(field, post)
+ field = field.reshape(shape)
+
+ field = field.transpose((3, 0, 1, 2))
+ return ndinterp.map_coordinates(image.data,
+ field,
+ order=order,
+ mode=mode,
+ cval=cval)
+
+
+def coefficientFieldToDeformationField(field, defType='relative', premat=True):
+ """Convert a :class:`CoefficientField` into a :class:`DeformationField`.
+
+ :arg field: :class:`CoefficientField` to convert
+
+ :arg defType: The type of deformation field - either ``'relative'`` (the
+ default) or ``'absolute'``.
+
+ :arg premat: If ``True`` (the default), the :meth:`srcToRefMat` is
+ encoded into the deformation field.
+
+ :return: :class:`DeformationField` calculated from ``field``.
+ """
+
+ # Generate coordinates for every
+ # voxel in the reference image
+ ix, iy, iz = field.ref.shape[:3]
+ x, y, z = np.meshgrid(np.arange(ix),
+ np.arange(iy),
+ np.arange(iz), indexing='ij')
+ x = x.flatten()
+ y = y.flatten()
+ z = z.flatten()
+ xyz = np.vstack((x, y, z)).T
+
+ # There are three spaces to consider here:
+ #
+ # - ref space: Reference image scaled voxels ("fsl" space)
+ #
+ # - aligned-src space: Source image scaled voxels, after the
+ # source image has been linearly aligned to
+ # the reference via the field.srcToRefMat
+ # This will typically be equivalent to ref
+ # space
+ #
+ # - orig-src space: Source image scaled voxels, in the coordinate
+ # system of the original source image, without
+ # linear alignment to the reference image
+
+ # The displacements method will
+ # return relative displacements
+ # from ref space to aligned-src
+ # space.
+ disps = field.displacements(xyz).reshape((ix, iy, iz, 3))
+ rdfield = DeformationField(disps,
+ header=field.ref.header,
+ src=field.src,
+ ref=field.ref,
+ srcSpace=field.srcSpace,
+ refSpace=field.refSpace,
+ defType='relative')
+
+ if (defType == 'relative') and (not premat):
+ return rdfield
+
+ # Convert to absolute - the
+ # deformations will now be
+ # absolute coordinates in
+ # aligned-src space
+ disps = convertDeformationType(rdfield)
+
+ # Apply the premat if requested -
+ # this will transform the coordinates
+ # from aligned-src to orig-src space.
+ if premat and field.srcToRefMat is not None:
+
+ # We apply the premat in the same way
+ # that fnirtfileutils does - applying
+ # the inverse affine to every ref space
+ # voxel coordinate, then adding it to
+ # the existing displacements.
+ shape = disps.shape
+ disps = disps.reshape(-1, 3)
+ premat = affine.concat(field.refToSrcMat - np.eye(4),
+ field.ref.getAffine('voxel', 'fsl'))
+ disps = disps + affine.transform(xyz, premat)
+ disps = disps.reshape(shape)
+
+ # note that convertwarp applies a premat
+ # differently - its method is equivalent
+ # to directly transforming the existing
+ # absolute displacements, i.e.:
+ #
+ # disps = affine.transform(disps, refToSrc)
+
+ adfield = DeformationField(disps,
+ header=field.ref.header,
+ src=field.src,
+ ref=field.ref,
+ srcSpace=field.srcSpace,
+ refSpace=field.refSpace,
+ defType='absolute')
+
+ # Not either return absolute displacements,
+ # or convert back to relative displacements
+ if defType == 'absolute':
+ return adfield
+ else:
+ return DeformationField(convertDeformationType(adfield),
+ src=field.src,
+ ref=field.ref,
+ srcSpace=field.srcSpace,
+ refSpace=field.refSpace,
+ header=field.ref.header,
+ defType='relative')
diff --git a/fsl/transform/x5.py b/fsl/transform/x5.py
new file mode 100644
index 0000000000000000000000000000000000000000..1ffe032821a5c3aeb231cf2835011fe5d67449bb
--- /dev/null
+++ b/fsl/transform/x5.py
@@ -0,0 +1,607 @@
+#!/usr/bin/env python
+#
+# x5.py - Functions for working with BIDS X5 files.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains functions for reading/writing linear/non-linear
+transformations from/to BIDS X5 files. The following functions are available:
+
+.. autosummary::
+ :nosignatures:
+
+ inferType
+ readLinearX5
+ writeLinearX5
+ readNonLinearX5
+ writeNonLinearX5
+
+
+An X5 file is a HDF5 container file which stores a linear or non-linear
+transformation from one NIfTI image to another.
+
+
+Several terms may be used to refer to these images, such as **source** and
+**reference**, **moving** and **fixed**, **from** and **to**, etc. In an X5
+file, the two images are simply referred to as **A** and **B**, where **A**
+refers to the starting point of the transformation, and **B** to the end
+point.
+
+
+X5 files enable a transformation from the **world coordinate system** of image
+**A** to the **world coordinate system** of image **B**. The **world
+coordinate system** of an image is defined by its ``sform`` or ``qform``
+(hereafter referred to as the ``sform``), which is contained in its NIfTI
+header.
+
+
+Custom HDF5 groups
+==================
+
+
+HDF5 files are composed primarily of *groups*, *attributes*, and
+*datasets*:
+
+
+ - *Groups* are simply containers for attributes, datasets, and other groups.
+ - *Datasets* are strongly-typed, structured N-dimensional arrays.
+ - *Attributes* are strongly-typed scalar values associated with a group or
+ dataset.
+
+
+To simplify the file format definitions below, we shall first define a few
+custom HDF5 groups. In the file format definitions, a HDF5 group which is
+listed as being of one of these custom types shall contain the
+attributes/datasets that are listed here.
+
+
+*affine*
+--------
+
+
+A HDF5 group which is listed as being of type *affine* contains an affine
+transformation, which can be used to transform coordinates from one space into
+another. Groups of type *affine* have the following fields:
+
+
++-------------+-----------+---------------------------------------------------+
+| **Name** | **Type** | **Value/Description** |
++-------------+-----------+---------------------------------------------------+
+| ``Type`` | attribute | ``'affine'`` |
++-------------+-----------+---------------------------------------------------+
+| ``Matrix`` | dataset | The affine transformation matrix - a ``float64`` |
+| | | array of shape ``(4, 4)`` |
++-------------+-----------+---------------------------------------------------+
+| ``Inverse`` | dataset | Optional pre-calculated inverse |
++-------------+-----------+---------------------------------------------------+
+
+
+*space*
+-------
+
+
+A HDF5 group which is listed as being of type *space* contains all of the
+information required to define the space of a NIfTI image, including its
+shape, dimensions, and voxel-to-world affine transformation.
+
+
+Groups of type *space* have the following fields:
+
+
++--------------+-----------+--------------------------------------------------+
+| **Name** | **Type** | **Value/Description** |
++--------------+-----------+--------------------------------------------------+
+| ``Type`` | attribute | ``'image'`` |
++--------------+-----------+--------------------------------------------------+
+| ``Size`` | attribute | ``uint64`` ``(X, Y, Z)`` voxel dimensions |
++--------------+-----------+--------------------------------------------------+
+| ``Scales`` | attribute | ``float64`` ``(X, Y, Z)`` voxel pixdims |
++--------------+-----------+--------------------------------------------------+
+| ``Mapping/`` | affine | The image voxel-to-world transformation (its |
+| | | ``sform``) |
++--------------+-----------+--------------------------------------------------+
+
+
+*deformation*
+-------------
+
+
+A HDF5 group which is listed as being of type *deformation* contains a
+non-linear transformation, which can be used to transform coordinates from
+one space (space **A**) into another (space **B**).
+
+
+The transformation is represented as a 3D **deformation field** which, at each
+voxel within the field, may contain:
+
+ - *relative displacements* from space **A** to space **B** (i.e. for a given
+ location in space **A**, you can add the displacement values to the
+ coordinates of that location to obtain the coordinates of the corresponding
+ location in space **B**).
+
+ - *absolute coordinates* in space **B**.
+
+
+The ``Mapping`` affine can be used to calculate a correspondence between the
+deformation field coordinate system and the coordinate system of space **A** -
+it is assumed that space **A** and the deformation field share a common world
+coordinate system.
+
+
+Groups of type *deformation* have the following fields:
+
+
++--------------+-----------+--------------------------------------------------+
+| **Name** | **Type** | **Value/Description** |
++--------------+-----------+--------------------------------------------------+
+| ``Type`` | attribute | ``'deformation'`` |
++--------------+-----------+--------------------------------------------------+
+| ``SubType`` | attribute | ``'absolute'`` or ``'relative'``. |
++--------------+-----------+--------------------------------------------------+
+| ``Matrix`` | dataset | The deformation field - a ``float64`` array of |
+| | | shape ``(X, Y, Z, 3)`` |
++--------------+-----------+--------------------------------------------------+
+| ``Mapping/`` | affine | The field voxel-to-world transformation (its |
+| | | ``sform``) |
++--------------+-----------+--------------------------------------------------+
+
+
+Linear X5 files
+===============
+
+
+Linear X5 transformation files contain an affine transformation matrix of
+shape ``(4, 4)``, which can be used to transform image **A** world
+coordinates into image **B** world coordinates.
+
+
+Linear X5 transformation files are assumed to adhere to the HDF5 structure
+defined in the table below. All fields are required unless otherwise noted.
+
+
++-----------------+-----------+-----------------------------------------------+
+| **Name** | **Type** | **Value/Description** |
++-----------------+-----------+-----------------------------------------------+
+| *Metadata* |
++-----------------+-----------+-----------------------------------------------+
+| ``/Format`` | attribute | ``'X5'`` |
++-----------------+-----------+-----------------------------------------------+
+| ``/Version`` | attribute | ``'0.0.1'`` |
++-----------------+-----------+-----------------------------------------------+
+| ``/Metadata`` | attribute | JSON string containing unstructured metadata. |
++-----------------+-----------+-----------------------------------------------+
+| *Transformation* |
++-----------------+-----------+-----------------------------------------------+
+| ``/Type`` | attribute | ``'linear'`` |
++-----------------+-----------+-----------------------------------------------+
+| ``/Transform/`` | affine | Affine transformation from image **A** world |
+| | | coordinates to image **B** world coordinates |
++-----------------+-----------+-----------------------------------------------+
+| ``/A/`` | space | Image **A** space |
++-----------------+-----------+-----------------------------------------------+
+| ``/B/`` | space | Image **B** space |
++-----------------+-----------+-----------------------------------------------+
+
+
+Storage of FSL FLIRT matrices in linear X5 files
+------------------------------------------------
+
+
+FLIRT outputs the result of a linear registration from a source image to a
+reference image as an affine matrix of shape ``(4, 4)``. This matrix encodes a
+transformation from source image **FSL coordinates** to reference image **FSL
+coordinates** [*]_.
+
+
+In contrast, X5 matrices encode a transformation in **world coordinates**,
+i.e. they can be used to transform coordinates from the source image to the
+reference image, after both images have been transformed into a common
+coordinate system via their respective ``sform`` affines.
+
+
+The :mod:`fsl.transform.flirt` module contains functions for converting
+between FLIRT-style matrices and X5 style matrices.
+
+
+.. [*] For a given image, FSL coordinates are voxel coordinates scaled by the
+ ``pixdim`` values in the NIFTI header, with an inversion along the X
+ axis if the voxel-to-world affine (the ``sform``) has a positive
+ determinant.
+
+
+Non-linear X5 files
+===================
+
+
+Non-linear X5 transformation files contain a non-linear transformation from
+image **A** world coordinates to image **B** world coordinates. The
+transformation is represented as a 3D **deformation field** which, at each
+voxel within the field, may contain:
+
+ - *relative displacements* from image **A** to image **B** (i.e. for a given
+ location in the image **A** world coordinate system, add the displacement
+ values to the coordinates to obtain the corresponding location in the
+ image **B** world coordinate system).
+
+ - *absolute coordinates* in the image **B** world coordinate system.
+
+
+File format specification
+-------------------------
+
+
+Non-linear X5 transformation files are assumed to adhere to the following
+HDF5 structure. All fields are required unless otherwise noted.
+
+
++---------------+-----------+-------------------------------------------------+
+| **Name** | **Type** | **Value/Description** |
++---------------+-----------+-------------------------------------------------+
+| *Metadata* |
++---------------+-----------+-------------------------------------------------+
+| ``/Format`` | attribute | ``'X5'`` |
++---------------+-----------+-------------------------------------------------+
+| ``/Version`` | attribute | ``'0.0.1'`` |
++---------------+-----------+-------------------------------------------------+
+| ``/Metadata`` | attribute | JSON string containing unstructured metadata. |
++---------------+-----------+-------------------------------------------------+
+| *Transformation* |
++-----------------+-------------+---------------------------------------------+
+| **Name** | **Type** | **Value/Description** |
++-----------------+-------------+---------------------------------------------+
+| ``/Type`` | attribute | ``'nonlinear'`` |
++-----------------+-------------+---------------------------------------------+
+| ``/Transform/`` | deformation | The deformation field, encoding a nonlinear |
+| | | transformation from image **A** to image |
+| | | **B** |
++-----------------+-------------+---------------------------------------------+
+| ``/Inverse/`` | deformation | Optional pre-calculated inverse, encoding a |
+| | | nonlinear transformation from image **B** |
+| | | to image **A** |
++-----------------+-------------+---------------------------------------------+
+| ``/A/`` | space | Image **A** space |
++-----------------+-------------+---------------------------------------------+
+| ``/B/`` | space | Image **B** space |
++-----------------+-------------+---------------------------------------------+
+
+
+Storage of FSL FNIRT warp fields in non-linear X5 files
+-------------------------------------------------------
+
+
+FLIRT outputs the result of a non-linear registration between a source image
+and a reference image as either a warp field, or a coefficient field which can
+be used to generate a warp field. A warp field is defined in terms of the
+reference image - the warp field has the same shape and FOV as the reference
+image, and contains either:
+
+ - relative displacements from the corresponding reference image location to
+ the unwarped source image location
+ - absolute unwarped source image coordinates
+
+
+The reference image for a FNIRT warp field thus corresponds to image **A** in
+a X5 non-linear transform, and the FNIRT source image to image **B**.
+
+
+FNIRT warp fields are defined in FSL coordinates - a relative warp contains
+displacements from reference image FSL coordinates to source image FSL
+coordinates, and an absolute warp contains source image FSL coordinates.
+
+
+When a FNIRT warp field is stored in an X5 file, the displacements/coordinates
+must be adjusted so that they encode a transformation from reference image
+world coordinates to source image world coordinates.
+
+
+Conversion of FNIRT coefficient fields
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+
+A FNIRT coefficient field can be used to generate a deformation field which
+contains relative displacements from reference image FSL coordinates to source
+image FSL coordinates. If an initial affine registration was used as the
+starting point for FNIRT, this generated displacement field contains relative
+displacements from the reference image to the *aligned* source image,
+i.e. after it has been transformed by the initial affine alignment.
+
+
+When a FNIRT coefficient field is stored in an X5 file, it must first be
+converted to a displacement field. The displacements must then be adjusted so
+that they take into account the initial affine alignment (if relevant), and
+so that they encode displacements from reference image world coordinates
+to source image world coordinates.
+
+
+The :mod:`fsl.transform.fnirt` module contains functions which can be used to
+perform all of the conversions and adjustments required to store FNIRT
+transformations as X5 files.
+
+"""
+
+
+import json
+
+import numpy as np
+import nibabel as nib
+import h5py
+
+import fsl.version as fslversion
+import fsl.data.image as fslimage
+from . import affine
+from . import nonlinear
+
+
+X5_FORMAT = 'X5'
+X5_VERSION = '0.1.0'
+
+
+class X5Error(Exception):
+ """Error raised if an invalid/incompatible file is detected. """
+ pass
+
+
+def inferType(fname):
+ """Return the type of the given X5 file - either ``'linear'`` or
+ ``'nonlinear'``.
+
+ :arg fname: Name of a X5 file
+ :returns: ``'linear'`` or ``'nonlinear'``
+ """
+ with h5py.File(fname, 'r') as f:
+
+ ftype = f.attrs.get('Type')
+
+ if ftype not in ('linear', 'nonlinear'):
+ raise X5Error('Unknown type: {}'.format(ftype))
+
+ return ftype
+
+
+def readLinearX5(fname):
+ """Read a linear X5 transformation file from ``fname``.
+
+ :arg fname: File name to read from
+ :returns: A tuple containing:
+
+ - A ``(4, 4)`` ``numpy`` array containing the affine
+ transformation
+ - A :class:`.Nifti` instance representing the source space
+ - A :class:`.Nifti` instance representing the reference space
+ """
+
+ with h5py.File(fname, 'r') as f:
+
+ if f.attrs['Type'] != 'linear':
+ raise X5Error('Not a linear transform')
+
+ _readMetadata( f['/'])
+ xform = _readAffine(f['/Transform'])
+ src = _readSpace( f['/A'])
+ ref = _readSpace( f['/B'])
+
+ return xform, src, ref
+
+
+def writeLinearX5(fname, xform, src, ref):
+ """Write a linear transformation to ``fname``.
+
+ :arg fname: File name to write to
+ :arg xform: ``(4, 4)`` ``numpy`` array containing the affine transformation
+ :arg src: :class:`.Nifti` image representing the source space
+ :arg ref: :class:`.Nifti` image representing the reference
+ """
+
+ with h5py.File(fname, 'w') as f:
+
+ f.attrs['Type'] = 'linear'
+
+ _writeMetadata(f)
+ _writeAffine( f.create_group('/Transform'), xform)
+ _writeSpace( f.create_group('/A'), src)
+ _writeSpace( f.create_group('/B'), ref)
+
+
+def readNonLinearX5(fname):
+ """Read a nonlinear X5 transformation file from ``fname``.
+
+ :arg fname: File name to read from
+ :returns: A :class:`.DeformationField`
+ """
+
+ with h5py.File(fname, 'r') as f:
+
+ if f.attrs.get('Type') != 'nonlinear':
+ raise X5Error('Not a nonlinear transform')
+
+ _readMetadata(f)
+
+ ref = _readSpace( f['/A'])
+ src = _readSpace( f['/B'])
+ field, xform, defType = _readDeformation(f['/Transform'])
+
+ return nonlinear.DeformationField(field,
+ xform=xform,
+ src=src,
+ ref=ref,
+ srcSpace='world',
+ refSpace='world',
+ defType=defType)
+
+
+def writeNonLinearX5(fname, field):
+ """Write a nonlinear X5 transformation to ``fname``.
+
+ :arg fname: File name to write to
+ :arg field: A :class:`.DeformationField`
+ """
+
+ with h5py.File(fname, 'w') as f:
+
+ f.attrs['Type'] = 'nonlinear'
+
+ _writeMetadata(f)
+ _writeSpace( f.create_group('/A'), field.ref)
+ _writeSpace( f.create_group('/B'), field.src)
+ _writeDeformation(f.create_group('/Transform'), field)
+
+
+def _readMetadata(group):
+ """Reads a metadata block from the given group, and raises a :exc:`X5Error`
+ if it does not look valid.
+
+ :arg group: A ``h5py.Group`` object.
+ :returns: A ``dict`` containing the metadata.
+ """
+
+ format = group.attrs.get('Format')
+ version = group.attrs.get('Version')
+ meta = group.attrs.get('Metadata')
+
+ parserver = fslversion.parseVersionString(X5_VERSION)
+ filever = fslversion.parseVersionString(version)
+
+ if (format != X5_FORMAT) or (filever[0] != parserver[0]):
+ raise X5Error('Incompatible format/version (required: {}/{}, '
+ 'present: {}/{})'.format(X5_FORMAT, X5_VERSION,
+ format, version))
+
+ meta = {}
+ meta['Format'] = format
+ meta['Version'] = version
+ meta['Metadata'] = meta
+
+ return meta
+
+
+def _writeMetadata(group):
+ """Writes a metadata block into the given group.
+
+ :arg group: A ``h5py.Group`` object.
+ """
+ group.attrs['Format'] = X5_FORMAT
+ group.attrs['Version'] = X5_VERSION
+ group.attrs['Metadata'] = json.dumps({'fslpy' : fslversion.__version__})
+
+
+def _readAffine(group):
+ """Reads an *affine* from the given group.
+
+ :arg group: A ``h5py.Group`` object.
+ :returns: ``numpy`` array containing a ``(4, 4)`` affine transformation.
+ """
+
+ if group.attrs.get('Type') != 'affine':
+ raise X5Error('Not an affine')
+
+ xform = group['Matrix']
+
+ if xform.shape != (4, 4):
+ raise X5Error('Not an affine')
+
+ return np.array(xform)
+
+
+def _writeAffine(group, xform):
+ """Writes the given affine transformation and its inverse to the given
+ group.
+
+ :arg group: A ``h5py.Group`` object.
+ :arg xform: ``numpy`` array containing a ``(4, 4)`` affine transformation.
+ """
+
+ xform = np.asarray(xform, dtype=np.float64)
+ inv = np.asarray(affine.invert(xform), dtype=np.float64)
+
+ group.attrs['Type'] = 'affine'
+ group.create_dataset('Matrix', data=xform)
+ group.create_dataset('Inverse', data=inv)
+
+
+def _readSpace(group):
+ """Reads a *space* group, defining a source or reference space.
+
+ :arg group: A ``h5py.Group`` object.
+ :returns: :class:`.Nifti` object. defining the mapping
+ """
+
+ if group.attrs.get('Type') != 'image':
+ raise X5Error('Not an image space')
+
+ shape = group.attrs['Size']
+ pixdim = group.attrs['Scales']
+ xform = _readAffine(group['Mapping'])
+
+ hdr = nib.Nifti2Header()
+ hdr.set_data_shape(shape)
+ hdr.set_zooms( pixdim)
+ hdr.set_sform( xform, 'aligned')
+ return fslimage.Nifti(hdr)
+
+
+def _writeSpace(group, img):
+ """Writes a space specified by ``img`` to the given group.
+
+ :arg group: A ``h5py.Group`` object.
+ :arg img: :class:`.Nifti` object. defining the mapping
+ """
+ group.attrs['Type'] = 'image'
+ group.attrs['Size'] = np.asarray(img.shape[ :3], np.uint32)
+ group.attrs['Scales'] = np.asarray(img.pixdim[:3], np.float32)
+
+ mapping = group.create_group('Mapping')
+ _writeAffine(mapping, img.getAffine('voxel', 'world'))
+
+
+def _readDeformation(group):
+ """Reads a *deformation* from the given group.
+
+ :arg group: A ``h5py.Group`` object
+ :returns: A tuple containing
+
+ - A ``numpy.array`` containing the deformation field
+
+ - A ``numpy.array`` of shape ``(4, 4)`` containing the
+ voxel to world affine for the deformation field
+
+ - The deformation type - either ``'absolute'`` or
+ ``'relative'``
+ """
+
+ type = group.attrs.get('Type')
+ subtype = group.attrs['SubType']
+
+ if type != 'deformation':
+ raise X5Error('Not a deformation')
+
+ if subtype not in ('absolute', 'relative'):
+ raise X5Error('Unknown deformation type: {}'.format(subtype))
+
+ mapping = _readAffine(group['Mapping'])
+ field = group['Matrix']
+
+ if len(field.shape) != 4 or field.shape[3] != 3:
+ raise X5Error('Invalid shape for deformation field')
+
+ return np.array(field), mapping, subtype
+
+
+def _writeDeformation(group, field):
+ """Write a deformation field to the given group.
+
+ :arg group: A ``h5py.Group`` object
+ :arg field: A :class:`.DeformationField` object
+ """
+
+ if field.srcSpace != 'world' or \
+ field.refSpace != 'world':
+ raise X5Error('Deformation field must encode a '
+ 'world<->world transformation')
+
+ group.attrs['Type'] = 'deformation'
+ group.attrs['SubType'] = field.deformationType
+
+ mapping = group.create_group('Mapping')
+
+ group.create_dataset('Matrix', data=field.data)
+ _writeAffine(mapping, field.getAffine('voxel', 'world'))
diff --git a/fsl/utils/deprecated.py b/fsl/utils/deprecated.py
index c79dddfeb9073420c4d5449edd42ed50d9783b40..64a0ad90885cd72fe550b61826983b66e171231e 100644
--- a/fsl/utils/deprecated.py
+++ b/fsl/utils/deprecated.py
@@ -6,8 +6,10 @@
#
"""This module provides the :func:`deprecated` function, a simple decorator
for deprecating functions and methods.
-"""
+The :func:`warn` function can also be called directly, to emit a
+``DeprecationWarning``
+"""
import functools as ft
@@ -16,14 +18,20 @@ import warnings
_warned_cache = set()
-"""Used by the :func:`deprecated` function to keep track of whether a warning
-has already been emitted for the use of a deprecated item.
+"""Used by to keep track of whether a warning has already been emitted for the
+use of a deprecated item.
"""
-def deprecated(vin=None, rin=None, msg=None):
- """Decorator to mark a function or method as deprecated. A
- ``DeprecationWarning`` is raised via the standard ``warnings`` module.
+def resetWarningCache():
+ """Clears the internal warning cache, so that the same line of code
+ may emit another deprecation warning.
+ """
+ _warned_cache.clear()
+
+
+def _buildMessageFormat(vin=None, rin=None, msg=None):
+ """Builds a deprecation warning message from the arguments.
:arg vin: Optional version - the warning message will mention that the
function is deprecated from this version.
@@ -32,9 +40,9 @@ def deprecated(vin=None, rin=None, msg=None):
function will be removed in this version.
:arg msg: Optional message to use in the warning.
- """
-
+ :returns: A format string which needs to be formatted with a ``{name}``.
+ """
if vin is not None and rin is not None:
msgfmt = '{{name}} is deprecated from version {vin} and will be ' \
'removed in {rin}.'.format(vin=vin, rin=rin)
@@ -49,21 +57,62 @@ def deprecated(vin=None, rin=None, msg=None):
if msg is not None:
msgfmt = msgfmt + ' ' + msg
- def wrapper(thing):
- name = thing.__name__
+ return msgfmt
- def decorator(*args, **kwargs):
- frame = inspect.stack()[1]
- ident = '{}:{}'.format(frame.filename, frame.lineno)
+def _buildWarningSourceIdentity(stacklevel=2):
+ """Creates a string to be used as an identifier for the calling code.
+
+ :arg stacklevel: How far up the calling stack the calling code is.
+ :returns: A string which can be used as an identifier for the
+ calling code.
+ """
+ frame = inspect.stack()[stacklevel]
+ ident = '{}:{}'.format(frame.filename, frame.lineno)
+ return ident
+
+
+def warn(name, vin=None, rin=None, msg=None, stacklevel=1):
+ """Emit a deprecation warning.
+
+ :arg name: Name of the thing (class, function, module, etc) that is
+ deprecated.
+
+ :arg vin: Optional version - the warning message will mention that
+ the function is deprecated from this version.
+
+ :arg rin: Optional version - the warning message will mention that
+ the function will be removed in this version.
- if ident not in _warned_cache:
- warnings.warn(msgfmt.format(name=name),
- category=DeprecationWarning,
- stacklevel=2)
- _warned_cache.add(ident)
+ :arg msg: Optional message to use in the warning.
+ :arg stacklevel: How far up the stack the calling code is.
+ """
+
+ msgfmt = _buildMessageFormat(vin=vin, rin=rin, msg=msg)
+ ident = _buildWarningSourceIdentity()
+ if ident not in _warned_cache:
+ warnings.warn(msgfmt.format(name=name),
+ category=DeprecationWarning,
+ stacklevel=stacklevel + 1)
+ _warned_cache.add(ident)
+
+
+def deprecated(vin=None, rin=None, msg=None):
+ """Decorator to mark a function or method as deprecated. A
+ ``DeprecationWarning`` is raised via the standard ``warnings`` module.
+
+ :arg vin: Optional version - the warning message will mention that the
+ function is deprecated from this version.
+
+ :arg rin: Optional version - the warning message will mention that the
+ function will be removed in this version.
+
+ :arg msg: Optional message to use in the warning.
+ """
+ def wrapper(thing):
+ def decorator(*args, **kwargs):
+ warn(thing.__name__, vin=vin, rin=rin, msg=msg, stacklevel=2)
return thing(*args, **kwargs)
return ft.update_wrapper(decorator, thing)
-
return wrapper
diff --git a/fsl/utils/idle.py b/fsl/utils/idle.py
index 1425ba191dd6d954c074369088c0a33d532dfeff..c1513708ebb401227dd1d0adb7cf9e0efd067631 100644
--- a/fsl/utils/idle.py
+++ b/fsl/utils/idle.py
@@ -278,7 +278,7 @@ class IdleTask(object):
def _wxIdleLoop(ev):
"""Function which is called on ``wx.EVT_IDLE`` events, and occasionally
- on ``wx.EVT_TIMER` events via the :attr:`_idleTimer`. If there
+ on ``wx.EVT_TIMER`` events via the :attr:`_idleTimer`. If there
is a function on the :attr:`_idleQueue`, it is popped and called.
.. note:: The ``wx.EVT_IDLE`` event is only triggered on user interaction
diff --git a/fsl/utils/image/__init__.py b/fsl/utils/image/__init__.py
index 36bde925ddd642203560caf5ac58ff99a364d386..54dde4b2a7c40e1c97f8b7fcfaa4f86d1b5781b5 100644
--- a/fsl/utils/image/__init__.py
+++ b/fsl/utils/image/__init__.py
@@ -13,4 +13,5 @@ The following modules are available:
:nosignature
.image.resample
+ .image.roi
"""
diff --git a/fsl/utils/image/resample.py b/fsl/utils/image/resample.py
index 91e79e6c0416cc01187bc4c636f59ec56aa23668..50dc0c4022cc5b32f763d8afd43f3d7fb31fba47 100644
--- a/fsl/utils/image/resample.py
+++ b/fsl/utils/image/resample.py
@@ -15,12 +15,12 @@ sub-functions of :func:`resample`.
"""
-import collections.abc as abc
+import collections.abc as abc
-import numpy as np
-import scipy.ndimage as ndimage
+import numpy as np
+import scipy.ndimage as ndimage
-import fsl.utils.transform as transform
+import fsl.transform.affine as affine
def resampleToPixdims(image, newPixdims, **kwargs):
@@ -39,7 +39,7 @@ def resampleToPixdims(image, newPixdims, **kwargs):
return resample(image, newShape, **kwargs)
-def resampleToReference(image, reference, **kwargs):
+def resampleToReference(image, reference, matrix=None, **kwargs):
"""Resample ``image`` into the space of the ``reference``.
This is a wrapper around :func:`resample` - refer to its documenttion
@@ -49,6 +49,7 @@ def resampleToReference(image, reference, **kwargs):
along the spatial (first three) dimensions.
:arg image: :class:`.Image` to resample
+ :arg matrix: Optional world-to-world affine alignment matrix
:arg reference: :class:`.Nifti` defining the space to resample ``image``
into
"""
@@ -56,6 +57,9 @@ def resampleToReference(image, reference, **kwargs):
oldShape = list(image.shape)
newShape = list(reference.shape[:3])
+ if matrix is None:
+ matrix = np.eye(4)
+
# If the input image is >3D, pad the
# new shape so that we only resample
# along the first 3 dimensions.
@@ -63,8 +67,12 @@ def resampleToReference(image, reference, **kwargs):
newShape = newShape + oldShape[len(newShape):]
# Align the two images together
- # via their vox-to-world affines.
- matrix = transform.concat(image.worldToVoxMat, reference.voxToWorldMat)
+ # via their vox-to-world affines,
+ # and the world-to-world affine
+ # if provided
+ matrix = affine.concat(image.worldToVoxMat,
+ affine.invert(matrix),
+ reference.voxToWorldMat)
# If the input image is >3D, we
# have to adjust the resampling
@@ -82,7 +90,12 @@ def resampleToReference(image, reference, **kwargs):
kwargs['newShape'] = newShape
kwargs['matrix'] = matrix
- return resample(image, **kwargs)
+ data = resample(image, **kwargs)[0]
+
+ # The image is now in the same space
+ # as the reference, so it inherits
+ # ref's voxel-to-world affine
+ return data, reference.voxToWorldMat
def resample(image,
@@ -91,9 +104,9 @@ def resample(image,
dtype=None,
order=1,
smooth=True,
- origin='centre',
+ origin=None,
matrix=None,
- mode='nearest',
+ mode=None,
cval=0):
"""Returns a copy of the data in the ``image``, resampled to the specified
``newShape``.
@@ -113,6 +126,12 @@ def resample(image,
particularly on the ``order``, ``matrix``, ``mode`` and
``cval`` arguments.
+ .. note:: If a custom resampling ``matrix`` is specified, the adjusted
+ voxel-to-world afffine cannot be calculated by this function,
+ so ``None`` will be returned instead.
+
+ :arg image: :class:`.Image` object to resample
+
:arg newShape: Desired shape. May containg floating point values, in which
case the resampled image will have shape
``round(newShape)``, but the voxel sizes will have scales
@@ -159,13 +178,18 @@ def resample(image,
- A ``numpy`` array of shape ``(4, 4)``, containing the
adjusted voxel-to-world transformation for the spatial
- dimensions of the resampled data.
+ dimensions of the resampled data, or ``None`` if a ``matrix``
+ was provided.
"""
if sliceobj is None: sliceobj = slice(None)
if dtype is None: dtype = image.dtype
+ if origin is None: origin = 'centre'
+ if mode is None: mode = 'nearest'
if origin == 'center': origin = 'centre'
+ ownMatrix = matrix is None
+
if origin not in ('centre', 'corner'):
raise ValueError('Invalid value for origin: {}'.format(origin))
@@ -216,7 +240,15 @@ def resample(image,
matrix[:3, :3] = rotmat
matrix[:3, -1] = offsets
- matrix = transform.concat(image.voxToWorldMat, matrix)
+ # We can only adjust the v2w affine if
+ # the input space and resampling space
+ # are aligned (e.g. if we're just
+ # resampling to different dimensions).
+ # We can't assume this when a custom
+ # matrix is specified, so we just give
+ # up and return None.
+ if ownMatrix: matrix = affine.concat(image.voxToWorldMat, matrix)
+ else: matrix = None
return data, matrix
@@ -239,7 +271,7 @@ def applySmoothing(data, matrix, newShape):
:returns: A smoothed copy of ``data``.
"""
- ratio = transform.decompose(matrix[:3, :3])[0]
+ ratio = affine.decompose(matrix[:3, :3])[0]
if len(newShape) > 3:
ratio = np.concatenate((
diff --git a/fsl/utils/image/roi.py b/fsl/utils/image/roi.py
index 84aee2434a30a6d16d19c3729466ca10845c4185..0b26d79d82fd93c3f2920da30b79a73ef4cb981f 100644
--- a/fsl/utils/image/roi.py
+++ b/fsl/utils/image/roi.py
@@ -11,8 +11,8 @@ a region-of-interest from, or expand the field-of-view of, an :class:`.Image`.
import numpy as np
-import fsl.data.image as fslimage
-import fsl.utils.transform as transform
+import fsl.data.image as fslimage
+import fsl.transform.affine as affine
def _normaliseBounds(shape, bounds):
@@ -97,8 +97,8 @@ def roi(image, bounds):
# each spatial dimension
oldaff = image.voxToWorldMat
offset = [lo for lo, hi in bounds[:3]]
- offset = transform.scaleOffsetXform([1, 1, 1], offset)
- newaff = transform.concat(oldaff, offset)
+ offset = affine.scaleOffsetXform([1, 1, 1], offset)
+ newaff = affine.concat(oldaff, offset)
return fslimage.Image(newdata,
xform=newaff,
diff --git a/fsl/utils/parse_data.py b/fsl/utils/parse_data.py
index 20e51d75218a87a850c3ce291ba75ef22a39c4ad..1c16e97dee1d8e3dc4b67618610461701d182a93 100644
--- a/fsl/utils/parse_data.py
+++ b/fsl/utils/parse_data.py
@@ -8,8 +8,9 @@
Argparse is the built-in python library for resolving command line arguments.
-The functions in this module can be passed on to the ``type`` argument in the ``ArgumentParser.add_command`` method
-to interpret command line arguments as neuroimageing objects (.e.g, NIFTI image files)
+The functions in this module can be passed on to the ``type`` argument in the
+``ArgumentParser.add_command`` method to interpret command line arguments as
+neuroimaging objects (.e.g, NIFTI image files)
.. autosummary::
@@ -27,18 +28,21 @@ from fsl.utils import path
import argparse
-def Image(filename):
+def Image(filename, *args, **kwargs):
"""
Reads in an image from a NIFTI or Analyze file.
:arg filename: filename provided by the user
:return: fsl.data.image.Image object
+
+ All other arguments are passed through to the :class:`.Image` upon
+ creation.
"""
try:
full_filename = image.addExt(filename)
except path.PathError as e:
raise argparse.ArgumentTypeError(*e.args)
- return image.Image(full_filename)
+ return image.Image(full_filename, *args, **kwargs)
def ImageOut(basename):
diff --git a/fsl/utils/transform.py b/fsl/utils/transform.py
index 8c9be2c440a139851371d9f1e86454915b11fb7a..b3d51afb21f0e5a5a861458d4fde9e2e62daebc2 100644
--- a/fsl/utils/transform.py
+++ b/fsl/utils/transform.py
@@ -1,647 +1,21 @@
#!/usr/bin/env python
#
-# transform.py - Functions for working with affine transformation matrices.
+# transforms.py - Deprecated
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
-"""This module provides functions related to 3D image transformations and
-spaces. The following functions are provided:
-
-.. autosummary::
- :nosignatures:
-
- transform
- transformNormal
- scaleOffsetXform
- invert
- concat
- compose
- decompose
- rotMatToAffine
- rotMatToAxisAngles
- axisAnglesToRotMat
- axisBounds
- flirtMatrixToSform
- sformToFlirtMatrix
- rmsdev
-
-And a few more functions are provided for working with vectors:
-
-.. autosummary::
- :nosignatures:
-
- veclength
- normalise
+"""The ``fsl.utils.transform`` module is deprecated - use the
+:mod:`fsl.transform` module instead.
"""
-import numpy as np
-import numpy.linalg as linalg
-import collections.abc as abc
-
-
-def invert(x):
- """Inverts the given matrix using ``numpy.linalg.inv``. """
- return linalg.inv(x)
-
-
-def concat(*xforms):
- """Combines the given matrices (returns the dot product)."""
-
- result = xforms[0]
-
- for i in range(1, len(xforms)):
- result = np.dot(result, xforms[i])
-
- return result
-
-
-def veclength(vec):
- """Returns the length of the given vector(s).
-
- Multiple vectors may be passed in, with a shape of ``(n, 3)``.
- """
- vec = np.array(vec, copy=False).reshape(-1, 3)
- return np.sqrt(np.einsum('ij,ij->i', vec, vec))
-
-
-def normalise(vec):
- """Normalises the given vector(s) to unit length.
-
- Multiple vectors may be passed in, with a shape of ``(n, 3)``.
- """
- vec = np.array(vec, copy=False).reshape(-1, 3)
- n = (vec.T / veclength(vec)).T
-
- if n.size == 3:
- n = n[0]
-
- return n
-
-
-def scaleOffsetXform(scales, offsets):
- """Creates and returns an affine transformation matrix which encodes
- the specified scale(s) and offset(s).
-
-
- :arg scales: A tuple of up to three values specifying the scale factors
- for each dimension. If less than length 3, is padded with
- ``1.0``.
-
- :arg offsets: A tuple of up to three values specifying the offsets for
- each dimension. If less than length 3, is padded with
- ``0.0``.
-
- :returns: A ``numpy.float32`` array of size :math:`4 \\times 4`.
- """
-
- oktypes = (abc.Sequence, np.ndarray)
-
- if not isinstance(scales, oktypes): scales = [scales]
- if not isinstance(offsets, oktypes): offsets = [offsets]
- if not isinstance(scales, list): scales = list(scales)
- if not isinstance(offsets, list): offsets = list(offsets)
-
- lens = len(scales)
- leno = len(offsets)
-
- if lens < 3: scales = scales + [1.0] * (3 - lens)
- if leno < 3: offsets = offsets + [0.0] * (3 - leno)
-
- xform = np.eye(4, dtype=np.float64)
-
- xform[0, 0] = scales[0]
- xform[1, 1] = scales[1]
- xform[2, 2] = scales[2]
-
- xform[0, 3] = offsets[0]
- xform[1, 3] = offsets[1]
- xform[2, 3] = offsets[2]
-
- return xform
-
-
-def compose(scales, offsets, rotations, origin=None):
- """Compose a transformation matrix out of the given scales, offsets
- and axis rotations.
-
- :arg scales: Sequence of three scale values.
-
- :arg offsets: Sequence of three offset values.
-
- :arg rotations: Sequence of three rotation values, in radians, or
- a rotation matrix of shape ``(3, 3)``.
-
- :arg origin: Origin of rotation - must be scaled by the ``scales``.
- If not provided, the rotation origin is ``(0, 0, 0)``.
- """
-
- preRotate = np.eye(4)
- postRotate = np.eye(4)
-
- rotations = np.array(rotations)
-
- if rotations.shape == (3,):
- rotations = axisAnglesToRotMat(*rotations)
-
- if origin is not None:
- preRotate[ 0, 3] = -origin[0]
- preRotate[ 1, 3] = -origin[1]
- preRotate[ 2, 3] = -origin[2]
- postRotate[0, 3] = origin[0]
- postRotate[1, 3] = origin[1]
- postRotate[2, 3] = origin[2]
-
- scale = np.eye(4, dtype=np.float64)
- offset = np.eye(4, dtype=np.float64)
- rotate = np.eye(4, dtype=np.float64)
-
- scale[ 0, 0] = scales[ 0]
- scale[ 1, 1] = scales[ 1]
- scale[ 2, 2] = scales[ 2]
- offset[ 0, 3] = offsets[0]
- offset[ 1, 3] = offsets[1]
- offset[ 2, 3] = offsets[2]
-
- rotate[:3, :3] = rotations
-
- return concat(offset, postRotate, rotate, preRotate, scale)
-
-
-def decompose(xform, angles=True):
- """Decomposes the given transformation matrix into separate offsets,
- scales, and rotations, according to the algorithm described in:
-
- Spencer W. Thomas, Decomposing a matrix into simple transformations, pp
- 320-323 in *Graphics Gems II*, James Arvo (editor), Academic Press, 1991,
- ISBN: 0120644819.
-
- It is assumed that the given transform has no perspective components. Any
- shears in the affine are discarded.
-
- :arg xform: A ``(3, 3)`` or ``(4, 4)`` affine transformation matrix.
-
- :arg angles: If ``True`` (the default), the rotations are returned
- as axis-angles, in radians. Otherwise, the rotation matrix
- is returned.
-
- :returns: The following:
-
- - A sequence of three scales
- - A sequence of three translations (all ``0`` if ``xform``
- was a ``(3, 3)`` matrix)
- - A sequence of three rotations, in radians. Or, if
- ``angles is False``, a rotation matrix.
- """
-
- # The inline comments in the code below are taken verbatim from
- # the referenced article, [except for notes in square brackets].
-
- # The next step is to extract the translations. This is trivial;
- # we find t_x = M_{4,1}, t_y = M_{4,2}, and t_z = M_{4,3}. At this
- # point we are left with a 3*3 matrix M' = M_{1..3,1..3}.
- xform = xform.T
-
- if xform.shape == (4, 4):
- translations = xform[ 3, :3]
- xform = xform[:3, :3]
- else:
- translations = np.array([0, 0, 0])
-
- M1 = xform[0]
- M2 = xform[1]
- M3 = xform[2]
-
- # The process of finding the scaling factors and shear parameters
- # is interleaved. First, find s_x = |M'_1|.
- sx = np.sqrt(np.dot(M1, M1))
-
- # Then, compute an initial value for the xy shear factor,
- # s_xy = M'_1 * M'_2. (this is too large by the y scaling factor).
- sxy = np.dot(M1, M2)
-
- # The second row of the matrix is made orthogonal to the first by
- # setting M'_2 = M'_2 - s_xy * M'_1.
- M2 = M2 - sxy * M1
-
- # Then the y scaling factor, s_y, is the length of the modified
- # second row.
- sy = np.sqrt(np.dot(M2, M2))
-
- # The second row is normalized, and s_xy is divided by s_y to
- # get its final value.
- M2 = M2 / sy
- sxy = sxy / sy
-
- # The xz and yz shear factors are computed as in the preceding,
- sxz = np.dot(M1, M3)
- syz = np.dot(M2, M3)
-
- # the third row is made orthogonal to the first two rows,
- M3 = M3 - sxz * M1 - syz * M2
-
- # the z scaling factor is computed,
- sz = np.sqrt(np.dot(M3, M3))
-
- # the third row is normalized, and the xz and yz shear factors are
- # rescaled.
- M3 = M3 / sz
- sxz = sxz / sz
- syz = sxz / sz
-
- # The resulting matrix now is a pure rotation matrix, except that it
- # might still include a scale factor of -1. If the determinant of the
- # matrix is -1, negate the matrix and all three scaling factors. Call
- # the resulting matrix R.
- #
- # [We do things different here - if the rotation matrix has negative
- # determinant, the flip is encoded in the x scaling factor.]
- R = np.array([M1, M2, M3])
- if linalg.det(R) < 0:
- R[0] = -R[0]
- sx = -sx
-
- # Finally, we need to decompose the rotation matrix into a sequence
- # of rotations about the x, y, and z axes. [This is done in the
- # rotMatToAxisAngles function]
- if angles: rotations = rotMatToAxisAngles(R.T)
- else: rotations = R.T
-
- return np.array([sx, sy, sz]), translations, rotations
-
-
-def rotMatToAffine(rotmat, origin=None):
- """Convenience function which encodes the given ``(3, 3)`` rotation
- matrix into a ``(4, 4)`` affine.
- """
- return compose([1, 1, 1], [0, 0, 0], rotmat, origin)
-
-
-def rotMatToAxisAngles(rotmat):
- """Given a ``(3, 3)`` rotation matrix, decomposes the rotations into
- an angle in radians about each axis.
- """
-
- yrot = np.sqrt(rotmat[0, 0] ** 2 + rotmat[1, 0] ** 2)
-
- if np.isclose(yrot, 0):
- xrot = np.arctan2(-rotmat[1, 2], rotmat[1, 1])
- yrot = np.arctan2(-rotmat[2, 0], yrot)
- zrot = 0
- else:
- xrot = np.arctan2( rotmat[2, 1], rotmat[2, 2])
- yrot = np.arctan2(-rotmat[2, 0], yrot)
- zrot = np.arctan2( rotmat[1, 0], rotmat[0, 0])
-
- return [xrot, yrot, zrot]
-
-
-def axisAnglesToRotMat(xrot, yrot, zrot):
- """Constructs a ``(3, 3)`` rotation matrix from the given angles, which
- must be specified in radians.
- """
-
- xmat = np.eye(3)
- ymat = np.eye(3)
- zmat = np.eye(3)
-
- xmat[1, 1] = np.cos(xrot)
- xmat[1, 2] = -np.sin(xrot)
- xmat[2, 1] = np.sin(xrot)
- xmat[2, 2] = np.cos(xrot)
-
- ymat[0, 0] = np.cos(yrot)
- ymat[0, 2] = np.sin(yrot)
- ymat[2, 0] = -np.sin(yrot)
- ymat[2, 2] = np.cos(yrot)
-
- zmat[0, 0] = np.cos(zrot)
- zmat[0, 1] = -np.sin(zrot)
- zmat[1, 0] = np.sin(zrot)
- zmat[1, 1] = np.cos(zrot)
-
- return concat(zmat, ymat, xmat)
-
-
-def axisBounds(shape,
- xform,
- axes=None,
- origin='centre',
- boundary='high',
- offset=1e-4):
- """Returns the ``(lo, hi)`` bounds of the specified axis/axes in the
- world coordinate system defined by ``xform``.
-
- If the ``origin`` parameter is set to ``centre`` (the default),
- this function assumes that voxel indices correspond to the voxel
- centre. For example, the voxel at ``(4, 5, 6)`` covers the space:
-
- ``[3.5 - 4.5, 4.5 - 5.5, 5.5 - 6.5]``
-
- So the bounds of the specified shape extends from the corner at
-
- ``(-0.5, -0.5, -0.5)``
-
- to the corner at
-
- ``(shape[0] - 0.5, shape[1] - 0.5, shape[1] - 0.5)``
-
- If the ``origin`` parameter is set to ``corner``, this function
- assumes that voxel indices correspond to the voxel corner. In this
- case, a voxel at ``(4, 5, 6)`` covers the space:
-
- ``[4 - 5, 5 - 6, 6 - 7]``
-
- So the bounds of the specified shape extends from the corner at
-
- ``(0, 0, 0)``
-
- to the corner at
-
- ``(shape[0], shape[1], shape[1])``.
-
-
- If the ``boundary`` parameter is set to ``high``, the high voxel bounds
- are reduced by a small amount (specified by the ``offset`` parameter)
- before they are transformed to the world coordinate system. If
- ``boundary`` is set to ``low``, the low bounds are increased by a small
- amount. The ``boundary`` parameter can also be set to ``'both'``, or
- ``None``. This option is provided so that you can ensure that the
- resulting bounds will always be contained within the image space.
-
- :arg shape: The ``(x, y, z)`` shape of the data.
-
- :arg xform: Transformation matrix which transforms voxel coordinates
- to the world coordinate system.
-
- :arg axes: The world coordinate system axis bounds to calculate.
-
- :arg origin: Either ``'centre'`` (the default) or ``'corner'``.
-
- :arg boundary: Either ``'high'`` (the default), ``'low'``, ''`both'``,
- or ``None``.
-
- :arg offset: Amount by which the boundary voxel coordinates should be
- offset. Defaults to ``1e-4``.
-
- :returns: A tuple containing the ``(low, high)`` bounds for each
- requested world coordinate system axis.
- """
-
- origin = origin.lower()
-
- # lousy US spelling
- if origin == 'center':
- origin = 'centre'
-
- if origin not in ('centre', 'corner'):
- raise ValueError('Invalid origin value: {}'.format(origin))
- if boundary not in ('low', 'high', 'both', None):
- raise ValueError('Invalid boundary value: {}'.format(boundary))
-
- scalar = False
-
- if axes is None:
- axes = [0, 1, 2]
-
- elif not isinstance(axes, abc.Iterable):
- scalar = True
- axes = [axes]
-
- x, y, z = shape[:3]
-
- points = np.zeros((8, 3), dtype=np.float32)
-
- if origin == 'centre':
- x0 = -0.5
- y0 = -0.5
- z0 = -0.5
- x -= 0.5
- y -= 0.5
- z -= 0.5
- else:
- x0 = 0
- y0 = 0
- z0 = 0
-
- if boundary in ('low', 'both'):
- x0 += offset
- y0 += offset
- z0 += offset
-
- if boundary in ('high', 'both'):
- x -= offset
- y -= offset
- z -= offset
-
- points[0, :] = [x0, y0, z0]
- points[1, :] = [x0, y0, z]
- points[2, :] = [x0, y, z0]
- points[3, :] = [x0, y, z]
- points[4, :] = [x, y0, z0]
- points[5, :] = [x, y0, z]
- points[6, :] = [x, y, z0]
- points[7, :] = [x, y, z]
-
- tx = transform(points, xform)
-
- lo = tx[:, axes].min(axis=0)
- hi = tx[:, axes].max(axis=0)
-
- if scalar: return (lo[0], hi[0])
- else: return (lo, hi)
-
-
-def transform(p, xform, axes=None, vector=False):
- """Transforms the given set of points ``p`` according to the given affine
- transformation ``xform``.
-
-
- :arg p: A sequence or array of points of shape :math:`N \\times 3`.
-
- :arg xform: A ``(4, 4)`` affine transformation matrix with which to
- transform the points in ``p``.
-
- :arg axes: If you are only interested in one or two axes, and the source
- axes are orthogonal to the target axes (see the note below),
- you may pass in a 1D, ``N*1``, or ``N*2`` array as ``p``, and
- use this argument to specify which axis/axes that the data in
- ``p`` correspond to.
-
- :arg vector: Defaults to ``False``. If ``True``, the points are treated
- as vectors - the translation component of the transformation
- is not applied. If you set this flag, you pass in a ``(3, 3)``
- transformation matrix.
-
- :returns: The points in ``p``, transformed by ``xform``, as a ``numpy``
- array with the same data type as the input.
-
-
- .. note:: The ``axes`` argument should only be used if the source
- coordinate system (the points in ``p``) axes are orthogonal
- to the target coordinate system (defined by the ``xform``).
-
- In other words, you can only use the ``axes`` argument if
- the ``xform`` matrix consists solely of translations and
- scalings.
- """
-
- p = _fillPoints(p, axes)
- t = np.dot(xform[:3, :3], p.T).T
-
- if not vector:
- t = t + xform[:3, 3]
-
- if axes is not None:
- t = t[:, axes]
-
- if t.size == 1: return t[0]
- else: return t
-
-
-def transformNormal(p, xform, axes=None):
- """Transforms the given point(s), under the assumption that they
- are normal vectors. In this case, the points are transformed by
- ``invert(xform[:3, :3]).T``.
- """
- return transform(p, invert(xform[:3, :3]).T, axes, vector=True)
-
-
-def _fillPoints(p, axes):
- """Used by the :func:`transform` function. Turns the given array p into
- a ``N*3`` array of ``x,y,z`` coordinates. The array p may be a 1D array,
- or an ``N*2`` or ``N*3`` array.
- """
-
- if not isinstance(p, abc.Iterable): p = [p]
-
- p = np.array(p)
-
- if axes is None: return p
-
- if not isinstance(axes, abc.Iterable): axes = [axes]
-
- if p.ndim == 1:
- p = p.reshape((len(p), 1))
-
- if p.ndim != 2:
- raise ValueError('Points array must be either one or two '
- 'dimensions')
-
- if len(axes) != p.shape[1]:
- raise ValueError('Points array shape does not match specified '
- 'number of axes')
-
- newp = np.zeros((len(p), 3), dtype=p.dtype)
-
- for i, ax in enumerate(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.voxToScaledVoxMat
- refInvScaledVoxelMat = refImage.scaledVoxToVoxMat
- refVoxToWorldMat = refImage.voxToWorldMat
-
- return concat(refVoxToWorldMat,
- refInvScaledVoxelMat,
- flirtMat,
- srcScaledVoxelMat)
-
-
-def sformToFlirtMatrix(srcImage, refImage, srcXform=None):
- """Under the assumption that the given ``srcImage`` and ``refImage`` share a
- common world coordinate system (defined by their
- :attr:`.Nifti.voxToWorldMat` attributes), this function will calculate and
- return a transformation matrix from the ``srcImage`` scaled voxel
- coordinate system to the ``refImage`` scaled voxel coordinate system, that
- can be saved to disk and used with FLIRT, to resample the source image to
- the reference image.
-
- :arg srcImage: Source :class:`.Image`
- :arg refImage: Reference :class:`.Image`
- :arg srcXform: Optionally used in place of the ``srcImage``
- :attr:`.Nifti.voxToWorldMat`
- """
-
- srcScaledVoxToVoxMat = srcImage.scaledVoxToVoxMat
- srcVoxToWorldMat = srcImage.voxToWorldMat
- refWorldToVoxMat = refImage.worldToVoxMat
- refVoxToScaledVoxMat = refImage.voxToScaledVoxMat
-
- if srcXform is not None:
- srcVoxToWorldMat = srcXform
-
- return concat(refVoxToScaledVoxMat,
- refWorldToVoxMat,
- srcVoxToWorldMat,
- srcScaledVoxToVoxMat)
-
-
-def rmsdev(T1, T2, R=None, xc=None):
- """Calculates the RMS deviation of the given affine transforms ``T1`` and
- ``T2``. This can be used as a measure of the 'distance' between two
- affines.
-
- The ``T1`` and ``T2`` arguments may be either full ``(4, 4)`` affines, or
- ``(3, 3)`` rotation matrices.
-
- See FMRIB technical report TR99MJ1, available at:
-
- https://www.fmrib.ox.ac.uk/datasets/techrep/
-
- :arg T1: First affine
- :arg T2: Second affine
- :arg R: Sphere radius
- :arg xc: Sphere centre
- :returns: The RMS deviation between ``T1`` and ``T2``.
- """
-
- if R is None:
- R = 1
-
- if xc is None:
- xc = np.zeros(3)
-
- # rotations only
- if T1.shape == (3, 3):
- M = np.dot(T2, invert(T1)) - np.eye(3)
- A = M[:3, :3]
- t = np.zeros(3)
-
- # full affine
- else:
- M = np.dot(T2, invert(T1)) - np.eye(4)
- A = M[:3, :3]
- t = M[:3, 3]
- Axc = np.dot(A, xc)
+import fsl.utils.deprecated as deprecated
+from fsl.transform.affine import * # noqa
+from fsl.transform.flirt import (flirtMatrixToSform, # noqa
+ sformToFlirtMatrix)
- erms = np.dot((t + Axc).T, t + Axc)
- erms = 0.2 * R ** 2 * np.dot(A.T, A).trace() + erms
- erms = np.sqrt(erms)
- return erms
+deprecated.warn('fsl.utils.transform',
+ vin='2.4.0',
+ rin='3.0.0',
+ msg='Use the fsl.transform module instead')
diff --git a/fsl/wrappers/wrapperutils.py b/fsl/wrappers/wrapperutils.py
index 373bcf08b1a0c53af416a074cc4c439cf534f792..c89afb2aadb002e1e6838235578ce124ce4caf41 100644
--- a/fsl/wrappers/wrapperutils.py
+++ b/fsl/wrappers/wrapperutils.py
@@ -446,9 +446,9 @@ class _FileOrThing(object):
Functions decorated with a ``_FileOrThing`` decorator will always return a
``dict``-like object, where the function's actual return value is
- accessible via an attribute called `output`. All output arguments with a
+ accessible via an attribute called ``output``. All output arguments with a
value of ``LOAD`` will be present as dictionary entries, with the keyword
- argument names used as keys. Any ``LOAD``ed output arguments which were not
+ argument names used as keys. Any ``LOAD`` output arguments which were not
generated by the function will not be present in the dictionary.
diff --git a/requirements.txt b/requirements.txt
index ef041ad9cd02c2d53af3875cb4cecf7d2d077d5f..f0abcb98be82381e843f9d303367b66ec73f898b 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,5 @@
-six==1.*
+h5py>=2.9
+nibabel>=2.3.1
numpy==1.*
scipy>=0.18
-nibabel>=2.3.1
+six==1.*
diff --git a/setup.py b/setup.py
index 8693d37c5eef988bc747e9b81c42afc5e852b81b..6b5ca7dd9d824851b44e1b4db37a91bb48abbf54 100644
--- a/setup.py
+++ b/setup.py
@@ -127,6 +127,8 @@ setup(
'atlasquery = fsl.scripts.atlasq:atlasquery_emulation',
'fsl_ents = fsl.scripts.fsl_ents:main',
'resample_image = fsl.scripts.resample_image:main',
+ 'fsl_convert_x5 = fsl.scripts.fsl_convert_x5:main',
+ 'fsl_apply_x5 = fsl.scripts.fsl_apply_x5:main'
]
}
)
diff --git a/tests/test_atlases.py b/tests/test_atlases.py
index b2dfec113a00349dea0601f781f565561fa18b2c..99e656e10f3a0e138645e92d2c2aa53c66d97a73 100644
--- a/tests/test_atlases.py
+++ b/tests/test_atlases.py
@@ -17,10 +17,10 @@ import mock
import pytest
import tests
-import fsl.utils.transform as transform
import fsl.utils.image.resample as resample
import fsl.data.atlases as atlases
import fsl.data.image as fslimage
+import fsl.transform.affine as affine
datadir = op.join(op.dirname(__file__), 'testdata')
@@ -326,7 +326,7 @@ def test_prepareMask():
np.random.random(list(ashape) + [2]))
wrongspace = fslimage.Image(
np.random.random((20, 20, 20)),
- xform=transform.concat(atlas.voxToWorldMat, np.diag([2, 2, 2, 1])))
+ xform=affine.concat(atlas.voxToWorldMat, np.diag([2, 2, 2, 1])))
with pytest.raises(atlases.MaskError):
atlas.prepareMask(wrongdims)
diff --git a/tests/test_atlases_query.py b/tests/test_atlases_query.py
index 693366a2e65000685196f4c7067304cd8a59fad9..051f9cd6b7bf5a22234c31ebd3273f6d924033db 100644
--- a/tests/test_atlases_query.py
+++ b/tests/test_atlases_query.py
@@ -11,11 +11,11 @@ import numpy as np
import pytest
-import fsl.data.atlases as fslatlases
-import fsl.data.image as fslimage
-import fsl.utils.transform as transform
+import fsl.data.atlases as fslatlases
+import fsl.data.image as fslimage
+import fsl.transform.affine as affine
import fsl.utils.image.resample as resample
-import fsl.utils.cache as cache
+import fsl.utils.cache as cache
from . import (testdir, make_random_mask)
@@ -155,7 +155,7 @@ def _gen_coord_voxel_query(atlas, qtype, qin, **kwargs):
dlo = (0, 0, 0)
dhi = atlas.shape
else:
- dlo, dhi = transform.axisBounds(atlas.shape, atlas.voxToWorldMat)
+ dlo, dhi = affine.axisBounds(atlas.shape, atlas.voxToWorldMat)
dlen = [hi - lo for lo, hi in zip(dlo, dhi)]
@@ -190,7 +190,7 @@ def _gen_coord_voxel_query(atlas, qtype, qin, **kwargs):
coords = np.array(coords, dtype=dtype)
if not voxel:
- coords = transform.transform(coords, atlas.voxToWorldMat)
+ coords = affine.transform(coords, atlas.voxToWorldMat)
return tuple([dtype(c) for c in coords])
@@ -200,7 +200,7 @@ def _eval_coord_voxel_query(atlas, query, qtype, qin):
voxel = qtype == 'voxel'
if voxel: vx, vy, vz = query
- else: vx, vy, vz = transform.transform(query, atlas.worldToVoxMat)
+ else: vx, vy, vz = affine.transform(query, atlas.worldToVoxMat)
vx, vy, vz = [int(round(v)) for v in [vx, vy, vz]]
diff --git a/tests/test_deprecated.py b/tests/test_deprecated.py
new file mode 100644
index 0000000000000000000000000000000000000000..307818960cbdac0f67b7c1f11c8de0c44bc9ebf1
--- /dev/null
+++ b/tests/test_deprecated.py
@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+
+
+import warnings
+import pytest
+
+import fsl.utils.deprecated as deprecated
+
+
+# the line number of the warning is hard coded in
+# the unit tests below. Don't change the line number!
+def emit_warning():
+ deprecated.warn('blag', vin='1.0.0', rin='2.0.0', msg='yo')
+
+WARNING_LINE_NUMBER = 13
+
+
+@deprecated.deprecated(vin='1.0.0', rin='2.0.0', msg='yo')
+def depfunc():
+ pass
+
+def call_dep_func():
+ depfunc()
+
+DEPRECATED_LINE_NUMBER = 23
+
+
+def _check_warning(w, name, lineno):
+ assert issubclass(w.category, DeprecationWarning)
+ assert '{} is deprecated'.format(name) in str(w.message)
+ assert 'test_deprecated.py' in str(w.filename)
+ assert w.lineno == lineno
+
+def test_warn():
+ deprecated.resetWarningCache()
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ emit_warning()
+ assert len(w) == 1
+ _check_warning(w[0], 'blag', WARNING_LINE_NUMBER)
+
+ # warning should only be emitted once
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ emit_warning()
+ assert len(w) == 0
+
+
+def test_deprecated():
+ deprecated.resetWarningCache()
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ call_dep_func()
+ assert len(w) == 1
+ _check_warning(w[0], 'depfunc', DEPRECATED_LINE_NUMBER)
+
+ # warning should only be emitted once
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ call_dep_func()
+ assert len(w) == 0
diff --git a/tests/test_image.py b/tests/test_image.py
index 5b1ba4c7ff15c2e7573fd2f577108fc7399fc23e..58d9fc625d5e107385c7280088b1c0c3d595af91 100644
--- a/tests/test_image.py
+++ b/tests/test_image.py
@@ -19,11 +19,10 @@ import nibabel as nib
from nibabel.spatialimages import ImageFileError
-import fsl.data.constants as constants
-import fsl.data.image as fslimage
-import fsl.data.imagewrapper as imagewrapper
-import fsl.utils.path as fslpath
-import fsl.utils.transform as transform
+import fsl.data.constants as constants
+import fsl.data.image as fslimage
+import fsl.utils.path as fslpath
+import fsl.transform.affine as affine
from fsl.utils.tempdir import tempdir
@@ -743,6 +742,48 @@ def _test_Image_changeXform(imgtype, sformcode=None, qformcode=None):
image = None
+def test_Image_changeIntent_analyze(): _test_Image_changeIntent(0)
+def test_Image_changeIntent_nifti1(): _test_Image_changeIntent(1)
+def test_Image_changeIntent_nifti2(): _test_Image_changeIntent(2)
+def _test_Image_changeIntent(imgtype):
+ """Test changing the Nifti.intent attribute. """
+
+ with tempdir() as testdir:
+ imagefile = op.join(testdir, 'image')
+
+ image = make_image(imagefile, imgtype)
+ if imgtype > 0:
+ image.header.set_intent(constants.NIFTI_INTENT_NONE)
+ nib.save(image, imagefile)
+
+ notified = {}
+ def onHdr( *a): notified['header'] = True
+ def onSave(*a): notified['save'] = True
+
+ img = fslimage.Image(imagefile)
+
+ img.register('name1', onHdr, 'header')
+ img.register('name2', onSave, 'saveState')
+
+ assert img.intent == constants.NIFTI_INTENT_NONE
+ img.intent = constants.NIFTI_INTENT_BETA
+
+ if imgtype == 0: exp = constants.NIFTI_INTENT_NONE
+ else: exp = constants.NIFTI_INTENT_BETA
+
+ assert img.intent == exp
+
+ if imgtype > 0:
+ assert img .header.get_intent('code')[0] == exp
+ assert img.nibImage.header.get_intent('code')[0] == exp
+
+ assert notified.get('header', False)
+ assert notified.get('save', False)
+
+
+
+
+
def test_Image_changeData_analyze(seed): _test_Image_changeData(0)
def test_Image_changeData_nifti1(seed): _test_Image_changeData(1)
def test_Image_changeData_nifti2(seed): _test_Image_changeData(2)
@@ -1072,12 +1113,12 @@ def _test_Image_init_xform(imgtype):
with tempdir() as td:
- sform = transform.compose(np.random.random(3),
- np.random.random(3),
- np.random.random(3))
- qform = transform.compose(np.random.random(3),
- np.random.random(3),
- np.random.random(3))
+ sform = affine.compose(np.random.random(3),
+ np.random.random(3),
+ np.random.random(3))
+ qform = affine.compose(np.random.random(3),
+ np.random.random(3),
+ np.random.random(3))
sform_code = 3
qform_code = 4
@@ -1128,9 +1169,9 @@ def _test_Image_init_xform(imgtype):
# to the xform. and its
# s/q form codes the same
# as what is in the header
- rxform = transform.compose(np.random.random(3),
- np.random.random(3),
- np.random.random(3))
+ rxform = affine.compose(np.random.random(3),
+ np.random.random(3),
+ np.random.random(3))
fimg = fslimage.Image(img.get_data(),
header=img.header,
xform=rxform)
@@ -1173,3 +1214,158 @@ def test_rgb_image():
assert img.nvals == 3
assert img.dataRange == (0, 255)
+
+
+def test_determineShape():
+ class MockHeader(object):
+ def __init__(self, shape, zooms):
+ self.shape = shape
+ self.zooms = zooms
+ def __getitem__(self, key):
+ return [len(self.zooms)] + self.zooms
+ def get_data_shape(self):
+ return self.shape
+ def get_zooms(self):
+ return self.zooms
+
+ # inshape, inzooms, outshape, outzooms)
+ tests = [
+ ([10], [2, 2, 2], [10, 1, 1], [2, 2, 2]),
+ ([10], [2], [10, 1, 1], [2, 1, 1]),
+ ([10], [2, 2, 2], [10, 1, 1], [2, 2, 2]),
+ ([10, 10], [2, 2], [10, 10, 1], [2, 2, 1]),
+ ([10, 10], [2, 2, 2], [10, 10, 1], [2, 2, 2]),
+ ([10, 10, 10], [2, 2, 2], [10, 10, 10], [2, 2, 2]),
+
+ ([10, 10, 10, 10], [2, 2, 2, 2],
+ [10, 10, 10, 10], [2, 2, 2, 2]),
+ ([10, 10, 10, 10, 10], [2, 2, 2, 2, 2],
+ [10, 10, 10, 10, 10], [2, 2, 2, 2, 2]),
+ ]
+
+ for inshape, inzooms, outshape, outzooms in tests:
+
+ hdr = MockHeader(inshape, inzooms)
+ origshape, gotshape, gotzooms = fslimage.Nifti.determineShape(hdr)
+
+ assert origshape == inshape
+ assert gotshape == outshape
+ assert gotzooms == outzooms
+
+
+def test_determineAffine():
+
+ # sformcode, qformcode, intent, expaff
+ tests = [
+ (constants.NIFTI_XFORM_ALIGNED_ANAT,
+ constants.NIFTI_XFORM_ALIGNED_ANAT,
+ constants.NIFTI_INTENT_NONE,
+ 'sform'),
+ (constants.NIFTI_XFORM_ALIGNED_ANAT,
+ constants.NIFTI_XFORM_UNKNOWN,
+ constants.NIFTI_INTENT_NONE,
+ 'sform'),
+ (constants.NIFTI_XFORM_UNKNOWN,
+ constants.NIFTI_XFORM_ALIGNED_ANAT,
+ constants.NIFTI_INTENT_NONE,
+ 'qform'),
+ (constants.NIFTI_XFORM_ALIGNED_ANAT,
+ constants.NIFTI_XFORM_ALIGNED_ANAT,
+ constants.FSL_FNIRT_DISPLACEMENT_FIELD,
+ 'sform'),
+ (constants.NIFTI_XFORM_ALIGNED_ANAT,
+ constants.NIFTI_XFORM_ALIGNED_ANAT,
+ constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
+ 'scaling'),
+ (constants.NIFTI_XFORM_UNKNOWN,
+ constants.NIFTI_XFORM_UNKNOWN,
+ constants.NIFTI_INTENT_NONE,
+ 'scaling'),
+ ]
+
+ for sformcode, qformcode, intent, exp in tests:
+
+ sform = affine.compose(np.random.random(3),
+ np.random.random(3),
+ np.random.random(3))
+ qform = affine.compose(np.random.random(3),
+ np.random.random(3),
+ np.random.random(3))
+ pixdims = np.random.randint(1, 10, 3)
+
+ hdr = nib.Nifti1Header()
+ hdr.set_data_shape((10, 10, 10))
+ hdr.set_sform(sform, sformcode)
+ hdr.set_qform(qform, qformcode)
+ hdr.set_intent(intent)
+ hdr.set_zooms(pixdims)
+
+ # the randomly generated qform might
+ # not be fully representable, so let
+ # nibabel fix it for us
+ sform = hdr.get_sform()
+ qform = hdr.get_qform()
+
+ got = fslimage.Nifti.determineAffine(hdr)
+
+ if exp == 'sform': exp = sform
+ elif exp == 'qform': exp = qform
+ elif exp == 'scaling': exp = affine.scaleOffsetXform(pixdims, 0)
+
+ assert np.all(np.isclose(got, exp))
+
+
+def test_generateAffines():
+
+ v2w = affine.compose(np.random.random(3),
+ np.random.random(3),
+ np.random.random(3))
+ shape = (10, 10, 10)
+ pixdim = (1, 1, 1)
+
+ got, isneuro = fslimage.Nifti.generateAffines(v2w, shape, pixdim)
+
+ w2v = npla.inv(v2w)
+
+ assert isneuro == (npla.det(v2w) > 0)
+
+ if not isneuro:
+ v2f = np.eye(4)
+ f2v = np.eye(4)
+ f2w = v2w
+ w2f = w2v
+ else:
+ v2f = affine.scaleOffsetXform([-1, 1, 1], [9, 0, 0])
+ f2v = npla.inv(v2f)
+ f2w = affine.concat(v2w, f2v)
+ w2f = affine.concat(v2f, w2v)
+
+ assert np.all(np.isclose(v2w, got['voxel', 'world']))
+ assert np.all(np.isclose(w2v, got['world', 'voxel']))
+ assert np.all(np.isclose(v2f, got['voxel', 'fsl']))
+ assert np.all(np.isclose(f2v, got['fsl', 'voxel']))
+ assert np.all(np.isclose(f2w, got['fsl' , 'world']))
+ assert np.all(np.isclose(w2f, got['world', 'fsl']))
+
+
+def test_identifyAffine():
+
+ identify = fslimage.Nifti.identifyAffine
+
+ assert identify(None, None, 'ho', 'hum') == ('ho', 'hum')
+
+ xform = affine.compose(0.1 + 5 * np.random.random(3),
+ -10 + 20 * np.random.random(3),
+ -np.pi / 2 + np.pi * np.random.random(3))
+
+ img = fslimage.Image(make_random_image(None, xform=xform))
+
+ for from_, to in it.permutations(('voxel', 'fsl', 'world'), 2):
+ assert identify(img, img.getAffine(from_, to)) == (from_, to)
+
+ assert identify(img, img.getAffine('voxel', 'world'), from_='voxel') == ('voxel', 'world')
+ assert identify(img, img.getAffine('voxel', 'world'), to='world') == ('voxel', 'world')
+
+ rubbish = np.random.random((4, 4))
+ with pytest.raises(ValueError):
+ identify(img, rubbish)
diff --git a/tests/test_image_resample.py b/tests/test_image_resample.py
index d9f474bf45330e5691535a545faa7ddcac9159a4..8fb488760ed74a14a0714c92d7815ac2a87f43f6 100644
--- a/tests/test_image_resample.py
+++ b/tests/test_image_resample.py
@@ -263,3 +263,28 @@ def test_resampleToReference3():
exp = np.array([refdata[..., 0]] * 15).transpose((1, 2, 3, 0))
resampled, xform = resample.resampleToReference(img, ref, order=0, mode='nearest')
assert np.all(resampled == exp)
+
+
+def test_resampleToReference4():
+
+ # the image and ref are out of
+ # alignment, but this affine
+ # will bring them into alignment
+ img2ref = transform.scaleOffsetXform([2, 2, 2], [10, 10, 10])
+
+ imgdata = np.random.randint(0, 65536, (5, 5, 5))
+ refdata = np.zeros((5, 5, 5))
+ img = fslimage.Image(imgdata)
+ ref = fslimage.Image(refdata, xform=img2ref)
+
+ # Without the affine, the image
+ # will be out of the FOV of the
+ # reference
+ resampled, xform = resample.resampleToReference(img, ref)
+ assert np.all(resampled == 0)
+
+ # But applying the affine will
+ # cause them to overlap
+ # perfectly in world coordinates
+ resampled, xform = resample.resampleToReference(img, ref, matrix=img2ref)
+ assert np.all(resampled == imgdata)
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 243944b2e3feb1b476792a3fa319a595a8143ae7..c5841df9d38338c601cc9e2a6ee4c5552d18f7b4 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -6,13 +6,13 @@
#
-import os.path as op
-import numpy as np
-import mock
-import pytest
+import os.path as op
+import numpy as np
+from unittest import mock
+import pytest
-import fsl.utils.transform as transform
-import fsl.data.mesh as fslmesh
+import fsl.transform.affine as affine
+import fsl.data.mesh as fslmesh
from . import tempdir
@@ -58,7 +58,7 @@ CUBE_CCW_VERTEX_NORMALS = np.zeros((8, 3))
for i in range(8):
faces = np.where(CUBE_TRIANGLES_CCW == i)[0]
CUBE_CCW_VERTEX_NORMALS[i] = CUBE_CCW_FACE_NORMALS[faces].sum(axis=0)
-CUBE_CCW_VERTEX_NORMALS = transform.normalise(CUBE_CCW_VERTEX_NORMALS)
+CUBE_CCW_VERTEX_NORMALS = affine.normalise(CUBE_CCW_VERTEX_NORMALS)
def test_mesh_create():
@@ -257,6 +257,13 @@ def test_needsFixing():
def test_trimesh_no_trimesh():
+ # Make sure trimesh and rtree
+ # are imported before messing
+ # with sys.modules, otherwise
+ # weird things can happen.
+ import trimesh
+ import rtree
+
mods = ['trimesh', 'rtree']
for mod in mods:
diff --git a/tests/test_mghimage.py b/tests/test_mghimage.py
index 2496ac4eb80cff687f6d03d522281dcd777eb45f..f34f635c17759eb9803d34c0405d4a2d93e4a9eb 100644
--- a/tests/test_mghimage.py
+++ b/tests/test_mghimage.py
@@ -12,10 +12,10 @@ import shutil
import numpy as np
import nibabel as nib
-import fsl.utils.tempdir as tempdir
-import fsl.utils.transform as transform
-import fsl.data.mghimage as fslmgh
-import fsl.data.image as fslimage
+import fsl.utils.tempdir as tempdir
+import fsl.transform.affine as affine
+import fsl.data.mghimage as fslmgh
+import fsl.data.image as fslimage
datadir = op.join(op.abspath(op.dirname(__file__)), 'testdata')
@@ -29,14 +29,14 @@ def test_MGHImage():
img = fslmgh.MGHImage(testfile)
nbimg = nib.load(testfile)
v2s = nbimg.header.get_vox2ras_tkr()
- w2s = transform.concat(v2s, transform.invert(nbimg.affine))
+ w2s = affine.concat(v2s, affine.invert(nbimg.affine))
assert np.all(np.isclose(img[:], nbimg.get_data()))
assert np.all(np.isclose(img.voxToWorldMat, nbimg.affine))
assert np.all(np.isclose(img.voxToSurfMat, v2s))
- assert np.all(np.isclose(img.surfToVoxMat, transform.invert(v2s)))
+ assert np.all(np.isclose(img.surfToVoxMat, affine.invert(v2s)))
assert np.all(np.isclose(img.worldToSurfMat, w2s))
- assert np.all(np.isclose(img.surfToWorldMat, transform.invert(w2s)))
+ assert np.all(np.isclose(img.surfToWorldMat, affine.invert(w2s)))
assert img.name == op.basename(testfile)
assert img.dataSource == testfile
diff --git a/tests/test_scripts/test_atlasq_query.py b/tests/test_scripts/test_atlasq_query.py
index a36e1eb5fa4f086b0e28bcd4440f08958ebb3d05..def92c956a08098afd779102a591044aa63740cf 100644
--- a/tests/test_scripts/test_atlasq_query.py
+++ b/tests/test_scripts/test_atlasq_query.py
@@ -15,11 +15,11 @@ import scipy.ndimage as ndi
import pytest
-import fsl.utils.transform as transform
+import fsl.transform.affine as affine
+import fsl.data.atlases as fslatlases
import fsl.utils.image.resample as resample
-import fsl.data.atlases as fslatlases
-import fsl.data.image as fslimage
-import fsl.scripts.atlasq as fslatlasq
+import fsl.data.image as fslimage
+import fsl.scripts.atlasq as fslatlasq
from .. import (tempdir,
make_random_mask,
@@ -152,7 +152,7 @@ def _gen_coord_voxel_query(atlas, use_label, q_type, q_in, res):
dlo = (0, 0, 0)
dhi = a_img.shape
else:
- dlo, dhi = transform.axisBounds(a_img.shape, a_img.voxToWorldMat)
+ dlo, dhi = affine.axisBounds(a_img.shape, a_img.voxToWorldMat)
dlen = [hi - lo for lo, hi in zip(dlo, dhi)]
@@ -187,7 +187,7 @@ def _gen_coord_voxel_query(atlas, use_label, q_type, q_in, res):
coords = np.array(coords, dtype=dtype)
if not voxel:
- coords = transform.transform(coords, a_img.voxToWorldMat)
+ coords = affine.transform(coords, a_img.voxToWorldMat)
return tuple([dtype(c) for c in coords])
@@ -453,8 +453,8 @@ def test_bad_mask(seed):
dtype=np.float32))
wrongspace = fslimage.Image(
np.random.random((20, 20, 20)),
- xform=transform.concat(atlas.voxToWorldMat,
- np.diag([2, 2, 2, 1])))
+ xform=affine.concat(atlas.voxToWorldMat,
+ np.diag([2, 2, 2, 1])))
print(wrongdims.shape)
print(wrongspace.shape)
diff --git a/tests/test_scripts/test_fsl_apply_x5.py b/tests/test_scripts/test_fsl_apply_x5.py
new file mode 100644
index 0000000000000000000000000000000000000000..dbe2bd3a7ec6cfae76fbb07e8a7e663fc5b135ab
--- /dev/null
+++ b/tests/test_scripts/test_fsl_apply_x5.py
@@ -0,0 +1,264 @@
+#!/usr/bin/env python
+
+
+import numpy as np
+
+import pytest
+
+import fsl.scripts.fsl_apply_x5 as fsl_apply_x5
+import fsl.data.image as fslimage
+import fsl.utils.image.resample as resample
+import fsl.utils.image.roi as roi
+import fsl.transform.x5 as x5
+import fsl.transform.affine as affine
+import fsl.utils.tempdir as tempdir
+
+from ..test_transform.test_nonlinear import _affine_field
+
+
+def _random_affine():
+ return affine.compose(
+ np.random.randint(2, 5, 3),
+ np.random.randint(1, 10, 3),
+ np.random.random(3))
+
+
+def _random_image(aff, shape=None):
+
+ if shape is None:
+ shape = np.random.randint(10, 50, 3)
+
+ data = (np.random.random(shape) - 0.5) * 10
+
+ return fslimage.Image(data, xform=aff)
+
+
+def test_help():
+ def run(args):
+ with pytest.raises(SystemExit) as e:
+ fsl_apply_x5.main(args)
+ assert e.value.code == 0
+
+ run([])
+ run(['-h'])
+
+
+
+def test_linear():
+ with tempdir.tempdir():
+
+ src2ref = _random_affine()
+ src = _random_image(np.eye(4))
+ ref = _random_image(src2ref)
+
+ x5.writeLinearX5('xform.x5', src2ref, src, ref)
+
+ src.save('src')
+ ref.save('ref')
+
+ fsl_apply_x5.main('src xform.x5 out'.split())
+
+ result = fslimage.Image('out')
+ expect = resample.resampleToReference(src, ref, matrix=src2ref)[0]
+
+ assert result.sameSpace(ref)
+ assert np.all(np.isclose(result.data, expect))
+
+
+def test_nonlinear():
+ with tempdir.tempdir():
+
+ src2ref = _random_affine()
+ ref2src = affine.invert(src2ref)
+ src = _random_image(np.eye(4))
+ ref = _random_image(src2ref)
+ field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+ x5.writeNonLinearX5('xform.x5', field)
+
+ src.save('src')
+
+ fsl_apply_x5.main('src xform.x5 out'.split())
+
+ result = fslimage.Image('out')
+ expect = resample.resampleToReference(src, ref, matrix=src2ref)[0]
+
+ assert result.sameSpace(ref)
+
+ # We might get truncation on the edges
+ result = result.data[1:-1, 1:-1, 1:-1]
+ expect = expect[ 1:-1, 1:-1, 1:-1]
+
+ assert np.all(np.isclose(result, expect))
+
+
+def test_linear_altref():
+ with tempdir.tempdir():
+
+ src2ref = affine.scaleOffsetXform([1, 1, 1], [5, 5, 5])
+ altv2w = affine.scaleOffsetXform([1, 1, 1], [10, 10, 10])
+
+ srcdata = np.random.randint(1, 65536, (10, 10, 10))
+ src = fslimage.Image(srcdata, xform=np.eye(4))
+ ref = fslimage.Image(src.data, xform=src2ref)
+ altref = fslimage.Image(src.data, xform=altv2w)
+
+ src .save('src')
+ ref .save('ref')
+ altref.save('altref')
+
+ x5.writeLinearX5('xform.x5', src2ref, src, ref)
+
+ fsl_apply_x5.main('src xform.x5 out -r altref'.split())
+
+ result = fslimage.Image('out')
+ expect = np.zeros(srcdata.shape)
+ expect[:5, :5, :5] = srcdata[5:, 5:, 5:]
+
+ assert result.sameSpace(altref)
+ assert np.all(result.data == expect)
+
+
+def test_nonlinear_altref():
+ with tempdir.tempdir():
+
+ src2ref = affine.scaleOffsetXform([1, 1, 1], [5, 5, 5])
+ ref2src = affine.invert(src2ref)
+ altv2w = affine.scaleOffsetXform([1, 1, 1], [10, 10, 10])
+
+ srcdata = np.random.randint(1, 65536, (10, 10, 10))
+ src = fslimage.Image(srcdata, xform=np.eye(4))
+ ref = fslimage.Image(src.data, xform=src2ref)
+ altref = fslimage.Image(src.data, xform=altv2w)
+
+ field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+ src .save('src')
+ ref .save('ref')
+ altref.save('altref')
+
+ x5.writeNonLinearX5('xform.x5', field)
+
+ fsl_apply_x5.main('src xform.x5 out -r altref'.split())
+
+ result = fslimage.Image('out')
+ expect = np.zeros(srcdata.shape)
+ expect[:5, :5, :5] = srcdata[5:, 5:, 5:]
+
+ assert result.sameSpace(altref)
+ assert np.all(result.data == expect)
+
+
+def test_linear_altsrc():
+ with tempdir.tempdir():
+
+ src2ref = _random_affine()
+ src = _random_image(np.eye(4), (20, 20, 20))
+ ref = _random_image(src2ref)
+
+ x5.writeLinearX5('xform.x5', src2ref, src, ref)
+
+ src.save('src')
+ ref.save('ref')
+
+ srclo, xf = resample.resample(src, (10, 10, 10))
+ srclo = fslimage.Image(srclo, xform=xf)
+ srchi, xf = resample.resample(src, (40, 40, 40))
+ srchi = fslimage.Image(srchi, xform=xf)
+
+ srcoff = roi.roi(src, [(-10, 10), (-10, 10), (-10, 10)])
+
+ srclo .save('srclo')
+ srchi .save('srchi')
+ srcoff.save('srcoff')
+
+ fsl_apply_x5.main('src xform.x5 out' .split())
+ fsl_apply_x5.main('srclo xform.x5 outlo' .split())
+ fsl_apply_x5.main('srchi xform.x5 outhi' .split())
+ fsl_apply_x5.main('srcoff xform.x5 outoff'.split())
+
+ out = fslimage.Image('out')
+ outlo = fslimage.Image('outlo')
+ outhi = fslimage.Image('outhi')
+ outoff = fslimage.Image('outoff')
+
+ exp = resample.resampleToReference(src, ref, matrix=src2ref)[0]
+ explo = resample.resampleToReference(srclo, ref, matrix=src2ref)[0]
+ exphi = resample.resampleToReference(srchi, ref, matrix=src2ref)[0]
+ expoff = resample.resampleToReference(srcoff, ref, matrix=src2ref)[0]
+
+ assert out .sameSpace(ref)
+ assert outlo .sameSpace(ref)
+ assert outhi .sameSpace(ref)
+ assert outoff.sameSpace(ref)
+
+ assert np.all(np.isclose(out .data, exp))
+ assert np.all(np.isclose(outlo .data, explo))
+ assert np.all(np.isclose(outhi .data, exphi))
+ assert np.all(np.isclose(outoff.data, expoff))
+
+
+def test_nonlinear_altsrc():
+ with tempdir.tempdir():
+
+ src2ref = _random_affine()
+ ref2src = affine.invert(src2ref)
+ src = _random_image(np.eye(4), (20, 20, 20))
+ ref = _random_image(src2ref, (20, 20, 20))
+
+ field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+ x5.writeNonLinearX5('xform.x5', field)
+
+ src.save('src')
+ ref.save('ref')
+
+ srclo, xf = resample.resample(src, (10, 10, 10), origin='corner')
+ srclo = fslimage.Image(srclo, xform=xf)
+ srchi, xf = resample.resample(src, (40, 40, 40), origin='corner')
+ srchi = fslimage.Image(srchi, xform=xf)
+
+ srcoff = roi.roi(src, [(-10, 10), (-10, 10), (-10, 10)])
+
+ srclo .save('srclo')
+ srchi .save('srchi')
+ srcoff.save('srcoff')
+
+ # The up-sampled source image is subject to
+ # an unavoidable interpolation, as the
+ # deformation field coordinates are affine-
+ # transformed to the space of the alternate
+ # source image. So in this test case we use
+ # nearest-neighbour interp, so that we can
+ # get the same output as a standard linear
+ # resample
+ fsl_apply_x5.main('src xform.x5 out' .split())
+ fsl_apply_x5.main('srclo xform.x5 outlo' .split())
+ fsl_apply_x5.main('srchi xform.x5 outhi -i nearest'.split())
+ fsl_apply_x5.main('srcoff xform.x5 outoff' .split())
+
+ out = fslimage.Image('out')
+ outlo = fslimage.Image('outlo')
+ outhi = fslimage.Image('outhi')
+ outoff = fslimage.Image('outoff')
+
+ exp, x1 = resample.resampleToReference(src, ref, matrix=src2ref, mode='constant')
+ explo, x2 = resample.resampleToReference(srclo, ref, matrix=src2ref, mode='constant')
+ exphi, x3 = resample.resampleToReference(srchi, ref, matrix=src2ref, mode='constant', order=0)
+ expoff, x4 = resample.resampleToReference(srcoff, ref, matrix=src2ref, mode='constant')
+
+ assert out .sameSpace(ref)
+ assert outlo .sameSpace(ref)
+ assert outhi .sameSpace(ref)
+ assert outoff.sameSpace(ref)
+
+ # We get boundary issues just at the first
+ # voxel, so I'm masking that voxel out
+ for img in (out, outlo, outhi, outoff,
+ exp, explo, exphi, expoff):
+ img[0, 0, 0] = 0
+
+ assert np.all(np.isclose(out .data, exp))
+ assert np.all(np.isclose(outlo .data, explo))
+ assert np.all(np.isclose(outhi .data, exphi))
+ assert np.all(np.isclose(outoff.data, expoff))
diff --git a/tests/test_scripts/test_fsl_convert_x5.py b/tests/test_scripts/test_fsl_convert_x5.py
new file mode 100644
index 0000000000000000000000000000000000000000..b27cffa872b3919f0da939d175ff306b81f3cd04
--- /dev/null
+++ b/tests/test_scripts/test_fsl_convert_x5.py
@@ -0,0 +1,209 @@
+#!/usr/bin/env python
+
+
+import os.path as op
+import hashlib
+
+import pytest
+
+import numpy as np
+
+import fsl.utils.tempdir as tempdir
+import fsl.transform.affine as affine
+import fsl.transform.flirt as flirt
+import fsl.transform.fnirt as fnirt
+import fsl.transform.x5 as x5
+import fsl.data.image as fslimage
+
+import fsl.scripts.fsl_convert_x5 as fsl_convert_x5
+
+
+
+def random_image():
+ vx, vy, vz = np.random.randint(10, 50, 3)
+ dx, dy, dz = np.random.randint( 1, 10, 3)
+ data = (np.random.random((vx, vy, vz)) - 0.5) * 10
+ aff = affine.compose(
+ (dx, dy, dz),
+ np.random.randint(1, 100, 3),
+ np.random.random(3) * np.pi / 2)
+
+ return fslimage.Image(data, xform=aff)
+
+
+def test_convert_flirt_help():
+
+ def run(args):
+ with pytest.raises(SystemExit) as e:
+ fsl_convert_x5.main(args)
+ assert e.value.code == 0
+
+ run([])
+ run(['-h'])
+ run(['flirt'])
+ run(['fnirt'])
+ run(['flirt', '-h'])
+ run(['fnirt', '-h'])
+
+
+def test_convert_flirt():
+ with tempdir.tempdir():
+ src = random_image()
+ ref = random_image()
+ src.save('src')
+ ref.save('ref')
+
+ xform = affine.compose(
+ np.random.randint(1, 10, 3),
+ np.random.randint(-100, 100, 3),
+ (np.random.random(3) - 0.5) * np.pi)
+
+ np.savetxt('src2ref.mat', xform)
+
+ fsl_convert_x5.main('flirt -s src -r ref '
+ 'src2ref.mat src2ref.x5'.split())
+ expxform = affine.concat(
+ ref.getAffine('fsl', 'world'),
+ xform,
+ src.getAffine('world', 'fsl'))
+ gotxform, gotsrc, gotref = x5.readLinearX5('src2ref.x5')
+ assert np.all(np.isclose(gotxform, expxform))
+ assert src.sameSpace(gotsrc)
+ assert ref.sameSpace(gotref)
+
+ fsl_convert_x5.main('flirt src2ref.x5 src2ref_copy.mat'.split())
+
+ gotxform = flirt.readFlirt('src2ref_copy.mat')
+ assert np.all(np.isclose(gotxform, xform))
+
+
+def test_convert_flirt_sameformat():
+ with tempdir.tempdir():
+ src = random_image()
+ ref = random_image()
+ src.save('src')
+ ref.save('ref')
+
+ xform = affine.compose(
+ np.random.randint(1, 10, 3),
+ np.random.randint(-100, 100, 3),
+ (np.random.random(3) - 0.5) * np.pi)
+
+ np.savetxt('src2ref.mat', xform)
+
+ # test both .mat and .x5
+ fsl_convert_x5.main('flirt -s src -r ref '
+ 'src2ref.mat src2ref.x5'.split())
+
+ # mat -> mat
+ fsl_convert_x5.main('flirt -s src -r ref '
+ 'src2ref.mat copy.mat'.split())
+
+ # x5 -> x5
+ fsl_convert_x5.main('flirt -s src -r ref '
+ 'src2ref.x5 copy.x5'.split())
+
+ with open('src2ref.mat', 'rb') as f: origmat = hashlib.md5(f.read()).digest()
+ with open('copy.mat', 'rb') as f: copymat = hashlib.md5(f.read()).digest()
+ with open('src2ref.x5', 'rb') as f: origx5 = hashlib.md5(f.read()).digest()
+ with open('copy.x5', 'rb') as f: copyx5 = hashlib.md5(f.read()).digest()
+
+ assert origmat == copymat
+ assert origx5 == copyx5
+
+
+
+def test_convert_fnirt_deformation_field():
+
+ datadir = op.join(op.dirname(__file__), '..',
+ 'test_transform', 'testdata', 'nonlinear')
+ srcfile = op.join(datadir, 'src.nii.gz')
+ reffile = op.join(datadir, 'ref.nii.gz')
+ dffile = op.join(datadir, 'displacementfield.nii.gz')
+
+ with tempdir.tempdir():
+
+ # nii -> x5
+ fsl_convert_x5.main('fnirt -s {} -r {} {} disp.x5'.format(
+ srcfile, reffile, dffile).split())
+
+ # x5 -> nii
+ fsl_convert_x5.main('fnirt disp.x5 disp.nii.gz'.split())
+
+ src = fslimage.Image(srcfile)
+ ref = fslimage.Image(reffile)
+ df = fnirt.readFnirt(dffile, src, ref)
+ dfnii = fnirt.readFnirt('disp.nii.gz', src, ref)
+
+ assert dfnii.src.sameSpace(src)
+ assert dfnii.ref.sameSpace(ref)
+ assert dfnii.srcSpace == df.srcSpace
+ assert dfnii.refSpace == df.refSpace
+ assert dfnii.deformationType == df.deformationType
+ assert np.all(np.isclose(dfnii.data, df.data))
+
+
+def test_convert_fnirt_coefficient_field():
+
+ datadir = op.join(op.dirname(__file__), '..',
+ 'test_transform', 'testdata', 'nonlinear')
+ srcfile = op.join(datadir, 'src.nii.gz')
+ reffile = op.join(datadir, 'ref.nii.gz')
+ cffile = op.join(datadir, 'coefficientfield.nii.gz')
+ dffile = op.join(datadir, 'displacementfield.nii.gz')
+
+ with tempdir.tempdir():
+
+ # nii -> x5
+ fsl_convert_x5.main('fnirt -s {} -r {} {} coef.x5'.format(
+ srcfile, reffile, cffile).split())
+
+ # x5 -> nii
+ fsl_convert_x5.main('fnirt coef.x5 coef.nii.gz'.split())
+
+ src = fslimage.Image(srcfile)
+ ref = fslimage.Image(reffile)
+ df = fnirt.readFnirt(dffile, src, ref)
+ dfnii = fnirt.readFnirt('coef.nii.gz', src, ref)
+
+ assert dfnii .sameSpace(df)
+ assert dfnii.src.sameSpace(src)
+ assert dfnii.ref.sameSpace(ref)
+
+ assert dfnii.srcSpace == df.srcSpace
+ assert dfnii.refSpace == df.refSpace
+ assert dfnii.deformationType == 'relative'
+
+ diff = np.abs(dfnii.data - df.data)
+ tols = {'rtol' : 1e-5, 'atol' : 1e-5}
+ assert np.all(np.isclose(dfnii.data, df.data, **tols))
+
+
+def test_convert_fnirt_sameformat():
+
+ datadir = op.join(op.dirname(__file__), '..',
+ 'test_transform', 'testdata', 'nonlinear')
+ srcfile = op.join(datadir, 'src.nii.gz')
+ reffile = op.join(datadir, 'ref.nii.gz')
+ dffile = op.join(datadir, 'displacementfield.nii.gz')
+
+ with tempdir.tempdir():
+
+ base = list('fnirt -s {} -r {}'.format(srcfile, reffile).split())
+
+ # test both .mat and .x5
+ fsl_convert_x5.main(base + [dffile, 'src2ref.x5'])
+
+ # nii -> nii
+ fsl_convert_x5.main(base + [dffile, 'copy.nii.gz'])
+
+ # x5 -> x5
+ fsl_convert_x5.main(base + ['src2ref.x5', 'copy.x5'])
+
+ with open(dffile, 'rb') as f: origdef = hashlib.md5(f.read()).digest()
+ with open('copy.nii.gz', 'rb') as f: copydef = hashlib.md5(f.read()).digest()
+ with open('src2ref.x5', 'rb') as f: origx5 = hashlib.md5(f.read()).digest()
+ with open('copy.x5', 'rb') as f: copyx5 = hashlib.md5(f.read()).digest()
+
+ assert origdef == copydef
+ assert origx5 == copyx5
diff --git a/tests/test_transform/__init__.py b/tests/test_transform/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/tests/test_transform.py b/tests/test_transform/test_affine.py
similarity index 72%
rename from tests/test_transform.py
rename to tests/test_transform/test_affine.py
index 70b86019aae9a7374a906dda830fd0a8f7c0966d..7a69b86e4d79531f53d52be991ba5afe994b248a 100644
--- a/tests/test_transform.py
+++ b/tests/test_transform/test_affine.py
@@ -19,8 +19,7 @@ import six
import pytest
-import fsl.utils.transform as transform
-import fsl.data.image as fslimage
+import fsl.transform.affine as affine
datadir = op.join(op.dirname(__file__), 'testdata')
@@ -57,7 +56,7 @@ def test_invert():
x = testdata[i * 4:i * 4 + 4, 0:4]
invx = testdata[i * 4:i * 4 + 4, 4:8]
- result = transform.invert(x)
+ result = affine.invert(x)
assert np.all(np.isclose(invx, result))
@@ -87,7 +86,7 @@ def test_concat():
for inputs, expected in tests:
- result = transform.concat(*inputs)
+ result = affine.concat(*inputs)
assert np.all(np.isclose(result, expected))
@@ -109,10 +108,10 @@ def test_veclength(seed):
vtype = random.choice((list, tuple, np.array))
v = vtype(v)
- assert np.isclose(transform.veclength(v), l(v))
+ assert np.isclose(affine.veclength(v), l(v))
# Multiple vectors in parallel
- result = transform.veclength(vectors)
+ result = affine.veclength(vectors)
expected = l(vectors)
assert np.all(np.isclose(result, expected))
@@ -122,8 +121,8 @@ def test_normalise(seed):
vectors = -100 + 200 * np.random.random((200, 3))
def parallel(v1, v2):
- v1 = v1 / transform.veclength(v1)
- v2 = v2 / transform.veclength(v2)
+ v1 = v1 / affine.veclength(v1)
+ v2 = v2 / affine.veclength(v2)
return np.isclose(np.dot(v1, v2), 1)
@@ -131,16 +130,16 @@ def test_normalise(seed):
vtype = random.choice((list, tuple, np.array))
v = vtype(v)
- vn = transform.normalise(v)
- vl = transform.veclength(vn)
+ vn = affine.normalise(v)
+ vl = affine.veclength(vn)
assert np.isclose(vl, 1.0)
assert parallel(v, vn)
# normalise should also be able
# to do multiple vectors at once
- results = transform.normalise(vectors)
- lengths = transform.veclength(results)
+ results = affine.normalise(vectors)
+ lengths = affine.veclength(results)
pars = np.zeros(200)
for i in range(200):
@@ -172,7 +171,7 @@ def test_scaleOffsetXform():
expected = [[float(v) for v in l.split()] for l in expected]
expected = np.array(expected)
- result = transform.scaleOffsetXform(scales, offsets)
+ result = affine.scaleOffsetXform(scales, offsets)
assert np.all(np.isclose(result, expected))
@@ -208,11 +207,11 @@ def test_scaleOffsetXform():
for (scale, expected) in stests:
expected = np.array(expected).reshape(4, 4)
- result = transform.scaleOffsetXform(scale, 0)
+ result = affine.scaleOffsetXform(scale, 0)
assert np.all(np.isclose(result, expected))
for (offset, expected) in otests:
expected = np.array(expected).reshape(4, 4)
- result = transform.scaleOffsetXform(1, offset)
+ result = affine.scaleOffsetXform(1, offset)
assert np.all(np.isclose(result, expected))
@@ -227,8 +226,8 @@ def test_compose_and_decompose():
xform = lines[i * 4: i * 4 + 4]
xform = np.genfromtxt(xform)
- scales, offsets, rotations = transform.decompose(xform)
- result = transform.compose(scales, offsets, rotations)
+ scales, offsets, rotations = affine.decompose(xform)
+ result = affine.compose(scales, offsets, rotations)
assert np.all(np.isclose(xform, result, atol=1e-5))
@@ -236,22 +235,22 @@ def test_compose_and_decompose():
# different rotation origin, but we test
# explicitly passing the origin for
# completeness
- scales, offsets, rotations = transform.decompose(xform)
- result = transform.compose(scales, offsets, rotations, [0, 0, 0])
+ scales, offsets, rotations = affine.decompose(xform)
+ result = affine.compose(scales, offsets, rotations, [0, 0, 0])
assert np.all(np.isclose(xform, result, atol=1e-5))
# compose should also accept a rotation matrix
rots = [np.pi / 5, np.pi / 4, np.pi / 3]
- rmat = transform.axisAnglesToRotMat(*rots)
- xform = transform.compose([1, 1, 1], [0, 0, 0], rmat)
+ rmat = affine.axisAnglesToRotMat(*rots)
+ xform = affine.compose([1, 1, 1], [0, 0, 0], rmat)
# And the angles flag should cause decompose
# to return the rotation matrix, instead of
# the axis angls
- sc, of, rot = transform.decompose(xform)
- scat, ofat, rotat = transform.decompose(xform, angles=True)
- scaf, ofaf, rotaf = transform.decompose(xform, angles=False)
+ sc, of, rot = affine.decompose(xform)
+ scat, ofat, rotat = affine.decompose(xform, angles=True)
+ scaf, ofaf, rotaf = affine.decompose(xform, angles=False)
sc, of, rot = np.array(sc), np.array(of), np.array(rot)
scat, ofat, rotat = np.array(scat), np.array(ofat), np.array(rotat)
@@ -270,8 +269,8 @@ def test_compose_and_decompose():
# decompose should accept a 3x3
# affine, and return translations of 0
- transform.decompose(xform[:3, :3])
- sc, of, rot = transform.decompose(xform[:3, :3])
+ affine.decompose(xform[:3, :3])
+ sc, of, rot = affine.decompose(xform[:3, :3])
sc, of, rot = np.array(sc), np.array(of), np.array(rot)
assert np.all(np.isclose(sc, [1, 1, 1]))
assert np.all(np.isclose(of, [0, 0, 0]))
@@ -291,8 +290,8 @@ def test_rotMatToAxisAngles(seed):
-pi2 + 2 * pi2 * np.random.random(),
-pi + 2 * pi * np.random.random()]
- rmat = transform.axisAnglesToRotMat(*rots)
- gotrots = transform.rotMatToAxisAngles(rmat)
+ rmat = affine.axisAnglesToRotMat(*rots)
+ gotrots = affine.rotMatToAxisAngles(rmat)
assert np.all(np.isclose(rots, gotrots))
@@ -311,9 +310,9 @@ def test_rotMatToAffine(seed):
if np.random.random() < 0.5: origin = None
else: origin = np.random.random(3)
- rmat = transform.axisAnglesToRotMat(*rots)
- mataff = transform.rotMatToAffine(rmat, origin)
- rotaff = transform.rotMatToAffine(rots, origin)
+ rmat = affine.axisAnglesToRotMat(*rots)
+ mataff = affine.rotMatToAffine(rmat, origin)
+ rotaff = affine.rotMatToAffine(rots, origin)
exp = np.eye(4)
exp[:3, :3] = rmat
@@ -350,7 +349,7 @@ def test_axisBounds():
shape, origin, boundary, xform, expected = readTest(i)
for axes in allAxes:
- result = transform.axisBounds(shape,
+ result = affine.axisBounds(shape,
xform,
axes=axes,
origin=origin,
@@ -372,14 +371,14 @@ def test_axisBounds():
# US-spelling
assert np.all(np.isclose(
expected,
- transform.axisBounds(
+ affine.axisBounds(
shape, xform, origin='center', boundary=boundary)))
# Bad origin/boundary values
with pytest.raises(ValueError):
- transform.axisBounds(shape, xform, origin='Blag', boundary=boundary)
+ affine.axisBounds(shape, xform, origin='Blag', boundary=boundary)
with pytest.raises(ValueError):
- transform.axisBounds(shape, xform, origin=origin, boundary='Blufu')
+ affine.axisBounds(shape, xform, origin=origin, boundary='Blufu')
def test_transform():
@@ -413,7 +412,7 @@ def test_transform():
lines = readlines(testfile)
xform = np.genfromtxt(lines[:4])
expected = np.genfromtxt(lines[ 4:])
- result = transform.transform(testcoords, xform)
+ result = affine.transform(testcoords, xform)
assert np.all(np.isclose(expected, result))
@@ -423,7 +422,7 @@ def test_transform():
for axes in allAxes:
atestcoords = testcoords[:, axes]
aexpected = expected[ :, axes]
- aresult = transform.transform(atestcoords, xform, axes=axes)
+ aresult = affine.transform(atestcoords, xform, axes=axes)
assert np.all(np.isclose(aexpected, aresult))
@@ -434,26 +433,26 @@ def test_transform():
coords = badcoords[:, :3]
with pytest.raises(IndexError):
- transform.transform(coords, badxform)
+ affine.transform(coords, badxform)
with pytest.raises(ValueError):
- transform.transform(badcoords, xform)
+ affine.transform(badcoords, xform)
with pytest.raises(ValueError):
- transform.transform(badcoords.reshape(5, 2, 4), xform)
+ affine.transform(badcoords.reshape(5, 2, 4), xform)
with pytest.raises(ValueError):
- transform.transform(badcoords.reshape(5, 2, 4), xform, axes=1)
+ affine.transform(badcoords.reshape(5, 2, 4), xform, axes=1)
with pytest.raises(ValueError):
- transform.transform(badcoords[:, (1, 2, 3)], xform, axes=[1, 2])
+ affine.transform(badcoords[:, (1, 2, 3)], xform, axes=[1, 2])
def test_transform_vector(seed):
# Some transform with a
# translation component
- xform = transform.compose([1, 2, 3],
+ xform = affine.compose([1, 2, 3],
[5, 10, 15],
[np.pi / 2, np.pi / 2, 0])
@@ -464,9 +463,9 @@ def test_transform_vector(seed):
vecExpected = np.dot(xform, list(v) + [0])[:3]
ptExpected = np.dot(xform, list(v) + [1])[:3]
- vecResult = transform.transform(v, xform, vector=True)
- vec33Result = transform.transform(v, xform[:3, :3], vector=True)
- ptResult = transform.transform(v, xform, vector=False)
+ vecResult = affine.transform(v, xform, vector=True)
+ vec33Result = affine.transform(v, xform[:3, :3], vector=True)
+ ptResult = affine.transform(v, xform, vector=False)
assert np.all(np.isclose(vecExpected, vecResult))
assert np.all(np.isclose(vecExpected, vec33Result))
@@ -489,83 +488,33 @@ def test_transformNormal(seed):
rotations = -np.pi + np.random.random(3) * 2 * np.pi
origin = -100 + np.random.random(3) * 200
- xform = transform.compose(scales,
+ xform = affine.compose(scales,
offsets,
rotations,
origin)
expected = tn(n, xform)
- result = transform.transformNormal(n, xform)
+ result = affine.transformNormal(n, xform)
assert np.all(np.isclose(expected, result))
-def test_flirtMatrixToSform():
-
- testfile = op.join(datadir, 'test_transform_test_flirtMatrixToSform.txt')
- lines = readlines(testfile)
- ntests = len(lines) // 18
-
- for i in range(ntests):
- tlines = lines[i * 18: i * 18 + 18]
- srcShape = [int( w) for w in tlines[0].split()]
- srcXform = np.genfromtxt(tlines[1:5])
- refShape = [int( w) for w in tlines[5].split()]
- refXform = np.genfromtxt(tlines[6:10])
- flirtMat = np.genfromtxt(tlines[10:14])
- expected = np.genfromtxt(tlines[14:18])
-
- srcImg = fslimage.Image(np.zeros(srcShape), xform=srcXform)
- refImg = fslimage.Image(np.zeros(refShape), xform=refXform)
-
- result = transform.flirtMatrixToSform(flirtMat, srcImg, refImg)
-
- assert np.all(np.isclose(result, expected))
-
-
-def test_sformToFlirtMatrix():
- testfile = op.join(datadir, 'test_transform_test_flirtMatrixToSform.txt')
- lines = readlines(testfile)
- ntests = len(lines) // 18
-
- for i in range(ntests):
- tlines = lines[i * 18: i * 18 + 18]
- srcShape = [int( w) for w in tlines[0].split()]
- srcXform = np.genfromtxt(tlines[1:5])
- refShape = [int( w) for w in tlines[5].split()]
- refXform = np.genfromtxt(tlines[6:10])
- expected = np.genfromtxt(tlines[10:14])
- srcXformOvr = np.genfromtxt(tlines[14:18])
-
- srcImg1 = fslimage.Image(np.zeros(srcShape), xform=srcXform)
- srcImg2 = fslimage.Image(np.zeros(srcShape), xform=srcXform)
- refImg = fslimage.Image(np.zeros(refShape), xform=refXform)
-
- srcImg2.voxToWorldMat = srcXformOvr
-
- result1 = transform.sformToFlirtMatrix(srcImg1, refImg, srcXformOvr)
- result2 = transform.sformToFlirtMatrix(srcImg2, refImg)
-
- assert np.all(np.isclose(result1, expected))
- assert np.all(np.isclose(result2, expected))
-
-
def test_rmsdev():
t1 = np.eye(4)
- t2 = transform.scaleOffsetXform([1, 1, 1], [2, 0, 0])
+ t2 = affine.scaleOffsetXform([1, 1, 1], [2, 0, 0])
- assert np.isclose(transform.rmsdev(t1, t2), 2)
- assert np.isclose(transform.rmsdev(t1, t2, R=2), 2)
- assert np.isclose(transform.rmsdev(t1, t2, R=2, xc=(1, 1, 1)), 2)
+ assert np.isclose(affine.rmsdev(t1, t2), 2)
+ assert np.isclose(affine.rmsdev(t1, t2, R=2), 2)
+ assert np.isclose(affine.rmsdev(t1, t2, R=2, xc=(1, 1, 1)), 2)
t1 = np.eye(3)
lastdist = 0
for i in range(1, 11):
rot = np.pi * i / 10.0
- t2 = transform.axisAnglesToRotMat(rot, 0, 0)
- result = transform.rmsdev(t1, t2)
+ t2 = affine.axisAnglesToRotMat(rot, 0, 0)
+ result = affine.rmsdev(t1, t2)
assert result > lastdist
@@ -573,8 +522,8 @@ def test_rmsdev():
for i in range(11, 20):
rot = np.pi * i / 10.0
- t2 = transform.axisAnglesToRotMat(rot, 0, 0)
- result = transform.rmsdev(t1, t2)
+ t2 = affine.axisAnglesToRotMat(rot, 0, 0)
+ result = affine.rmsdev(t1, t2)
assert result < lastdist
diff --git a/tests/test_transform/test_flirt.py b/tests/test_transform/test_flirt.py
new file mode 100644
index 0000000000000000000000000000000000000000..d2dd68e29b114ca9813ebd1eb8511df825d8edaf
--- /dev/null
+++ b/tests/test_transform/test_flirt.py
@@ -0,0 +1,132 @@
+#!/usr/bin/env python
+#
+# test_transform_flirt.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import itertools as it
+import os.path as op
+
+import numpy as np
+
+import fsl.data.image as fslimage
+import fsl.transform.flirt as flirt
+import fsl.transform.affine as affine
+import fsl.utils.tempdir as tempdir
+
+from .test_affine import readlines
+
+from .. import make_random_image
+
+
+datadir = op.join(op.dirname(__file__), 'testdata')
+
+
+def test_read_write():
+ with tempdir.tempdir():
+ aff = np.random.random((4, 4))
+ flirt.writeFlirt(aff, 'aff.mat')
+ got = flirt.readFlirt('aff.mat')
+ assert np.all(np.isclose(aff, got))
+
+
+def test_fromFlirt():
+
+ src = affine.compose(np.random.randint( 1, 5, 3),
+ np.random.randint(-20, 20, 3),
+ np.random.random(3) - 0.5)
+ ref = affine.compose(np.random.randint( 1, 5, 3),
+ np.random.randint(-20, 20, 3),
+ np.random.random(3) - 0.5)
+
+ src = fslimage.Image(make_random_image(xform=src))
+ ref = fslimage.Image(make_random_image(xform=ref))
+ flirtmat = affine.concat(ref.getAffine('world', 'fsl'),
+ src.getAffine('fsl', 'world'))
+
+
+ spaces = it.permutations(('voxel', 'fsl', 'world'), 2)
+
+ for from_, to in spaces:
+ got = flirt.fromFlirt(flirtmat, src, ref, from_, to)
+ exp = affine.concat(ref.getAffine('fsl', to),
+ flirtmat,
+ src.getAffine(from_, 'fsl'))
+
+ assert np.all(np.isclose(got, exp))
+
+
+def test_toFlirt():
+
+ src = affine.compose(np.random.randint( 1, 5, 3),
+ np.random.randint(-20, 20, 3),
+ np.random.random(3) - 0.5)
+ ref = affine.compose(np.random.randint( 1, 5, 3),
+ np.random.randint(-20, 20, 3),
+ np.random.random(3) - 0.5)
+
+ src = fslimage.Image(make_random_image(xform=src))
+ ref = fslimage.Image(make_random_image(xform=ref))
+ flirtmat = affine.concat(ref.getAffine('world', 'fsl'),
+ src.getAffine('fsl', 'world'))
+
+
+ spaces = it.permutations(('voxel', 'fsl', 'world'), 2)
+
+ for from_, to in spaces:
+ xform = affine.concat(ref.getAffine('world', to),
+ src.getAffine(from_, 'world'))
+ got = flirt.toFlirt(xform, src, ref, from_, to)
+
+ assert np.all(np.isclose(got, flirtmat))
+
+
+def test_flirtMatrixToSform():
+
+ testfile = op.join(datadir, 'test_transform_test_flirtMatrixToSform.txt')
+ lines = readlines(testfile)
+ ntests = len(lines) // 18
+
+ for i in range(ntests):
+ tlines = lines[i * 18: i * 18 + 18]
+ srcShape = [int( w) for w in tlines[0].split()]
+ srcXform = np.genfromtxt(tlines[1:5])
+ refShape = [int( w) for w in tlines[5].split()]
+ refXform = np.genfromtxt(tlines[6:10])
+ flirtMat = np.genfromtxt(tlines[10:14])
+ expected = np.genfromtxt(tlines[14:18])
+
+ srcImg = fslimage.Image(np.zeros(srcShape), xform=srcXform)
+ refImg = fslimage.Image(np.zeros(refShape), xform=refXform)
+
+ result = flirt.flirtMatrixToSform(flirtMat, srcImg, refImg)
+
+ assert np.all(np.isclose(result, expected))
+
+
+def test_sformToFlirtMatrix():
+ testfile = op.join(datadir, 'test_transform_test_flirtMatrixToSform.txt')
+ lines = readlines(testfile)
+ ntests = len(lines) // 18
+
+ for i in range(ntests):
+ tlines = lines[i * 18: i * 18 + 18]
+ srcShape = [int( w) for w in tlines[0].split()]
+ srcXform = np.genfromtxt(tlines[1:5])
+ refShape = [int( w) for w in tlines[5].split()]
+ refXform = np.genfromtxt(tlines[6:10])
+ expected = np.genfromtxt(tlines[10:14])
+ srcXformOvr = np.genfromtxt(tlines[14:18])
+
+ srcImg1 = fslimage.Image(np.zeros(srcShape), xform=srcXform)
+ srcImg2 = fslimage.Image(np.zeros(srcShape), xform=srcXform)
+ refImg = fslimage.Image(np.zeros(refShape), xform=refXform)
+
+ srcImg2.voxToWorldMat = srcXformOvr
+
+ result1 = flirt.sformToFlirtMatrix(srcImg1, refImg, srcXformOvr)
+ result2 = flirt.sformToFlirtMatrix(srcImg2, refImg)
+
+ assert np.all(np.isclose(result1, expected))
+ assert np.all(np.isclose(result2, expected))
diff --git a/tests/test_transform/test_fnirt.py b/tests/test_transform/test_fnirt.py
new file mode 100644
index 0000000000000000000000000000000000000000..eefc68b64e471dc052ff4fb92b637f56a599a8fa
--- /dev/null
+++ b/tests/test_transform/test_fnirt.py
@@ -0,0 +1,127 @@
+#!/usr/bin/env python
+#
+# test_fnirt.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import os.path as op
+import itertools as it
+
+import numpy as np
+
+import pytest
+
+import fsl.data.image as fslimage
+import fsl.transform.affine as affine
+import fsl.transform.nonlinear as nonlinear
+import fsl.transform.fnirt as fnirt
+
+from .test_nonlinear import _random_affine_field
+
+
+datadir = op.join(op.dirname(__file__), 'testdata', 'nonlinear')
+
+
+def test_readFnirt():
+
+ src = op.join(datadir, 'src')
+ ref = op.join(datadir, 'ref')
+ coef = op.join(datadir, 'coefficientfield')
+ disp = op.join(datadir, 'displacementfield')
+
+ src = fslimage.Image(src)
+ ref = fslimage.Image(ref)
+ coef = fnirt.readFnirt(coef, src, ref)
+ disp = fnirt.readFnirt(disp, src, ref)
+
+ with pytest.raises(ValueError):
+ fnirt.readFnirt(src, src, ref)
+
+ assert isinstance(coef, nonlinear.CoefficientField)
+ assert isinstance(disp, nonlinear.DeformationField)
+
+ assert coef.src.sameSpace(src)
+ assert coef.ref.sameSpace(ref)
+ assert disp.src.sameSpace(src)
+ assert disp.ref.sameSpace(ref)
+ assert coef.srcSpace == 'fsl'
+ assert coef.refSpace == 'fsl'
+ assert disp.srcSpace == 'fsl'
+ assert disp.refSpace == 'fsl'
+
+
+def test_toFnirt():
+
+ def check(got, exp):
+ tol = dict(atol=1e-5, rtol=1e-5)
+ assert np.all(np.isclose(got.data, exp.data, **tol))
+ assert got.src.sameSpace(exp.src)
+ assert got.ref.sameSpace(exp.ref)
+ assert got.srcSpace == 'fsl'
+ assert got.refSpace == 'fsl'
+
+ basefield, xform = _random_affine_field()
+ src = basefield.src
+ ref = basefield.ref
+
+ spaces = it.permutations(('voxel', 'fsl', 'world'), 2)
+
+ for from_, to in spaces:
+ field = nonlinear.convertDeformationSpace(basefield, from_, to)
+ got = fnirt.toFnirt(field)
+ check(got, basefield)
+
+ src = fslimage.Image(op.join(datadir, 'src'))
+ ref = fslimage.Image(op.join(datadir, 'ref'))
+ coef = fnirt.readFnirt(op.join(datadir, 'coefficientfield'), src, ref)
+ got = fnirt.toFnirt(coef)
+ check(got, coef)
+
+
+def test_fromFnirt():
+
+ basefield, basexform = _random_affine_field()
+ src = basefield.src
+ ref = basefield.ref
+ spaces = list(it.permutations(('voxel', 'fsl', 'world'), 2))
+
+ for from_, to in spaces:
+
+ got = fnirt.fromFnirt(basefield, from_, to)
+
+ assert got.srcSpace == to
+ assert got.refSpace == from_
+
+ coords = [np.random.randint(0, basefield.shape[0], 5),
+ np.random.randint(0, basefield.shape[1], 5),
+ np.random.randint(0, basefield.shape[2], 5)]
+ coords = np.array(coords).T
+
+ coords = affine.transform(coords, ref.getAffine('voxel', from_))
+
+ aff = affine.concat(src.getAffine('fsl', to),
+ basexform,
+ ref.getAffine(from_, 'fsl'))
+
+ got = got.transform(coords)
+ exp = affine.transform(coords, aff)
+
+ enan = np.isnan(exp)
+ gnan = np.isnan(got)
+
+ assert np.all(np.isclose(enan, gnan))
+ assert np.all(np.isclose(exp[~enan], got[~gnan]))
+
+ # Converting from a FNIRT coefficient field
+ src = fslimage.Image(op.join(datadir, 'src'))
+ ref = fslimage.Image(op.join(datadir, 'ref'))
+ coef = fnirt.readFnirt(op.join(datadir, 'coefficientfield'), src, ref)
+ disp = fnirt.readFnirt(op.join(datadir, 'displacementfield'), src, ref)
+
+ for from_, to in spaces:
+
+ cnv = fnirt.fromFnirt(coef, from_, to)
+ exp = nonlinear.convertDeformationSpace(disp, from_, to)
+ tol = dict(atol=1e-5, rtol=1e-5)
+ assert np.all(np.isclose(cnv.data, exp.data, **tol))
diff --git a/tests/test_transform/test_nonlinear.py b/tests/test_transform/test_nonlinear.py
new file mode 100644
index 0000000000000000000000000000000000000000..694162c3b3496cec470c68d57eb3142506e25cf0
--- /dev/null
+++ b/tests/test_transform/test_nonlinear.py
@@ -0,0 +1,481 @@
+#!/usr/bin/env python
+
+import itertools as it
+import os.path as op
+
+import numpy as np
+
+import fsl.data.image as fslimage
+import fsl.utils.image.resample as resample
+import fsl.utils.image.roi as roi
+import fsl.transform.affine as affine
+import fsl.transform.nonlinear as nonlinear
+import fsl.transform.fnirt as fnirt
+
+
+datadir = op.join(op.dirname(__file__), 'testdata')
+
+
+def _random_image():
+ vx, vy, vz = np.random.randint(10, 50, 3)
+ dx, dy, dz = np.random.randint( 1, 10, 3)
+ data = (np.random.random((vx, vy, vz)) - 0.5) * 10
+ aff = affine.compose(
+ (dx, dy, dz),
+ np.random.randint(1, 100, 3),
+ np.random.random(3) * np.pi / 2)
+
+ return fslimage.Image(data, xform=aff)
+
+
+def _random_field():
+
+ src = _random_image()
+ vx, vy, vz = np.random.randint(10, 50, 3)
+ dx, dy, dz = np.random.randint( 1, 10, 3)
+
+ field = (np.random.random((vx, vy, vz, 3)) - 0.5) * 10
+ aff = affine.compose(
+ (dx, dy, dz),
+ np.random.randint(1, 100, 3),
+ np.random.random(3) * np.pi / 2)
+
+ return nonlinear.DeformationField(field, src=src, xform=aff)
+
+
+def _affine_field(src, ref, xform, srcSpace, refSpace, shape=None, fv2w=None):
+
+ if shape is None: shape = ref.shape[:3]
+ if fv2w is None: fv2w = ref.getAffine('voxel', 'world')
+
+ rx, ry, rz = np.meshgrid(np.arange(shape[0]),
+ np.arange(shape[1]),
+ np.arange(shape[2]), indexing='ij')
+
+ rvoxels = np.vstack((rx.flatten(), ry.flatten(), rz.flatten())).T
+ f2r = affine.concat(ref.getAffine('world', refSpace), fv2w)
+ rcoords = affine.transform(rvoxels, f2r)
+ scoords = affine.transform(rcoords, xform)
+
+ field = np.zeros(list(shape[:3]) + [3])
+ field[:] = (scoords - rcoords).reshape(*it.chain(shape, [3]))
+ field = nonlinear.DeformationField(field, src, ref,
+ srcSpace=srcSpace,
+ refSpace=refSpace,
+ xform=fv2w,
+ header=ref.header,
+ defType='relative')
+ return field
+
+
+def _random_affine_field():
+
+ src = _random_image()
+ ref = _random_image()
+
+ # our test field just encodes an affine
+ xform = affine.compose(
+ np.random.randint(2, 5, 3),
+ np.random.randint(1, 10, 3),
+ np.random.random(3))
+
+ rx, ry, rz = np.meshgrid(np.arange(ref.shape[0]),
+ np.arange(ref.shape[1]),
+ np.arange(ref.shape[2]), indexing='ij')
+
+ rvoxels = np.vstack((rx.flatten(), ry.flatten(), rz.flatten())).T
+ rcoords = affine.transform(rvoxels, ref.voxToScaledVoxMat)
+ scoords = affine.transform(rcoords, xform)
+
+ field = np.zeros(list(ref.shape[:3]) + [3])
+ field[:] = (scoords - rcoords).reshape(*it.chain(ref.shape, [3]))
+ field = nonlinear.DeformationField(field, src, ref,
+ header=ref.header,
+ defType='relative')
+ return field, xform
+
+
+def _field_coords(field):
+ vx, vy, vz = field.shape[ :3]
+ coords = np.meshgrid(np.arange(vx),
+ np.arange(vy),
+ np.arange(vz), indexing='ij')
+ coords = np.array(coords).transpose((1, 2, 3, 0))
+ return affine.transform(
+ coords.reshape(-1, 3),
+ field.getAffine('voxel', 'fsl')).reshape(field.shape)
+
+
+def test_detectDeformationType():
+ relfield = _random_field()
+ coords = _field_coords(relfield)
+ absfield = nonlinear.DeformationField(
+ relfield.data + coords,
+ src=relfield.src,
+ xform=relfield.voxToWorldMat)
+
+ assert nonlinear.detectDeformationType(relfield) == 'relative'
+ assert nonlinear.detectDeformationType(absfield) == 'absolute'
+
+
+def test_convertDeformationType():
+
+ relfield = _random_field()
+ coords = _field_coords(relfield)
+ absfield = nonlinear.DeformationField(
+ relfield.data + coords,
+ src=relfield.src,
+ xform=relfield.voxToWorldMat)
+
+ gotconvrel1 = nonlinear.convertDeformationType(relfield)
+ gotconvabs1 = nonlinear.convertDeformationType(absfield)
+ gotconvrel2 = nonlinear.convertDeformationType(relfield, 'absolute')
+ gotconvabs2 = nonlinear.convertDeformationType(absfield, 'relative')
+
+ tol = dict(atol=1e-5, rtol=1e-5)
+
+ assert np.all(np.isclose(gotconvrel1, absfield.data, **tol))
+ assert np.all(np.isclose(gotconvabs1, relfield.data, **tol))
+ assert np.all(np.isclose(gotconvrel2, absfield.data, **tol))
+ assert np.all(np.isclose(gotconvabs2, relfield.data, **tol))
+
+
+def test_convertDeformationSpace():
+
+ basefield, xform = _random_affine_field()
+ src = basefield.src
+ ref = basefield.ref
+
+ # generate reference fsl->fsl coordinate mappings
+
+ # For each combination of srcspace->tospace
+ # Generate random coordinates, check that
+ # displacements are correct
+ spaces = ['fsl', 'voxel', 'world']
+ spaces = list(it.combinations_with_replacement(spaces, 2))
+ spaces = spaces + [(r, s) for s, r in spaces]
+ spaces = list(set(spaces))
+
+ for from_, to in spaces:
+
+ refcoords = [np.random.randint(0, basefield.shape[0], 5),
+ np.random.randint(0, basefield.shape[1], 5),
+ np.random.randint(0, basefield.shape[2], 5)]
+ refcoords = np.array(refcoords, dtype=np.int).T
+ refcoords = affine.transform(refcoords, ref.voxToScaledVoxMat)
+ srccoords = basefield.transform(refcoords)
+
+ field = nonlinear.convertDeformationSpace(basefield, from_, to)
+ premat = ref.getAffine('fsl', from_)
+ postmat = src.getAffine('fsl', to)
+
+ input = affine.transform(refcoords, premat)
+ expect = affine.transform(srccoords, postmat)
+
+ got = field.transform(input)
+ enan = np.isnan(expect)
+ gnan = np.isnan(got)
+
+ assert np.all(np.isclose(enan, gnan))
+ assert np.all(np.isclose(expect[~enan], got[~gnan]))
+
+
+def test_DeformationField_transform():
+
+ relfield, xform = _random_affine_field()
+ src = relfield.src
+ ref = relfield.ref
+
+ rx, ry, rz = np.meshgrid(np.arange(ref.shape[0]),
+ np.arange(ref.shape[1]),
+ np.arange(ref.shape[2]), indexing='ij')
+ rvoxels = np.vstack((rx.flatten(), ry.flatten(), rz.flatten())).T
+ rcoords = affine.transform(rvoxels, ref.voxToScaledVoxMat)
+ scoords = affine.transform(rcoords, xform)
+ svoxels = affine.transform(scoords, src.scaledVoxToVoxMat)
+
+ absfield = np.zeros(list(ref.shape[:3]) + [3])
+ absfield[:] = scoords.reshape(*it.chain(ref.shape, [3]))
+ absfield = nonlinear.DeformationField(absfield, src, ref,
+ header=ref.header,
+ defType='absolute')
+
+ got = relfield.transform(rcoords)
+ assert np.all(np.isclose(got, scoords))
+ got = absfield.transform(rcoords)
+ assert np.all(np.isclose(got, scoords))
+
+ got = relfield.transform(rvoxels, from_='voxel', to='voxel')
+ assert np.all(np.isclose(got, svoxels))
+ got = absfield.transform(rvoxels, from_='voxel', to='voxel')
+ assert np.all(np.isclose(got, svoxels))
+
+ # test out of bounds are returned as nan
+ rvoxels = np.array([[-1, -1, -1],
+ [ 0, 0, 0]])
+ rcoords = affine.transform(rvoxels, ref.voxToScaledVoxMat)
+ scoords = affine.transform(rcoords, xform)
+ svoxels = affine.transform(scoords, src.scaledVoxToVoxMat)
+
+ got = relfield.transform(rcoords)
+ assert np.all(np.isnan(got[0, :]))
+ assert np.all(np.isclose(got[1, :], scoords[1, :]))
+ got = absfield.transform(rcoords)
+ assert np.all(np.isnan(got[0, :]))
+ assert np.all(np.isclose(got[1, :], scoords[1, :]))
+
+
+def test_CoefficientField_displacements():
+
+ 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')
+ df = op.join(nldir, 'displacementfield_no_premat.nii.gz')
+
+ src = fslimage.Image(src)
+ ref = fslimage.Image(ref)
+ cf = fnirt.readFnirt(cf, src, ref)
+ df = fnirt.readFnirt(df, src, ref)
+
+ ix, iy, iz = ref.shape[:3]
+ x, y, z = np.meshgrid(np.arange(ix),
+ np.arange(iy),
+ np.arange(iz), indexing='ij')
+ x = x.flatten()
+ y = y.flatten()
+ z = z.flatten()
+ xyz = np.vstack((x, y, z)).T
+
+ disps = cf.displacements(xyz)
+ disps = disps.reshape(df.shape)
+
+ tol = dict(atol=1e-5, rtol=1e-5)
+ assert np.all(np.isclose(disps, df.data, **tol))
+
+
+def test_CoefficientField_transform():
+ 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')
+ df = op.join(nldir, 'displacementfield.nii.gz')
+ dfnp = op.join(nldir, 'displacementfield_no_premat.nii.gz')
+
+ src = fslimage.Image(src)
+ ref = fslimage.Image(ref)
+ cf = fnirt.readFnirt(cf, src, ref)
+ df = fnirt.readFnirt(df, src, ref)
+ dfnp = fnirt.readFnirt(dfnp, src, ref)
+
+ spaces = ['fsl', 'voxel', 'world']
+ spaces = list(it.combinations_with_replacement(spaces, 2))
+ spaces = spaces + [(r, s) for s, r in spaces]
+ spaces = list(set(spaces))
+
+ rx, ry, rz = np.meshgrid(np.arange(ref.shape[0]),
+ np.arange(ref.shape[1]),
+ np.arange(ref.shape[2]), indexing='ij')
+ rvoxels = np.vstack((rx.flatten(), ry.flatten(), rz.flatten())).T
+
+ refcoords = {
+ 'voxel' : rvoxels,
+ 'fsl' : affine.transform(rvoxels, ref.getAffine('voxel', 'fsl')),
+ 'world' : affine.transform(rvoxels, ref.getAffine('voxel', 'world'))
+ }
+
+ srccoords = refcoords['fsl'] + df.data.reshape(-1, 3)
+ srccoords = {
+ 'voxel' : affine.transform(srccoords, src.getAffine('fsl', 'voxel')),
+ 'fsl' : srccoords,
+ 'world' : affine.transform(srccoords, src.getAffine('fsl', 'world'))
+ }
+
+ srccoordsnp = refcoords['fsl'] + dfnp.data.reshape(-1, 3)
+ srccoordsnp = {
+ 'voxel' : affine.transform(srccoordsnp, src.getAffine('fsl', 'voxel')),
+ 'fsl' : srccoordsnp,
+ 'world' : affine.transform(srccoordsnp, src.getAffine('fsl', 'world'))
+ }
+
+ tol = dict(atol=1e-5, rtol=1e-5)
+ for srcspace, refspace in spaces:
+ got = cf.transform(refcoords[refspace], refspace, srcspace)
+ gotnp = cf.transform(refcoords[refspace], refspace, srcspace,
+ premat=False)
+ assert np.all(np.isclose(got, srccoords[ srcspace], **tol))
+ assert np.all(np.isclose(gotnp, srccoordsnp[srcspace], **tol))
+
+
+def test_coefficientFieldToDeformationField():
+
+ 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')
+ df = op.join(nldir, 'displacementfield.nii.gz')
+ dfnp = op.join(nldir, 'displacementfield_no_premat.nii.gz')
+
+ src = fslimage.Image(src)
+ ref = fslimage.Image(ref)
+ cf = fnirt.readFnirt(cf, src, ref)
+ rdf = fnirt.readFnirt(df, src, ref)
+ rdfnp = fnirt.readFnirt(dfnp, src, ref)
+ adf = nonlinear.convertDeformationType(rdf)
+ adfnp = nonlinear.convertDeformationType(rdfnp)
+
+ rcnv = nonlinear.coefficientFieldToDeformationField(cf)
+ acnv = nonlinear.coefficientFieldToDeformationField(cf,
+ defType='absolute')
+ acnvnp = nonlinear.coefficientFieldToDeformationField(cf,
+ defType='absolute',
+ premat=False)
+ rcnvnp = nonlinear.coefficientFieldToDeformationField(cf,
+ premat=False)
+
+ tol = dict(atol=1e-5, rtol=1e-5)
+ assert np.all(np.isclose(rcnv.data, rdf .data, **tol))
+ assert np.all(np.isclose(acnv.data, adf .data, **tol))
+ assert np.all(np.isclose(rcnvnp.data, rdfnp.data, **tol))
+ assert np.all(np.isclose(acnvnp.data, adfnp.data, **tol))
+
+
+def test_applyDeformation():
+
+ src2ref = affine.compose(
+ np.random.randint(2, 5, 3),
+ np.random.randint(1, 10, 3),
+ np.random.random(3))
+ ref2src = affine.invert(src2ref)
+
+ srcdata = np.random.randint(1, 65536, (10, 10, 10))
+ refdata = np.random.randint(1, 65536, (10, 10, 10))
+
+ src = fslimage.Image(srcdata)
+ ref = fslimage.Image(refdata, xform=src2ref)
+ field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+ expect, xf = resample.resampleToReference(
+ src, ref, matrix=src2ref, order=1, mode='nearest')
+ result = nonlinear.applyDeformation(
+ src, field, order=1, mode='nearest')
+
+ assert np.all(np.isclose(expect, result))
+
+
+def test_applyDeformation_altsrc():
+
+ src2ref = affine.compose(
+ np.random.randint(2, 5, 3),
+ np.random.randint(1, 10, 3),
+ [0, 0, 0])
+ ref2src = affine.invert(src2ref)
+
+ srcdata = np.random.randint(1, 65536, (10, 10, 10))
+ refdata = np.random.randint(1, 65536, (10, 10, 10))
+
+ src = fslimage.Image(srcdata)
+ ref = fslimage.Image(refdata, xform=src2ref)
+ field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+ # First try a down-sampled version
+ # of the original source image
+ altsrc, xf = resample.resample(src, (5, 5, 5), origin='corner')
+ altsrc = fslimage.Image(altsrc, xform=xf, header=src.header)
+ expect, xf = resample.resampleToReference(
+ altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+ result = nonlinear.applyDeformation(
+ altsrc, field, order=1, mode='nearest')
+ assert np.all(np.isclose(expect, result))
+
+ # Now try a down-sampled ROI
+ # of the original source image
+ altsrc = roi.roi(src, [(2, 9), (2, 9), (2, 9)])
+ altsrc, xf = resample.resample(altsrc, (4, 4, 4))
+ altsrc = fslimage.Image(altsrc, xform=xf, header=src.header)
+ expect, xf = resample.resampleToReference(
+ altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+ result = nonlinear.applyDeformation(
+ altsrc, field, order=1, mode='nearest')
+ assert np.all(np.isclose(expect, result))
+
+ # down-sampled and offset ROI
+ # of the original source image
+ altsrc = roi.roi(src, [(-5, 8), (-5, 8), (-5, 8)])
+ altsrc, xf = resample.resample(altsrc, (6, 6, 6))
+ altsrc = fslimage.Image(altsrc, xform=xf, header=src.header)
+ expect, xf = resample.resampleToReference(
+ altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+ result = nonlinear.applyDeformation(
+ altsrc, field, order=1, mode='nearest')
+ assert np.all(np.isclose(expect, result))
+
+
+def test_applyDeformation_altref():
+ src2ref = affine.compose(
+ np.random.randint(2, 5, 3),
+ np.random.randint(1, 10, 3),
+ np.random.random(3))
+ ref2src = affine.invert(src2ref)
+
+ srcdata = np.random.randint(1, 65536, (10, 10, 10))
+ refdata = np.random.randint(1, 65536, (10, 10, 10))
+
+ src = fslimage.Image(srcdata)
+ ref = fslimage.Image(refdata, xform=src2ref)
+ field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+ altrefxform = affine.concat(
+ src2ref,
+ affine.scaleOffsetXform([1, 1, 1], [5, 0, 0]))
+
+ altref = fslimage.Image(refdata, xform=altrefxform)
+
+ expect, xf = resample.resampleToReference(
+ src, altref, matrix=src2ref, order=1, mode='constant', cval=0)
+ result = nonlinear.applyDeformation(
+ src, field, ref=altref, order=1, mode='constant', cval=0)
+
+ # boundary voxels can get truncated
+ # (4 is the altref-ref overlap boundary)
+ expect[4, :, :] = 0
+ result[4, :, :] = 0
+ expect = expect[1:-1, 1:-1, 1:-1]
+ result = result[1:-1, 1:-1, 1:-1]
+
+ assert np.all(np.isclose(expect, result))
+
+
+# test when reference/field
+# are not voxel-aligned
+def test_applyDeformation_worldAligned():
+ refv2w = affine.scaleOffsetXform([1, 1, 1], [10, 10, 10])
+ fieldv2w = affine.scaleOffsetXform([2, 2, 2], [10.5, 10.5, 10.5])
+ src2ref = refv2w
+ ref2src = affine.invert(src2ref)
+
+ srcdata = np.random.randint(1, 65536, (10, 10, 10))
+
+ src = fslimage.Image(srcdata)
+ ref = fslimage.Image(srcdata, xform=src2ref)
+ field = _affine_field(src, ref, ref2src, 'world', 'world',
+ shape=(5, 5, 5), fv2w=fieldv2w)
+
+ field = nonlinear.DeformationField(
+ nonlinear.convertDeformationType(field, 'absolute'),
+ header=field.header,
+ src=src,
+ ref=ref,
+ srcSpace='world',
+ refSpace='world',
+ defType='absolute')
+
+ expect, xf = resample.resampleToReference(
+ src, ref, matrix=src2ref, order=1, mode='constant', cval=0)
+ result = nonlinear.applyDeformation(
+ src, field, order=1, mode='constant', cval=0)
+
+ expect = expect[1:-1, 1:-1, 1:-1]
+ result = result[1:-1, 1:-1, 1:-1]
+
+ assert np.all(np.isclose(expect, result))
diff --git a/tests/test_transform/test_x5.py b/tests/test_transform/test_x5.py
new file mode 100644
index 0000000000000000000000000000000000000000..3a92ece655197f3fa77051c858b9fcd03bc7bb57
--- /dev/null
+++ b/tests/test_transform/test_x5.py
@@ -0,0 +1,112 @@
+#!/usr/bin/env python
+#
+# test_x5.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+import os.path as op
+import numpy as np
+
+import pytest
+
+import h5py
+
+import fsl.data.image as fslimage
+import fsl.utils.tempdir as tempdir
+import fsl.transform.affine as affine
+import fsl.transform.fnirt as fnirt
+import fsl.transform.nonlinear as nonlinear
+import fsl.transform.x5 as x5
+
+from .. import make_random_image
+
+
+def _check_metadata(group):
+ assert group.attrs['Format'] == x5.X5_FORMAT
+ assert group.attrs['Version'] == x5.X5_VERSION
+
+
+def _check_affine(group, xform):
+ assert group.attrs['Type'] == 'affine'
+ gotxform = np.array(group['Matrix'])
+ assert np.all(np.isclose(gotxform, xform))
+
+
+def _check_space(group, img):
+ assert group.attrs['Type'] == 'image'
+ assert np.all(np.isclose(group.attrs['Size'], img.shape[ :3]))
+ assert np.all(np.isclose(group.attrs['Scales'], img.pixdim[:3]))
+ _check_affine(group['Mapping'], img.voxToWorldMat)
+
+
+def _check_deformation(group, field):
+ assert group.attrs['Type'] == 'deformation'
+ assert group.attrs['SubType'] == field.deformationType
+ xform = np.array(group['Matrix'])
+ assert np.all(np.isclose(xform, field.data))
+ _check_affine(group['Mapping'], field.voxToWorldMat)
+
+
+def test_readWriteLinearX5():
+ with tempdir.tempdir():
+ make_random_image('src.nii')
+ make_random_image('ref.nii')
+ xform = affine.compose(
+ np.random.randint(1, 5, 3),
+ np.random.randint(-10, 10, 3),
+ -np.pi / 4 + np.random.random(3) * np.pi / 2)
+
+ src = fslimage.Image('src.nii')
+ ref = fslimage.Image('ref.nii')
+
+ x5.writeLinearX5('linear.x5', xform, src, ref)
+
+ gotxform, gotsrc, gotref = x5.readLinearX5('linear.x5')
+ assert np.all(np.isclose(gotxform, xform))
+ assert gotsrc.sameSpace(src)
+ assert gotref.sameSpace(ref)
+
+ with h5py.File('linear.x5', 'r') as f:
+ _check_metadata(f)
+ assert f.attrs['Type'] == 'linear'
+ _check_affine(f['/Transform'], xform)
+ _check_space( f['/A'], src)
+ _check_space( f['/B'], ref)
+
+
+def test_readWriteNonLinearX5():
+ datadir = op.join(op.dirname(__file__), 'testdata', 'nonlinear')
+ dffile = op.join(datadir, 'displacementfield.nii.gz')
+ srcfile = op.join(datadir, 'src.nii.gz')
+ reffile = op.join(datadir, 'ref.nii.gz')
+
+ src = fslimage.Image(srcfile)
+ ref = fslimage.Image(reffile)
+ dfield = fnirt.readFnirt(dffile, src, ref)
+ wdfield = nonlinear.convertDeformationSpace(dfield, 'world', 'world')
+
+ with tempdir.tempdir():
+
+ # field must be world->world
+ with pytest.raises(x5.X5Error):
+ x5.writeNonLinearX5('nonlinear.x5', dfield)
+
+ x5.writeNonLinearX5('nonlinear.x5', wdfield)
+
+ gotdfield = x5.readNonLinearX5('nonlinear.x5')
+
+ assert gotdfield.src.sameSpace(src)
+ assert gotdfield.ref.sameSpace(ref)
+ assert gotdfield.srcSpace == wdfield.srcSpace
+ assert gotdfield.refSpace == wdfield.refSpace
+ assert gotdfield.deformationType == wdfield.deformationType
+ assert np.all(np.isclose(gotdfield.data, wdfield.data))
+
+ with h5py.File('nonlinear.x5', 'r') as f:
+ assert f.attrs['Type'] == 'nonlinear'
+ _check_metadata(f)
+ _check_deformation(f['/Transform'], wdfield)
+ _check_space( f['/A'], ref)
+ _check_space( f['/B'], src)
diff --git a/tests/test_transform/testdata/nonlinear/coefficientfield.nii.gz b/tests/test_transform/testdata/nonlinear/coefficientfield.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..8d4b58668d2f100819ccf6263a1e2c77286b1bd0
Binary files /dev/null and b/tests/test_transform/testdata/nonlinear/coefficientfield.nii.gz differ
diff --git a/tests/test_transform/testdata/nonlinear/displacementfield.nii.gz b/tests/test_transform/testdata/nonlinear/displacementfield.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..f1e28a1bd1c072de4720add4401699b8dcb3a62e
Binary files /dev/null and b/tests/test_transform/testdata/nonlinear/displacementfield.nii.gz differ
diff --git a/tests/test_transform/testdata/nonlinear/displacementfield_no_premat.nii.gz b/tests/test_transform/testdata/nonlinear/displacementfield_no_premat.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..26578f9dba2911634789e2d8d2ebac8ad8d5224f
Binary files /dev/null and b/tests/test_transform/testdata/nonlinear/displacementfield_no_premat.nii.gz differ
diff --git a/tests/test_transform/testdata/nonlinear/ref.nii.gz b/tests/test_transform/testdata/nonlinear/ref.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..ee4f0a98a8cd4cc12a5d1039b36a73419fdcffdc
Binary files /dev/null and b/tests/test_transform/testdata/nonlinear/ref.nii.gz differ
diff --git a/tests/test_transform/testdata/nonlinear/src.nii.gz b/tests/test_transform/testdata/nonlinear/src.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..01b6290f590a8d5110e71978e4dbbfe90fc475ed
Binary files /dev/null and b/tests/test_transform/testdata/nonlinear/src.nii.gz differ
diff --git a/tests/test_transform/testdata/nonlinear/src2ref.mat b/tests/test_transform/testdata/nonlinear/src2ref.mat
new file mode 100644
index 0000000000000000000000000000000000000000..521dd0b7411a1fe958d9d4feb910be9442647f68
--- /dev/null
+++ b/tests/test_transform/testdata/nonlinear/src2ref.mat
@@ -0,0 +1,4 @@
+1.130654884 -0.02090760163 0.05273442572 -0.523881946
+0.01777234539 0.9990179929 0.2039857328 5.806280637
+-0.003319001082 -0.2740303429 1.071017063 0.4915523243
+0 0 0 1
diff --git a/tests/testdata/test_transform_test_axisBounds.txt b/tests/test_transform/testdata/test_transform_test_axisBounds.txt
similarity index 100%
rename from tests/testdata/test_transform_test_axisBounds.txt
rename to tests/test_transform/testdata/test_transform_test_axisBounds.txt
diff --git a/tests/testdata/test_transform_test_compose.txt b/tests/test_transform/testdata/test_transform_test_compose.txt
similarity index 100%
rename from tests/testdata/test_transform_test_compose.txt
rename to tests/test_transform/testdata/test_transform_test_compose.txt
diff --git a/tests/testdata/test_transform_test_concat.txt b/tests/test_transform/testdata/test_transform_test_concat.txt
similarity index 100%
rename from tests/testdata/test_transform_test_concat.txt
rename to tests/test_transform/testdata/test_transform_test_concat.txt
diff --git a/tests/testdata/test_transform_test_flirtMatrixToSform.txt b/tests/test_transform/testdata/test_transform_test_flirtMatrixToSform.txt
similarity index 100%
rename from tests/testdata/test_transform_test_flirtMatrixToSform.txt
rename to tests/test_transform/testdata/test_transform_test_flirtMatrixToSform.txt
diff --git a/tests/testdata/test_transform_test_invert.txt b/tests/test_transform/testdata/test_transform_test_invert.txt
similarity index 100%
rename from tests/testdata/test_transform_test_invert.txt
rename to tests/test_transform/testdata/test_transform_test_invert.txt
diff --git a/tests/testdata/test_transform_test_scaleoffsetxform.txt b/tests/test_transform/testdata/test_transform_test_scaleoffsetxform.txt
similarity index 100%
rename from tests/testdata/test_transform_test_scaleoffsetxform.txt
rename to tests/test_transform/testdata/test_transform_test_scaleoffsetxform.txt
diff --git a/tests/testdata/test_transform_test_transform_00.txt b/tests/test_transform/testdata/test_transform_test_transform_00.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_00.txt
rename to tests/test_transform/testdata/test_transform_test_transform_00.txt
diff --git a/tests/testdata/test_transform_test_transform_01.txt b/tests/test_transform/testdata/test_transform_test_transform_01.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_01.txt
rename to tests/test_transform/testdata/test_transform_test_transform_01.txt
diff --git a/tests/testdata/test_transform_test_transform_02.txt b/tests/test_transform/testdata/test_transform_test_transform_02.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_02.txt
rename to tests/test_transform/testdata/test_transform_test_transform_02.txt
diff --git a/tests/testdata/test_transform_test_transform_03.txt b/tests/test_transform/testdata/test_transform_test_transform_03.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_03.txt
rename to tests/test_transform/testdata/test_transform_test_transform_03.txt
diff --git a/tests/testdata/test_transform_test_transform_04.txt b/tests/test_transform/testdata/test_transform_test_transform_04.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_04.txt
rename to tests/test_transform/testdata/test_transform_test_transform_04.txt
diff --git a/tests/testdata/test_transform_test_transform_05.txt b/tests/test_transform/testdata/test_transform_test_transform_05.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_05.txt
rename to tests/test_transform/testdata/test_transform_test_transform_05.txt
diff --git a/tests/testdata/test_transform_test_transform_06.txt b/tests/test_transform/testdata/test_transform_test_transform_06.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_06.txt
rename to tests/test_transform/testdata/test_transform_test_transform_06.txt
diff --git a/tests/testdata/test_transform_test_transform_07.txt b/tests/test_transform/testdata/test_transform_test_transform_07.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_07.txt
rename to tests/test_transform/testdata/test_transform_test_transform_07.txt
diff --git a/tests/testdata/test_transform_test_transform_08.txt b/tests/test_transform/testdata/test_transform_test_transform_08.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_08.txt
rename to tests/test_transform/testdata/test_transform_test_transform_08.txt
diff --git a/tests/testdata/test_transform_test_transform_09.txt b/tests/test_transform/testdata/test_transform_test_transform_09.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_09.txt
rename to tests/test_transform/testdata/test_transform_test_transform_09.txt
diff --git a/tests/testdata/test_transform_test_transform_10.txt b/tests/test_transform/testdata/test_transform_test_transform_10.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_10.txt
rename to tests/test_transform/testdata/test_transform_test_transform_10.txt
diff --git a/tests/testdata/test_transform_test_transform_11.txt b/tests/test_transform/testdata/test_transform_test_transform_11.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_11.txt
rename to tests/test_transform/testdata/test_transform_test_transform_11.txt
diff --git a/tests/testdata/test_transform_test_transform_12.txt b/tests/test_transform/testdata/test_transform_test_transform_12.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_12.txt
rename to tests/test_transform/testdata/test_transform_test_transform_12.txt
diff --git a/tests/testdata/test_transform_test_transform_13.txt b/tests/test_transform/testdata/test_transform_test_transform_13.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_13.txt
rename to tests/test_transform/testdata/test_transform_test_transform_13.txt
diff --git a/tests/testdata/test_transform_test_transform_14.txt b/tests/test_transform/testdata/test_transform_test_transform_14.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_14.txt
rename to tests/test_transform/testdata/test_transform_test_transform_14.txt
diff --git a/tests/testdata/test_transform_test_transform_15.txt b/tests/test_transform/testdata/test_transform_test_transform_15.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_15.txt
rename to tests/test_transform/testdata/test_transform_test_transform_15.txt
diff --git a/tests/testdata/test_transform_test_transform_16.txt b/tests/test_transform/testdata/test_transform_test_transform_16.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_16.txt
rename to tests/test_transform/testdata/test_transform_test_transform_16.txt
diff --git a/tests/testdata/test_transform_test_transform_17.txt b/tests/test_transform/testdata/test_transform_test_transform_17.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_17.txt
rename to tests/test_transform/testdata/test_transform_test_transform_17.txt
diff --git a/tests/testdata/test_transform_test_transform_18.txt b/tests/test_transform/testdata/test_transform_test_transform_18.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_18.txt
rename to tests/test_transform/testdata/test_transform_test_transform_18.txt
diff --git a/tests/testdata/test_transform_test_transform_19.txt b/tests/test_transform/testdata/test_transform_test_transform_19.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_19.txt
rename to tests/test_transform/testdata/test_transform_test_transform_19.txt
diff --git a/tests/testdata/test_transform_test_transform_20.txt b/tests/test_transform/testdata/test_transform_test_transform_20.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_20.txt
rename to tests/test_transform/testdata/test_transform_test_transform_20.txt
diff --git a/tests/testdata/test_transform_test_transform_coords.txt b/tests/test_transform/testdata/test_transform_test_transform_coords.txt
similarity index 100%
rename from tests/testdata/test_transform_test_transform_coords.txt
rename to tests/test_transform/testdata/test_transform_test_transform_coords.txt