nonlinear.py 31.2 KB
Newer Older
1
2
#!/usr/bin/env python
#
3
# nonlinear.py - Functions/classes for non-linear transformations.
4
5
6
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
7
"""This module contains data structures and functions for working with
8
9
10
FNIRT-style nonlinear transformations.


11
The :class:`DeformationField` and :class:`CoefficientField` can be used to
12
13
load and interact with FNIRT-style transformation images. The following
utility functions are also available:
14

15

16
17
18
.. autosummary::
   :nosignatures:

19
20
21
   detectDeformationType
   convertDeformationType
   convertDeformationSpace
22
   applyDeformation
23
   coefficientFieldToDeformationField
24
"""
25
26


27
28
import                                logging
import itertools                   as it
29

30
31
import numpy                       as np
import scipy.ndimage.interpolation as ndinterp
32

33
34
35
36
import fsl.utils.memoize           as memoize
import fsl.data.image              as fslimage
import fsl.utils.image.resample    as resample
from . import                         affine
37

38

39
40
41
log = logging.getLogger(__name__)


42
class NonLinearTransform(fslimage.Image):
43
    """Class which represents a nonlinear transformation. This is just a base
44
    class for the :class:`DeformationField` and :class:`CoefficientField`
45
46
47
48
    classes.


    A nonlinear transformation is an :class:`.Image` which contains
49
    some mapping between a source image coordinate system and a reference image
50
    coordinate system.
51

52

53
54
55
56
    In FSL, non-linear transformations are defined in terms of the reference
    image coordinate system.  At a given location in the reference image
    space, the non-linear mapping at that location can be used to calculate
    the corresponding location in the source image space. Therefore, these
57
    non-linear transformation effectively encode a transformation *from* the
58
    reference image *to* the (unwarped) source image.
59
60
    """

61
62
63
64

    def __init__(self,
                 image,
                 src,
65
                 ref,
66
67
68
                 srcSpace=None,
                 refSpace=None,
                 **kwargs):
69
70
71
        """Create a ``NonLinearTransform``. See the :meth:`.Nifti.getAffine`
        method for an overview of the values that ``srcSpace`` and ``refSpace``
        may take.
72

73
74
75
        :arg image:       A string containing the name of an image file to
                          load, or a :mod:`numpy` array, or a :mod:`nibabel`
                          image object.
76

77
        :arg src:         :class:`.Nifti` representing the source image.
78

79
        :arg ref:         :class:`.Nifti` representing the reference image.
80

81
82
83
        :arg srcSpace:    Coordinate system in the source image that this
                          ``NonLinearTransform`` maps to. Defaults to
                          ``'fsl'``.
84

85
86
87
88
        :arg refSpace:    Coordinate system in the reference image that this
                          ``NonLinearTransform`` maps from. Defaults to
                          ``'fsl'``.

89
        All other arguments are passed through to :meth:`.Image.__init__`.
90
91
        """

92
93
94
95
96
97
        # 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

98
99
        if srcSpace is None: srcSpace = 'fsl'
        if refSpace is None: refSpace = 'fsl'
100
101
102
103
104

        if srcSpace not in ('fsl', 'voxel', 'world') or \
           refSpace not in ('fsl', 'voxel', 'world'):
            raise ValueError('Invalid source/reference space: {} -> {}'.format(
                srcSpace, refSpace))
105

106
107
        fslimage.Image.__init__(self, image, **kwargs)

108
109
110
111
        self.__src      = fslimage.Nifti(src.header.copy())
        self.__ref      = fslimage.Nifti(ref.header.copy())
        self.__srcSpace = srcSpace
        self.__refSpace = refSpace
112
113
114


    @property
115
    def src(self):
116
117
118
        """Return a reference to the :class:`.Nifti` instance representing
        the source image.
        """
119
        return self.__src
120
121
122


    @property
123
    def ref(self):
124
125
126
        """Return a reference to the :class:`.Nifti` instance representing
        the reference image.
        """
127
        return self.__ref
