image.py 59.4 KB
Newer Older
1
#!/usr/bin/env python
2
#
3
4
# image.py - Provides the :class:`Image` class, for representing 3D/4D NIFTI
#            images.
5
6
7
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
8
"""This module provides the :class:`Nifti` and :class:`Image` classes, for
9
representing NIFTI1 and NIFTI2 images. The ``nibabel`` package is used
Paul McCarthy's avatar
Paul McCarthy committed
10
for file I/O.
11
12
13
14
15
16
17
18


It is very easy to load a NIFTI image::

    from fsl.data.image import Image
    myimg = Image('MNI152_T1_2mm.nii.gz')


19
20
A handful of other functions are also provided for working with image files
and file names:
21
22
23
24

.. autosummary::
   :nosignatures:

25
   canonicalShape
26
   looksLikeImage
27
   addExt
28
   splitExt
29
30
   getExt
   removeExt
31
   defaultExt
32
"""
33

34

35
36
import                      os
import os.path           as op
37
import itertools         as it
38
import                      json
39
import                      string
40
import                      logging
41
import                      tempfile
42
import                      warnings
43

44
import                      six
45
46
47
48
import numpy             as np

import nibabel           as nib
import nibabel.fileslice as fileslice
49

50
import fsl.utils.meta        as meta
51
import fsl.transform.affine  as affine
52
import fsl.utils.notifier    as notifier
53
import fsl.utils.memoize     as memoize
54
import fsl.utils.path        as fslpath
55
import fsl.utils.deprecated  as deprecated
56
import fsl.utils.bids        as fslbids
57
58
import fsl.data.constants    as constants
import fsl.data.imagewrapper as imagewrapper
59

60

61
62
log = logging.getLogger(__name__)

63

Paul McCarthy's avatar
Paul McCarthy committed
64
65
ALLOWED_EXTENSIONS = ['.nii.gz', '.nii', '.img', '.hdr', '.img.gz', '.hdr.gz']
"""The file extensions which we understand. This list is used as the default
Paul McCarthy's avatar
Paul McCarthy committed
66
67
if the ``allowedExts`` parameter is not passed to any of the ``*Ext``
functions, or the :func:`looksLikeImage` function.
Paul McCarthy's avatar
Paul McCarthy committed
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
"""


EXTENSION_DESCRIPTIONS = ['Compressed NIFTI images',
                          'NIFTI images',
                          'ANALYZE75 images',
                          'NIFTI/ANALYZE75 headers',
                          'Compressed NIFTI/ANALYZE75 images',
                          'Compressed NIFTI/ANALYZE75 headers']
"""Descriptions for each of the extensions in :data:`ALLOWED_EXTENSIONS`. """


FILE_GROUPS = [('.hdr',    '.img'),
               ('.hdr.gz', '.img.gz')]
"""File suffix groups used by :func:`addExt` to resolve file path
ambiguities - see :func:`fsl.utils.path.addExt`.
"""


PathError = fslpath.PathError
"""Error raised by :mod:`fsl.utils.path` functions when an error occurs.
Made available in this module for convenience.
"""


