fnirt.py 13.3 KB
Newer Older
1
2
3
4
5
6
7
#!/usr/bin/env python
#
# fnirt.py - Functions for working with FNIRT non-linear transformations.
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module contains functions for working with FNIRT non-linear
8
transformations. The following functions are available:
9
10
11
12
13

.. autosummary::
   :nosignatures:

   readFnirt
14
15
   toFnirt
   fromFnirt
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


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.

85
86
87
88
89
90
.. 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.
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

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
128
129
130
131
132
"""


import logging

133
import nibabel as nib
134
import numpy   as np
135

136
import fsl.data.constants as constants
137
138
139
import fsl.data.image     as fslimage
from . import                affine
from . import                nonlinear
140
141
142
143
144


log = logging.getLogger(__name__)


145
146
def _readFnirtDeformationField(fname, img, src, ref, defType=None):
    """Loads ``fname``, assumed to be a FNIRT deformation field.
147

148
149
150
151
152
153
154
155
156
157
158
159
    :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
160
    """
161
162
163
164
165
166
    return nonlinear.DeformationField(fname,
                                      src,
                                      ref,
                                      srcSpace='fsl',
                                      refSpace='fsl',
                                      defType=defType)
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212


def _readFnirtCoefficientField(fname, img, src, ref):
    """Loads ``fname``, assumed to be a FNIRT coefficient field.

    :arg fname: File name of FNIRT coefficient field
    :arg img:   ``fname`` loaded as an :class:`.Image`
    :arg src:   Source image
    :arg ref:   Reference image
    :return:    A :class:`.CoefficientField`
    """

    # FNIRT uses NIFTI header fields in
    # non-standard ways to store some
    # additional information about the
    # coefficient field. See
    # $FSLDIR/src/fnirt/fnirt_file_writer.cpp
    # for more details.

    # The field type (quadratic, cubic,
    # or discrete-cosine-transform) is
    # inferred from the intent. There is
    # no support in this implementation
    # for DCT fields
    cubics = (constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
              constants.FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS)
    quads  = (constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS,
              constants.FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS)

    if img.intent in cubics:
        fieldType = 'cubic'
    elif img.intent in quads:
        fieldType = 'quadratic'
    else:
        fieldType = 'cubic'
        log.warning('Unrecognised/unsupported coefficient '
                    'field type (assuming cubic b-spline): '
                    '{}'.format(img.intent))

    # Knot spacing (in voxels) is
    # stored in the pixdims
    knotSpacing = img.pixdim[:3]

    # The sform contains an initial
    # global src-to-ref affine
    # (the starting point for the
213
214
215
216
217
    # non-linear registration). This
    # is encoded as a flirt matrix,
    # i.e. it transforms from
    # source-scaled-voxels to
    # ref-scaled-voxels
218
219
    srcToRefMat = img.header.get_sform()

220
221
222
223
224
225
    # The fieldToRefMat affine tells
    # the CoefficientField class how
    # to transform coefficient field
    # voxel coordinates into
    # displacement field/reference
    # image voxel coordinates.
226
227
228
229
230
231
232
233
234
235
236
237
238
    fieldToRefMat = affine.scaleOffsetXform(knotSpacing, 0)

    return nonlinear.CoefficientField(fname,
                                      src,
                                      ref,
                                      srcSpace='fsl',
                                      refSpace='fsl',
                                      fieldType=fieldType,
                                      knotSpacing=knotSpacing,
                                      srcToRefMat=srcToRefMat,
                                      fieldToRefMat=fieldToRefMat)


239
def readFnirt(fname, src, ref, defType=None):
240
    """Reads a non-linear FNIRT transformation image, returning
241
    a :class:`.DeformationField` or :class:`.CoefficientField` depending
242
243
    on the file type.

Paul McCarthy's avatar
Paul McCarthy committed
244
245
246
    :arg fname:    File name of FNIRT transformation
    :arg src:      Source image
    :arg ref:      Reference image
247
248
249
    :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.
