image.py 58.9 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.bids        as fslbids
56
57
import fsl.data.constants    as constants
import fsl.data.imagewrapper as imagewrapper
58

59

60
61
log = logging.getLogger(__name__)

62

Paul McCarthy's avatar
Paul McCarthy committed
63
64
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
65
66
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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
"""


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


92
class Nifti(notifier.Notifier, meta.Meta):
93
94
95
96
    """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
97
    a file, use the :class:`Image` class instead.
Paul McCarthy's avatar
Paul McCarthy committed
98

99

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

103

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

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

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

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

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

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

127

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

133

134
135
136
137
138
    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.

139

140
    **Affine transformations**
141

142

143
144
145
146
147
148
149
150
    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
151
      - The ``fsl`` coordinate system, where voxel coordinates are scaled by
152
153
154
155
156
157
158
159
160
161
162
163
164
165
        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
166
167
    http://nipy.org/nibabel/nifti_images.html#the-nifti-affines):

168

169
170
171
172
     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.

173

174
175
176
177
178
179
180
181
182
183
184
    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
185

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

188
189
190
     - 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.
191

192
193
194
195

    **ANALYZE support**


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

      - 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
204
        the ANALYZE ``origin`` fields. Construction of the affine is handled
205
206
207
208
209
210
211
212
        by ``nibabel``.

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

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


213
214
215
    **Metadata**


216
217
    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.
218
219


220
221
222
    **Notification**


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

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

234

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

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

Paul McCarthy's avatar
Paul McCarthy committed
244
        # Nifti2Header is a sub-class of Nifti1Header,
245
246
247
        # 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
248
249
            raise ValueError('Unrecognised header: {}'.format(header))

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

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

263

264
265
266
267
268
    def __del__(self):
        """Clears the reference to the ``nibabel`` header object. """
        self.__header = None


269
270
    @staticmethod
    def determineShape(header):
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
        """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)]

297
298
299
300
301
302
        # 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))

303
304
305
        return origShape, shape, pixdims


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

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

321
322
323
324
325
326
327
328
        # 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.
329
330
331
332
333
        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):
334
335
336

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

Paul McCarthy's avatar
Paul McCarthy committed
337
338
339
340
341
342
343
            # 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
344
            knotpix       =  header.get_zooms()[:3]
345
346
347
            refpix        = (header.get('intent_p1', 1) or 1,
                             header.get('intent_p2', 1) or 1,
                             header.get('intent_p3', 1) or 1)
348
349
350
            voxToWorldMat = affine.concat(
                affine.scaleOffsetXform(refpix,  0),
                affine.scaleOffsetXform(knotpix, 0))
351

352
353
354
355
        # 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.
356
        #
357
358
359
        # n.b. For images like this, nibabel returns
        # a scaling matrix where the centre voxel
        # corresponds to world location (0, 0, 0).
360
        # This goes against the NIFTI spec - it
361
        # should just be a straight scaling matrix.
362
        elif qform == 0 and sform == 0:
363
            pixdims       = header.get_zooms()
364
            voxToWorldMat = affine.scaleOffsetXform(pixdims, 0)
365

366
        # Otherwise we let nibabel decide
367
        # which transform to use.
368
        else:
369
            voxToWorldMat = np.array(header.get_best_affine())
370

371
        return voxToWorldMat
372

373

374
375
    @staticmethod
    def generateAffines(voxToWorldMat, shape, pixdim):
376
377
378
379
380
        """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.
381

382
        :arg voxToWorldMat: The voxel-to-world affine transformation