93
class Nifti(notifier.Notifier, meta.Meta):
94
95
96
97
    """The ``Nifti`` class is intended to be used as a base class for
    things which either are, or are associated with, a NIFTI image.
    The ``Nifti`` class is intended to represent information stored in
    the header of a NIFTI file - if you want to load the data from
98
    a file, use the :class:`Image` class instead.
Paul McCarthy's avatar
Paul McCarthy committed
99

100

101
    When a ``Nifti`` instance is created, it adds the following attributes
102
    to itself:
103

104

105
    ================= ====================================================
106
107
    ``header``        The :mod:`nibabel` NIFTI1/NIFTI2/Analyze header
                      object.
108

109
    ``shape``         A list/tuple containing the number of voxels along
110
                      each image dimension.
111
112

    ``pixdim``        A list/tuple containing the length of one voxel
113
                      along each image dimension.
114

115
116
117
    ``voxToWorldMat`` A 4*4 array specifying the affine transformation
                      for transforming voxel coordinates into real world
                      coordinates.
118

119
120
121
    ``worldToVoxMat`` A 4*4 array specifying the affine transformation
                      for transforming real world coordinates into voxel
                      coordinates.
122

123
124
125
    ``intent``        The NIFTI intent code specified in the header (or
                      :attr:`.constants.NIFTI_INTENT_NONE` for Analyze
                      images).
126
    ================= ====================================================
127

128

129
130
131
132
133
    The ``header`` field may either be a ``nifti1``, ``nifti2``, or
    ``analyze`` header object. Make sure to take this into account if you are
    writing code that should work with all three. Use the :meth:`niftiVersion`
    property if you need to know what type of image you are dealing with.

134

135
136
137
138
139
    The ``shape`` attribute may not precisely match the image shape as
    reported in the NIFTI header, because trailing dimensions of size 1 are
    squeezed out. See the :meth:`__determineShape` and :meth:`mapIndices`
    methods.

140

141
    **Affine transformations**
142

143

144
145
146
147
148
149
150
151
    The ``Nifti`` class is aware of three coordinate systems:

      - The ``voxel`` coordinate system, used to access image data

      - The ``world`` coordinate system, where voxel coordinates are
        transformed into a millimetre coordinate system, defined by the
        ``sform`` and/or ``qform`` elements of the NIFTI header.

Paul McCarthy's avatar
Paul McCarthy committed
152
      - The ``fsl`` coordinate system, where voxel coordinates are scaled by
153
154
155
156
157
158
159
160
161
162
163
164
165
166
        the ``pixdim`` values in the NIFTI header, and the X axis is inverted
        if the voxel-to-world affine has a positive determinant.


    The :meth:`getAffine` method is a simple means of acquiring an affine
    which will transform between any of these coordinate systems.


    See `here <http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT/FAQ#What_is_the_format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to_the_transformation_parameters.3F>`_
    for more details on the ``fsl`` coordinate system..


    The ``Nifti`` class follows the same process as ``nibabel`` in determining
    the ``voxel`` to ``world`` affine (see
167
168
    http://nipy.org/nibabel/nifti_images.html#the-nifti-affines):

169

170
171
172
173
     1. If ``sform_code != 0`` ("unknown") use the sform affine; else
     2. If ``qform_code != 0`` ("unknown") use the qform affine; else
     3. Use the fall-back affine.

174

175
176
177
178
179
180
181
182
183
184
185
    However, the *fall-back* affine used by the ``Nifti`` class differs to
    that used by ``nibabel``. In ``nibabel``, the origin (world coordinates
    (0, 0, 0)) is set to the centre of the image. Here in the ``Nifti``
    class, we set the world coordinate orign to be the corner of the image,
    i.e. the corner of voxel (0, 0, 0).


    You may change the ``voxToWorldMat`` of a ``Nifti`` instance (unless it
    is an Analyze image). When you do so:

     - Only the ``sform`` of the underlying ``Nifti1Header`` object is changed
186

187
     - The ``qform`` is not modified.
188

189
190
191
     - If the ``sform_code`` was previously set to ``NIFTI_XFORM_UNKNOWN``,
       it is changed to ``NIFTI_XFORM_ALIGNED_ANAT``. Otherwise, the
       ``sform_code`` is not modified.
192

193
194
195
196

    **ANALYZE support**


197
198
    A ``Nifti`` instance expects to be passed either a
    ``nibabel.nifti1.Nifti1Header`` or a ``nibabel.nifti2.Nifti2Header``, but
199
    can also encapsulate a ``nibabel.analyze.AnalyzeHeader``. In this case:
200
201
202
203
204

      - The image voxel orientation is assumed to be R->L, P->A, I->S.

      - The affine will be set to a diagonal matrix with the header pixdims as
        its elements (with the X pixdim negated), and an offset specified by
205
        the ANALYZE ``origin`` fields. Construction of the affine is handled
206
207
208
209
210
211
212
213
        by ``nibabel``.

      - The :meth:`niftiVersion` method will return ``0``.

      - The :meth:`getXFormCode` method will return
        :attr:`.constants.NIFTI_XFORM_ANALYZE`.


214
215
216
    **Metadata**


217
218
    The ``Image`` class inherits from the :class:`.Meta` class - its methods
    can be used to store and query any meta-data associated with the image.
219
220


221
222
223
    **Notification**


224
    The ``Nifti`` class implements the :class:`.Notifier` interface -
225
226
227
228
    listeners may register to be notified on the following topics:

    =============== ========================================================
    ``'transform'`` The affine transformation matrix has changed. This topic
229
                    will occur when the :meth:`voxToWorldMat` is changed.
Paul McCarthy's avatar
Paul McCarthy committed
230
    ``'header'``    A header field has changed. This will occur when the
231
                    :meth:`intent` is changed.
232
    =============== ========================================================
233
    """  # noqa
234

235

236
    def __init__(self, header):
237
        """Create a ``Nifti`` object.
Paul McCarthy's avatar
Paul McCarthy committed
238

239
        :arg header: A :class:`nibabel.nifti1.Nifti1Header`,
240
241
                       :class:`nibabel.nifti2.Nifti2Header`, or
                       ``nibabel.analyze.AnalyzeHeader`` to be used as the
242
                       image header.
243
        """
244

Paul McCarthy's avatar
Paul McCarthy committed
245
        # Nifti2Header is a sub-class of Nifti1Header,
246
247
248
        # and Nifti1Header a sub-class of AnalyzeHeader,
        # so we only need to test for the latter.
        if not isinstance(header, nib.analyze.AnalyzeHeader):
Paul McCarthy's avatar
Paul McCarthy committed
249
250
            raise ValueError('Unrecognised header: {}'.format(header))