128
129
130


    @property
131
    def srcSpace(self):
132
133
134
        """Return the source image coordinate system this
        ``NonLinearTransform`` maps from - see :meth:`.Nifti.getAffine`.
        """
135
        return self.__srcSpace
136
137
138


    @property
139
    def refSpace(self):
140
141
142
        """Return the reference image coordinate system this
        ``NonLinearTransform`` maps to - see :meth:`.Nifti.getAffine`.
        """
143
144
145
        return self.__refSpace


146
    def transform(self, coords, from_=None, to=None):
147
148
149
        """Transform coordinates from the reference image space to the source
        image space. Implemented by sub-classes.

150
151
152
        See the :meth:`.Nifti.getAffine` method for an overview of the values
        that ``from_`` and ``to`` may take.

153
154
155
        :arg coords: A sequence of XYZ coordinates, or ``numpy`` array of shape
                     ``(n, 3)`` containing ``n`` sets of coordinates in the
                     reference space.
156

157
        :arg from_:  Reference image space that ``coords`` are defined in
158

159
        :arg to:     Source image space to transform ``coords`` into
160

161
        :returns:    The corresponding coordinates in the source image space.
162
163
164
165
        """
        raise NotImplementedError()


166
167
168
169
class DeformationField(NonLinearTransform):
    """Class which represents a deformation (a.k.a. warp) field which, at each
    voxel, contains either:

170
171
      - a relative displacement from the reference image space to the source
        image space, or
172
      - absolute coordinates in the source space
173
174
175
176


    It is assumed that the a ``DeformationField`` is aligned with the
    reference image in their world coordinate systems (i.e. their ``sform``
Paul McCarthy's avatar
Paul McCarthy committed
177
178
    affines project the reference image and the deformation field into
    alignment).
179
180
    """

181

182
    def __init__(self, image, src, ref=None, **kwargs):
183
184
        """Create a ``DisplacementField``.

185
186
        :arg ref:     Optional. If not provided, it is assumed that the
                      reference is defined in the same space as ``image``.
187

188
189
190
191
        :arg defType: Either ``'absolute'`` or ``'relative'``, indicating
                      the type of this displacement field. If not provided,
                      will be inferred via the :func:`detectDeformationType`
                      function.
192
193
194

        All other arguments are passed through to
        :meth:`NonLinearTransform.__init__`.
195
196
        """

197
198
199
        if ref is None:
            ref = self

200
        defType = kwargs.pop('defType', None)
201

202
203
        if defType not in (None, 'relative', 'absolute'):
            raise ValueError('Invalid value for defType: {}'.format(defType))
204

205
206
        NonLinearTransform.__init__(self, image, src, ref, **kwargs)

207
        self.__defType = defType
208
209
210


    @property
211
212
    def deformationType(self):
        """The type of this ``DeformationField`` - ``'absolute'`` or
213
214
        ``'relative'``.
        """
215
216
217
        if self.__defType is None:
            self.__defType = detectDeformationType(self)
        return self.__defType
218
219
220
221


    @property
    def absolute(self):
222
223
        """``True`` if this ``DeformationField`` contains absolute
        coordinates.
224
        """
225
        return self.deformationType == 'absolute'
226
227
228
229


    @property
    def relative(self):
230
        """``True`` if this ``DeformationField`` contains relative
231
232
        displacements.
        """
233
        return self.deformationType == 'relative'
234
235


236
    def transform(self, coords, from_=None, to=None):
237
238
        """Transform the given XYZ coordinates from the reference image space
        to the source image space.
239
240
241
242

        :arg coords: A sequence of XYZ coordinates, or ``numpy`` array of shape
                     ``(n, 3)`` containing ``n`` sets of coordinates in the
                     reference space.
243

244
        :arg from_:  Reference image space that ``coords`` are defined in
245

246
        :arg to:     Source image space to transform ``coords`` into
247

Paul McCarthy's avatar
Paul McCarthy committed
248
        :returns:    ``coords``, transformed into the source image space
249
250
251
252
253
        """

        if from_ is None: from_ = self.refSpace
        if to    is None: to    = self.srcSpace