383
384
385
        :arg shape:         Image shape (number of voxels along each dimension
        :arg pixdim:        Image pixdims (size of one voxel along each
                            dimension)
386
387
388
389
390
391
392
        :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.
393
394
        """

395
        import numpy.linalg as npla
396

397
        affines = {}
398
399
        shape             = list(shape[ :3])
        pixdim            = list(pixdim[:3])
400
401
        voxToScaledVoxMat = np.diag(pixdim + [1.0])
        isneuro           = npla.det(voxToWorldMat) > 0
402

403
404
        if isneuro:
            x                 = (shape[0] - 1) * pixdim[0]
405
406
407
            flip              = affine.scaleOffsetXform([-1, 1, 1],
                                                        [ x, 0, 0])
            voxToScaledVoxMat = affine.concat(flip, voxToScaledVoxMat)
408

409
410
411
        affines['fsl',   'fsl']   = np.eye(4)
        affines['voxel', 'voxel'] = np.eye(4)
        affines['world', 'world'] = np.eye(4)
412
        affines['voxel', 'world'] = voxToWorldMat
413
        affines['world', 'voxel'] = affine.invert(voxToWorldMat)
414
        affines['voxel', 'fsl']   = voxToScaledVoxMat
415
416
417
418
419
        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'])
420
421

        return affines, isneuro
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
468
469
470
471
    @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')


472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
    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
491
492
        try:              val = bytes(val).partition(b'\0')[0]
        except Exception: val = bytes(val)
493
494
495

        val = val.decode('ascii')

496
        return ''.join([c for c in val if c in string.printable]).strip()
497
498


499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
    @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')
514
        self.__header = header
515
516


517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
    @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)

541

542
543
544
545
546
547
548
549
550
551
    @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)


552
553
554
555
556
557
558
559
    @property
    def pixdim(self):
        """Returns a tuple containing the image pixdims (voxel sizes)."""
        return tuple(self.__pixdim)


    @property
    def intent(self):
560
561
562
563
564
565
566
567
        """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
568
        if (self.niftiVersion > 0) and (val != self.intent):
569
570
            self.header.set_intent(val, allow_unknown=True)
            self.notify(topic='header')
571
572


573
574
575
576
    @property
    def xyzUnits(self):
        """Returns the NIFTI XYZ dimension unit code. """

577
578
579
580
        # analyze images have no unit field
        if self.niftiVersion == 0:
            return constants.NIFTI_UNITS_MM

581
582
583
584
        # 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
585
586
587
        units = self.header.get_xyzt_units()[0]
        units = nib.nifti1.unit_codes[units]
        return units
588
589


Paul McCarthy's avatar
Paul McCarthy committed
590
    @property
591
592
593
    def timeUnits(self):
        """Returns the NIFTI time dimension unit code. """

594
595
596
597
        # analyze images have no unit field
        if self.niftiVersion == 0:
            return constants.NIFTI_UNITS_SEC

598
        # See xyzUnits
Paul McCarthy's avatar
Paul McCarthy committed
599
600
601
        units = self.header.get_xyzt_units()[1]
        units = nib.nifti1.unit_codes[units]
        return units
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
627
628
629
630
    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])


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


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


    @voxToWorldMat.setter
    def voxToWorldMat(self, xform):
649
650
651
        """Update the ``voxToWorldMat``. The ``worldToVoxMat`` value is also
        updated. This will result in notification on the ``'transform'``
        topic.
652
653
        """

654
655
656
657
658
659
660
        # 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
661
        sformCode = int(header['sform_code'])
662

663
664
665
666
        if sformCode == constants.NIFTI_XFORM_UNKNOWN:
            sformCode = constants.NIFTI_XFORM_ALIGNED_ANAT

        header.set_sform(xform, code=sformCode)
667

668
669
670
        affines, isneuro = Nifti.generateAffines(xform,
                                                 self.shape,
                                                 self.pixdim)
671
672
673

        self.__affines        = affines
        self.__isNeurological = isneuro
674

Paul McCarthy's avatar
Paul McCarthy committed
675
        log.debug('Affine changed:\npixdims: '
676
677
678
679
                  '%s\nsform: %s\nqform: %s',
                  header.get_zooms(),
                  header.get_sform(),
                  header.get_qform())
680
681
682
683

        self.notify(topic='transform')


684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
    @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')


705
706
    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
707
        underlying ``nibabel`` NIFTI image object.
708
709

        See the :meth:`__determineShape` method.
710

711
712
713
        :arg sliceobj: Something that can be used to slice a
                       multi-dimensional array, e.g. ``arr[sliceobj]``.
        """
714

715
716
717
        # How convenient - nibabel has a function
        # that does the dirty work for us.
        return fileslice.canonical_slicers(sliceobj, self.__origShape)
718
719


720
    def getXFormCode(self, code=None):
721
        """This method returns the code contained in the NIFTI header,
722
        indicating the space to which the (transformed) image is oriented.
723
724

        The ``code`` parameter may be either ``sform`` (the default) or