251
        origShape, shape, pixdim = Nifti.determineShape(header)
Paul McCarthy's avatar
Paul McCarthy committed
252
253
254
255
        voxToWorldMat            = Nifti.determineAffine(header)
        affines, isneuro         = Nifti.generateAffines(voxToWorldMat,
                                                         shape,
                                                         pixdim)
256

257
        self.__header         = header
Paul McCarthy's avatar
Paul McCarthy committed
258
259
260
        self.__shape          = shape
        self.__origShape      = origShape
        self.__pixdim         = pixdim
261
262
        self.__affines        = affines
        self.__isNeurological = isneuro
Paul McCarthy's avatar
Paul McCarthy committed
263

264

265
266
    @staticmethod
    def determineShape(header):
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
        """This method is called by :meth:`__init__`. It figures out the actual
        shape of the image data, and the zooms/pixdims for each data axis. Any
        empty trailing dimensions are squeezed, but the returned shape is
        guaranteed to be at least 3 dimensions. Returns:

         - A sequence/tuple containing the image shape, as reported in the
           header.
         - A sequence/tuple containing the effective image shape.
         - A sequence/tuple containing the zooms/pixdims.
        """

        # The canonicalShape function figures out
        # the data shape that we should use.
        origShape = list(header.get_data_shape())
        shape     = canonicalShape(origShape)
        pixdims   = list(header.get_zooms())

        # if get_zooms() doesn't return at
        # least len(shape) values, we'll
        # fall back to the pixdim field in
        # the header.
        if len(pixdims) < len(shape):
            pixdims = header['pixdim'][1:]

        pixdims = pixdims[:len(shape)]

293
294
295
296
297
298
        # should never happen, but if we only
        # have zoom values for the original
        # (< 3D) shape, pad them with 1s.
        if len(pixdims) < len(shape):
            pixdims = pixdims + [1] * (len(shape) - len(pixdims))

299
300
301
        return origShape, shape, pixdims


302
303
    @staticmethod
    def determineAffine(header):
304
305
        """Called by :meth:`__init__`. Figures out the voxel-to-world
        coordinate transformation matrix that is associated with this
306
        ``Nifti`` instance.
307
        """
308

309
310
311
        # We have to treat FSL/FNIRT images
        # specially, as FNIRT clobbers the
        # sform section of the NIFTI header
312
        # to store other data.
313
314
315
        intent = header.get('intent_code', -1)
        qform  = header.get('qform_code',  -1)
        sform  = header.get('sform_code',  -1)
316

317
318
319
320
321
322
323
324
        # FNIRT non-linear coefficient files
        # clobber the sform/qform/intent
        # and pixdims of the nifti header,
        # so we can't correctly place it in
        # the world coordinate system. See
        # $FSLDIR/src/fnirt/fnirt_file_writer.cpp
        # and fsl.transform.nonlinear for more
        # details.
325
326
327
328
329
        if intent in (constants.FSL_DCT_COEFFICIENTS,
                      constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
                      constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS,
                      constants.FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS,
                      constants.FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS):
330
331
332

            log.debug('FNIRT coefficient field detected - generating affine')

Paul McCarthy's avatar
Paul McCarthy committed
333
334
335
336
337
338
339
            # Knot spacing is stored in the pixdims
            # (specified in terms of reference image
            # voxels), and reference image pixdims
            # are stored as intent code parameters.
            # If we combine the two, we can at least
            # get the shape/size of the coefficient
            # field about right
340
            knotpix       =  header.get_zooms()[:3]
341
342
343
            refpix        = (header.get('intent_p1', 1) or 1,
                             header.get('intent_p2', 1) or 1,
                             header.get('intent_p3', 1) or 1)
344
345
346
            voxToWorldMat = affine.concat(
                affine.scaleOffsetXform(refpix,  0),
                affine.scaleOffsetXform(knotpix, 0))
347

348
349
350
351
        # If the qform or sform codes are unknown,
        # then we can't assume that the transform
        # matrices are valid. So we fall back to a
        # pixdim scaling matrix.
352
        #
353
354
355
        # n.b. For images like this, nibabel returns
        # a scaling matrix where the centre voxel
        # corresponds to world location (0, 0, 0).
356
        # This goes against the NIFTI spec - it
357
        # should just be a straight scaling matrix.
358
        elif qform == 0 and sform == 0:
359
            pixdims       = header.get_zooms()
360
            voxToWorldMat = affine.scaleOffsetXform(pixdims, 0)
361

362
        # Otherwise we let nibabel decide
363
        # which transform to use.
364
        else:
365
            voxToWorldMat = np.array(header.get_best_affine())
366

367
        return voxToWorldMat
368

369

370
371
    @staticmethod
    def generateAffines(voxToWorldMat, shape, pixdim):