254
255
256
        coords   = np.asanyarray(coords)
        outshape = coords.shape
        coords   = coords.reshape((-1, 3))
257
258
259
260
261
262

        # We may need to pre-transform the
        # coordinates so they are in the
        # same reference image space as the
        # displacements
        if from_ != self.refSpace:
263
264
            xform  = self.ref.getAffine(from_, self.refSpace)
            coords = affine.transform(coords, xform)
265

266
267
268
269
270
271
        # 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
272
273
        xform = affine.concat(self    .getAffine('world',       'voxel'),
                              self.ref.getAffine(self.refSpace, 'world'))
274
275

        if np.all(np.isclose(xform, np.eye(4))):
276
            voxels = coords
277
278
        else:
            voxels = affine.transform(coords, xform)
279
280

        # Mask out the coordinates
281
282
        # that are out of bounds of
        # the deformation field
283
        voxels  = np.round(voxels).astype(np.int32)
284
285
286
287
288
289
290
        voxmask = (voxels >= [0, 0, 0]) & (voxels < self.shape[:3])
        voxmask = voxmask.all(axis=1)
        voxels  = voxels[voxmask]

        xs, ys, zs = voxels.T

        if self.absolute: disps = self.data[xs, ys, zs, :]
291
        else:             disps = self.data[xs, ys, zs, :] + coords[voxmask]
292

293
294
        # Make sure the coordinates are in
        # the requested source image space
295
        if to != self.srcSpace:
296
            xform = self.src.getAffine(self.srcSpace, to)
297
            disps = affine.transform(disps, xform)
298

299
300
        # Nans for input coordinates which
        # were outside of the field
301
302
303
        outcoords          = np.full(coords.shape, np.nan)
        outcoords[voxmask] = disps

304
        return outcoords.reshape(outshape)
305

306

307
class CoefficientField(NonLinearTransform):
308
    """Class which represents a B-spline coefficient field generated by FNIRT.
309
310
311

    The :meth:`displacements` method can be used to calculate relative
    displacements for a set of reference space voxel coordinates.
312
313
314
315
316
317
318
319
320
321
322


    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.
323
324
    """

325

326
327
328
329
    def __init__(self,
                 image,
                 src,
                 ref,
330
331
332
333
334
                 srcSpace=None,
                 refSpace=None,
                 fieldType='cubic',
                 knotSpacing=None,
                 fieldToRefMat=None,
335
                 srcToRefMat=None,
336
                 **kwargs):
337
338
        """Create a ``CoefficientField``.

339
        :arg fieldType:     Must currently be ``'cubic'``
340
341
342
343
344
345
346
347

        :arg knotSpacing:   A tuple containing the spline knot spacings along
                            each axis.

        :arg fieldToRefMat: Affine transformation which can transform reference
                            image voxel coordinates into coefficient field
                            voxel coordinates.

348
349
350
351
352
353
354
355
        :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).

356
357
        See the :class:`NonLinearTransform` class for details on the other
        arguments.
358
359
        """

360
        if fieldType not in ('cubic',):
361
362
            raise ValueError('Unsupported field type: {}'.format(fieldType))

363
364
        if srcToRefMat   is not None: srcToRefMat   = np.copy(srcToRefMat)
        if fieldToRefMat is     None: fieldToRefMat = np.eye(4)
365

366
367
368
369
370
371
372
373
374
375
        NonLinearTransform.__init__(self,
                                    image,
                                    src,
                                    ref,
                                    srcSpace,
                                    refSpace,
                                    **kwargs)

        self.__fieldType     = fieldType
        self.__knotSpacing   = tuple(knotSpacing)
376
377
        self.__refToSrcMat   = None
        self.__srcToRefMat   = srcToRefMat
378
        self.__fieldToRefMat = np.copy(fieldToRefMat)
379
380
        self.__refToFieldMat = affine.invert(self.__fieldToRefMat)