725
726
727
728
729
730
731
732
        ``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`
733
                    - :data:`~.constants.NIFTI_XFORM_ANALYZE`
734
735
        """

736
737
738
        if self.niftiVersion == 0:
            return constants.NIFTI_XFORM_ANALYZE

739
        if   code == 'sform' : code = 'sform_code'
740
        elif code == 'qform' : code = 'qform_code'
741
742
743
744
745
746
        elif code is not None:
            raise ValueError('code must be None, sform, or qform')

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

747
748
749
750
751
        # 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.
752
        else:
753

754
755
            sform_code = self.header['sform_code']
            qform_code = self.header['qform_code']
756

757
758
            if   sform_code != constants.NIFTI_XFORM_UNKNOWN: code = sform_code
            elif qform_code != constants.NIFTI_XFORM_UNKNOWN: code = qform_code
759

760
        # Invalid values
761
762
        if code not in range(5):
            code = constants.NIFTI_XFORM_UNKNOWN
763

764
        return int(code)
765

766

767
768
769
770
    # 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.
771
    @memoize.Instanceify(memoize.memoizeMD5)
772
773
774
775
776
777
778
779
780
781
782
783
    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)


784
    def isNeurological(self):
785
786
787
788
789
        """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
790
791
        radiological orientation.

Paul McCarthy's avatar
Paul McCarthy committed
792
793
794
795
796
797
        ..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
798
799
800
        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
801
        """
802
        return self.__isNeurological
803
804


805
806
807
808
809
    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
810
811
812
813
814
815
        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))
816
817


818
    def getOrientation(self, axis, xform):
819
        """Returns a code representing the orientation of the specified
820
        axis in the input coordinate system of the given transformation
821
822
823
        matrix.

        :arg xform: A transformation matrix which is assumed to transform
824
825
826
827
828
829
830
                    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.
831
832
833

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

835
836
837
838
839
840
841
          - :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
842
843

        The returned value is dictated by the XForm code contained in the
844
845
846
847
848
849
        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).
850
851
        """

852
        if self.getXFormCode() == constants.NIFTI_XFORM_UNKNOWN:
853
854
            return constants.ORIENT_UNKNOWN

855
        code = nib.orientations.aff2axcodes(
856
            xform,
857
858
859
            ((constants.ORIENT_R2L, constants.ORIENT_L2R),
             (constants.ORIENT_A2P, constants.ORIENT_P2A),
             (constants.ORIENT_S2I, constants.ORIENT_I2S)))[axis]
860

861
        return code
862
863


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

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

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

Paul McCarthy's avatar
Paul McCarthy committed
872
873
874
875
        Only the spatial dimensions may be adjusted - use the functions in
        the :mod:`.image.resample` module if you need to adjust non-spatial
        dimensions.

876
877
878
879
880
881
882
883
884
885
886
887
888
        :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
889
890
891
892
893
894
895
896
897
898
        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])
899
900
901
902
        newShape  = shape
        newPixdim = pixdim

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

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


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

937

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


943
    ============== ===========================================================
944
945
    ``name``       The name of this ``Image`` - defaults to the image
                   file name, sans-suffix.
946

947
948
949
    ``dataSource`` The data source of this ``Image`` - the name of the
                   file from where it was loaded, or some other string
                   describing its origin.
950

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

953
954
955
    ``saveState``  A boolean value which is ``True`` if this image is
                   saved to disk, ``False`` if it is in-memory, or has
                   been edited.
956

957
958
959
    ``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``
960
                   is created, or may be incrementally updated as more image
961
                   data is loaded from disk.
962
    ============== ===========================================================
963

964

965
966
967
968
    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):
969

970

971
972
    =============== ======================================================
    ``'data'``      This topic is notified whenever the image data changes
973
974
975
976
                    (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`).
977

978
    ``'saveState'`` This topic is notified whenever the saved state of the
979
980
                    image changes (i.e. data or ``voxToWorldMat`` is
                    edited, or the image saved to disk).
981

982
983
984
    ``'dataRange'`` This topic is notified whenever the image data range
                    is changed/adjusted.
    =============== ======================================================
985
986
987
    """


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

1001
1002
        :arg image:      A string containing the name of an image file to load,
                         or a :mod:`numpy` array, or a :mod:`nibabel` image
Paul McCarthy's avatar
Paul McCarthy committed
1003
                         object, or an ``Image``object.
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042

        :arg name:       A name for the image.

        :arg header:     If not ``None``, assumed to be a
                         :class:`nibabel.nifti1.Nifti1Header` or
                         :class:`nibabel.nifti2.Nifti2Header` to be used as the
                         image header. Not applied to images loaded from file,
                         or existing :mod:`nibabel` images.

        :arg xform:      A :math:`4\\times 4` affine transformation matrix
                         which transforms voxel coordinates into real world
                         coordinates. If not provided, and a ``header`` is
                         provided, the transformation in the header is used.
                         If neither a ``xform`` nor a ``header`` are provided,
                         an identity matrix is used. If both a ``xform`` and a
                         ``header`` are provided, the ``xform`` is used in
                         preference to the header transformation.

        :arg loadData:   If ``True`` (the default) the image data is loaded
                         in to memory.  Otherwise, only the image header
                         information is read, and the image data is kept
                         from disk. In either case, the image data is
                         accessed through an :class:`.ImageWrapper` instance.
                         The data may be loaded into memory later on via the
                         :meth:`loadData` method.

        :arg calcRange:  If ``True`` (the default), the image range is
                         calculated immediately (vi a call to
                         :meth:`calcRange`). Otherwise, the image range is
                         incrementally updated as more data is read from memory
                         or disk.

        :arg threaded:   If ``True``, the :class:`.ImageWrapper` will use a
                         separate thread for data range calculation. Defaults
                         to ``False``. Ignored if ``loadData`` is ``True``.

        :arg dataSource: If ``image`` is not a file name, this argument may be
                         used to specify the file from which the image was
                         loaded.
1043

1044
1045
1046
1047
1048
1049
        :arg loadMeta:   If ``True``, any metadata contained in JSON sidecar
                         files is loaded and attached to this ``Image`` via
                         the :class:`.Meta` interface. if ``False``, metadata
                         can be loaded at a later stage via the
                         :func:`loadMeta` function. Defaults to ``False``.

1050
1051
        All other arguments are passed through to the ``nibabel.load`` function
        (if it is called).
1052
        """