372
373
374
375
376
        """Called by :meth:`__init__`, and the :meth:`voxToWorldMat` setter.
        Generates and returns a dictionary containing affine transformations
        between the ``voxel``, ``fsl``, and ``world`` coordinate
        systems. These affines are accessible via the :meth:`getAffine`
        method.
377

378
        :arg voxToWorldMat: The voxel-to-world affine transformation
379
380
381
        :arg shape:         Image shape (number of voxels along each dimension
        :arg pixdim:        Image pixdims (size of one voxel along each
                            dimension)
382
383
384
385
386
387
388
        :returns:           A tuple containing:

                             - a dictionary of affine transformations between
                               each pair of coordinate systems
                             - ``True`` if the image is to be considered
                               "neurological", ``False`` otherwise - see the
                               :meth:`isNeurological` method.
389
390
        """

391
        import numpy.linalg as npla
392

393
        affines = {}
394
395
        shape             = list(shape[ :3])
        pixdim            = list(pixdim[:3])
396
397
        voxToScaledVoxMat = np.diag(pixdim + [1.0])
        isneuro           = npla.det(voxToWorldMat) > 0
398

399
400
        if isneuro:
            x                 = (shape[0] - 1) * pixdim[0]
401
402
403
            flip              = affine.scaleOffsetXform([-1, 1, 1],
                                                        [ x, 0, 0])
            voxToScaledVoxMat = affine.concat(flip, voxToScaledVoxMat)
404

405
406
407
        affines['fsl',   'fsl']   = np.eye(4)
        affines['voxel', 'voxel'] = np.eye(4)
        affines['world', 'world'] = np.eye(4)
408
        affines['voxel', 'world'] = voxToWorldMat
409
        affines['world', 'voxel'] = affine.invert(voxToWorldMat)
410
        affines['voxel', 'fsl']   = voxToScaledVoxMat
411
412
413
414
415
        affines['fsl',   'voxel'] = affine.invert(voxToScaledVoxMat)
        affines['fsl',   'world'] = affine.concat(affines['voxel', 'world'],
                                                  affines['fsl',   'voxel'])
        affines['world', 'fsl']   = affine.concat(affines['voxel',   'fsl'],
                                                  affines['world', 'voxel'])
416
417

        return affines, isneuro
418
419


420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
    @staticmethod
    def identifyAffine(image, xform, from_=None, to=None):
        """Attempt to identify the source or destination space for the given
        affine.

        ``xform`` is assumed to be an affine transformation which can be used
        to transform coordinates between two coordinate systems associated with
        ``image``.

        If one of ``from_`` or ``to`` is provided, the other will be derived.
        If neither are provided, both will be derived. See the
        :meth:`.Nifti.getAffine` method for details on the valild values that
        ``from_`` and ``to`` may take.

        :arg image: :class:`.Nifti` instance associated with the affine.

        :arg xform: ``(4, 4)`` ``numpy`` array encoding an affine
                    transformation

        :arg from_: Label specifying the coordinate system which ``xform``
                    takes as input

        :arg to:    Label specifying the coordinate system which ``xform``
                    produces as output

        :returns:   A tuple containing:
                      - A label for the ``from_`` coordinate system
                      - A label for the ``to`` coordinate system
        """

        if (from_ is not None) and (to is not None):
            return from_, to

        if from_ is not None: froms = [from_]
        else:                 froms = ['voxel', 'fsl', 'world']
        if to    is not None: tos   = [to]
        else:                 tos   = ['voxel', 'fsl', 'world']

        for from_, to in it.product(froms, tos):

            candidate = image.getAffine(from_, to)

            if np.all(np.isclose(candidate, xform)):
                return from_, to

        raise ValueError('Could not identify affine')


468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
    def strval(self, key):
        """Returns the specified NIFTI header field, converted to a python
        string, correctly null-terminated, and with non-printable characters
        removed.

        This method is used to sanitise some NIFTI header fields. The default
        Python behaviour for converting a sequence of bytes to a string is to
        strip all termination characters (bytes with value of ``0x00``) from
        the end of the sequence.

        This default behaviour does not handle the case where a sequence of
        bytes which did contain a long string is subsequently overwritten with
        a shorter string - the short string will be terminated, but that
        termination character will be followed by the remainder of the
        original string.
        """

        val = self.header[key]

Paul McCarthy's avatar
Paul McCarthy committed
487
488
        try:              val = bytes(val).partition(b'\0')[0]
        except Exception: val = bytes(val)
489
490
491

        val = val.decode('ascii')

492
        return ''.join([c for c in val if c in string.printable]).strip()
493
494


495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
    @property
    def header(self):
        """Return a reference to the ``nibabel`` header object. """
        return self.__header


    @header.setter
    def header(self, header):
        """Replace the ``nibabel`` header object managed by this ``Nifti``
        with a new header. The new header must have the same dimensions,
        voxel size, and orientation as the old one.
        """
        new = Nifti(header)
        if not (self.sameSpace(new) and self.ndim == new.ndim):
            raise ValueError('Incompatible header')
        self.__header = new