381
382
383
        if srcToRefMat is not None:
            self.__refToSrcMat = affine.invert(srcToRefMat)

384
385
386

    @property
    def fieldType(self):
387
388
        """Return the type of the coefficient field, currently always
        ``'cubic'``.
389
390
391
392
393
394
395
396
397
398
399
        """
        return self.__fieldType


    @property
    def knotSpacing(self):
        """Return a tuple containing spline knot spacings along the x, y, and
        z axes.
        """
        return self.__knotSpacing

400

401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
    @property
    def fieldToRefMat(self):
        """Return an affine transformation which can transform coefficient
        field voxel coordinates into reference image voxel coordinates.
        """
        return np.copy(self.__fieldToRefMat)


    @property
    def refToFieldMat(self):
        """Return an affine transformation which can transform reference
        image voxel coordinates into coefficient field voxel coordinates.
        """
        return np.copy(self.__refToFieldMat)


417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
    @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


433
    @memoize.Instanceify(memoize.memoize)
434
435
    def asDeformationField(self, defType='relative', premat=True):
        """Convert this ``CoefficientField`` to a :class:`DeformationField`.
436
437
438

        This method is a wrapper around
        :func:`coefficientFieldToDeformationField`
439
        """
440
        return coefficientFieldToDeformationField(self, defType, premat)
441
442


443
444
445
446
447
448
449
450
    def transform(self, coords, from_=None, to=None, premat=True):
        """Overrides :meth:`NonLinearTransform.transform`. Transforms the
        given ``coords`` from the reference image space into the source image
        space.

        :arg coords: A sequence of XYZ coordinates, or ``numpy`` array of shape
                     ``(n, 3)`` containing ``n`` sets of coordinates in the
                     reference space.
451

452
        :arg from_:  Reference image space that ``coords`` are defined in
453

454
        :arg to:     Source image space to transform ``coords`` into
455

456
        :arg premat: If ``True``, the inverse :meth:`srcToRefMat` is applied
457
                     to the coordinates after the displacements have been
Paul McCarthy's avatar
Paul McCarthy committed
458
                     added.
Paul McCarthy's avatar
Paul McCarthy committed
459
460

        :returns:    ``coords``, transformed into the source image space
461
        """
462
        df = self.asDeformationField(premat=premat)
463
        return df.transform(coords, from_, to)
464
465
466


    def displacements(self, coords):
Paul McCarthy's avatar
Paul McCarthy committed
467
        """Calculate the relative displacements for the given coordinates.
468
469
470
471
472
473
474
475
476
477

        :arg coords: ``(N, 3)`` array of reference image voxel coordinates.
        :return:      A ``(N, 3)`` array  of relative displacements to the
                      source image for ``coords``
        """

        if self.fieldType != 'cubic':
            raise NotImplementedError()

        # See
478
479
        #   https://www.cs.jhu.edu/~cis/cista/746/papers/\
        #     RueckertFreeFormBreastMRI.pdf
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
        #   https://www.fmrib.ox.ac.uk/datasets/techrep/tr07ja2/tr07ja2.pdf

        # Cubic b-spline basis functions
        def b0(u):
            return ((1 - u) ** 3) / 6

        def b1(u):
            return (3 * (u ** 3) - 6 * (u ** 2) + 4) / 6

        def b2(u):
            return (-3 * (u ** 3) + 3 * (u ** 2)  + 3 * u + 1) / 6

        def b3(u):
            return (u ** 3) / 6

        b = [b0, b1, b2, b3]

        fdata      = self.data
        nx, ny, nz = self.shape[:3]

        # Convert the given voxel coordinates
        # into the corresponding coefficient
        # field voxel coordinates
Paul McCarthy's avatar
Paul McCarthy committed
503
        i, j, k = affine.transform(coords, self.refToFieldMat).T
504
505
506
507
508
509
510

        # i, j, k: coefficient field indices
        # u, v, w: position of the ref voxel
        #          on the current spline
        u = np.remainder(i, 1)
        v = np.remainder(j, 1)
        w = np.remainder(k, 1)