1053

1054
        nibImage = None
1055
        saved    = False
1056

1057
1058
1059
        if loadData:
            threaded = False

1060
1061
1062
1063
1064
1065
1066
1067
1068
        # Take a copy of the header if one has
        # been provided
        #
        # NOTE: Nifti extensions are copied by
        # reference, which may cause issues in
        # the future.
        if header is not None:
            header = header.copy()

1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
        # if a header and xform are provided,
        # make sure the xform gets used. Does
        # not apply to ANALYZE images,
        if header is not None and \
           xform  is not None and \
           isinstance(header, nib.nifti1.Nifti1Header):
            sform = int(header.get_sform(True)[1])
            qform = int(header.get_qform(True)[1])
            header.set_sform(xform, code=sform)
            header.set_qform(xform, code=qform)

1080
1081
        # The image parameter may be the name of an image file
        if isinstance(image, six.string_types):
1082
1083
            image      = op.abspath(addExt(image))
            nibImage   = nib.load(image, **kwargs)
1084
            dataSource = image
1085
            saved      = True
1086

1087
1088
1089
1090
1091
        # Or a numpy array - we wrap it in a nibabel image,
        # with an identity transformation (each voxel maps
        # to 1mm^3 in real world space)
        elif isinstance(image, np.ndarray):

1092
1093
1094
            if xform is None:
                if header is not None: xform = header.get_best_affine()
                else:                  xform = np.identity(4)
Paul McCarthy's avatar
Paul McCarthy committed
1095
1096
1097

            # We default to NIFTI1 and not
            # NIFTI2, because the rest of
1098
            # FSL is not yet NIFTI2 compatible.
1099
1100
1101
            if header is None:
                ctr = nib.nifti1.Nifti1Image

1102
1103
1104
1105
1106
1107
            # make sure that the data type is correct,
            # in case this header was passed in from
            # a different image
            if header is not None:
                header.set_data_dtype(image.dtype)

1108
1109
1110
1111
1112
            # But if a nibabel header has been provided,
            # we use the corresponding image type
            if isinstance(header, nib.nifti2.Nifti2Header):
                ctr = nib.nifti2.Nifti2Image
            elif isinstance(header, nib.nifti1.Nifti1Header):
1113
                ctr = nib.nifti1.Nifti1Image
1114
1115
1116
1117
            elif isinstance(header, nib.analyze.AnalyzeHeader):
                ctr = nib.analyze.AnalyzeImage

            nibImage = ctr(image, xform, header=header)
1118

1119
1120
1121
1122
1123
1124
1125
        # If it's an Image object, we
        # just take the nibabel image
        elif isinstance(image, Image):
            nibImage = image.nibImage

        # otherwise, we assume that
        # it is a nibabel image
1126
        else:
1127
            nibImage = image
1128

1129
        # Figure out the name of this image, if
1130
1131
        # it has not beenbeen explicitly passed in
        if name is None:
1132

1133
1134
1135
1136
            # If this image was loaded
            # from disk, use the file name.
            if isinstance(image, six.string_types):
                name = removeExt(op.basename(image))
1137

1138
1139
1140
            # Or the image was created from a numpy array
            elif isinstance(image, np.ndarray):
                name = 'Numpy array'
1141

1142
1143
1144
            # Or image from a nibabel image
            else:
                name = 'Nibabel image'
1145

1146
        Nifti.__init__(self, nibImage.header)
1147

1148
1149
1150
1151
1152
        self.name           = name
        self.__lName        = '{}_{}'.format(id(self), self.name)
        self.__dataSource   = dataSource
        self.__threaded     = threaded
        self.__nibImage     = nibImage
1153
        self.__saveState    = saved
1154
1155
1156
1157
        self.__imageWrapper = imagewrapper.ImageWrapper(self.nibImage,
                                                        self.name,
                                                        loadData=loadData,
                                                        threaded=threaded)
1158

1159
        # Listen to ourself for changes
1160
        # to header attributse so we
1161
        # can update the saveState.
1162
1163
        self.register(self.name, self.__headerChanged, topic='transform')
        self.register(self.name, self.__headerChanged, topic='header')
1164

1165
1166
        # calculate min/max
        # of image data
1167
1168
1169