513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
    @property
    def niftiVersion(self):
        """Returns the NIFTI file version:

           - ``0`` for ANALYZE
           - ``1`` for NIFTI1
           - ``2`` for NIFTI2
        """

        # nib.Nifti2 is a subclass of Nifti1,
        # and Nifti1 a subclass of Analyze,
        # so we have to check in order
        if   isinstance(self.header, nib.nifti2.Nifti2Header):   return 2
        elif isinstance(self.header, nib.nifti1.Nifti1Header):   return 1
        elif isinstance(self.header, nib.analyze.AnalyzeHeader): return 0

        else: raise RuntimeError('Unrecognised header: {}'.format(self.header))


    @property
    def shape(self):
        """Returns a tuple containing the image data shape. """
        return tuple(self.__shape)

537

538
539
540
541
542
543
544
545
546
547
    @property
    def ndim(self):
        """Returns the number of dimensions in this image. This number may not
        match the number of dimensions specified in the NIFTI header, as
        trailing dimensions of length 1 are ignored. But it is guaranteed to be
        at least 3.
        """
        return len(self.__shape)


548
549
550
551
552
553
554
555
    @property
    def pixdim(self):
        """Returns a tuple containing the image pixdims (voxel sizes)."""
        return tuple(self.__pixdim)


    @property
    def intent(self):
556
557
558
559
560
561
562
563
        """Returns the NIFTI intent code of this image. """
        return self.header.get('intent_code', constants.NIFTI_INTENT_NONE)


    @intent.setter
    def intent(self, val):
        """Sets the NIFTI intent code of this image. """
        # analyze has no intent
564
        if (self.niftiVersion > 0) and (val != self.intent):
565
566
            self.header.set_intent(val, allow_unknown=True)
            self.notify(topic='header')
567
568


569
570
571
572
    @property
    def xyzUnits(self):
        """Returns the NIFTI XYZ dimension unit code. """

573
574
575
576
        # analyze images have no unit field
        if self.niftiVersion == 0:
            return constants.NIFTI_UNITS_MM

577
578
579
580
        # The nibabel get_xyzt_units returns labels,
        # but we want the NIFTI codes. So we use
        # the (undocumented) nifti1.unit_codes field
        # to convert back to the raw codes.
Paul McCarthy's avatar
Paul McCarthy committed
581
582
583
        units = self.header.get_xyzt_units()[0]
        units = nib.nifti1.unit_codes[units]
        return units
584
585


Paul McCarthy's avatar
Paul McCarthy committed
586
    @property
587
588
589
    def timeUnits(self):
        """Returns the NIFTI time dimension unit code. """

590
591
592
593
        # analyze images have no unit field
        if self.niftiVersion == 0:
            return constants.NIFTI_UNITS_SEC

594
        # See xyzUnits
Paul McCarthy's avatar
Paul McCarthy committed
595
596
597
        units = self.header.get_xyzt_units()[1]
        units = nib.nifti1.unit_codes[units]
        return units
598
599


600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
    def getAffine(self, from_, to):
        """Return an affine transformation which can be used to transform
        coordinates from ``from_`` to ``to``.

        Valid values for the ``from_`` and ``to`` arguments are:

         - ``'voxel'``: The voxel coordinate system
         - ``'world'``: The world coordinate system, as defined by the image
           sform/qform
         - ``'fsl'``: The FSL coordinate system (scaled voxels, with a
           left-right flip if the sform/qform has a positive determinant)

        :arg from_: Source coordinate system
        :arg to:    Destination coordinate system
        :returns:   A ``numpy`` array of shape ``(4, 4)``
        """
        from_ = from_.lower()
        to    = to   .lower()

        if from_ not in ('voxel', 'fsl', 'world') or \
           to    not in ('voxel', 'fsl', 'world'):
            raise ValueError('Invalid source/reference spaces: '
                             '{} -> {}'.format(from_, to))

        return np.copy(self.__affines[from_, to])


627
628
629
630
631
    @property
    def worldToVoxMat(self):
        """Returns a ``numpy`` array of shape ``(4, 4)`` containing an
        affine transformation from world coordinates to voxel coordinates.
        """
632
        return self.getAffine('world', 'voxel')
633
634
635
636
637
638
639


    @property
    def voxToWorldMat(self):
        """Returns a ``numpy`` array of shape ``(4, 4)`` containing an
        affine transformation from voxel coordinates to world coordinates.
        """
640
        return self.getAffine('voxel', 'world')
641
642
643
644


    @voxToWorldMat.setter
    def voxToWorldMat(self, xform):
645
646
647
        """Update the ``voxToWorldMat``. The ``worldToVoxMat`` value is also
        updated. This will result in notification on the ``'transform'``
        topic.
648
649
        """