511
512
513
        i = np.floor(i).astype(np.int32)
        j = np.floor(j).astype(np.int32)
        k = np.floor(k).astype(np.int32)
514
515
516
517
518
519
520

        disps = np.zeros(coords.shape)

        for l, m, n in it.product(range(4), range(4), range(4)):
            il   = i + l
            jm   = j + m
            kn   = k + n
521
522
523
524
525
526
            mask = (il >= 0)  & \
                   (il <  nx) & \
                   (jm >= 0)  & \
                   (jm <  ny) & \
                   (kn >= 0)  & \
                   (kn <  nz)
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544

            il = il[mask]
            jm = jm[mask]
            kn = kn[mask]
            uu = u[ mask]
            vv = v[ mask]
            ww = w[ mask]

            cx, cy, cz = fdata[il, jm, kn, :].T
            c          = b[l](uu) * b[m](vv) * b[n](ww)

            disps[mask, 0] += c * cx
            disps[mask, 1] += c * cy
            disps[mask, 2] += c * cz

        return disps


545
546
def detectDeformationType(field):
    """Attempt to automatically determine whether a deformation field is
547
548
    specified in absolute or relative coordinates.

549
    :arg field: A :class:`DeformationField`
550
551

    :returns:   ``'absolute'`` if it looks like ``field`` contains absolute
552
                coordinates, ``'relative'`` otherwise.
553
554
555
    """

    # This test is based on the assumption
556
    # that a deformation field containing
557
558
559
    # absolute coordinates will have a
    # greater standard deviation than one
    # which contains relative coordinates.
560
    absdata = field[:]
561
    reldata = convertDeformationType(field, 'relative')
562
563
564
565
566
567
568
    stdabs  = absdata.std(axis=(0, 1, 2)).sum()
    stdrel  = reldata.std(axis=(0, 1, 2)).sum()

    if stdabs > stdrel: return 'absolute'
    else:               return 'relative'


569
570
571
def convertDeformationType(field, defType=None):
    """Convert a deformation field between storing absolute coordinates or
    relative displacements.
572

573
574
575
576
    :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.
577
578
    """

579
580
581
    if defType is None:
        if field.deformationType == 'absolute': defType = 'relative'
        else:                                   defType = 'absolute'
582

Paul McCarthy's avatar
Paul McCarthy committed
583
584
585
586
    if defType not in ('absolute', 'relative'):
        raise ValueError('defType must be "absolute" or "relative" '
                         '("{}" passed)'.format(defType))

587
588
    # Regardless of the conversion direction,
    # we need the coordinates of every voxel
589
    # in the reference coordinate system.
590
    dx, dy, dz = field.shape[:3]
591
592
    xform      = affine.concat(field.ref.getAffine('world', field.refSpace),
                               field    .getAffine('voxel', 'world'))
593

594
595
596
597
    coords     = np.meshgrid(np.arange(dx),
                             np.arange(dy),
                             np.arange(dz), indexing='ij')
    coords     = np.array(coords).transpose((1, 2, 3, 0))
598
    coords     = affine.transform(coords.reshape((-1, 3)), xform)
599
600
601
    coords     = coords.reshape((dx, dy, dz, 3))

    # If converting from relative to absolute,
602
603
604
    # we just add the coordinates to (what is
    # assumed to be) the relative shift. Or,
    # to convert from absolute to relative,
605
    # we subtract the reference image voxels.
Paul McCarthy's avatar
Paul McCarthy committed
606
607
    if defType == 'absolute': return field.data + coords
    else:                     return field.data - coords
608
609


610
611
def convertDeformationSpace(field, from_, to):
    """Adjust the source and/or reference spaces of the given deformation
612
613
614
    field. See the :meth:`.Nifti.getAffine` method for the valid values for
    the ``from_`` and ``to`` arguments.

615
    :arg field: A :class:`DeformationField` instance
616
617
    :arg from_: New reference image coordinate system
    :arg to:    New source image coordinate system
618

619
    :returns:   A new :class:`DeformationField` which transforms between
620
                the reference ``from_`` coordinate system and the source ``to``
621
622
                coordinate system.
    """
623

624
625
626
    if field.srcSpace == to and field.refSpace == from_:
        return field

627
628
629
    # Get the field in absolute coordinates
    # if necessary - these are our source
    # coordinates in the original "to" space.
630
    fieldcoords = field.data
631
    if field.relative: srccoords = convertDeformationType(field)
632
633
634
    else:              srccoords = fieldcoords

    srccoords = srccoords.reshape((-1, 3))
635
636
637
638
639
640
641
642

    # Now transform those source coordinates
    # from the original source space to the
    # source space specified by "to"
    if to != field.srcSpace:

        srcmat    = field.src.getAffine(field.srcSpace, to)
        srccoords = affine.transform(srccoords, srcmat)
643
644

    # If we have been asked to return
645
    # absolute coordinates, the
646
647
    # reference "from_" coordinate
    # system is irrelevant - we're done.
648
    if field.absolute:
649
650
        fieldcoords = srccoords

651
    # Otherwise our deformation field
652
    # will contain relative displacements
653
    # between the reference image "from_"
654
    # coordinate system and the source
655
656
657
658
659
    # image "to" coordinate system. We
    # need to re-calculate the relative
    # displacements between the new
    # reference "from_" space and source
    # "to" space.
660
    else:
661
662
663
664
665
666
        refcoords = np.meshgrid(np.arange(field.shape[0]),
                                np.arange(field.shape[1]),
                                np.arange(field.shape[2]), indexing='ij')
        refcoords = np.array(refcoords)
        refcoords = refcoords.transpose((1, 2, 3, 0)).reshape((-1, 3))

667
668
669
670
671
672
        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)
673

674
675
        fieldcoords = srccoords - refcoords

676
    return DeformationField(
677
678
        fieldcoords.reshape(field.shape),
        header=field.header,
679
680
        src=field.src,
        ref=field.ref,
681
682
        srcSpace=to,
        refSpace=from_,
683
        defType=field.deformationType)
684
685


686
687
688
689
690
691
692
def applyDeformation(image,
                     field,
                     ref=None,
                     order=1,
                     mode=None,
                     cval=None,
                     premat=None):
693
694
695
696
697
698
699
700
    """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
Paul McCarthy's avatar
Paul McCarthy committed
701
    the input image. It is therefore assumed that an alternate ``ref`` is
702
703
    aligned in world coordinates with the field's actual reference image.

704
    :arg image:  :class:`.Image` to be transformed
705

706
    :arg field:  :class:`DeformationField` to use
707

708
709
    :arg ref:    Alternate reference image - if not provided, ``field.ref``
                 is used
710

711
712
713
714
715
    :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.
716

717
718
    :arg mode:   How to handle regions which are outside of the image FOV.
                 Defaults to `''nearest'``.
719

720
    :arg cval:   Constant value to use when ``mode='constant'``.
721

722
723
724
725
726
727
    :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:     ``numpy.array`` containing the transformed image data.
728
729
730
731
732
    """

    if order is None: order = 1
    if mode  is None: mode  = 'nearest'
    if cval  is None: cval  = 0
733
    if ref   is None: ref   = field.ref
734
735
736
737
738
739
740
741
742
743
744
745
746
747

    # 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')

748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
    # 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.
763
764
    src = field.src

765
    if not field.sameSpace(ref):
766
767
        field = resample.resampleToReference(field,
                                             ref,
768
                                             order=order,
769
770
                                             mode='constant',
                                             cval=-1)[0]
771
772
773
    else:
        field = field.data

774
775
776
777
778
779
780
    # 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
781
782
783
784
785
786
    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)
787
788
        shape = field.shape
        field = field.reshape((-1, 3))
789
        field = affine.transform(field, premat)
790
791
        field = field.reshape(shape)