250
251
252
    """

    # Figure out whether the file
253
    # is a deformation field or
254
    # a coefficient field
255
256
257
258
259
260
261
262
263
264
265

    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:
266
        return _readFnirtDeformationField(fname, img, src, ref, defType)
267
268
269
270
271
    elif img.intent in coefs:
        return _readFnirtCoefficientField(fname, img, src, ref)
    else:
        raise ValueError('Cannot determine type of nonlinear '
                         'file {}'.format(fname))
272
273


274
275
def toFnirt(field):
    """Convert a :class:`.NonLinearTransform` to a FNIRT-compatible
276
    :class:`.DeformationField` or :class:`.CoefficientField`.
277

278
    :arg field: :class:`.NonLinearTransform` to convert
279
    :return:    A FNIRT-compatible :class:`.DeformationField` or
280
                :class:`.CoefficientField`.
281
    """
282

283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
    # If we have a coefficient field
    # which transforms between fsl
    # space, we can just create a copy.
    if isinstance(field, nonlinear.CoefficientField) and \
       (field.srcSpace == 'fsl' and field.refSpace == 'fsl'):

        # We start with a nibabel image,
        # as we need to mess with the header
        # fields to store all of the FNIRT
        # coefficient field information
        fieldBase = nib.nifti1.Nifti1Image(field.data, None)

        # Set the intent
        if field.fieldType == 'cubic':
            intent = constants.FSL_CUBIC_SPLINE_COEFFICIENTS
        elif field.fieldType == 'quadratic':
            intent = constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS
        fieldBase.header['intent_code'] = intent

        # Reference image pixdims are
        # stored in the intent code
        # parameters.
        fieldBase.header['intent_p1'] = field.ref.pixdim[0]
        fieldBase.header['intent_p2'] = field.ref.pixdim[1]
        fieldBase.header['intent_p3'] = field.ref.pixdim[2]

        # Pixdims are used to
        # store the knot spacing,
        pixdims      = list(field.knotSpacing) + [1]
        qform        = np.diag(pixdims)

        # The sform is used to store the
        # initial src-to-ref affine
        if field.srcToRefMat is not None: sform = field.srcToRefMat
        else:                             sform = np.eye(4)

        # The qform offsets are
        # used to store the
        # reference image shape
        qform[:3, 3] = field.ref.shape[:3]

        fieldBase.header.set_zooms(pixdims)
        fieldBase.set_sform(sform, 1)
        fieldBase.set_qform(qform, 1)
        fieldBase.update_header()

        field = nonlinear.CoefficientField(
            fieldBase,
            src=field.src,
            ref=field.ref,
            srcSpace='fsl',
            refSpace='fsl',
            fieldType=field.fieldType,
            knotSpacing=field.knotSpacing,
            fieldToRefMat=field.fieldToRefMat,
            srcToRefMat=field.srcToRefMat)

    # Otherwise we have a non-FSL coefficient
341
    # field, or a deformation field.
342
    else:
343

344
345
346
347
348
349
350
351
        # 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.
352
        if isinstance(field, nonlinear.CoefficientField):
353
            field = nonlinear.coefficientFieldToDeformationField(field)
354
355
356
357
358

        # 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':
359
            field = nonlinear.DeformationField(
360
                field.data,
361
                header=field.header,
362
363
                src=field.src,
                ref=field.ref,
364
                defType=field.deformationType)
365
366
367
368

        # Otherwise we have to adjust the
        # displacements so they transform
        # between fsl coordinates.
369
        field = nonlinear.convertDeformationSpace(
370
371
372
            field, from_='fsl', to='fsl')

        field.header['intent_code'] = constants.FSL_FNIRT_DISPLACEMENT_FIELD
373

374
    return field
375
376


377
378
def fromFnirt(field, from_='world', to='world'):
    """Convert a FNIRT-style :class:`.NonLinearTransform` to a generic
379
    :class:`.DeformationField`.
380

381
    :arg field: A FNIRT-style :class:`.CoefficientField` or
382
                :class:`.DeformationField`
383

384
    :arg from_: Desired reference image coordinate system
385

386
    :arg to:    Desired source image coordinate system
387

388
    :return:    A :class:`.DeformationField` which contains displacements
389
390
391
392
                from the reference image ``from_`` cordinate system to the
                source image ``to`` coordinate syste.
    """
    if isinstance(field, nonlinear.CoefficientField):
393
394
        field = nonlinear.coefficientFieldToDeformationField(field)
    return nonlinear.convertDeformationSpace(field, from_=from_, to=to)