650
651
652
653
654
655
656
        # Can't do much with
        # an analyze image
        if self.niftiVersion == 0:
            raise Exception('voxToWorldMat cannot be '
                            'changed for an ANALYZE image')

        header    = self.header
657
        sformCode = int(header['sform_code'])
658

659
660
661
662
        if sformCode == constants.NIFTI_XFORM_UNKNOWN:
            sformCode = constants.NIFTI_XFORM_ALIGNED_ANAT

        header.set_sform(xform, code=sformCode)
663

664
665
666
        affines, isneuro = Nifti.generateAffines(xform,
                                                 self.shape,
                                                 self.pixdim)
667
668
669

        self.__affines        = affines
        self.__isNeurological = isneuro
670

Paul McCarthy's avatar
Paul McCarthy committed
671
        log.debug('Affine changed:\npixdims: '
672
673
674
675
                  '%s\nsform: %s\nqform: %s',
                  header.get_zooms(),
                  header.get_sform(),
                  header.get_qform())
676
677
678
679

        self.notify(topic='transform')


680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
    @property
    def voxToScaledVoxMat(self):
        """Returns a transformation matrix which transforms from voxel
        coordinates into scaled voxel coordinates, with a left-right flip
        if the image appears to be stored in neurological order.

        See http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT/FAQ#What_is_the\
        _format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to\
        _the_transformation_parameters.3F
        """
        return self.getAffine('voxel', 'fsl')


    @property
    def scaledVoxToVoxMat(self):
        """Returns a transformation matrix which transforms from scaled voxels
        into voxels, the inverse of the :meth:`voxToScaledVoxMat` transform.
        """
        return self.getAffine('fsl', 'voxel')


701
702
    def mapIndices(self, sliceobj):
        """Adjusts the given slice object so that it may be used to index the
Paul McCarthy's avatar
Paul McCarthy committed
703
        underlying ``nibabel`` NIFTI image object.
704
705

        See the :meth:`__determineShape` method.
706

707
708
709
        :arg sliceobj: Something that can be used to slice a
                       multi-dimensional array, e.g. ``arr[sliceobj]``.
        """
710

711
712
713
        # How convenient - nibabel has a function
        # that does the dirty work for us.
        return fileslice.canonical_slicers(sliceobj, self.__origShape)
714
715


716
    def getXFormCode(self, code=None):
717
        """This method returns the code contained in the NIFTI header,
718
        indicating the space to which the (transformed) image is oriented.
719
720

        The ``code`` parameter may be either ``sform`` (the default) or
721
722
723
724
725
726
727
728
        ``qform`` in which case the corresponding matrix is used.

        :returns: one of the following codes:
                    - :data:`~.constants.NIFTI_XFORM_UNKNOWN`
                    - :data:`~.constants.NIFTI_XFORM_SCANNER_ANAT`
                    - :data:`~.constants.NIFTI_XFORM_ALIGNED_ANAT`
                    - :data:`~.constants.NIFTI_XFORM_TALAIRACH`
                    - :data:`~.constants.NIFTI_XFORM_MNI_152`
729
                    - :data:`~.constants.NIFTI_XFORM_ANALYZE`
730
731
        """

732
733
734
        if self.niftiVersion == 0:
            return constants.NIFTI_XFORM_ANALYZE

735
        if   code == 'sform' : code = 'sform_code'
736
        elif code == 'qform' : code = 'qform_code'
737
738
739
740
741
742
        elif code is not None:
            raise ValueError('code must be None, sform, or qform')

        if code is not None:
            code = self.header[code]

743
744
745
746
747
        # If the caller did not specify
        # a code, we check both. If the
        # sform is present, we return it.
        # Otherwise, if the qform is
        # present, we return that.
748
        else:
749

750
751
            sform_code = self.header['sform_code']
            qform_code = self.header['qform_code']
752

753
754
            if   sform_code != constants.NIFTI_XFORM_UNKNOWN: code = sform_code
            elif qform_code != constants.NIFTI_XFORM_UNKNOWN: code = qform_code
755

756
        # Invalid values
757
758
        if code not in range(5):
            code = constants.NIFTI_XFORM_UNKNOWN
759

760
        return int(code)
761

762

763
764
765
766
    # TODO Check what has worse performance - hashing
    #      a 4x4 array (via memoizeMD5), or the call
    #      to aff2axcodes. I'm guessing the latter,
    #      but am not 100% sure.
767
    @memoize.Instanceify(memoize.memoizeMD5)
768
769
770
771
772
773
774
775
776
777
778
779
    def axisMapping(self, xform):
        """Returns the (approximate) correspondence of each axis in the source
        coordinate system to the axes in the destination coordinate system,
        where the source and destinations are defined by the given affine
        transformation matrix.
        """

        inaxes = [[-1, 1], [-2, 2], [-3, 3]]

        return nib.orientations.aff2axcodes(xform, inaxes)


780
    def isNeurological(self):
781
782
783
784
785
        """Returns ``True`` if it looks like this ``Nifti`` object has a
        neurological voxel orientation, ``False`` otherwise. This test is
        purely based on the determinant of the voxel-to-mm transformation
        matrix - if it has a positive determinant, the image is assumed to
        be in neurological orientation, otherwise it is assumed to be in
786
787
        radiological orientation.

Paul McCarthy's avatar
Paul McCarthy committed
788
789
790
791
792
793
        ..warning:: This method will return ``True`` for images with an
                    unknown orientation (e.g. the ``sform_code`` and
                    ``qform_code`` are both set to ``0``). Therefore, you
                    must check the orientation via the :meth:`getXFormCode`
                    before trusting the result of this method.

Paul McCarthy's avatar
Paul McCarthy committed
794
795
796
        See http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT/FAQ#What_is_the\
        _format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to\
        _the_transformation_parameters.3F
797
        """
798
        return self.__isNeurological
799
800


801
802
803
804
805
    def sameSpace(self, other):
        """Returns ``True`` if the ``other`` image (assumed to be a
        :class:`Nifti` instance) has the same dimensions and is in the
        same space as this image.
        """
Paul McCarthy's avatar
Paul McCarthy committed
806
807
808
809
810
811
        return np.all(np.isclose(self .shape[:3],
                                 other.shape[:3]))  and \
               np.all(np.isclose(self .pixdim[:3],
                                 other.pixdim[:3])) and \
               np.all(np.isclose(self .voxToWorldMat,
                                 other.voxToWorldMat))
812
813


814
    def getOrientation(self, axis, xform):
815
        """Returns a code representing the orientation of the specified
816
        axis in the input coordinate system of the given transformation
817
818
819
        matrix.

        :arg xform: A transformation matrix which is assumed to transform
820
821
822
823
824
825
826
                    coordinates from some coordinate system (the one
                    which you want an orientation for) into the image
                    world coordinate system.

                    For example, if you pass in the voxel-to-world
                    transformation matrix, you will get an orientation
                    for axes in the voxel coordinate system.
827
828
829

        This method returns one of the following values, indicating the
        direction in which coordinates along the specified axis increase:
830

831
832
833
834
835
836
837
          - :attr:`~.constants.ORIENT_L2R`:     Left to right
          - :attr:`~.constants.ORIENT_R2L`:     Right to left
          - :attr:`~.constants.ORIENT_A2P`:     Anterior to posterior
          - :attr:`~.constants.ORIENT_P2A`:     Posterior to anterior
          - :attr:`~.constants.ORIENT_I2S`:     Inferior to superior
          - :attr:`~.constants.ORIENT_S2I`:     Superior to inferior
          - :attr:`~.constants.ORIENT_UNKNOWN`: Orientation is unknown
838
839

        The returned value is dictated by the XForm code contained in the
840
841
842
843
844
845
        image file header (see the :meth:`getXFormCode` method). Basically, if
        the XForm code is *unknown*, this method will return
        ``ORIENT_UNKNOWN`` for all axes. Otherwise, it is assumed that the
        image is in RAS orientation (i.e. the X axis increases from left to
        right, the Y axis increases from posterior to anterior, and the Z axis
        increases from inferior to superior).
846
847
        """

848
        if self.getXFormCode() == constants.NIFTI_XFORM_UNKNOWN:
849
850
            return constants.ORIENT_UNKNOWN

851
        code = nib.orientations.aff2axcodes(
852
            xform,
853
854
855
            ((constants.ORIENT_R2L, constants.ORIENT_L2R),
             (constants.ORIENT_A2P, constants.ORIENT_P2A),
             (constants.ORIENT_S2I, constants.ORIENT_I2S)))[axis]
856

857
        return code
858
859


860
    def adjust(self, pixdim=None, shape=None, origin=None):
Paul McCarthy's avatar
Paul McCarthy committed
861
862
        """Return a new ``Nifti`` object with the specified ``pixdim`` or
        ``shape``. The affine of the new ``Nifti`` is adjusted accordingly.
863
864
865

        Only one of ``pixdim`` or ``shape`` can be specified.

Paul McCarthy's avatar
Paul McCarthy committed
866
        See :func:`.affine.rescale` for the meaning of the ``origin`` argument.
867

Paul McCarthy's avatar
Paul McCarthy committed
868
869
870
871
        Only the spatial dimensions may be adjusted - use the functions in
        the :mod:`.image.resample` module if you need to adjust non-spatial
        dimensions.

872
873
874
875
876
877
878
879
880
881
882
883
884
        :arg pixdim: New voxel dimensions
        :arg shape:  New image shape
        :arg origin: Voxel grid alignment - either ``'centre'`` (the default)
                     or ``'corner'``
        :returns:    A new ``Nifti`` object based on this one, with adjusted
                     pixdims, shape and affine.
        """

        if ((pixdim is not None) and (shape is not None)) or \
           ((pixdim is     None) and (shape is     None)):
            raise ValueError('Exactly one of pixdim or '
                             'shape must be specified')

Paul McCarthy's avatar
Paul McCarthy committed
885
886
887
888
889
890
891
892
893
894
        if shape is not None: ndim = len(shape)
        else:                 ndim = len(pixdim)

        # We only allow adjustment of
        # the spatial dimensions
        if ndim != 3:
            raise ValueError('Three dimensions must be specified')

        oldShape  = np.array(self.shape[ :ndim])
        oldPixdim = np.array(self.pixdim[:ndim])
895
896
897
898
        newShape  = shape
        newPixdim = pixdim

        # if pixdims were specified,
Paul McCarthy's avatar
Paul McCarthy committed
899
900
        # convert them into a shape,
        # and vice versa
901
        if newPixdim is not None:
Paul McCarthy's avatar
Paul McCarthy committed
902
903
904
            newShape = oldShape * (oldPixdim / newPixdim)
        else:
            newPixdim = oldPixdim * (oldShape / newShape)
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926

        # Rescale the voxel-to-world affine
        xform = affine.rescale(oldShape, newShape, origin)
        xform = affine.concat(self.getAffine('voxel', 'world'), xform)

        # Now that we've got our spatial
        # scaling/offset matrix, pad the
        # new shape/pixdims with those
        # from any non-spatial dimensions
        newShape  = list(newShape)  + list(self.shape[ 3:])
        newPixdim = list(newPixdim) + list(self.pixdim[3:])

        # And create the new header
        # and we're away
        header = self.header.copy()
        header.set_data_shape(newShape)
        header.set_zooms(     newPixdim)
        header.set_sform(     xform)
        header.set_qform(     xform)
        return Nifti(header)


927
class Image(Nifti):
928
    """Class which represents a NIFTI image. Internally, the image is
Paul McCarthy's avatar
Paul McCarthy committed
929
930
931
    loaded/stored using a :mod:`nibabel.nifti1.Nifti1Image` or
    :mod:`nibabel.nifti2.Nifti2Image`, and data access managed by a
    :class:`.ImageWrapper`.
932

933

934
    In addition to the attributes added by the :meth:`Nifti.__init__` method,
935
    the following attributes/properties are present on an ``Image`` instance
936
    as properties (https://docs.python.org/2/library/functions.html#property):
937
938


939
    ============== ===========================================================
940
941
    ``name``       The name of this ``Image`` - defaults to the image
                   file name, sans-suffix.
942

943
944
945
    ``dataSource`` The data source of this ``Image`` - the name of the
                   file from where it was loaded, or some other string
                   describing its origin.
946

Paul McCarthy's avatar
Paul McCarthy committed
947
    ``nibImage``   A reference to the ``nibabel`` NIFTI image object.
948

949
950
951
    ``saveState``  A boolean value which is ``True`` if this image is
                   saved to disk, ``False`` if it is in-memory, or has
                   been edited.
952

953
954
955
    ``dataRange``  The minimum/maximum values in the image. Depending upon
                   the value of the ``calcRange`` parameter to
                   :meth:`__init__`, this may be calculated when the ``Image``
956
                   is created, or may be incrementally updated as more image
957
                   data is loaded from disk.
958
    ============== ===========================================================
959

960

961
962
963
964
    The ``Image`` class adds some :class:`.Notifier` topics to those which are
    already provided by the :class:`Nifti` class - listeners may register to
    be notified of changes to the above properties, by registering on the
    following _topic_ names (see the :class:`.Notifier` class documentation):
965

966

967
968
    =============== ======================================================
    ``'data'``      This topic is notified whenever the image data changes
969
970
971
972
                    (via the :meth:`__setitem__` method). The indices/
                    slices of the portion of data that was modified is
                    passed to registered listeners as the notification
                    value (see :meth:`.Notifier.notify`).
973

974
    ``'saveState'`` This topic is notified whenever the saved state of the
975
976
                    image changes (i.e. data or ``voxToWorldMat`` is
                    edited, or the image saved to disk).
977

978
979
980
    ``'dataRange'`` This topic is notified whenever the image data range
                    is changed/adjusted.
    =============== ======================================================
981
982
983
    """


984
985
986
987
988
    def __init__(self,
                 image,
                 name=None,
                 header=None,
                 xform=None,
989
                 loadData=True,
990
                 calcRange=True,
991
                 indexed=False,
992
                 threaded=False,
993
                 dataSource=None,
994
                 loadMeta=False,
995
                 **kwargs):
996
997
        """Create an ``Image`` object with the given image data or file name.

998
999
1000
        :arg image:      A string containing the name of an image file to load,
                         or a :mod:`numpy` array, or a :mod:`nibabel` image
                         object.
For faster browsing, not all history is shown. View entire blame