792
793
794
795
796
797
798
799
    field = field.transpose((3, 0, 1, 2))
    return ndinterp.map_coordinates(image.data,
                                    field,
                                    order=order,
                                    mode=mode,
                                    cval=cval)


800
def coefficientFieldToDeformationField(field, defType='relative', premat=True):
801
    """Convert a :class:`CoefficientField` into a :class:`DeformationField`.
802

803
    :arg field:   :class:`CoefficientField` to convert
804

805
806
    :arg defType: The type of deformation field - either ``'relative'`` (the
                  default) or ``'absolute'``.
807

808
809
    :arg premat:  If ``True`` (the default), the :meth:`srcToRefMat` is
                  encoded into the deformation field.
810

811
    :return:      :class:`DeformationField` calculated from ``field``.
812
    """
813

814
815
    # Generate coordinates for every
    # voxel in the reference image
816
817
818
819
820
821
822
823
    ix, iy, iz = field.ref.shape[:3]
    x,  y,  z  = np.meshgrid(np.arange(ix),
                             np.arange(iy),
                             np.arange(iz), indexing='ij')
    x          = x.flatten()
    y          = y.flatten()
    z          = z.flatten()
    xyz        = np.vstack((x, y, z)).T
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843

    # There are three spaces to consider here:
    #
    #  - ref space:         Reference image scaled voxels ("fsl" space)
    #
    #  - aligned-src space: Source image scaled voxels, after the
    #                       source image has been linearly aligned to
    #                       the reference via the field.srcToRefMat
    #                       This will typically be equivalent to ref
    #                       space
    #
    #  - orig-src space:    Source image scaled voxels, in the coordinate
    #                       system of the original source image, without
    #                       linear alignment to the reference image

    # The displacements method will
    # return relative displacements
    # from ref space to aligned-src
    # space.
    disps   = field.displacements(xyz).reshape((ix, iy, iz, 3))
844
    rdfield = DeformationField(disps,
845
                               header=field.ref.header,
846
847
848
849
850
851
852
                               src=field.src,
                               ref=field.ref,
                               srcSpace=field.srcSpace,
                               refSpace=field.refSpace,
                               defType='relative')

    if (defType == 'relative') and (not premat):
853
854
855
        return rdfield

    # Convert to absolute - the
856
    # deformations will now be
857
858
    # absolute coordinates in
    # aligned-src space
859
    disps = convertDeformationType(rdfield)
860
861
862
863

    # Apply the premat if requested -
    # this will transform the coordinates
    # from aligned-src to orig-src space.
864
    if premat and field.srcToRefMat is not None:
865
866
867
868
869
870

        # We apply the premat in the same way
        # that fnirtfileutils does - applying
        # the inverse affine to every ref space
        # voxel coordinate, then adding it to
        # the existing displacements.
871
872
873
874
875
876
        shape  = disps.shape
        disps  = disps.reshape(-1, 3)
        premat = affine.concat(field.refToSrcMat - np.eye(4),
                               field.ref.getAffine('voxel', 'fsl'))
        disps  = disps + affine.transform(xyz, premat)
        disps  = disps.reshape(shape)
877
878
879
880
881
882
883
884

        # note that convertwarp applies a premat
        # differently - its method is equivalent
        # to directly transforming the existing
        # absolute displacements, i.e.:
        #
        #   disps = affine.transform(disps, refToSrc)

885
    adfield = DeformationField(disps,
886
                               header=field.ref.header,
887
888
889
890
891
                               src=field.src,
                               ref=field.ref,
                               srcSpace=field.srcSpace,
                               refSpace=field.refSpace,
                               defType='absolute')
892
893
894

    # Not either return absolute displacements,
    # or convert back to relative displacements
895
    if defType == 'absolute':
896
897
        return adfield
    else:
898
899
900
901
902
903
904
        return DeformationField(convertDeformationType(adfield),
                                src=field.src,
                                ref=field.ref,
                                srcSpace=field.srcSpace,
                                refSpace=field.refSpace,
                                header=field.ref.header,
                                defType='relative')