Skip to content
Snippets Groups Projects
Commit ed5835ab authored by Paul McCarthy's avatar Paul McCarthy :mountain_bicyclist:
Browse files

RF: Update x5 i/o code to adhere to new spec; fix bug in 

nonlinear.DeformationField
parent 4e6e03ca
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fsl/transform/nonlinear.py
+ 2
− 2
View file @ ed5835ab
...@@ -267,8 +267,8 @@ class DeformationField(NonLinearTransform): ...@@ -267,8 +267,8 @@ class DeformationField(NonLinearTransform):
# can get this through the assumption # can get this through the assumption
# that field and ref are aligned in # that field and ref are aligned in
# the world coordinate system # the world coordinate system
xform = affine.concat(self .getAffine('world', 'voxel'), xform = affine.concat(self .getAffine('world', 'voxel'),
self.ref.getAffine(from_, 'world')) self.ref.getAffine(self.refSpace, 'world'))
if np.all(np.isclose(xform, np.eye(4))): if np.all(np.isclose(xform, np.eye(4))):
voxels = coords voxels = coords
......
fsl/transform/x5.py
+ 91
− 227
View file @ ed5835ab
...@@ -350,11 +350,16 @@ def readLinearX5(fname): ...@@ -350,11 +350,16 @@ def readLinearX5(fname):
- A :class:`.Nifti` instance representing the source space - A :class:`.Nifti` instance representing the source space
- A :class:`.Nifti` instance representing the reference space - A :class:`.Nifti` instance representing the reference space
""" """
with h5py.File(fname, 'r') as f: with h5py.File(fname, 'r') as f:
_validateMetadata( f['/'])
xform = _readAffine( f['/']) if f.attrs['Type'] != 'linear':
src = _readSpace( f['/From']) raise X5Error('Not a linear transform')
ref = _readSpace( f['/To'])
_readMetadata( f['/'])
xform = _readAffine(f['/Transform'])
src = _readSpace( f['/A'])
ref = _readSpace( f['/B'])
return xform, src, ref return xform, src, ref
...@@ -369,14 +374,13 @@ def writeLinearX5(fname, xform, src, ref): ...@@ -369,14 +374,13 @@ def writeLinearX5(fname, xform, src, ref):
""" """
with h5py.File(fname, 'w') as f: with h5py.File(fname, 'w') as f:
_writeMetadata(f)
_writeAffine(f, xform)
from_ = f.create_group('/From') f.attrs['Type'] = 'linear'
to = f.create_group('/To')
_writeSpace(from_, src) _writeMetadata(f)
_writeSpace(to, ref) _writeAffine( f.create_group('/Transform'), xform)
_writeSpace( f.create_group('/A'), src)
_writeSpace( f.create_group('/B'), ref)
def readNonLinearX5(fname): def readNonLinearX5(fname):
...@@ -388,15 +392,21 @@ def readNonLinearX5(fname): ...@@ -388,15 +392,21 @@ def readNonLinearX5(fname):
with h5py.File(fname, 'r') as f: with h5py.File(fname, 'r') as f:
root = f['/'] if f.attrs.get('Type') != 'nonlinear':
_validateNonLinearTransform(root) raise X5Error('Not a nonlinear transform')
subtype = root.attrs['SubType'] _readMetadata(f)
if subtype == 'displacement': field = _readDisplacementField(root) ref = _readSpace( f['/A'])
elif subtype == 'coefficient': field = _readCoefficientField(root) src = _readSpace( f['/B'])
field, space = _readDeformation(f['/Transform'])
return field return nonlinear.DeformationField(field,
header=space.header,
src=src,
ref=ref,
srcSpace='world',
refSpace='world')
def writeNonLinearX5(fname, field): def writeNonLinearX5(fname, field):
...@@ -407,60 +417,64 @@ def writeNonLinearX5(fname, field): ...@@ -407,60 +417,64 @@ def writeNonLinearX5(fname, field):
""" """
with h5py.File(fname, 'w') as f: with h5py.File(fname, 'w') as f:
_writeMetadata(f)
f.attrs['Type'] = 'nonlinear' f.attrs['Type'] = 'nonlinear'
if isinstance(field, nonlinear.DisplacementField): _writeMetadata(f)
_writeDisplacementField(f, field) _writeSpace( f.create_group('/A'), f.ref)
elif isinstance(field, nonlinear.CoefficientField): _writeSpace( f.create_group('/B'), f.src)
_writeCoefficientField(f, field) _writeDeformation(f.create_group('/Transform'), 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__})
def _validateMetadata(group): def _readMetadata(group):
"""Reads a metadata block from the given group, and raises a :exc:`X5Error` """Reads a metadata block from the given group, and raises a :exc:`X5Error`
if it does not look valid. if it does not look valid.
:arg group: A ``h5py.Group`` object. :arg group: A ``h5py.Group`` object.
:returns: A ``dict`` containing the metadata.
""" """
try: format = group.attrs.get('Format')
format = group.attrs['Format'] version = group.attrs.get('Version')
version = group.attrs['Version'] meta = group.attrs.get('Metadata')
except Exception:
raise X5Error('File does not appear to be an X5 file')
if (format != X5_FORMAT) or (version != X5_VERSION): if (format != X5_FORMAT) or (version != X5_VERSION):
raise X5Error('Incompatible format/version (required: {}/{}, ' raise X5Error('Incompatible format/version (required: {}/{}, '
'present: {}/{})'.format(X5_FORMAT, X5_VERSION, 'present: {}/{})'.format(X5_FORMAT, X5_VERSION,
format, 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' : version.__version__})
def _readAffine(group): def _readAffine(group):
"""Reads an affine from the given group. """Reads an *affine* from the given group.
:arg group: A ``h5py.Group`` object. :arg group: A ``h5py.Group`` object.
:returns: ``numpy`` array containing a ``(4, 4)`` affine transformation. :returns: ``numpy`` array containing a ``(4, 4)`` affine transformation.
""" """
if group.attrs['Type'] != 'linear': if group.attrs.get('Type') != 'affine':
raise X5Error('Not a linear transform') raise X5Error('Not an affine')
xform = group['Transform'] xform = group['Matrix']
if xform.shape != (4, 4): if xform.shape != (4, 4):
raise X5Error('Not a linear transform') raise X5Error('Not an affine')
return np.array(xform) return np.array(xform)
...@@ -476,21 +490,19 @@ def _writeAffine(group, xform): ...@@ -476,21 +490,19 @@ def _writeAffine(group, xform):
xform = np.asarray(xform, dtype=np.float64) xform = np.asarray(xform, dtype=np.float64)
inv = np.asarray(affine.invert(xform), dtype=np.float64) inv = np.asarray(affine.invert(xform), dtype=np.float64)
group.attrs['Type'] = 'linear' group.attrs['Type'] = 'affine'
group.create_dataset('Transform', data=xform) group.create_dataset('Matrix', data=xform)
group.create_dataset('Inverse', data=inv) group.create_dataset('Inverse', data=inv)
def _readSpace(group): 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. :arg group: A ``h5py.Group`` object.
:returns: :class:`.Nifti` object. defining the mapping :returns: :class:`.Nifti` object. defining the mapping
""" """
import fsl.data.image as fslimage if group.attrs.get('Type') != 'image':
if group.attrs['Type'] != 'image':
raise X5Error('Not an image space') raise X5Error('Not an image space')
shape = group.attrs['Size'] shape = group.attrs['Size']
...@@ -518,198 +530,50 @@ def _writeSpace(group, img): ...@@ -518,198 +530,50 @@ def _writeSpace(group, img):
_writeAffine(mapping, img.getAffine('voxel', 'world')) _writeAffine(mapping, img.getAffine('voxel', 'world'))
def _validateNonLinearTransform(group): def _readDeformation(group):
"""Checks that the attributes of the given group, assumed to contain a """Reads a *deformation* from the given group.
nonlinear transform, are valid. Raises a :exc:`X5Error` if not.
:arg group: A ``h5py.Group`` object.
"""
type = group.attrs['Type']
subtype = group.attrs['SubType']
repr = group.attrs['Representation']
if type != 'nonlinear':
raise X5Error('Not a nonlinear transform')
if subtype not in ('displacement', 'coefficient'):
raise X5Error('Unrecognised nonlinear subtype: {}'.format(subtype))
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. :arg group: A ``h5py.Group`` object
:returns: A tuple containing
: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`
"""
src = _readSpace(group['From'])
ref = _readSpace(group['To'])
pre = group.get('Pre')
post = group.get('Post')
init = group.get('InitialAlignment')
if pre is not None: pre = _readAffine(pre)
if post is not None: post = _readAffine(post)
if init is not None: init = _readAffine(init)
refSpace = 'world' - A ``numpy.arrayThe`` containing the deformation field
srcSpace = 'world'
try: - A :class:`.Nifti` object representing the deformation
if pre is not None: field space
refSpace = fslimage.Nifti.identifyAffine(
ref, pre, from_='world')[1]
if post is not None:
srcSpace = fslimage.Nifti.identifyAffine(
src, post, to='world')[ 0]
except ValueError: - The deformation type - either ``'absolute'`` or
raise X5Error('Invalid pre/post affine') ``'relative'``
return src, ref, pre, post, init, srcSpace, refSpace
def _writeNonLinearCommon(group, field):
"""Writes the spaces and affines for the given nonlinear transform to the
given group.
:arg group: A ``h5py.Group`` object.
:arg field: A :class:`.NonLinearTransform` object
""" """
_writeSpace(group.create_group('From'), field.src) type = group.attrs.get('Type')
_writeSpace(group.create_group('To'), field.ref) subtype = group.attrs['SubType']
pre = field.ref.getAffine('world', field.refSpace) if type != 'deformation':
post = field.src.getAffine(field.srcSpace, 'world') raise X5Error('Not a deformation')
_writeAffine(group.create_group('Pre'), pre) if subtype not in ('absolute', 'relative'):
_writeAffine(group.create_group('Post'), post) raise X5Error('Unknown deformation type: {}'.format(subtype))
if field.srcToRefMat is not None: mapping = _readSpace(group['Mapping'])
_writeAffine(group.create_group('InitialAlignment'), field.srcToRefMat) field = group['Matrix']
if len(field.shape) != 4 or field.shape[3] != 3:
raise X5Error('Invalid shape for deformation field')
def _readDisplacementField(group): return np.array(field), mapping, subtype
"""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 _writeDeformation(group, field):
"""Write a deformation fieldto the given group.
def _writeDisplacementField(group, field): :arg group: A ``h5py.Group`` object
"""Writes a nonlinear displacement field to the given group. :arg field: A :class:`.DeformationField` object
:arg group: A ``h5py.Group`` object.
:arg field: A :class:`.DisplacementField` object
""" """
_writeNonLinearCommon(group, field) group.attrs['Type'] = 'deformation'
group.attrs['SubType'] = field.deformationType
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. mapping = group.create_group('Mapping')
: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')
params.attrs['Spacing'] = field.knotSpacing group.create_dataset('Matrix', data=field.data)
_writeAffine(refToField, field.refToFieldMat) _writeSpace(mapping, field)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment