Newer
Older
#!/usr/bin/env python
#
# x5.py - Functions for working with BIDS X5 files.
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
"""This module contains functions for reading/writing linear/non-linear FSL
transformations from/to BIDS X5 files.
An X5 file is a HDF5 container file which stores a linear or non-linear
transformation between a **source** NIfTI image and a **refernece** image. In
an X5 file, the source image is referred to as the ``From`` space, and the
reference image the ``To`` space.
Custom HDF5 groups
------------------
HDF5 files are composed of *groups*, *attributes*, and *datasets*. Groups are
simply containers for attributes, datasets, and other groups. Attributes are
strongly-typed scalar values. Datasets are strongly-typed, structured
N-dimensional arrays.
To simplify the file format efinitions below, we shall first define a few
custom HDF5 groups. In the file format definitions, a HDF5 group which is
listed as being of one of these custom types shall contain the
attributes/datasets that are listed here.
*affine*
^^^^^^^^
A HDF5 group which is listed as being of type "affine" contains an affine
transformation, which can be used to transform coordinates from one space into
another. Groups of type "affine" have the following fields:
+---------------+-----------+-------------------------------------------------+
| Name | Type | Value/Description |
+---------------+-----------+-------------------------------------------------+
| ``Type`` | attribute | ``'linear'`` |
| ``Transform`` | dataset | The transformation itself - a 4x4 ``float64`` |
| | | affine transformation |
| ``Inverse`` | dataset | Optional pre-calculated inverse |
+---------------+-----------+-------------------------------------------------+
*mapping*
^^^^^^^^^
A HDF5 group which is listed as being of type "mapping" contains all of the
information required to define the space of a NIfTI image, including its
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 "mapping" have the following fields:
+-------------+-----------+---------------------------------------------------+
| Name | Type | Value/Description |
+-------------+-----------+---------------------------------------------------+
| ``Type`` | attribute | ``'image'`` |
| ``Size`` | attribute | ``uint64`` ``(X, Y, Z)`` voxel dimensions |
| ``Scales`` | attribute | ``float64`` ``(X, Y, Z)`` voxel pixdims |
| ``Mapping`` | affine | The image voxel-to-world transformation (its |
| | | "sform") |
+-------------+-----------+---------------------------------------------------+
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**.
Linear X5 transformation files are assumed to adhere to the HDF5 structure
defined in the table below. All fields are required.
+-----------------------------+-----------+-----------------------------------+
| Name | Type | Value/Description |
+-----------------------------+-----------+-----------------------------------+
| *Metadata* |
+-----------------------------+-----------+-----------------------------------+
| ``/Format`` | attribute | ``'X5'`` |
| ``/Version`` | attribute | ``'0.0.1'`` |
| ``/Metadata`` | attribute | JSON string containing |
| | | unstructured metadata. |
+-----------------------------+-----------+-----------------------------------+
| *Transformation* |
+-----------------------------+-----------+-----------------------------------+
| ``/From/`` | mapping | Source image mapping |
| ``/To/`` | mapping | Reference image mapping |
| ``/`` | affine | Affine transformation from source |
| | | image world coordinates to |
| | | reference image world coordinates |
+-----------------------------+-----------+-----------------------------------+
Non-linear X5 transformation files
----------------------------------
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:
- A *displacement field*, which is defined in the same space as the reference
image, and which contains displacements.
- A quadratic or cubic B-spline *coefficient field*, which contains
coefficients defined in a coarse grid which overlays the same space as the
reference image, and from which a displacement field can be calculated.
Displacement fields may contain either:
- *relative displacements*, which contains displacements *from* the
reference image coordinates *to* the source image coordinates (i.e. add the
displacements to the reference image coordinates to get the corresponding
source image coordinates)
- *absolute displacement* which simply contain the source image coordinates
at each voxel in the reference image space.
Displacement field transformations define displacements between the refernece
image *world* coordinate system, and the source image *world* coordinate
system.
The displacement field generated from a coefficient field defines relative
displacements from a reference image space to a source image space. These
spaces are undefined; however, the ``/Pre`` affine transformation may be used
to transform reference image world coordinates into the reference image space,
and the ``/Post`` affine transformation may be used to transform the resulting
source image coordinates into the source image world coordinate system.
Non-linear X5 transformation files are assumed to adhere to the following
HDF5 structure::
+-----------------------------+-----------+-----------------------------------+
| Name | Type | Value/Description |
+-----------------------------+-----------+-----------------------------------+
| *Metadata* - same as for linear transformation files |
+-----------------------------------------------------------------------------+
| *Transformation* |
+-----------------------------+-----------+-----------------------------------+
| ``/Type`` | attribute | ``'nonlinear'`` |
| ``/SubType`` | attribute | ``'displacement'`` or |
| | | ``'coefficient'`` |
| ``/Representation`` | attribute | If ``/SubType`` is |
| | | ``'displacement'``, its |
| | | ``/Representation`` may be either |
| | | ``'absolute'`` or ``'relative'``. |
| | | If ``/SubType`` is |
| | | ``'coefficient'``, its |
| | | ``/Representation`` may be either |
| | | ``'quadratic bspline'`` or |
| | | ``'cubic bspline'``. |
| ``/Transform`` | dataset | The transformation itself - see |
| | | above. |
| ``/Inverse`` | dataset | Optional pre-calculated inverse |
| ``/From/`` | mapping | Source image mapping |
| ``/To/`` | mapping | Reference image mapping |
|-----------------------------+-----------+-----------------------------------+
| *Coefficient field parameters* (required for ``'coefficient'`` files) |
|-----------------------------+-----------+-----------------------------------+
| ``/Parameters/Spacing`` | attribute |
| ``/Parameters/ReferenceToField`` | affine | | attribute |
+-----------------------------+-----------+-----------------------------------+
"""
import json
import numpy as np
import nibabel as nib
import h5py
import fsl.version as version
from . import affine
X5_FORMAT = 'X5'
X5_VERSION = '0.0.1'
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 _readLinearTransform(group):
"""Reads a linear transformation 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 ValueError('Not a linear transform')
xform = np.array(group['Transform'])
if xform.shape != (4, 4):
raise ValueError('Not a linear transform')
return xform
def _writeLinearTransform(group, xform):
"""Writes the given linear transformation and its inverse to the given
group.
:arg group: A ``h5py.Group`` object.
:arg xform: ``numpy`` array containing a ``(4, 4)`` affine transformation.
"""
xform = np.asarray(xform, dtype=np.float64)
inv = np.asarray(affine.inverse(xform), dtype=np.float64)
group.attrs['Type'] = 'linear'
group.create_dataset('Transform', data=xform)
group.create_dataset('Inverse', data=inv)
def _readLinearMapping(group):
"""Reads a linear mapping, defining a source or reference space, from the
given group.
:arg group: A ``h5py.Group`` object.
:returns: :class:`.Nifti` object. defining the mapping
"""
import fsl.data.image as fslimage
if group.attrs['Type'] != 'image':
raise ValueError('Not an image mapping')
shape = group.attrs['Size']
pixdim = group.attrs['Scales']
xform = _readLinearTransform(group['Mapping'])
hdr = nib.Nifti2Header()
hdr.set_data_shape(shape)
hdr.set_zooms( pixdim)
hdr.set_sform( xform, 'aligned')
return fslimage.Nifti(hdr)
def _writeLinearMapping(group, img):
"""Writes a linear mapping specified by ``img`` to the given group.
:arg group: A ``h5py.Group`` object.
:arg img: :class:`.Nifti` object. defining the mapping
"""
group.attrs['Type'] = 'image'
group.attrs['Size'] = np.asarray(img.shape[ :3], np.uint32)
group.attrs['Scales'] = np.asarray(img.pixdim[:3], np.float32)
mapping = group.create_group('Mapping')
_writeLinearTransform(mapping, img.getAffine('voxel', 'world'))
def _readNonLinearTransform(group):
if group.attrs['Type'] != 'nonlinear':
raise ValueError('Not a nonlinear transform')
return np.array(group['Transform'])
def _writeNonLinearTransform(group, field):
"""
"""
group.attrs['Type'] = 'nonlinear'
group.create_dataset('Transform', data=field, dtype=np.float32)
def readLinearX5(fname):
"""
"""
with h5py.File(fname, 'r') as f:
xform = _readLinearTransform(f['/'])
src = _readLinearMapping( f['/From'])
ref = _readLinearMapping( f['/To'])
return xform, src, ref
def writeLinearX5(fname, xform, src, ref):
"""
with h5py.File(fname, 'w') as f:
_writeMetadata(f)
_writeLinearTransform(f, xform)
from_ = f.create_group('/From')
to = f.create_group('/To')
_writeLinearMapping(from_, src)
_writeLinearMapping(to, ref)
def readNonLinearX5(fname):
"""
"""
from . import nonlinear
with h5py.File(fname, 'r') as f:
field = _readNonLinearTransform(f['/'])
src = _readLinearMapping(f['/From'])
ref = _readLinearMapping(f['/To'])
# TODO coefficient fields
return nonlinear.DisplacementField(field,
src=src,
ref=ref,
srcSpace='world',
refSpace='world')
def writeNonLinearX5(fname, field):
"""
::
/Format # "X5"
/Version # "0.0.1"
/Metadata # json string containing unstructured metadata
/Transform # the displacement/coefficient field itself
Paul McCarthy
committed
/Type # "nonlinear"
/SubType # "displacement" / "deformation"
/Representation # "cubic bspline" / "quadratic bspline"
/Inverse # optional pre-calculated inverse
Paul McCarthy
committed
/Pre/Type # "linear"
/Pre/Transform # ref world-to-[somewhere], to prepare ref
# world coordinates as inputs to the nonlinear
# transform
/Pre/Inverse # optional pre-calculated inverse
/Post/Type # "linear"
/Post/Transform # source [somewhere]-to-world, to transform
# source coordinates produced by the nonlinear
# transform into world coordinates
/Post/Inverse # optional pre-calculated inverse
/From/Type # "image"
/From/Size # voxel dimensions
/From/Scales # voxel pixdims
/From/Mapping/Transform # source voxel-to-world sform
Paul McCarthy
committed
/From/Mapping/Type # "linear"
/From/Mapping/Inverse # optional inverse
/To/Type # "image"
/To/Size # voxel dimensions
/To/Scales # voxel pixdims
/To/Mapping/Transform # reference voxel-to-world sform
Paul McCarthy
committed
/To/Mapping/Type # "linear"
/To/Mapping/Inverse # optional inverse
Paul McCarthy
committed
""" # noqa
# TODO coefficient fields
with h5py.File(fname, 'w') as f:
_writeMetadata(f)
_writeNonLinearTransform(f, field.data)
from_ = f.create_group('/From')
to = f.create_group('/To')
_writeLinearMapping(from_, field.src)
_writeLinearMapping(to, field.ref)