Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • paulmc/fslpy
  • ndcn0236/fslpy
  • seanf/fslpy
3 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1809 additions and 986 deletions
fsl/transform/__init__.py
View file @ 3b2a186e
......@@ -11,9 +11,9 @@ transformations. Functionality is split across the following modules:
.. autosummary::
:nosignatures:
.affine
.flirt
.fnirt
.nonlinear
.x5
~fsl.transform.affine
~fsl.transform.flirt
~fsl.transform.fnirt
~fsl.transform.nonlinear
~fsl.transform.x5
"""
fsl/transform/affine.py
View file @ 3b2a186e
......@@ -13,6 +13,7 @@ transformations. The following functions are available:
transform
scaleOffsetXform
invert
flip
concat
compose
decompose
......@@ -20,7 +21,9 @@ transformations. The following functions are available:
rotMatToAxisAngles
axisAnglesToRotMat
axisBounds
mergeBounds
rmsdev
rescale
And a few more functions are provided for working with vectors:
......@@ -59,7 +62,7 @@ def veclength(vec):
Multiple vectors may be passed in, with a shape of ``(n, 3)``.
"""
vec = np.array(vec, copy=False).reshape(-1, 3)
vec = np.asarray(vec).reshape(-1, 3)
return np.sqrt(np.einsum('ij,ij->i', vec, vec))
......@@ -68,7 +71,7 @@ def normalise(vec):
Multiple vectors may be passed in, with a shape of ``(n, 3)``.
"""
vec = np.array(vec, copy=False).reshape(-1, 3)
vec = np.asarray(vec).reshape(-1, 3)
n = (vec.T / veclength(vec)).T
if n.size == 3:
......@@ -77,7 +80,41 @@ def normalise(vec):
return n
def scaleOffsetXform(scales, offsets):
def flip(shape, xform, *axes):
"""Applies a flip/inversion to an affine transform along specified axes.
:arg shape: The ``(x, y, z)`` shape of the data.
:arg xform: Transformation matrix which transforms voxel coordinates
to a world coordinate system.
:arg axes: Indices of axes to flip
"""
if len(axes) == 0:
return xform
axes = sorted(set(axes))
if any(a not in (0, 1, 2) for a in axes):
raise ValueError(f'Invalid axis: {axes}')
los, his = axisBounds(shape, xform, axes, boundary=None)
scales = [1, 1, 1]
pre = [0, 0, 0]
post = [0, 0, 0]
for ax, lo, hi in zip(axes, los, his):
mid = lo + (hi - lo) / 2
scales[ax] = -1
pre[ ax] = -mid
post[ ax] = mid
pre = scaleOffsetXform(None, pre)
post = scaleOffsetXform(None, post)
scales = scaleOffsetXform(scales)
return concat(post, scales, pre, xform)
def scaleOffsetXform(scales=None, offsets=None):
"""Creates and returns an affine transformation matrix which encodes
the specified scale(s) and offset(s).
......@@ -93,6 +130,9 @@ def scaleOffsetXform(scales, offsets):
:returns: A ``numpy.float32`` array of size :math:`4 \\times 4`.
"""
if scales is None: scales = [1, 1, 1]
if offsets is None: offsets = [0, 0, 0]
oktypes = (abc.Sequence, np.ndarray)
if not isinstance(scales, oktypes): scales = [scales]
......@@ -119,19 +159,28 @@ def scaleOffsetXform(scales, offsets):
return xform
def compose(scales, offsets, rotations, origin=None):
def compose(scales, offsets, rotations, origin=None, shears=None,
scaleAtOrigin=False):
"""Compose a transformation matrix out of the given scales, offsets
and axis rotations.
:arg scales: Sequence of three scale values.
:arg scales: Sequence of three scale values.
:arg offsets: Sequence of three offset 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 rotations: Sequence of three rotation values, in radians, or
a rotation matrix of shape ``(3, 3)``.
:arg origin: Origin of rotation. If not provided, the rotation
origin is ``(0, 0, 0)``.
:arg origin: Origin of rotation - must be scaled by the ``scales``.
If not provided, the rotation origin is ``(0, 0, 0)``.
:arg shears: Sequence of three shear values
:arg scaleAtOrigin: If ``True``, the scaling parameters are applied with
respect to the ``origin``, i.e. so that the location
of ``origin`` is unchanged after scaling. If
``False`` (the default), the scaling parameters are
applied with respect to location ``(0, 0, 0)``.
"""
preRotate = np.eye(4)
......@@ -153,6 +202,7 @@ def compose(scales, offsets, rotations, origin=None):
scale = np.eye(4, dtype=np.float64)
offset = np.eye(4, dtype=np.float64)
rotate = np.eye(4, dtype=np.float64)
shear = np.eye(4, dtype=np.float64)
scale[ 0, 0] = scales[ 0]
scale[ 1, 1] = scales[ 1]
......@@ -163,10 +213,18 @@ def compose(scales, offsets, rotations, origin=None):
rotate[:3, :3] = rotations
return concat(offset, postRotate, rotate, preRotate, scale)
if shears is not None:
shear[0, 1] = shears[0]
shear[0, 2] = shears[1]
shear[1, 2] = shears[2]
if scaleAtOrigin:
return concat(offset, postRotate, rotate, scale, preRotate, shear)
else:
return concat(offset, postRotate, rotate, preRotate, scale, shear)
def decompose(xform, angles=True):
def decompose(xform, angles=True, shears=False):
"""Decomposes the given transformation matrix into separate offsets,
scales, and rotations, according to the algorithm described in:
......@@ -174,8 +232,7 @@ def decompose(xform, angles=True):
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.
It is assumed that the given transform has no perspective components.
:arg xform: A ``(3, 3)`` or ``(4, 4)`` affine transformation matrix.
......@@ -183,6 +240,8 @@ def decompose(xform, angles=True):
as axis-angles, in radians. Otherwise, the rotation matrix
is returned.
:arg shears: Defaults to ``False``. If ``True``, shears are returned.
:returns: The following:
- A sequence of three scales
......@@ -190,6 +249,7 @@ def decompose(xform, angles=True):
was a ``(3, 3)`` matrix)
- A sequence of three rotations, in radians. Or, if
``angles is False``, a rotation matrix.
- If ``shears is True``, a sequence of three shears.
"""
# The inline comments in the code below are taken verbatim from
......@@ -198,7 +258,7 @@ def decompose(xform, angles=True):
# 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
xform = np.array(xform).T
if xform.shape == (4, 4):
translations = xform[ 3, :3]
......@@ -213,6 +273,7 @@ def decompose(xform, angles=True):
# 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))
M1 = M1 / sx
# 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).
......@@ -229,7 +290,7 @@ def decompose(xform, angles=True):
# The second row is normalized, and s_xy is divided by s_y to
# get its final value.
M2 = M2 / sy
sxy = sxy / sy
sxy = sxy / sx
# The xz and yz shear factors are computed as in the preceding,
sxz = np.dot(M1, M3)
......@@ -244,8 +305,8 @@ def decompose(xform, angles=True):
# the third row is normalized, and the xz and yz shear factors are
# rescaled.
M3 = M3 / sz
sxz = sxz / sz
syz = sxz / sz
sxz = sxz / sx
syz = syz / sy
# 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
......@@ -265,7 +326,13 @@ def decompose(xform, angles=True):
if angles: rotations = rotMatToAxisAngles(R.T)
else: rotations = R.T
return np.array([sx, sy, sz]), translations, rotations
retval = [np.array([sx, sy, sz]), translations, rotations]
if shears:
retval.append(np.array((sxy, sxz, syz)))
return tuple(retval)
def rotMatToAffine(rotmat, origin=None):
......@@ -450,11 +517,45 @@ def axisBounds(shape,
else: return (lo, hi)
def mergeBounds(*bounds):
"""Merge multiple ``(lo, hi)`` bounds together. Returns the union of
all bounds, i.e. a new set of bounds which encompasses all given bounds.
Each set of bounds can be provided in one of the following formats:
- ``((lo, lo, lo), (hi, hi, hi))``
- ``((lo, hi), (lo, hi), (lo, hi))``
The return value will be a tuple of tuples, otherwise having the same
format as the inputs.
"""
nbounds = len(bounds)
bounds = np.array(bounds)
if bounds.shape not in ((nbounds, 3, 2), (nbounds, 2, 3)):
raise ValueError('Unsupported bounds format')
if bounds.shape == (nbounds, 2, 3):
return ((bounds[:, 0, 0].min(),
bounds[:, 0, 1].min(),
bounds[:, 0, 2].min()),
(bounds[:, 1, 0].max(),
bounds[:, 1, 1].max(),
bounds[:, 1, 2].max()))
else:
return ((bounds[:, 0, 0].min(),
bounds[:, 0, 1].max()),
(bounds[:, 1, 0].min(),
bounds[:, 1, 1].max()),
(bounds[:, 2, 0].min(),
bounds[:, 2, 1].max()))
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
......@@ -582,3 +683,59 @@ def rmsdev(T1, T2, R=None, xc=None):
erms = np.sqrt(erms)
return erms
def rescale(oldShape, newShape, origin=None):
"""Calculates an affine matrix to use for resampling.
This function generates an affine transformation matrix that can be used
to resample an N-D array from ``oldShape`` to ``newShape`` using, for
example, ``scipy.ndimage.affine_transform``.
The matrix will contain scaling factors derived from the ``oldShape /
newShape`` ratio, and an offset determined by the ``origin``.
The default value for ``origin`` (``'centre'``) causes the corner voxel of
the output to have the same centre as the corner voxel of the input. If
the origin is ``'corner'``, we apply an offset which effectively causes
the voxel grid corners of the input and output to be aligned.
:arg oldShape: Shape of input data
:arg newShape: Shape to resample data to
:arg origin: Voxel grid alignment - either ``'centre'`` (the default) or
``'corner'``
:returns: An affine resampling matrix
"""
if origin is None:
origin = 'centre'
oldShape = np.array(oldShape, dtype=float)
newShape = np.array(newShape, dtype=float)
ndim = len(oldShape)
if len(oldShape) != len(newShape):
raise ValueError('Shape mismatch')
# shapes are the same - no rescaling needed
if np.all(np.isclose(oldShape, newShape)):
return np.eye(ndim + 1)
# Otherwise we calculate a scaling
# matrix from the old/new shape
# ratio, and specify an offset
# according to the origin
ratio = oldShape / newShape
scale = np.diag(ratio)
# Calculate an offset from the origin
if origin == 'centre': offset = [0] * ndim
elif origin == 'corner': offset = (ratio - 1) / 2
# combine the scales and translations
# to form thte final affine
xform = np.eye(ndim + 1)
xform[:ndim, :ndim] = scale
xform[:ndim, -1] = offset
return xform
fsl/transform/flirt.py
View file @ 3b2a186e
......@@ -4,7 +4,7 @@
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module contains functions for working with FLIRT affime transformation
"""This module contains functions for working with FLIRT affine transformation
matrices. The following functions are available:
.. autosummary::
......@@ -16,6 +16,18 @@ matrices. The following functions are available:
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.
"""
......@@ -51,9 +63,8 @@ def fromFlirt(xform, src, ref, from_='voxel', to='world'):
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)
- ``fsl``: The FSL coordinate 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
......@@ -64,7 +75,7 @@ def fromFlirt(xform, src, ref, from_='voxel', to='world'):
: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 target 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.
......@@ -97,13 +108,8 @@ def flirtMatrixToSform(flirtMat, srcImage, refImage):
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:
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)
......@@ -121,10 +127,10 @@ 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.
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`
......@@ -132,15 +138,15 @@ def sformToFlirtMatrix(srcImage, refImage, srcXform=None):
:attr:`.Nifti.voxToWorldMat`
"""
srcScaledVoxToVoxMat = srcImage.scaledVoxToVoxMat
srcVoxToWorldMat = srcImage.voxToWorldMat
refWorldToVoxMat = refImage.worldToVoxMat
refVoxToScaledVoxMat = refImage.voxToScaledVoxMat
srcFSLToVoxMat = srcImage.getAffine('fsl', 'voxel')
srcVoxToWorldMat = srcImage.getAffine('voxel', 'world')
refWorldToVoxMat = refImage.getAffine('world', 'voxel')
refVoxToFSLMat = refImage.getAffine('voxel', 'fsl')
if srcXform is not None:
srcVoxToWorldMat = srcXform
return concat(refVoxToScaledVoxMat,
return concat(refVoxToFSLMat,
refWorldToVoxMat,
srcVoxToWorldMat,
srcScaledVoxToVoxMat)
srcFSLToVoxMat)
fsl/transform/fnirt.py
View file @ 3b2a186e
......@@ -5,7 +5,7 @@
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module contains functions for working with FNIRT non-linear
transformation matrices. The following functions are available:
transformations. The following functions are available:
.. autosummary::
:nosignatures:
......@@ -13,38 +13,157 @@ transformation matrices. The following functions are available:
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.
.. note:: FNIRT deformation field files give no indication as to whether they
contain relative displacements or absolute coordinates, so heuristics
must be used to infer what is stored in a a particular file. The
:func:`.nonlinear.detectDeformationType` function can be used to
determine whether a file contains relative displacements or absolute
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 numpy as np
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 _readFnirtDisplacementField(fname, img, src, ref, dispType=None):
"""Loads ``fname``, assumed to be a FNIRT displacement field.
def _readFnirtDeformationField(fname, img, src, ref, defType=None):
"""Loads ``fname``, assumed to be a FNIRT deformation field.
:arg fname: File name of FNIRT displacement field
:arg img: ``fname`` loaded as an :class:`.Image`
:arg src: Source image
:arg ref: Reference image
:arg dispType: Displacement type - either ``'absolute'`` or ``'relative'``.
If not provided, is automatically inferred from the data.
:return: A :class:`.DisplacementField`
: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
"""
from . import nonlinear
return nonlinear.DisplacementField(fname,
src,
ref,
srcSpace='fsl',
refSpace='fsl',
dispType=dispType)
return nonlinear.DeformationField(fname,
src,
ref,
srcSpace='fsl',
refSpace='fsl',
defType=defType)
def _readFnirtCoefficientField(fname, img, src, ref):
......@@ -57,9 +176,6 @@ def _readFnirtCoefficientField(fname, img, src, ref):
:return: A :class:`.CoefficientField`
"""
from . import affine
from . import nonlinear
# FNIRT uses NIFTI header fields in
# non-standard ways to store some
# additional information about the
......@@ -109,6 +225,27 @@ def _readFnirtCoefficientField(fname, img, src, ref):
# image voxel coordinates.
fieldToRefMat = affine.scaleOffsetXform(knotSpacing, 0)
# But if the provided reference has
# different resolution to the
# reference that was originally
# used to generate the warp field,
# we need to adjust the field
# accordingly. We assume that the
# references are aligned in the FSL
# coordinate system, so simply apply
# a scaling factor calculated by
# dividing the original reference
# pixdims by the provided reference
# pixdims.
refPixdims = np.array([img.header['intent_p1'],
img.header['intent_p2'],
img.header['intent_p3']])
if not np.all(np.isclose(ref.pixdim[:3], refPixdims)):
fieldToRefMat = affine.concat(
affine.scaleOffsetXform(refPixdims / ref.pixdim[:3], 0),
fieldToRefMat)
return nonlinear.CoefficientField(fname,
src,
ref,
......@@ -120,52 +257,60 @@ def _readFnirtCoefficientField(fname, img, src, ref):
fieldToRefMat=fieldToRefMat)
def readFnirt(fname, src, ref, dispType=None):
def readFnirt(fname, src, ref, defType=None, intent=None):
"""Reads a non-linear FNIRT transformation image, returning
a :class:`.DisplacementField` or :class:`.CoefficientField` depending
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 dispType: Displacement type - either ``'absolute'`` or ``'relative'``.
If not provided, is automatically inferred from the data.
:arg defType: Deformation type - either ``'absolute'`` or ``'relative'``.
Only used if the file is a deformation field. If not
provided, is automatically inferred from the data.
:arg intent: NIFTI intent code of ``fname``. e.g.
:attr:`.constants.FSL_FNIRT_DISPLACEMENT_FIELD`. If not
provided, the intent is read from the image header.
"""
# Figure out whether the file
# is a displacement field or
# a coefficient field
import fsl.data.image as fslimage
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 _readFnirtDisplacementField(fname, img, src, ref, dispType)
elif img.intent in coefs:
if defType not in (None, 'absolute', 'relative'):
raise ValueError('defType must be None, "absolute" or "relative" '
'(passed in as {})'.format(defType))
# Figure out whether the file is a
# deformation field or a coefficient
# field by checking the intent code.
# If the intent is provided, assume
# that the caller knows the type of
# the field.
img = fslimage.Image(fname)
intent = intent or img.intent
disps = (constants.FSL_FNIRT_DISPLACEMENT_FIELD,
constants.FSL_TOPUP_FIELD)
coefs = (constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
constants.FSL_DCT_COEFFICIENTS,
constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS,
constants.FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS,
constants.FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS)
if intent in disps:
return _readFnirtDeformationField(fname, img, src, ref, defType)
elif intent in coefs:
return _readFnirtCoefficientField(fname, img, src, ref)
else:
raise ValueError('Cannot determine type of nonlinear '
'file {}'.format(fname))
raise ValueError('Cannot determine type of nonlinear warp field '
'{} (intent code: {})'.format(fname, intent))
def toFnirt(field):
"""Convert a :class:`.NonLinearTransform` to a FNIRT-compatible
:class:`.DisplacementField` or:class:`.CoefficientField`.
:class:`.DeformationField` or :class:`.CoefficientField`.
:arg field: :class:`.NonLinearTransform` to convert
:return: A FNIRT-compatible :class:`.DisplacementField` or
:return: A FNIRT-compatible :class:`.DeformationField` or
:class:`.CoefficientField`.
"""
from . import nonlinear
# If we have a coefficient field
# which transforms between fsl
# space, we can just create a copy.
......@@ -224,36 +369,35 @@ def toFnirt(field):
srcToRefMat=field.srcToRefMat)
# Otherwise we have a non-FSL coefficient
# field, or a displacement field.
#
# 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.
# 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.coefficientFieldToDisplacementField(field)
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.DisplacementField(
field = nonlinear.DeformationField(
field.data,
header=field.header,
src=field.src,
ref=field.ref,
xform=field.voxToWorldMat,
dispType=field.displacementType)
defType=field.deformationType)
# Otherwise we have to adjust the
# displacements so they transform
# between fsl coordinates.
field = nonlinear.convertDisplacementSpace(
field = nonlinear.convertDeformationSpace(
field, from_='fsl', to='fsl')
field.header['intent_code'] = constants.FSL_FNIRT_DISPLACEMENT_FIELD
......@@ -263,23 +407,19 @@ def toFnirt(field):
def fromFnirt(field, from_='world', to='world'):
"""Convert a FNIRT-style :class:`.NonLinearTransform` to a generic
:class:`.DisplacementField`.
:class:`.DeformationField`.
:arg field: A FNIRT-style :class:`.CoefficientField` or
:class:`.DisplacementField`
:class:`.DeformationField`
:arg from_: Desired reference image coordinate system
:arg to: Desired source image coordinate system
:return: A :class:`.DisplacementField` which contains displacements
:return: A :class:`.DeformationField` which contains displacements
from the reference image ``from_`` cordinate system to the
source image ``to`` coordinate syste.
"""
from . import nonlinear
# see comments in toFnirt
if isinstance(field, nonlinear.CoefficientField):
field = nonlinear.coefficientFieldToDisplacementField(field)
return nonlinear.convertDisplacementSpace(field, from_=from_, to=to)
field = nonlinear.coefficientFieldToDeformationField(field)
return nonlinear.convertDeformationSpace(field, from_=from_, to=to)
fsl/transform/nonlinear.py
View file @ 3b2a186e
......@@ -8,29 +8,33 @@
FNIRT-style nonlinear transformations.
The :class:`DisplacementField` and :class:`CoefficientField` can be used to
load and interact with FNIRT transformation images. The following utility
functions are also available:
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:
detectDisplacementType
convertDisplacementType
convertDisplacementSpace
coefficientFieldToDisplacementField
detectDeformationType
convertDeformationType
convertDeformationSpace
applyDeformation
coefficientFieldToDeformationField
"""
import logging
import itertools as it
import numpy as np
import logging
import itertools as it
import fsl.utils.memoize as memoize
import fsl.data.image as fslimage
import numpy as np
import scipy.ndimage as ndimage
from . import affine
import fsl.utils.memoize as memoize
import fsl.data.image as fslimage
import fsl.data.constants as fslconstants
import fsl.utils.image.resample as resample
from . import affine
log = logging.getLogger(__name__)
......@@ -38,7 +42,7 @@ log = logging.getLogger(__name__)
class NonLinearTransform(fslimage.Image):
"""Class which represents a nonlinear transformation. This is just a base
class for the :class:`DisplacementField` and :class:`CoefficientField`
class for the :class:`DeformationField` and :class:`CoefficientField`
classes.
......@@ -52,18 +56,7 @@ class NonLinearTransform(fslimage.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 source image.
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 ``NonLinearTransform`` as the
``srcToRefMat`` argument to :meth:`__init__`, and is subsequently accessed
via the :meth:`srcToRefMat` attribute.
reference image *to* the (unwarped) source image.
"""
......@@ -73,9 +66,10 @@ class NonLinearTransform(fslimage.Image):
ref,
srcSpace=None,
refSpace=None,
srcToRefMat=None,
**kwargs):
"""Create a ``NonLinearTransform``.
"""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`
......@@ -93,20 +87,17 @@ class NonLinearTransform(fslimage.Image):
``NonLinearTransform`` maps from. Defaults to
``'fsl'``.
: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).
All other arguments are passed through to :meth:`.Image.__init__`.
"""
if srcSpace is None: srcSpace = 'fsl'
if refSpace is None: refSpace = 'fsl'
if srcToRefMat is not None: srcToRefMat = np.copy(srcToRefMat)
# 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'):
......@@ -115,15 +106,10 @@ class NonLinearTransform(fslimage.Image):
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
self.__srcToRefMat = srcToRefMat
self.__refToSrcMat = None
if srcToRefMat is not None:
self.__refToSrcMat = affine.invert(srcToRefMat)
self.__src = fslimage.Nifti(src.header.copy())
self.__ref = fslimage.Nifti(ref.header.copy())
self.__srcSpace = srcSpace
self.__refSpace = refSpace
@property
......@@ -158,26 +144,13 @@ class NonLinearTransform(fslimage.Image):
return self.__refSpace
@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
def transform(self, coords, from_=None, to=None, premat=True):
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.
......@@ -186,31 +159,37 @@ class NonLinearTransform(fslimage.Image):
:arg to: Source image space to transform ``coords`` into
:arg premat: If ``True``, the inverse :meth:`srcToRefMat` is applied
to the coordinates after they have been calculated.
:returns: ``coords``, transformed into the source image space
:returns: The corresponding coordinates in the source image space.
"""
raise NotImplementedError()
class DisplacementField(NonLinearTransform):
"""Class which represents a displacement field which, at each voxel,
contains an absolute or relative displacement between a source space and a
reference space.
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 ref: Optional. If not provided, it is assumed that the
reference is defined in the same space as ``image``.
:arg dispType: Either ``'absolute'`` or ``'relative'``, indicating
the type of this displacement field. If not provided,
will be inferred via the :func:`detectDisplacementType`
function.
: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__`.
......@@ -219,46 +198,62 @@ class DisplacementField(NonLinearTransform):
if ref is None:
ref = self
dispType = kwargs.pop('dispType', None)
defType = kwargs.pop('defType', None)
if dispType not in (None, 'relative', 'absolute'):
raise ValueError('Invalid value for dispType: {}'.format(dispType))
if defType not in (None, 'relative', 'absolute'):
raise ValueError(f'Invalid value for defType: {defType}')
NonLinearTransform.__init__(self, image, src, ref, **kwargs)
if not self.sameSpace(self.ref):
raise ValueError('Invalid reference image: {}'.format(self.ref))
# Displacement fields generated by the FSL
# convertwarp tool prior to FSL 6.0.5 would
# often not have their sform/qform affines
# set correctly. Here we check the sform/
# qform and, if necessary, construct a new
# affine from the reference image.
if self.getXFormCode() == fslconstants.NIFTI_XFORM_UNKNOWN:
xform = self.ref.voxToWorldMat
self.__dispType = dispType
# It is possible for a displacement field
# to have different resolution to its
# reference space - they just have to be
# aligned in the world coordinate system.
if self.shape[:3] != self.ref.shape[:3]:
rescale = affine.rescale(self.shape[:3], self.ref.shape[:3])
xform = affine.concat(xform, rescale)
self.voxToWorldMat = xform
self.__defType = defType
@property
def displacementType(self):
"""The type of this ``DisplacementField`` - ``'absolute'`` or
def deformationType(self):
"""The type of this ``DeformationField`` - ``'absolute'`` or
``'relative'``.
"""
if self.__dispType is None:
self.__dispType = detectDisplacementType(self)
return self.__dispType
if self.__defType is None:
self.__defType = detectDeformationType(self)
return self.__defType
@property
def absolute(self):
"""``True`` if this ``DisplacementField`` contains absolute
displacements.
"""``True`` if this ``DeformationField`` contains absolute
coordinates.
"""
return self.displacementType == 'absolute'
return self.deformationType == 'absolute'
@property
def relative(self):
"""``True`` if this ``DisplacementField`` contains relative
"""``True`` if this ``DeformationField`` contains relative
displacements.
"""
return self.displacementType == 'relative'
return self.deformationType == 'relative'
def transform(self, coords, from_=None, to=None, premat=True):
def transform(self, coords, from_=None, to=None):
"""Transform the given XYZ coordinates from the reference image space
to the source image space.
......@@ -270,16 +265,15 @@ class DisplacementField(NonLinearTransform):
:arg to: Source image space to transform ``coords`` into
:arg premat: If ``True``, the inverse :meth:`srcToRefMat` is applied
to the coordinates after they have been calculated.
: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)
coords = np.asanyarray(coords)
outshape = coords.shape
coords = coords.reshape((-1, 3))
# We may need to pre-transform the
# coordinates so they are in the
......@@ -289,18 +283,24 @@ class DisplacementField(NonLinearTransform):
xform = self.ref.getAffine(from_, self.refSpace)
coords = affine.transform(coords, xform)
# We also need to make sure that the
# coordinates are in voxels, so we
# can look up the displacements
if self.refSpace != 'voxel':
xform = self.ref.getAffine(self.refSpace, 'voxel')
voxels = affine.transform(coords, xform)
else:
# 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
voxels = np.round(voxels).astype(np.int)
# that are out of bounds of
# the deformation field
voxels = np.round(voxels).astype(np.int32)
voxmask = (voxels >= [0, 0, 0]) & (voxels < self.shape[:3])
voxmask = voxmask.all(axis=1)
voxels = voxels[voxmask]
......@@ -310,32 +310,18 @@ class DisplacementField(NonLinearTransform):
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, and apply the initial
# inv(srcToRefMat) if it there
# is one
# And here the premat becomes
# a postmat; how confusing
if premat: postmat = self.refToSrcMat
else: postmat = None
# Make sure the coordinates are in
# the requested source image space
if to != self.srcSpace:
xform = self.src.getAffine(self.srcSpace, to)
if postmat is not None: postmat = affine.concat(xform, postmat)
else: postmat = xform
if postmat is not None:
disps = affine.transform(disps, postmat)
disps = affine.transform(disps, xform)
# Nans for input coordinates
# which were outside of the
# field
# Nans for input coordinates which
# were outside of the field
outcoords = np.full(coords.shape, np.nan)
outcoords[voxmask] = disps
return outcoords
return outcoords.reshape(outshape)
class CoefficientField(NonLinearTransform):
......@@ -343,6 +329,17 @@ class CoefficientField(NonLinearTransform):
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.
"""
......@@ -350,15 +347,16 @@ class CoefficientField(NonLinearTransform):
image,
src,
ref,
srcSpace,
refSpace,
fieldType,
knotSpacing,
fieldToRefMat,
srcSpace=None,
refSpace=None,
fieldType='cubic',
knotSpacing=None,
fieldToRefMat=None,
srcToRefMat=None,
**kwargs):
"""Create a ``CoefficientField``.
:arg fieldType: Must be ``'cubic'``
:arg fieldType: Must currently be ``'cubic'``
:arg knotSpacing: A tuple containing the spline knot spacings along
each axis.
......@@ -367,6 +365,14 @@ class CoefficientField(NonLinearTransform):
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.
"""
......@@ -374,6 +380,9 @@ class CoefficientField(NonLinearTransform):
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,
......@@ -384,9 +393,14 @@ class CoefficientField(NonLinearTransform):
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):
......@@ -420,11 +434,30 @@ class CoefficientField(NonLinearTransform):
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 asDisplacementField(self, dispType='relative', premat=True):
"""Convert this ``CoefficientField`` to a :class:`DisplacementField`.
def asDeformationField(self, defType='relative', premat=True):
"""Convert this ``CoefficientField`` to a :class:`DeformationField`.
This method is a wrapper around
:func:`coefficientFieldToDeformationField`
"""
return coefficientFieldToDisplacementField(self, dispType, premat)
return coefficientFieldToDeformationField(self, defType, premat)
def transform(self, coords, from_=None, to=None, premat=True):
......@@ -442,16 +475,16 @@ class CoefficientField(NonLinearTransform):
:arg premat: If ``True``, the inverse :meth:`srcToRefMat` is applied
to the coordinates after the displacements have been
addd.
added.
:returns: ``coords``, transformed into the source image space
"""
df = self.asDisplacementField(premat=premat)
df = self.asDeformationField(premat=premat)
return df.transform(coords, from_, to)
def displacements(self, coords):
"""Calculate the relative displacemenets for the given coordinates.
"""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
......@@ -483,13 +516,11 @@ class CoefficientField(NonLinearTransform):
fdata = self.data
nx, ny, nz = self.shape[:3]
ix, iy, iz = self.ref.shape[:3]
# Convert the given voxel coordinates
# into the corresponding coefficient
# field voxel coordinates
x, y, z = coords.T
i, j, k = affine.transform(coords, self.refToFieldMat).T
i, j, k = affine.transform(coords, self.refToFieldMat).T
# i, j, k: coefficient field indices
# u, v, w: position of the ref voxel
......@@ -497,9 +528,9 @@ class CoefficientField(NonLinearTransform):
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)
i = np.floor(i).astype(np.int32)
j = np.floor(j).astype(np.int32)
k = np.floor(k).astype(np.int32)
disps = np.zeros(coords.shape)
......@@ -531,23 +562,23 @@ class CoefficientField(NonLinearTransform):
return disps
def detectDisplacementType(field):
"""Attempt to automatically determine whether a displacement field is
def detectDeformationType(field):
"""Attempt to automatically determine whether a deformation field is
specified in absolute or relative coordinates.
:arg field: A :class:`DisplacementField`
:arg field: A :class:`DeformationField`
:returns: ``'absolute'`` if it looks like ``field`` contains absolute
displacements, ``'relative'`` otherwise.
coordinates, ``'relative'`` otherwise.
"""
# This test is based on the assumption
# that a displacement field containing
# that a deformation field containing
# absolute coordinates will have a
# greater standard deviation than one
# which contains relative coordinates.
absdata = field[:]
reldata = convertDisplacementType(field, 'relative')
reldata = convertDeformationType(field, 'relative')
stdabs = absdata.std(axis=(0, 1, 2)).sum()
stdrel = reldata.std(axis=(0, 1, 2)).sum()
......@@ -555,26 +586,30 @@ def detectDisplacementType(field):
else: return 'relative'
def convertDisplacementType(field, dispType=None):
"""Convert a displacement field between storing absolute and relative
displacements.
def convertDeformationType(field, defType=None):
"""Convert a deformation field between storing absolute coordinates or
relative displacements.
:arg field: A :class:`DisplacementField` instance
:arg dispType: Either ``'absolute'`` or ``'relative'``. If not provided,
the opposite type to ``field.displacementType`` is used.
:returns: A ``numpy.array`` containing the adjusted displacement
field.
: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 dispType is None:
if field.displacementType == 'absolute': dispType = 'relative'
else: dispType = 'absolute'
if defType is None:
if field.deformationType == 'absolute': defType = 'relative'
else: defType = 'absolute'
if defType not in ('absolute', 'relative'):
raise ValueError('defType must be "absolute" or "relative" '
'("{}" passed)'.format(defType))
# Regardless of the conversion direction,
# we need the coordinates of every voxel
# in the reference FSL coordinate system.
# in the reference coordinate system.
dx, dy, dz = field.shape[:3]
xform = field.getAffine('voxel', field.refSpace)
xform = affine.concat(field.ref.getAffine('world', field.refSpace),
field .getAffine('voxel', 'world'))
coords = np.meshgrid(np.arange(dx),
np.arange(dy),
......@@ -584,40 +619,36 @@ def convertDisplacementType(field, dispType=None):
coords = coords.reshape((dx, dy, dz, 3))
# If converting from relative to absolute,
# we just add the voxel coordinates to
# (what is assumed to be) the relative shift.
# Or, to convert from absolute to relative,
# 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 dispType == 'absolute': return field.data + coords
elif dispType == 'relative': return field.data - coords
if defType == 'absolute': return field.data + coords
else: return field.data - coords
def convertDisplacementSpace(field, from_, to):
"""Adjust the source and/or reference spaces of the given displacement
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:`DisplacementField` instance
:arg field: A :class:`DeformationField` instance
:arg from_: New reference image coordinate system
:arg to: New source image coordinate system
:returns: A new :class:`DisplacementField` which transforms between
:returns: A new :class:`DeformationField` which transforms between
the reference ``from_`` coordinate system and the source ``to``
coordinate system.
"""
# We can't adjust the displacements
# for fields with an initial
# source-to-ref affine
if field.srcToRefMat is not None:
raise ValueError('Cannot adjust displacements of '
'fields with an initial affine')
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 = convertDisplacementType(field)
if field.relative: srccoords = convertDeformationType(field)
else: srccoords = fieldcoords
srccoords = srccoords.reshape((-1, 3))
......@@ -631,13 +662,13 @@ def convertDisplacementSpace(field, from_, to):
srccoords = affine.transform(srccoords, srcmat)
# If we have been asked to return
# an absolute displacement, the
# absolute coordinates, the
# reference "from_" coordinate
# system is irrelevant - we're done.
if field.absolute:
fieldcoords = srccoords
# Otherwise our displacement field
# Otherwise our deformation field
# will contain relative displacements
# between the reference image "from_"
# coordinate system and the source
......@@ -653,36 +684,151 @@ def convertDisplacementSpace(field, from_, to):
refcoords = np.array(refcoords)
refcoords = refcoords.transpose((1, 2, 3, 0)).reshape((-1, 3))
if from_ != 'voxel':
refmat = field.ref.getAffine('voxel', from_)
refcoords = affine.transform(refcoords, refmat)
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 DisplacementField(
return DeformationField(
fieldcoords.reshape(field.shape),
header=field.header,
src=field.src,
ref=field.ref,
srcSpace=to,
refSpace=from_,
dispType=field.displacementType)
defType=field.deformationType)
def applyDeformation(image,
field,
ref=None,
order=1,
mode=None,
cval=None,
premat=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
def coefficientFieldToDisplacementField(field,
dispType='relative',
premat=True):
"""Convert a :class:`CoefficientField` into a :class:`DisplacementField`.
: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 field: :class:`CoefficientField` to convert
:arg mode: How to handle regions which are outside of the image FOV.
Defaults to `''nearest'``.
:arg dispType: The type of displacement field - either ``'relative'`` (the
default) or ``'absolute'``.
:arg cval: Constant value to use when ``mode='constant'``.
:arg premat: If ``True`` (the default), the :meth:`srcToRefMat` is
encoded into the displacements.
:arg premat: Optional affine transform which can be used if ``image`` is
not in the same space as ``field.src``. Assumed to transform
from ``image`` **voxel** coordinates into ``field.src``
**voxel** coordinates.
:return: :class:`DisplacementField` calculated from ``field``.
: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 (premat is not None) or (not image.sameSpace(src)):
if premat is None:
premat = affine.concat(image.getAffine('world', 'voxel'),
src .getAffine('voxel', 'world'))
else:
premat = affine.invert(premat)
shape = field.shape
field = field.reshape((-1, 3))
field = affine.transform(field, premat)
field = field.reshape(shape)
field = field.transpose((3, 0, 1, 2))
return ndimage.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
......@@ -715,22 +861,22 @@ def coefficientFieldToDisplacementField(field,
# from ref space to aligned-src
# space.
disps = field.displacements(xyz).reshape((ix, iy, iz, 3))
rdfield = DisplacementField(disps,
src=field.src,
ref=field.ref,
srcSpace=field.srcSpace,
refSpace=field.refSpace,
header=field.ref.header,
dispType='relative')
if (dispType == 'relative') and (not premat):
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
# displacements will now be
# deformations will now be
# absolute coordinates in
# aligned-src space
disps = convertDisplacementType(rdfield)
disps = convertDeformationType(rdfield)
# Apply the premat if requested -
# this will transform the coordinates
......@@ -756,23 +902,23 @@ def coefficientFieldToDisplacementField(field,
#
# disps = affine.transform(disps, refToSrc)
adfield = DisplacementField(disps,
src=field.src,
ref=field.ref,
srcSpace=field.srcSpace,
refSpace=field.refSpace,
header=field.ref.header,
dispType='absolute')
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 dispType == 'absolute':
if defType == 'absolute':
return adfield
else:
return DisplacementField(convertDisplacementType(adfield),
src=field.src,
ref=field.ref,
srcSpace=field.srcSpace,
refSpace=field.refSpace,
header=field.ref.header,
dispType='relative')
return DeformationField(convertDeformationType(adfield),
src=field.src,
ref=field.ref,
srcSpace=field.srcSpace,
refSpace=field.refSpace,
header=field.ref.header,
defType='relative')
fsl/transform/x5.py
View file @ 3b2a186e
......@@ -4,22 +4,38 @@
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module contains functions for reading/writing linear/non-linear FSL
"""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
.. warning:: This is a development release, and is subject to change.
An X5 file is a HDF5 container file which stores a linear or non-linear
transformation between a **source** NIfTI image and a **reference** image. In
an X5 file, the source image is referred to as the ``From`` space, and the
reference image the ``To`` space.
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
......@@ -31,8 +47,9 @@ HDF5 files are composed primarily of *groups*, *attributes*, and
- *Groups* are simply containers for attributes, datasets, and other groups.
- *Attributes* are strongly-typed scalar values.
- *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
......@@ -45,334 +62,264 @@ attributes/datasets that are listed here.
--------
A HDF5 group which is listed as being of type "affine" contains an 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 | ``'linear'`` |
+---------------+-----------+-------------------------------------------------+
| ``Transform`` | dataset | The affine transformation - 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. The *world*
coordinate system of an image is defined by its ``sform`` or ``qform``
(hereafter referred to as the ``sform``), which is contained in the NIfTI
header.
Groups of type "space" have the following fields:
another. Groups of type *affine* have the following fields:
+-------------+-----------+---------------------------------------------------+
| **Name** | **Type** | **Value/Description** |
+-------------+-----------+---------------------------------------------------+
| ``Type`` | attribute | ``'image'`` |
| ``Type`` | attribute | ``'affine'`` |
+-------------+-----------+---------------------------------------------------+
| ``Size`` | attribute | ``uint64`` ``(X, Y, Z)`` voxel dimensions |
| ``Matrix`` | dataset | The affine transformation matrix - a ``float64`` |
| | | array of shape ``(4, 4)`` |
+-------------+-----------+---------------------------------------------------+
| ``Scales`` | attribute | ``float64`` ``(X, Y, Z)`` voxel pixdims |
| ``Inverse`` | dataset | Optional pre-calculated inverse |
+-------------+-----------+---------------------------------------------------+
| ``Mapping`` | affine | The image 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 **source image world
coordinates** into **reference image world coordinates**.
*space*
-------
Linear X5 transformation files are assumed to adhere to the HDF5 structure
defined in the table below. All fields are required.
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.
+-----------------------------+-----------+-----------------------------------+
| **Name** | **Type** | **Value/Description** |
+-----------------------------+-----------+-----------------------------------+
| *Metadata* |
+-----------------------------+-----------+-----------------------------------+
| ``/Format`` | attribute | ``'X5'`` |
+-----------------------------+-----------+-----------------------------------+
| ``/Version`` | attribute | ``'0.0.1'`` |
+-----------------------------+-----------+-----------------------------------+
| ``/Metadata`` | attribute | JSON string containing |
| | | unstructured metadata. |
+-----------------------------+-----------+-----------------------------------+
| *Transformation* |
+-----------------------------+-----------+-----------------------------------+
| ``/`` | affine | Affine transformation from source |
| | | image world coordinates to |
| | | reference image world coordinates |
+-----------------------------+-----------+-----------------------------------+
| ``/From/`` | space | Source image definition |
+-----------------------------+-----------+-----------------------------------+
| ``/To/`` | space | Reference image definition |
+-----------------------------+-----------+-----------------------------------+
Groups of type *space* have the following fields:
Storage of FSL FLIRT matrices in linear X5 files
------------------------------------------------
+--------------+-----------+--------------------------------------------------+
| **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``) |
+--------------+-----------+--------------------------------------------------+
.. image:: images/x5_linear_transform_file.png
:width: 80%
:align: center
*deformation*
-------------
Non-linear X5 transformation files
==================================
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**).
Non-linear X5 transformation files contain a non-linear transformation between
a source image coordinate system and a reference image coordinate system. The
transformation is represented as either:
The transformation is represented as a 3D **deformation field** which, at each
voxel within the field, may contain:
- A *displacement field*, which is defined in the same space as the reference
image, and which contains displacements from reference image coordinates to
source image coordinates.
- *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**).
- A quadratic or cubic B-spline *coefficient field*, which contains
coefficients defined in a coarse grid overlaid onto the same space as the
reference image, and from which a displacement field can be calculated.
- *absolute coordinates* in space **B**.
Displacement fields
-------------------
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.
A displacement field is a ``float64`` array of shape ``(X, Y, Z, 3)``, defined
in the same space as the reference image (i.e. the reference image is assumed
to have shape ``(X, Y, Z)``. A displacement field may contain either:
Groups of type *deformation* have the following fields:
- *relative* displacements, where each voxel in the displacement 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* displacements, where each voxel in the displacement field simply
contains the source image coordinates which correspond to those reference
image coordinates.
+--------------+-----------+--------------------------------------------------+
| **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``) |
+--------------+-----------+--------------------------------------------------+
Coefficient fields
------------------
Linear X5 files
===============
A coefficient field is a ``float64`` array of shape ``(X, Y, Z, 3)`` which
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*.
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.
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.
Linear X5 transformation files are assumed to adhere to the HDF5 structure
defined in the table below. All fields are required unless otherwise noted.
The shape of this coefficient grid is not necessarily the same as the shape of
the reference image grid. For this reason, some additional parameters are
stored in coefficient field files, in a sub-group called ``/Parameters/``:
+-----------------+-----------+-----------------------------------------------+
| **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 |
+-----------------+-----------+-----------------------------------------------+
- The distance between control points, defined in terms of reference image
voxels.
- An affine transformation which can be used to transform reference image
voxel coordinates into coefficient field voxel coordinates.
Storage of FSL FLIRT matrices in linear X5 files
------------------------------------------------
Non-linear coordinate systems
-----------------------------
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** [*]_.
The coordinate systems used in a displacement field, or in a displacement
field that has been generated from a coefficient field, defines relative
displacements from a reference image coordinate system to a source image
coordinate system. The coordinate systems used are not defined in this
specification - they may be voxels, world coordinates, FSL coordinates, or
some other coordinate system.
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.
Howewer, if the transformation does not transform between source and reference
image **world** coordinates, the ``/Pre/`` and ``/Post/`` affine
transformations must be provided.
The :mod:`fsl.transform.flirt` module contains functions for converting
between FLIRT-style matrices and X5 style matrices.
The ``/Pre/`` affine transformation will be used to transform reference image
world coordinates into the reference image coordinate system required for use
with the displacement values (or required for the spline coefficient
evaluation). Similarly, the ``/Post/`` affine transformation will be used to
transform the resulting source image coordinates into the source image world
coordinate system.
.. [*] 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.
Initial affine alignment
------------------------
Non-linear X5 files
===================
Non-linear transformations are often accompanied by an initial affine
transformation, which provides a coarse global initial alignment that is used
as the starting point for the non-linear registration process. The non-linear
transformation will then typically encode displacements between the
registered, or aligned source image, and the reference image.
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:
Now we have three spaces, and three sets of coordinate systems, to consider:
- *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).
1. Source image space - the source image, before initial linear registration
to the reference image
- *absolute coordinates* in the image **B** world coordinate system.
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 encodes a transformation between spaces 1 and
2, and the non-linear registration encodes a transformation between spaces 2
and 3. Note that the fields-of-view for spaces 2 and 3 are typically
equivalent.
File format specification
-------------------------
The initial affine transformation may be included in an X5 file, in the
``/InitialAlignment/`` group. If provided, this initial transformation is
assumed to provide a transformation:
Non-linear X5 transformation files are assumed to adhere to the following
HDF5 structure. All fields are required unless otherwise noted.
- *From* the **source image** world coordinate system (or the coordinate
system used as input to the ``/Post/`` affine, if provided).
- *To* the **aligned source** image coordinate system used within the
non-linear transformation.
+---------------+-----------+-------------------------------------------------+
| **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 |
+-----------------+-------------+---------------------------------------------+
File format specification
-------------------------
Storage of FSL FNIRT warp fields in non-linear X5 files
-------------------------------------------------------
Non-linear X5 transformation files are assumed to adhere to the following
HDF5 structure. All fields are required unless otherwise noted:
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
+---------------------+-----------+-------------------------------------------+
| **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'`` |
+------------------------+-----------+----------------------------------------+
| ``/SubType`` | attribute | ``'displacement'`` or |
| | | ``'coefficient'`` |
+------------------------+-----------+----------------------------------------+
| ``/Representation`` | attribute | If ``/SubType`` is ``'displacement'``, |
| | | ``/Representation`` may be either |
| | | ``'absolute'`` or ``'relative'``. |
| | | If ``/SubType`` is ``'coefficient'``, |
| | | ``/Representation`` may be either |
| | | ``'quadratic bspline'`` or |
| | | ``'cubic bspline'``. |
+------------------------+-----------+----------------------------------------+
| ``/Transform`` | dataset | The displacement/coefficient field - |
| | | see above description. |
+------------------------+-----------+----------------------------------------+
| ``/Inverse`` | dataset | Optional pre-calculated inverse |
+------------------------+-----------+----------------------------------------+
| ``/From/`` | space | Source image definition |
+------------------------+-----------+----------------------------------------+
| ``/To/`` | space | Reference image definition |
+------------------------+-----------+----------------------------------------+
| ``/Pre/`` | affine | Optional affine transformation from |
| | | reference image world coordinates to |
| | | the reference image coordinate system |
| | | required by the displacement/ |
| | | coefficient field. |
+------------------------+-----------+----------------------------------------+
| ``/Post/`` | affine | Optional affine transformation from |
| | | the source image coordinate system |
| | | produced by the displacement/ |
| | | coefficient field to the source image |
| | | world coordinate system. |
+------------------------+-----------+----------------------------------------+
| ``/InitialAlignment/`` | affine | Optional initial affine registration |
| | | from the source image to the reference |
| | | image. |
+------------------------+-----------+----------------------------------------+
| *Coefficient field parameters* (required for ``'coefficient'`` files) |
+-----------------------------------+-----------+-----------------------------+
| **Name** | **Type** | **Value/Description** |
+-----------------------------------+-----------+-----------------------------+
| ``/Parameters/Spacing`` | attribute | ``uint64`` ``(X, Y, Z)`` |
| | | knot spacing (defined in |
| | | reference image voxels) |
+-----------------------------------+-----------+-----------------------------+
| ``/Parameters/ReferenceToField/`` | affine | Reference image voxel to |
| | | coefficient field voxel |
| | | affine transformation. |
+-----------------------------------+-----------+-----------------------------+
Storage of FSL FNIRT files in nonlinear X5 files
------------------------------------------------
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**.
.. image:: images/nonlinear_registration_process.png
:width: 80%
:align: center
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.
Displacement fields
^^^^^^^^^^^^^^^^^^^
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.
.. image:: images/fnirt_displacement_field.png
:width: 80%
:align: center
Conversion of FNIRT coefficient fields
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: images/x5_nonlinear_displacement_field_file.png
:width: 95%
:align: center
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.
Coefficient fields
^^^^^^^^^^^^^^^^^^
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.
.. image:: images/fnirt_coefficient_field.png
:width: 80%
:align: center
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.
.. image:: images/x5_nonlinear_coefficient_field_file.png
:width: 95%
:align: center
"""
......@@ -382,14 +329,14 @@ import numpy as np
import nibabel as nib
import h5py
import fsl.version as version
import fsl.version as fslversion
import fsl.data.image as fslimage
from . import affine
from . import nonlinear
X5_FORMAT = 'X5'
X5_VERSION = '0.0.1'
X5_VERSION = '0.1.0'
class X5Error(Exception):
......@@ -397,6 +344,23 @@ class X5Error(Exception):
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``.
......@@ -408,17 +372,22 @@ def readLinearX5(fname):
- A :class:`.Nifti` instance representing the source space
- A :class:`.Nifti` instance representing the reference space
"""
with h5py.File(fname, 'r') as f:
_validateMetadata( f['/'])
xform = _readAffine( f['/'])
src = _readSpace( f['/From'])
ref = _readSpace( f['/To'])
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 the specified file.
"""Write a linear transformation to ``fname``.
:arg fname: File name to write to
:arg xform: ``(4, 4)`` ``numpy`` array containing the affine transformation
......@@ -427,98 +396,111 @@ def writeLinearX5(fname, xform, src, ref):
"""
with h5py.File(fname, 'w') as f:
_writeMetadata(f)
_writeAffine(f, xform)
from_ = f.create_group('/From')
to = f.create_group('/To')
f.attrs['Type'] = 'linear'
_writeSpace(from_, src)
_writeSpace(to, ref)
_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:`.DisplacementField` or :class:`.CoefficientField`
:returns: A :class:`.DeformationField`
"""
with h5py.File(fname, 'r') as f:
root = f['/']
_validateNonLinearTransform(root)
if f.attrs.get('Type') != 'nonlinear':
raise X5Error('Not a nonlinear transform')
subtype = root.attrs['SubType']
_readMetadata(f)
if subtype == 'displacement': field = _readDisplacementField(root)
elif subtype == 'coefficient': field = _readCoefficientField(root)
ref = _readSpace( f['/A'])
src = _readSpace( f['/B'])
field, xform, defType = _readDeformation(f['/Transform'])
return field
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 file to ``fname``.
"""Write a nonlinear X5 transformation to ``fname``.
:arg fname: File name to write to
:arg field: A :class:`.DisplacementField` or :class:`.CoefficientField`
:arg field: A :class:`.DeformationField`
"""
with h5py.File(fname, 'w') as f:
_writeMetadata(f)
f.attrs['Type'] = 'nonlinear'
if isinstance(field, nonlinear.DisplacementField):
_writeDisplacementField(f, field)
elif isinstance(field, nonlinear.CoefficientField):
_writeCoefficientField(f, field)
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' : version.__version__})
_writeMetadata(f)
_writeSpace( f.create_group('/A'), field.ref)
_writeSpace( f.create_group('/B'), field.src)
_writeDeformation(f.create_group('/Transform'), field)
def _validateMetadata(group):
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.
"""
try:
format = group.attrs['Format']
version = group.attrs['Version']
format = group.attrs.get('Format')
version = group.attrs.get('Version')
meta = group.attrs.get('Metadata')
except Exception:
raise X5Error('File does not appear to be an X5 file')
parserver = fslversion.parseVersionString(X5_VERSION)
filever = fslversion.parseVersionString(version)
if (format != X5_FORMAT) or (version != X5_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.
"""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['Type'] != 'linear':
raise X5Error('Not a linear transform')
if group.attrs.get('Type') != 'affine':
raise X5Error('Not an affine')
xform = group['Transform']
xform = group['Matrix']
if xform.shape != (4, 4):
raise X5Error('Not a linear transform')
raise X5Error('Not an affine')
return np.array(xform)
......@@ -534,21 +516,19 @@ def _writeAffine(group, xform):
xform = np.asarray(xform, dtype=np.float64)
inv = np.asarray(affine.invert(xform), dtype=np.float64)
group.attrs['Type'] = 'linear'
group.create_dataset('Transform', data=xform)
group.create_dataset('Inverse', data=inv)
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.
"""Reads a *space* group, defining a source or reference space.
:arg group: A ``h5py.Group`` object.
:returns: :class:`.Nifti` object. defining the mapping
"""
import fsl.data.image as fslimage
if group.attrs['Type'] != 'image':
if group.attrs.get('Type') != 'image':
raise X5Error('Not an image space')
shape = group.attrs['Size']
......@@ -576,198 +556,55 @@ def _writeSpace(group, img):
_writeAffine(mapping, img.getAffine('voxel', 'world'))
def _validateNonLinearTransform(group):
"""Checks that the attributes of the given group, assumed to contain a
nonlinear transform, are valid. Raises a :exc:`X5Error` if not.
:arg group: A ``h5py.Group`` object.
"""
def _readDeformation(group):
"""Reads a *deformation* from the given group.
type = group.attrs['Type']
subtype = group.attrs['SubType']
repr = group.attrs['Representation']
:arg group: A ``h5py.Group`` object
:returns: A tuple containing
if type != 'nonlinear':
raise X5Error('Not a nonlinear transform')
- A ``numpy.array`` containing the deformation field
if subtype not in ('displacement', 'coefficient'):
raise X5Error('Unrecognised nonlinear subtype: {}'.format(subtype))
- A ``numpy.array`` of shape ``(4, 4)`` containing the
voxel to world affine for the deformation field
if (subtype == 'displacement') and \
(repr not in ('absolute', 'relative')):
raise X5Error('Unrecognised nonlinear representation: '
'{}'.format(repr))
if (subtype == 'coefficient') and \
(repr not in ('quadratic bspline', 'cubic bspline')):
raise X5Error('Unrecognised nonlinear representation: '
'{}'.format(repr))
def _readNonLinearCommon(group):
"""Reads the spaces and affines from the given group, assumed to contain a
nonlinear transform.
:arg group: A ``h5py.Group`` object.
:returns: A tuple containing:
- A :class:`.Nifti` representing the source
- A :class:`.Nifti` representing the reference
- A ``(4, 4)`` ``numpy`` array containing the pre affine, or
``None`` if there is not one.
- A ``(4, 4)`` ``numpy`` array containing the post affine, or
``None`` if there is not one.
- A ``(4, 4)`` ``numpy`` array containing the initial
alignment affine, or ``None`` if there is not one.
- A string describing the source space - see
:meth:`.Nifti.getAffine`
- A string describing the reference space - see
:meth:`.Nifti.getAffine`
- The deformation type - either ``'absolute'`` or
``'relative'``
"""
src = _readSpace(group['From'])
ref = _readSpace(group['To'])
pre = group.get('Pre')
post = group.get('Post')
init = group.get('InitialAlignment')
type = group.attrs.get('Type')
subtype = group.attrs['SubType']
if pre is not None: pre = _readAffine(pre)
if post is not None: post = _readAffine(post)
if init is not None: init = _readAffine(init)
if type != 'deformation':
raise X5Error('Not a deformation')
refSpace = 'world'
srcSpace = 'world'
if subtype not in ('absolute', 'relative'):
raise X5Error('Unknown deformation type: {}'.format(subtype))
try:
if pre is not None:
refSpace = fslimage.Nifti.identifyAffine(
ref, pre, from_='world')[1]
if post is not None:
srcSpace = fslimage.Nifti.identifyAffine(
src, post, to='world')[ 0]
mapping = _readAffine(group['Mapping'])
field = group['Matrix']
except ValueError:
raise X5Error('Invalid pre/post affine')
if len(field.shape) != 4 or field.shape[3] != 3:
raise X5Error('Invalid shape for deformation field')
return src, ref, pre, post, init, srcSpace, refSpace
return np.array(field), mapping, subtype
def _writeNonLinearCommon(group, field):
"""Writes the spaces and affines for the given nonlinear transform to the
given group.
def _writeDeformation(group, field):
"""Write a deformation field to the given group.
:arg group: A ``h5py.Group`` object.
:arg field: A :class:`.NonLinearTransform` object
:arg group: A ``h5py.Group`` object
:arg field: A :class:`.DeformationField` object
"""
_writeSpace(group.create_group('From'), field.src)
_writeSpace(group.create_group('To'), field.ref)
pre = field.ref.getAffine('world', field.refSpace)
post = field.src.getAffine(field.srcSpace, 'world')
if field.srcSpace != 'world' or \
field.refSpace != 'world':
raise X5Error('Deformation field must encode a '
'world<->world transformation')
_writeAffine(group.create_group('Pre'), pre)
_writeAffine(group.create_group('Post'), post)
group.attrs['Type'] = 'deformation'
group.attrs['SubType'] = field.deformationType
if field.srcToRefMat is not None:
_writeAffine(group.create_group('InitialAlignment'), field.srcToRefMat)
def _readDisplacementField(group):
"""Reads a nonlinear displacement field from the given group.
:arg group: A ``h5py.Group`` object.
:returns: A :class:`.DisplacementField` object
"""
src, ref, pre, post, init, srcSpace, refSpace = _readNonLinearCommon(group)
field = np.array(group['Transform'])
dtype = group.attrs['Representation']
field = nonlinear.DisplacementField(field,
src=src,
ref=ref,
srcSpace=srcSpace,
refSpace=refSpace,
dispType=dtype,
srcToRefMat=init,
xform=ref.voxToWorldMat)
return field
def _writeDisplacementField(group, field):
"""Writes a nonlinear displacement field to the given group.
:arg group: A ``h5py.Group`` object.
:arg field: A :class:`.DisplacementField` object
"""
_writeNonLinearCommon(group, field)
group.attrs['SubType'] = 'displacement'
group.attrs['Representation'] = field.displacementType
xform = field.data.astype(np.float64)
group.create_dataset('Transform', data=xform)
def _readCoefficientField(group):
"""Reads a nonlinear coefficient field from the given group.
:arg group: A ``h5py.Group`` object.
:returns: A :class:`.CoefficientField` object
"""
src, ref, pre, post, init, srcSpace, refSpace = _readNonLinearCommon(group)
field = np.array(group['Transform'])
ftype = group.attrs['Representation']
spacing = group['Parameters'].attrs['Spacing']
refToField = _readAffine(group['Parameters/ReferenceToField'])
fieldToRef = affine.invert(refToField)
if ftype == 'quadratic bspline': ftype = 'quadratic'
elif ftype == 'cubic bspline': ftype = 'cubic'
if spacing.shape != (3,):
raise X5Error('Invalid spacing: {}'.format(spacing))
field = nonlinear.CoefficientField(field,
src=src,
ref=ref,
srcSpace=srcSpace,
refSpace=refSpace,
fieldType=ftype,
knotSpacing=spacing,
fieldToRefMat=fieldToRef,
srcToRefMat=init)
return field
def _writeCoefficientField(group, field):
"""Writes a nonlinear coefficient field to the given group.
:arg group: A ``h5py.Group`` object.
:arg field: A :class:`.CoefficientField` object
"""
_writeNonLinearCommon(group, field)
group.attrs['SubType'] = 'coefficient'
if field.fieldType == 'cubic':
group.attrs['Representation'] = 'cubic bspline'
elif field.fieldType == 'quadratic':
group.attrs['Representation'] = 'quadratic bspline'
xform = field.data.astype(np.float64)
group.create_dataset('Transform', data=xform)
params = group.create_group('Parameters')
refToField = params.create_group('ReferenceToField')
mapping = group.create_group('Mapping')
params.attrs['Spacing'] = field.knotSpacing
_writeAffine(refToField, field.refToFieldMat)
group.create_dataset('Matrix', data=field.data)
_writeAffine(mapping, field.getAffine('voxel', 'world'))
fsl/utils/assertions.py
View file @ 3b2a186e
......@@ -33,25 +33,32 @@ import fsl.utils.ensure as ensure
import fsl.data.melodicanalysis as fslma
_DISABLE_ASSERTIONS = False
"""
"""
_DISABLE_ASSERTIONS = 0
"""Semaphore used by the :func:`disabled` context manager. """
@contextlib.contextmanager
def disabled():
def disabled(disable=True):
"""Context manager which allows assertion checks to be temporarily
disabled.
If calls to this function are nested, only one of the calls need to be made
with ``disable=True`` for assertions to be disabled; any other calls which
are part of the call stack which set ``disable=False`` will have no effect.
:arg disable: Set to ``True`` (the default) to disable assertions,
or ``False`` to enable them.
"""
global _DISABLE_ASSERTIONS
oldval = _DISABLE_ASSERTIONS
_DISABLE_ASSERTIONS = True
if disable:
_DISABLE_ASSERTIONS += 1
try:
yield
finally:
_DISABLE_ASSERTIONS = oldval
if disable:
_DISABLE_ASSERTIONS -= 1
def _canDisable(func):
......@@ -59,7 +66,7 @@ def _canDisable(func):
via the :func:`disabled` context manager.
"""
def wrapper(*args, **kwargs):
if not _DISABLE_ASSERTIONS:
if _DISABLE_ASSERTIONS == 0:
return func(*args, **kwargs)
return wrapper
......
fsl/utils/bids.py 0 → 100644
View file @ 3b2a186e
#!/usr/bin/env python
#
# bids.py - Simple BIDS metadata reader.
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module provides a few functions for working with BIDS data sets.
.. autosummary::
:nosignatures:
BIDSFile
isBIDSDir
inBIDSDir
isBIDSFile
loadMetadata
All of the other functions in this module should not be considered part of the
public API.
See https://bids-specification.readthedocs.io/en/stable/ for more information
about BIDS.
.. note:: The `pybids <https://bids-standard.github.io/pybids/>`_ library is
a more suitable choice if you are after a more robust and featured
interface for working with BIDS datasets.
"""
import os.path as op
import itertools as it
import re
import glob
import json
import fsl.utils.memoize as memoize
import fsl.utils.path as fslpath
class BIDSFile(object):
"""The ``BIDSFile`` class parses and stores the entities and suffix
contained in a BIDS file. See the :func:`parseFilename` function.
The :meth:`match` method can be used to compare two ``BIDSFile`` instances.
The following attributes are available on a ``BIDSFile`` instance:
- ``filename``: Absolute path to the file
- ``entities``: Dict of ``key : value`` pairs, the entities that are
present in the file name (e.g. ``{'sub' : '01}``)
- ``suffix``: File suffix (e.g. ``T1w``, ``bold``, etc.)
"""
def __init__(self, filename):
"""Create a ``BIDSFile``. """
entities, suffix = parseFilename(filename)
self.filename = op.abspath(filename)
self.entities = entities
self.suffix = suffix
def __str__(self):
"""Return a strimg representation of this ``BIDSFile``. """
return 'BIDSFile({})'.format(self.filename)
def __repr__(self):
"""Return a strimg representation of this ``BIDSFile``. """
return str(self)
def match(self, other, suffix=True):
"""Compare this ``BIDSFile`` to ``other``.
:arg other: ``BIDSFile`` to compare
:arg suffix: Defaults to ``True``. If ``False``, the comparison
is made solely on the entity values.
:returns: ``True`` if ``self.suffix == other.suffix`` (unless
``suffix`` is ``False``) and if all of the entities in
``other`` are present in ``self``, ``False`` otherwise.
"""
suffix = (not suffix) or (self.suffix == other.suffix)
entities = True
for key, value in other.entities.items():
entities = entities and (self.entities.get(key, None) == value)
return suffix and entities
def parseFilename(filename):
"""Parses a BIDS-like file name, returning the entities and suffix encoded
in the name. See the :func:`isBIDSFile` function for an explanation of
what is considered to be a valid BIDS file name.
.. note:: This function assumes that no period (``.``) characters occur in
the body of a BIDS filename.
:returns: A tuple containing:
- A dict containing the entities
- The suffix
"""
if not isBIDSFile(filename, strict=False):
raise ValueError('Does not look like a BIDS '
'file: {}'.format(filename))
suffix = None
entities = []
filename = op.basename(filename)
filename = fslpath.removeExt(filename, firstDot=True)
parts = filename.split('_')
for part in parts[:-1]:
entities.append(part.split('-'))
suffix = parts[-1]
entities = dict(entities)
return entities, suffix
def isBIDSDir(dirname):
"""Returns ``True`` if ``dirname`` is the root directory of a BIDS dataset.
"""
return op.exists(op.join(dirname, 'dataset_description.json'))
def inBIDSDir(filename):
"""Returns ``True`` if ``filename`` looks like it is within a BIDS dataset
directory, ``False`` otherwise.
"""
dirname = op.abspath(op.dirname(filename))
inBIDS = False
while True:
if isBIDSDir(dirname):
inBIDS = True
break
prevdir = dirname
dirname = op.dirname(dirname)
# at filesystem root
if prevdir == dirname:
break
return inBIDS
def isBIDSFile(filename, strict=True):
"""Returns ``True`` if ``filename`` looks like a BIDS image or JSON file.
A BIDS file name must consist of zero or more "entities" (alpha-numeric
``name-value`` pairs), a "suffix", all separated by underscores, and a
regular file extension. For example, the following file::
sub-01_ses-01_task-stim_bold.nii.gz
has suffix ``bold``, entities ``sub=01``, ``ses=01`` and ``task=stim``, and
extension ``.nii.gz``.
:arg filename: Name of file to check
:arg strict: If ``True`` (the default), the file must be within a BIDS
dataset directory, as defined by :func:`inBIDSDir`.
"""
# Zero or more entities because sidecar files
# do not necessarily need to contain any
# entities (e.g. ``T1w.json`` is valid)
name = op.basename(filename)
pattern = r'([a-z0-9]+-[a-z0-9]+_)*([a-z0-9])+\.(.+)'
flags = re.ASCII | re.IGNORECASE
match = re.fullmatch(pattern, name, flags) is not None
return ((not strict) or inBIDSDir(filename)) and match
@memoize.memoize
def loadMetadataFile(filename):
"""Load ``filename`` (assumed to be JSON), returning its contents. """
with open(filename, 'rt') as f:
return json.load(f, strict=False)
def loadMetadata(filename):
"""Load all of the metadata associated with ``filename``.
:arg filename: Path to a data file in a BIDS dataset.
:returns: A dict containing all of the metadata associated with
``filename``
"""
filename = op.abspath(filename)
bfile = BIDSFile(filename)
dirname = op.dirname(filename)
prevdir = filename
metafiles = []
metadata = {}
# Walk up the directory tree until
# we hit the BIDS dataset root, or
# the filesystem root
while True:
# Gather all json files in this
# directory with matching entities
# and suffix, sorted alphabetically
# and reversed, so that earlier
# ones take precedence
files = reversed(sorted(glob.glob(op.join(dirname, '*.json'))))
files = [BIDSFile(f) for f in files if isBIDSFile(f)]
files = [f.filename for f in files if bfile.match(f)]
# build a list of all files
metafiles.append(files)
# move to the next dir up
prevdir = dirname
dirname = op.dirname(dirname)
# stop when we hit the dataset or filesystem root
if isBIDSDir(prevdir) or dirname == prevdir:
break
# Load in each json file, from
# shallowest to deepest, so entries
# in deeper files take precedence
# over shallower ones.
for f in it.chain(*reversed(metafiles)):
# assuming here that every file contains a dict
metadata.update(loadMetadataFile(f))
return metadata
fsl/utils/cache.py
View file @ 3b2a186e
......@@ -16,10 +16,9 @@ class Expired(Exception):
"""``Exception`` raised by the :meth:`Cache.get` metho when an attempt is
made to access a cache item that has expired.
"""
pass
class CacheItem(object):
class CacheItem:
"""Internal container class used to store :class:`Cache` items. """
def __init__(self, key, value, expiry=0):
......@@ -29,7 +28,7 @@ class CacheItem(object):
self.storetime = time.time()
class Cache(object):
class Cache:
"""The ``Cache`` is a simple in-memory cache built on a
``collections.OrderedDict``. The ``Cache`` class has the following
features:
......@@ -42,14 +41,20 @@ class Cache(object):
raised.
"""
def __init__(self, maxsize=100):
def __init__(self, maxsize=100, lru=False):
"""Create a ``Cache``.
:arg maxsize: Maximum number of items allowed in the ``Cache`` before
it starts dropping old items
:arg lru: (least recently used) If ``False`` (the default), items
are dropped according to their insertion time. Otherwise,
items are dropped according to their most recent access
time.
"""
self.__cache = collections.OrderedDict()
self.__maxsize = maxsize
self.__lru = lru
def put(self, key, value, expiry=0):
......@@ -94,14 +99,26 @@ class Cache(object):
else:
entry = self.__cache[key]
# Check to see if the entry
# has expired
now = time.time()
if entry.expiry > 0:
if time.time() - entry.storetime > entry.expiry:
if now - entry.storetime > entry.expiry:
self.__cache.pop(key)
if defaultSpecified: return default
else: raise Expired(key)
# If we are an lru cache, update
# this entry's expiry, and update
# its order in the cache dict
if self.__lru:
entry.storetime = now
self.__cache.pop(key)
self.__cache[key] = entry
return entry.value
......@@ -125,6 +142,30 @@ class Cache(object):
return self.put(key, value)
def __contains__(self, key):
"""Check whether an item is in the cache. Note that the item may
be in the cache, but it may be expired.
"""
return key in self.__cache
def keys(self):
"""Return all keys in the cache. """
return self.__cache.keys()
def values(self):
"""Return all values in the cache. """
for item in self.__cache.values():
yield item.value
def items(self):
"""Return all (key, value) pairs in the cache. """
for key, item in self.__cache.items():
yield key, item.value
def __parseDefault(self, *args, **kwargs):
"""Used by the :meth:`get` method. Parses the ``default`` argument,
which may be specified as either a positional or keyword argumnet.
......@@ -134,7 +175,7 @@ class Cache(object):
- ``True`` if a default argument was specified, ``False``
otherwise.
- The specifeid default value, or ``None`` if it wasn't
- The specified default value, or ``None`` if it wasn't
specified.
"""
......
fsl/utils/ensure.py
View file @ 3b2a186e
......@@ -13,8 +13,7 @@ that some condition is met.
ensureIsImage
"""
import six
import pathlib
import nibabel as nib
......@@ -24,7 +23,7 @@ import fsl.data.image as fslimage
def ensureIsImage(img):
"""Ensures that the given ``img`` is an in-memory ``nibabel`` object.
"""
if isinstance(img, six.string_types):
if isinstance(img, (str, pathlib.Path)):
img = fslimage.addExt(img)
img = nib.load(img)
return img
fsl/utils/filetree/__init__.py
View file @ 3b2a186e
"""
Easy format to define input/output files in a python pipeline.
"""Easy format to define input/output files in a python pipeline.
.. warning::
File-tree is now an independent Python library, and will eventually be
removed from ``fslpy`` - visit the `file-tree` :ref:`API documentation
<https://open.win.ox.ac.uk/pages/fsl/file-tree/>`_ for more details.
The goal is to separate the definition of the input/output filenames from the actual code
by defining a directory tree (i.e., FileTree) in a separate file from the code.
......@@ -177,7 +181,7 @@ which amongst others refers to
Example pipeline
----------------
A very simple pipeline to run BET on every subject can start with a simply FileTree like
A very simple pipeline to run BET on every subject can start with a FileTree like
::
{subject}
......@@ -200,6 +204,12 @@ Assuming that the input T1w's already exist, we can then simply run BET for ever
# make_dir=True ensures that the output directory containing the "bet_output" actually exists
bet(input=T1w_tree.get('T1w'), output=T1w_tree.get('bet_output', make_dir=True), mask=True)
Useful tips
-----------
Changing directory structure
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
If later on in our input files change, because for some subjects we added a second session, we could keep our script
and simply update the FileTree:
::
......@@ -207,8 +217,8 @@ and simply update the FileTree:
{subject}
[ses-{session}]
T1w.nii.gz
T1w_brain.nii.gz (bed_output)
T1w_brain_mask.nii.gz (bed_mask)
T1w_brain.nii.gz (bet_output)
T1w_brain_mask.nii.gz (bet_mask)
Note the square brackets around the session sub-directory. This indicates that this sub-directory is optional and
will only be present if the "session" variable is defined (see `Optional variables`_).
......@@ -230,17 +240,24 @@ altering this behaviour is again as simple as altering the FileTree to something
::
raw_data
{subject}
[ses-{session}]
{subject} (input_subject_dir)
[ses-{session}] (input_session_dir)
T1w.nii.gz
processed_data
{subject}
[ses-{session}]
{subject} (output_subject_dir)
[ses-{session}] (output_session_dir)
bet
{subject}[_{session}]_T1w_brain.nii.gz (bet_output)
{subject}[_{session}]_T1w_brain_mask.nii.gz (bet_mask)
Note that we also encoded the subject and session ID in the output filename.
We also have to explicitly assign short names to the subject and session directories,
even though we don't explicitly reference these in the script.
The reason for this is that each directory and filename template must have a unique short name and
in this case the default short names (respectively, "{subject}" and "[ses-{session}]") would not have been unique.
Output "basenames"
^^^^^^^^^^^^^^^^^^
Some tools like FSL's FAST produce many output files. Rather than entering all
of these files in our FileTree by hand you can include them all at once by including `Sub-trees`_:
......@@ -248,12 +265,12 @@ of these files in our FileTree by hand you can include them all at once by inclu
::
raw_data
{subject}
[ses-{session}]
{subject} (input_subject_dir)
[ses-{session}] (input_session_dir)
T1w.nii.gz
processed_data
{subject}
[ses-{session}]
{subject} (output_subject_dir)
[ses-{session}] (output_session_dir)
bet
{subject}[_{session}]_T1w_brain.nii.gz (bet_output)
{subject}[_{session}]_T1w_brain_mask.nii.gz (bet_mask)
......@@ -272,6 +289,35 @@ Within the script we can generate the fast output by running
The output files will be available as `T1w_tree.get('segment/<variable name>')`, where `<variable name>` is one
of the short variable names defined in the
`FAST FileTree <https://git.fmrib.ox.ac.uk/fsl/fslpy/blob/master/fsl/utils/filetree/trees/fast.tree>`_.
Running a pipeline on a subset of participants/sessions/runs
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Suppose you want to run your pipeline on a subset of your data while testing.
You may want to do this if your data has a a hierarchy of variables (e.g. participant, session, run) as in the example below.
::
sub-001
ses-01
sub-001_ses-01_run-1.feat
sub-001_ses-01_run-2.feat
ses-02
sub-{participant}_ses-{session}_run-{run}.feat (feat_dir)
...
sub-002
sub-003
...
You can update the FileTree with one or more variables before calling `get_all_trees` as follows:
.. code-block:: python
for participant in ("001", "002"):
for t in tree.update(participant=participant, run="1").get_all_trees("feat_dir", glob_vars="all"):
my_pipeline(t)
This code will iterate over all sessions that have a run="1" for participants "001" and "002".
"""
__author__ = 'Michiel Cottaar <Michiel.Cottaar@ndcn.ox.ac.uk>'
......@@ -279,3 +325,13 @@ __author__ = 'Michiel Cottaar <Michiel.Cottaar@ndcn.ox.ac.uk>'
from .filetree import FileTree, register_tree, MissingVariable
from .parse import tree_directories, list_all_trees
from .query import FileTreeQuery
import fsl.utils.deprecated as deprecated
deprecated.warn('fsl.utils.filetree',
stacklevel=2,
vin='3.6.0',
rin='4.0.0',
msg='The filetree package is now released as a separate '
'Python library ("file-tree" on PyPi), and will be '
'removed in a future version of fslpy.')
fsl/utils/filetree/filetree.py
View file @ 3b2a186e
from pathlib import Path, PurePath
from typing import Tuple, Optional, Dict, Any, Set
from copy import deepcopy
from . import parse
import pickle
import json
import os.path as op
from . import utils
from fsl.utils.deprecated import deprecated
class MissingVariable(KeyError):
......@@ -25,12 +24,14 @@ class FileTree(object):
- ``variables``: dictionary mapping variables in the templates to specific values (variables set to None are explicitly unset)
- ``sub_trees``: filename trees describing specific sub-directories
- ``parent``: parent FileTree, of which this sub-tree is a sub-directory
- ``name``: descriptive name of the tree
"""
def __init__(self,
templates: Dict[str, str],
variables: Dict[str, Any],
sub_trees: Dict[str, "FileTree"]=None,
parent: Optional["FileTree"]=None):
templates: Dict[str, str],
variables: Dict[str, Any],
sub_trees: Dict[str, "FileTree"] = None,
parent: Optional["FileTree"] = None,
name: str = None):
"""
Creates a new filename tree.
"""
......@@ -40,6 +41,7 @@ class FileTree(object):
sub_trees = {}
self.sub_trees = sub_trees
self._parent = parent
self._name = name
@property
def parent(self, ):
......@@ -48,6 +50,13 @@ class FileTree(object):
"""
return self._parent
@property
def name(self, ):
"""
Name of this ``FileTree``, or ``None`` if it has no name.
"""
return self._name
@property
def all_variables(self, ):
"""
......@@ -169,10 +178,10 @@ class FileTree(object):
:param glob_vars: sequence of undefined variables that can take any possible values when looking for matches on the disk.
Any defined variables in `glob_vars` will be ignored.
If glob_vars is set to 'all', all undefined variables will be used to look up matches.
:return: sorted sequence of paths
:return: sequence of paths
"""
text, variables = self.get_template(short_name)
return tuple(str(Path(fn)) for fn in utils.get_all(text, variables, glob_vars=glob_vars))
return tuple([self.update(**vars).get(short_name)
for vars in self.get_all_vars(short_name, glob_vars=glob_vars)])
def get_all_vars(self, short_name: str, glob_vars=()) -> Tuple[Dict[str, str]]:
"""
......@@ -184,7 +193,8 @@ class FileTree(object):
If glob_vars is set to 'all', all undefined variables will be used to look up matches.
:return: sequence of dictionaries with the variables settings used to generate each filename
"""
return tuple(self.extract_variables(short_name, fn) for fn in self.get_all(short_name, glob_vars=glob_vars))
text, variables = self.get_template(short_name)
return utils.get_all(text, variables, glob_vars=glob_vars)
def get_all_trees(self, short_name: str, glob_vars=(), set_parent=True) -> Tuple["FileTree"]:
"""
......@@ -213,7 +223,7 @@ class FileTree(object):
Setting a variable to None will cause the variable to be unset
:return: New FileTree with same templates for directory names and filenames, but updated variables
"""
new_tree = deepcopy(self)
new_tree = self.copy()
set_tree = new_tree
while set_parent and set_tree.parent is not None:
set_tree = set_tree.parent
......@@ -244,12 +254,28 @@ class FileTree(object):
with open(filename, 'wb') as f:
pickle.dump(self, f)
def save_json(self, filename):
"""
Saves the Filetree to a JSON file
:param filename: filename to store the file tree in
"""
def default(obj):
if isinstance(obj, FileTree):
res = dict(obj.__dict__)
del res['_parent']
return res
return obj
with open(filename, 'w') as f:
json.dump(self, f, default=default, indent=2)
@classmethod
def load_pickle(cls, filename):
"""
Loads the Filetree from a pickle file
:param filename: filename produced from Filetree.save
:param filename: filename produced from Filetree.save_pickle
:return: stored Filetree
"""
with open(filename, 'rb') as f:
......@@ -258,6 +284,29 @@ class FileTree(object):
raise IOError("Pickle file did not contain %s object" % cls)
return res
@classmethod
def load_json(cls, filename):
"""
Loads the FileTree from a JSON file
:param filename: filename produced by FileTree.save_json
:return: stored FileTree
"""
def from_dict(input_dict):
res_tree = FileTree(
templates=input_dict['templates'],
variables=input_dict['variables'],
sub_trees={name: from_dict(value) for name, value in input_dict['sub_trees'].items()},
name=input_dict['_name'],
)
for sub_tree in res_tree.sub_trees.values():
sub_tree._parent = res_tree
return res_tree
with open(filename, 'r') as f:
as_dict = json.load(f)
return from_dict(as_dict)
def defines(self, short_names, error=False):
"""
Checks whether templates are defined for all the `short_names`
......@@ -315,19 +364,68 @@ class FileTree(object):
return False
return True
@deprecated(rin='2.4', msg='Use FileTree.defines or FileTree.on_disk instead')
def exists(self, short_names, on_disk=False, error=False, glob_vars=()):
def partial_fill(self, ) -> "FileTree":
"""
Deprecated in favor of :meth:`on_disk` and :meth:`defines`.
Fills in known variables into the templates
:return: The resulting tree will have empty `variables` dictionaries and updated templates
"""
if on_disk:
return self.on_disk(short_names, error=error, glob_vars=glob_vars)
else:
return self.defines(short_names, error=error, glob_vars=glob_vars)
new_tree = self.copy()
to_update = new_tree
while to_update.parent is not None:
to_update = to_update.parent
to_update._update_partial_fill()
return new_tree
def _update_partial_fill(self, ):
"""
Helper function for `partial_fill` that updates the templates in place
"""
new_templates = {}
for short_name in self.templates:
template, variables = self.get_template(short_name)
new_templates[short_name] = str(utils.Template.parse(template).fill_known(variables))
self.templates = new_templates
for tree in self.sub_trees.values():
tree._update_partial_fill()
self.variables = {}
def copy(self, ):
"""
Copies the FileTree
Copies the templates, variables, sub_trees, and parent
:return: a copy of the FileTree
"""
return self._copy()
def _copy(self, new_parent=None, new_sub_tree=None):
"""
Helper function for copying a FileTree
"""
if new_sub_tree is None:
new_sub_tree = (None, None)
new_copy = type(self)(
templates=self.templates.copy(),
variables=self.variables.copy(),
name=self.name,
parent=new_parent
)
new_copy.sub_trees = {name: new_sub_tree[1] if new_sub_tree[0] == name else tree._copy(new_parent=new_copy)
for name, tree in self.sub_trees.items()}
if self.parent is not None and new_parent is None:
for my_key, ref_tree in self.parent.sub_trees.items():
if self is ref_tree:
break
else:
raise ValueError(f"Sub-tree {self} not found in parent tree")
new_copy._parent = self.parent._copy(new_sub_tree=(my_key, new_copy))
return new_copy
@classmethod
def read(cls, tree_name: str, directory='.', **variables) -> "FileTree":
def read(cls, tree_name: str, directory='.', partial_fill=False, **variables) -> "FileTree":
"""
Reads a FileTree from a specific file
......@@ -337,6 +435,7 @@ class FileTree(object):
:param tree_name: file containing the filename tree.
Can provide the filename of a tree file or the name for a tree in the ``filetree.tree_directories``.
:param directory: parent directory of the full tree (defaults to current directory)
:param partial_fill: By default any known `variables` are filled into the `template` immediately
:param variables: variable settings
:return: dictionary from specifier to filename
"""
......@@ -346,7 +445,9 @@ class FileTree(object):
filename = tree_name + '.tree'
else:
filename = parse.search_tree(tree_name)
tree_name = op.splitext(op.basename(filename))[0]
filename = Path(filename)
dirname = str(filename.parent)
templates = {}
nspaces_level = []
......@@ -379,11 +480,13 @@ class FileTree(object):
else:
sub_dir = current_filename.parents[0]
_, sub_tree, short_name = parse.read_subtree_line(line, sub_dir)
with parse.extra_tree_dirs([dirname]):
_, sub_tree, short_name = parse.read_subtree_line(line, sub_dir)
if short_name in sub_trees:
raise ValueError("Name of sub_tree {short_name} used multiple times in {tree_name}.tree".format(**locals()))
sub_trees[short_name] = sub_tree
sub_tree._name = short_name
elif '=' in line:
key, value = line.split('=')
if len(key.split()) != 1:
......@@ -413,9 +516,11 @@ class FileTree(object):
templates[short_name] = str(current_filename)
file_variables.update(variables)
res = get_registered(tree_name, cls)(templates, variables=file_variables, sub_trees=sub_trees)
res = get_registered(tree_name, cls)(templates, variables=file_variables, sub_trees=sub_trees, name=tree_name)
for tree in sub_trees.values():
tree._parent = res
if partial_fill:
res = res.partial_fill()
return res
......
fsl/utils/filetree/parse.py
View file @ 3b2a186e
import glob
import os.path as op
from . import filetree
from contextlib import contextmanager
from pathlib import PurePath
from typing import Tuple, List
import re
......@@ -9,6 +10,25 @@ import re
tree_directories = ['.', op.join(op.split(__file__)[0], 'trees')]
@contextmanager
def extra_tree_dirs(extra_dirs):
"""Temporarily insert ``extra_dirs`` to the beginning of :attr:`tree_directories`.
:arg extra_dirs: Sequence of additional tree file directories to search.
"""
global tree_directories
old_tree_directories = list(tree_directories)
tree_directories = list(extra_dirs) + list(tree_directories)
try:
yield
finally:
tree_directories = old_tree_directories
def search_tree(name: str) -> str:
"""
Searches for the file defining the specific tree
......@@ -57,7 +77,7 @@ def read_line(line: str) -> Tuple[int, PurePath, str]:
"""
Parses line from the tree file
:param line: input line from a \*.tree file
:param line: input line from a ``*.tree`` file
:return: Tuple with:
- number of spaces in front of the name
......@@ -85,7 +105,7 @@ def read_subtree_line(line: str, directory: str) -> Tuple[int, "filetree.FileTre
"""
Parses the line defining a sub_tree
:param line: input line from a \*.tree file
:param line: input line from a ``*.tree`` file
:param directory: containing directory
:return: Tuple with
......
fsl/utils/filetree/query.py
View file @ 3b2a186e
......@@ -6,7 +6,7 @@
# Author: Michiel Cottaar <michiel.cottaar@.ndcn.ox.ac.uk>
#
"""This module contains the :class:`FileTreeQuery` class, which can be used to
search for files in a directory described by a `.FileTree`. A
search for files in a directory described by a :class:`.FileTree`. A
``FileTreeQuery`` object returns :class:`Match` objects which each represent a
file that is described by the ``FileTree``, and which is present in the
directory.
......@@ -22,8 +22,9 @@ defined in this module:
"""
import logging
import collections
import logging
import collections
import functools as ft
import os.path as op
from typing import Dict, List, Tuple
......@@ -44,11 +45,10 @@ class FileTreeQuery(object):
by a :class:`.FileTree`, and identifies all file types (a.k.a. *templates*
or *short names*) that are present, and the values of variables within each
short name that are present. The :meth:`query` method can be used to
retrieve files which match a specific short name, and variable values.
retrieve files which match a specific template, and variable values.
The :meth:`query` method returns a multi-dimensional ``numpy.array``
which contains :class:`Match` objects, where each dimension one
represents variable for the short name in question.
The :meth:`query` method returns a collection of :class:`Match` objects,
each of which represents one file which matches the query.
Example usage::
......@@ -71,15 +71,13 @@ class FileTreeQuery(object):
'session': [None]}
>>> query.query('anat_image', participant='01')
array([[[[[[[Match(./my_bids_data/sub-01/anat/sub-01_T1w.nii.gz)],
[nan],
[nan],
[nan]]]],
[[[[Match(./my_bids_data/sub-01/anat/sub-01_T2w.nii.gz)],
[nan],
[nan],
[nan]]]]]]], dtype=object)
[Match(./my_bids_data/sub-01/anat/sub-01_T1w.nii.gz),
Match(./my_bids_data/sub-01/anat/sub-01_T2w.nii.gz)]
Matches for templates contained within sub-trees are referred to by
constructing a hierarchical path from the sub-tree template name(s),
and the template name - see the :meth:`Match.full_name` method.
"""
......@@ -90,118 +88,137 @@ class FileTreeQuery(object):
:arg tree: The :class:`.FileTree` object
"""
# Hard-code into the templates any pre-defined variables
tree = tree.partial_fill()
# Find all files present in the directory
# (as Match objects), and find all variables,
# plus their values, and all short names,
# plus their values, and all templates,
# that are present in the directory.
matches = scan(tree)
allvars, shortnamevars = allVariables(tree, matches)
matches = scan(tree)
allvars, templatevars = allVariables(tree, matches)
# Now we are going to build a series of ND
# arrays to store Match objects. We create
# one array for each short name. Each axis
# one array for each template. Each axis
# in an array corresponds to a variable
# present in files of that short name type,
# present in files of that template type,
# and each position along an axis corresponds
# to one value of that variable.
#
# These arrays will be used to store and
# retrieve Match objects - given a short
# name and a set of variable values, we
# can quickly find the corresponding Match
# retrieve Match objects - given a template
# and a set of variable values, we can
# quickly find the corresponding Match
# object (or objects).
# matcharrays contains {shortname : ndarray}
# matcharrays contains {template : ndarray}
# mappings, and varidxs contains
# {shortname : {varvalue : index}} mappings
# {template : {varvalue : index}} mappings
matcharrays = {}
varidxs = {}
for shortname in shortnamevars.keys():
for template, tvars in templatevars.items():
snvars = shortnamevars[shortname]
snvarlens = [len(allvars[v]) for v in snvars]
tvarlens = [len(allvars[v]) for v in tvars]
# "Scalar" match objects - templates
# which have no variables, and for
# which zero or one file is present
if len(tvarlens) == 0:
tvarlens = 1
# An ND array for this short
# name. Each element is a
# Match object, or nan.
matcharray = np.zeros(snvarlens, dtype=np.object)
matcharray = np.zeros(tvarlens, dtype=object)
matcharray[:] = np.nan
# indices into the match array
# for each variable value
snvaridxs = {}
for v in snvars:
snvaridxs[v] = {n : i for i, n in enumerate(allvars[v])}
tvaridxs = {}
for v in tvars:
tvaridxs[v] = {n : i for i, n in enumerate(allvars[v])}
matcharrays[shortname] = matcharray
varidxs[ shortname] = snvaridxs
matcharrays[template] = matcharray
varidxs[ template] = tvaridxs
# Populate the match arrays
for match in matches:
snvars = shortnamevars[match.short_name]
snvaridxs = varidxs[ match.short_name]
snarr = matcharrays[ match.short_name]
idx = []
for var in snvars:
tvars = templatevars[match.full_name]
tvaridxs = varidxs[ match.full_name]
tarr = matcharrays[ match.full_name]
idx = []
val = match.variables[var]
idx.append(snvaridxs[var][val])
if len(match.variables) == 0:
idx = [0]
else:
for var in tvars:
val = match.variables[var]
idx.append(tvaridxs[var][val])
snarr[tuple(idx)] = match
tarr[tuple(idx)] = match
self.__tree = tree
self.__allvars = allvars
self.__shortnamevars = shortnamevars
self.__templatevars = templatevars
self.__matches = matches
self.__matcharrays = matcharrays
self.__varidxs = varidxs
def axes(self, short_name) -> List[str]:
def axes(self, template) -> List[str]:
"""Returns a list containing the names of variables present in files
of the given ``short_name`` type, in the same order of the axes of
of the given ``template`` type, in the same order of the axes of
:class:`Match` arrays that are returned by the :meth:`query` method.
"""
return self.__shortnamevars[short_name]
return self.__templatevars[template]
def variables(self, short_name=None) -> Dict[str, List]:
def variables(self, template=None) -> Dict[str, List]:
"""Return a dict of ``{variable : [values]}`` mappings.
This dict describes all variables and their possible values in
the tree.
If a ``short_name`` is specified, only variables which are present in
files of that ``short_name`` type are returned.
If a ``template`` is specified, only variables which are present in
files of that ``template`` type are returned.
"""
if short_name is None:
if template is None:
return {var : list(vals) for var, vals in self.__allvars.items()}
else:
varnames = self.__shortnamevars[short_name]
varnames = self.__templatevars[template]
return {var : list(self.__allvars[var]) for var in varnames}
@property
def short_names(self) -> List[str]:
"""Returns a list containing all short names of the ``FileTree`` that
def tree(self):
"""Returns the :class:`.FileTree` associated with this
``FileTreeQuery``.
"""
return self.__tree
@property
def templates(self) -> List[str]:
"""Returns a list containing all templates of the ``FileTree`` that
are present in the directory.
"""
return list(self.__shortnamevars.keys())
return list(self.__templatevars.keys())
def query(self, short_name, asarray=False, **variables):
"""Search for files of the given ``short_name``, which match
def query(self, template, asarray=False, **variables):
"""Search for files of the given ``template``, which match
the specified ``variables``. All hits are returned for variables
that are unspecified.
:arg short_name: Short name of files to search for.
:arg template: Template of files to search for.
:arg asarray: If ``True``, the relevant :class:`Match` objects are
returned in a in a ND ``numpy.array`` where each
dimension corresponds to a variable for the
``short_name`` in question (as returned by
:meth:`axes`). Otherwise (the default), they are
returned in a list.
:arg asarray: If ``True``, the relevant :class:`Match` objects are
returned in a in a ND ``numpy.array`` where each
dimension corresponds to a variable for the
``templates`` in question (as returned by
:meth:`axes`). Otherwise (the default), they are
returned in a list.
All other arguments are assumed to be ``variable=value`` pairs,
used to restrict which matches are returned. All values are returned
......@@ -213,9 +230,9 @@ class FileTreeQuery(object):
"""
varnames = list(variables.keys())
allvarnames = self.__shortnamevars[short_name]
varidxs = self.__varidxs[ short_name]
matcharray = self.__matcharrays[short_name]
allvarnames = self.__templatevars[template]
varidxs = self.__varidxs[ template]
matcharray = self.__matcharrays[ template]
slc = []
for var in allvarnames:
......@@ -237,6 +254,7 @@ class FileTreeQuery(object):
else: return [m for m in result.flat if isinstance(m, Match)]
@ft.total_ordering
class Match(object):
"""A ``Match`` object represents a file with a name matching a template in
a ``FileTree``. The :func:`scan` function and :meth:`FileTree.query`
......@@ -244,16 +262,18 @@ class Match(object):
"""
def __init__(self, filename, short_name, variables):
def __init__(self, filename, template, tree, variables):
"""Create a ``Match`` object. All arguments are added as attributes.
:arg filename: name of existing file
:arg short_name: template identifier
:arg template: template identifier
:arg tree: :class:`.FileTree` which contains this ``Match``
:arg variables: Dictionary of ``{variable : value}`` mappings
containing all variables present in the file name.
"""
self.__filename = filename
self.__short_name = short_name
self.__template = template
self.__tree = tree
self.__variables = dict(variables)
......@@ -263,8 +283,38 @@ class Match(object):
@property
def short_name(self):
return self.__short_name
def template(self):
return self.__template
@property
def full_name(self):
"""The ``full_name`` of a ``Match`` is a combination of the
``template`` (i.e. the matched template), and the name(s) of
the relevant ``FileTree`` objects.
It allows one to unamiguously identify the location of a ``Match``
in a ``FileTree`` hierarchy, where the same ``short_name`` may be
used in different sub-trees.
"""
def parents(tree):
if tree.parent is None:
return []
else:
return [tree.parent] + parents(tree.parent)
trees = [self.tree] + parents(self.tree)
# Drop the root tree
trees = list(reversed(trees))[1:]
return '/'.join([t.name for t in trees] + [self.template])
@property
def tree(self):
return self.__tree
@property
......@@ -275,7 +325,8 @@ class Match(object):
def __eq__(self, other):
return (isinstance(other, Match) and
self.filename == other.filename and
self.short_name == other.short_name and
self.template == other.template and
self.tree is other.tree and
self.variables == other.variables)
......@@ -283,13 +334,9 @@ class Match(object):
return isinstance(other, Match) and self.filename < other.filename
def __le__(self, other):
return isinstance(other, Match) and self.filename <= other.filename
def __repr__(self):
"""Returns a string representation of this ``Match``. """
return 'Match({})'.format(self.filename)
return 'Match({}: {})'.format(self.full_name, self.filename)
def __str__(self):
......@@ -301,19 +348,21 @@ def scan(tree : FileTree) -> List[Match]:
"""Scans the directory of the given ``FileTree`` to find all files which
match a tree template.
:return: list of :class:`Match` objects
:arg tree: :class:`.FileTree` to scan
:returns: list of :class:`Match` objects
"""
matches = []
for template in tree.templates:
for filename in tree.get_all(template, glob_vars='all'):
for variables in tree.get_all_vars(template, glob_vars='all'):
filename = tree.update(**variables).get(template)
if not op.isfile(filename):
continue
variables = dict(tree.extract_variables(template, filename))
matches.append(Match(filename, template, variables))
matches.append(Match(filename, template, tree, variables))
for tree_name, sub_tree in tree.sub_trees.items():
matches.extend(scan(sub_tree))
......@@ -336,26 +385,29 @@ def allVariables(
variables and their possible values present in the given list
of ``Match`` objects.
- A dict of ``{ short_name : [variables] }`` mappings,
containing the variables which are relevant to each short
name.
- A dict of ``{ full_name : [variables] }`` mappings,
containing the variables which are relevant to each template.
"""
allvars = collections.defaultdict(set)
allshortnames = collections.defaultdict(set)
allvars = collections.defaultdict(set)
alltemplates = {}
for m in matches:
if m.full_name not in alltemplates:
alltemplates[m.full_name] = set()
for var, val in m.variables.items():
allvars[ var] .add(val)
allshortnames[m.short_name].add(var)
allvars[ var] .add(val)
alltemplates[m.full_name].add(var)
# allow us to compare None with strings
def key(v):
if v is None: return ''
else: return v
allvars = {var : list(sorted(vals, key=key))
for var, vals in allvars.items()}
allshortnames = {sn : list(sorted(vars))
for sn, vars in allshortnames.items()}
allvars = {var : list(sorted(vals, key=key))
for var, vals in allvars.items()}
alltemplates = {tn : list(sorted(vars))
for tn, vars in alltemplates.items()}
return allvars, allshortnames
return allvars, alltemplates
fsl/utils/filetree/trees/bids_raw.tree
View file @ 3b2a186e
......@@ -2,19 +2,22 @@ ext=.nii.gz
dataset_description.json
participants.tsv
README
CHANGES
README (readme)
CHANGES (changes)
LICENSE (license)
genetic_info.json
sub-{participant}
[ses-{session}]
sub-{participant}_sessions.tsv (sessions_tsv)
anat (anat_dir)
sub-{participant}[_ses-{session}][_acq-{acq}][_rec-{rec}][_run-{run_index}]_{modality}{ext} (anat_image)
sub-{participant}[_ses-{session}][_acq-{acq}][_ce-{ce}][_rec-{rec}][_run-{run_index}]_{modality}{ext} (anat_image)
sub-{participant}[_ses-{session}][_acq-{acq}][_ce-{ce}][_rec-{rec}][_run-{run_index}][_mod-{modality}]_defacemask{ext} (anat_deface)
func (func_dir)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_rec-{rec}][_run-{run_index}]_bold{ext} (task_image)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_rec-{rec}][_run-{run_index}]_sbref{ext} (task_sbref)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_rec-{rec}][_run-{run_index}]_events.tsv (task_events)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_rec-{rec}][_run-{run_index}][_recording-{recording}]_physio{ext} (task_physio)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_rec-{rec}][_run-{run_index}][_recording-{recording}]_stim{ext} (task_stim)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_ce-{ce}][_dir-{dir}][_rec-{rec}][_run-{run_index}][_echo-{echo}]_bold{ext} (task_image)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_ce-{ce}][_dir-{dir}][_rec-{rec}][_run-{run_index}][_echo-{echo}]_sbref{ext} (task_sbref)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_ce-{ce}][_dir-{dir}][_rec-{rec}][_run-{run_index}][_echo-{echo}]_events.tsv (task_events)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_ce-{ce}][_dir-{dir}][_rec-{rec}][_run-{run_index}][_echo-{echo}][_recording-{recording}]_physio.tsv.gz (task_physio)
sub-{participant}[_ses-{session}]_task-{task}[_acq-{acq}][_ce-{ce}][_dir-{dir}][_rec-{rec}][_run-{run_index}][_echo-{echo}][_recording-{recording}]_stim.tsv.gz (task_stim)
dwi (dwi_dir)
sub-{participant}[_ses-{session}][_acq-{acq}][_run-{run_index}]_dwi{ext} (dwi_image)
sub-{participant}[_ses-{session}][_acq-{acq}][_run-{run_index}]_dwi.bval (bval)
......@@ -28,3 +31,15 @@ sub-{participant}
sub-{participant}[_ses-{session}][_acq-{acq}][_run-{run_index}]_phase2{ext} (fmap_phase2)
sub-{participant}[_ses-{session}][_acq-{acq}][_run-{run_index}]_fieldmap{ext} (fmap)
sub-{participant}[_ses-{session}][_acq-{acq}]_dir-{dir}[_run-{run_index}]_epi{ext} (fmap_epi)
meg (meg_dir)
sub-{participant}[_ses-{session}]_task-{task}[_run-{run}][_proc-{proc}]_meg.{meg_ext} (meg)
eeg (eeg_dir)
sub-{participant}[_ses-{session}]_task-{task}[_run-{run}][_proc-{proc}]_eeg.{eeg_ext} (eeg)
ieeg (ieeg_dir)
sub-{participant}[_ses-{session}]_task-{task}[_run-{run}][_proc-{proc}]_ieeg.{ieeg_ext} (ieeg)
beh (behavioral_dir)
sub-{participant}[_ses-{session}]_task-{task}_events.tsv (behavioural_events)
sub-{participant}[_ses-{session}]_task-{task}_beh.tsv (behavioural)
sub-{participant}[_ses-{session}]_task-{task}_physio.tsv.gz (behavioural_physio)
sub-{participant}[_ses-{session}]_task-{task}_stim.tsv.gz (behavioral_stim)