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

43
import                      six
44
45
46
47
import numpy             as np

import nibabel           as nib
import nibabel.fileslice as fileslice
48

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

58

59
60
log = logging.getLogger(__name__)

61

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


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


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

98

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

102

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

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

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

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

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

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

126

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

132

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

138

139
    **Affine transformations**
140

141

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

167

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

172

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

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

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

191
192
193
194

    **ANALYZE support**


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

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

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

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


212
213
214
    **Metadata**


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


219
220
221
    **Notification**


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

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

233

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

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

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

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

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

262

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


268
269
    @staticmethod
    def determineShape(header):
270
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
        """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)]

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

302
303
304
        return origShape, shape, pixdims


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

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

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

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

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

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

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

370
        return voxToWorldMat
371

372

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

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

394
        import numpy.linalg as npla
395

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

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

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

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


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

        val = val.decode('ascii')

495
        return ''.join([c for c in val if c in string.printable]).strip()
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')
513
        self.__header = header
514
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
    @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)

540

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


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


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

562
563
564
565
566
    @property
    def niftiDataType(self):
        """Returns the NIFTI data type code of this image. """
        return self.header.get('datatype', constants.NIFTI_DT_UNKNOWN)

567
568
569
570
571

    @intent.setter
    def intent(self, val):
        """Sets the NIFTI intent code of this image. """
        # analyze has no intent
572
        if (self.niftiVersion > 0) and (val != self.intent):
573
574
            self.header.set_intent(val, allow_unknown=True)
            self.notify(topic='header')
575
576


577
578
579
580
    @property
    def xyzUnits(self):
        """Returns the NIFTI XYZ dimension unit code. """

581
582
583
584
        # analyze images have no unit field
        if self.niftiVersion == 0:
            return constants.NIFTI_UNITS_MM

585
586
587
588
        # 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
589
590
591
        units = self.header.get_xyzt_units()[0]
        units = nib.nifti1.unit_codes[units]
        return units
592
593


Paul McCarthy's avatar
Paul McCarthy committed
594
    @property
595
596
597
    def timeUnits(self):
        """Returns the NIFTI time dimension unit code. """

598
599
600
601
        # analyze images have no unit field
        if self.niftiVersion == 0:
            return constants.NIFTI_UNITS_SEC

602
        # See xyzUnits
Paul McCarthy's avatar
Paul McCarthy committed
603
604
605
        units = self.header.get_xyzt_units()[1]
        units = nib.nifti1.unit_codes[units]
        return units
606
607


608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
    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'):
629
630
631
            raise ValueError('Invalid source/reference spaces: "{}" -> "{}".'
                             'Recognised spaces are "voxel", "fsl", and '
                             '"world"'.format(from_, to))
632
633
634
635

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


636
637
638
639
640
    @property
    def worldToVoxMat(self):
        """Returns a ``numpy`` array of shape ``(4, 4)`` containing an
        affine transformation from world coordinates to voxel coordinates.
        """
641
        return self.getAffine('world', 'voxel')
642
643
644
645
646
647
648


    @property
    def voxToWorldMat(self):
        """Returns a ``numpy`` array of shape ``(4, 4)`` containing an
        affine transformation from voxel coordinates to world coordinates.
        """
649
        return self.getAffine('voxel', 'world')
650
651
652
653


    @voxToWorldMat.setter
    def voxToWorldMat(self, xform):
654
655
656
        """Update the ``voxToWorldMat``. The ``worldToVoxMat`` value is also
        updated. This will result in notification on the ``'transform'``
        topic.
657
658
        """

659
660
661
662
663
664
665
        # 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
666
        sformCode = int(header['sform_code'])
667

668
669
670
671
        if sformCode == constants.NIFTI_XFORM_UNKNOWN:
            sformCode = constants.NIFTI_XFORM_ALIGNED_ANAT

        header.set_sform(xform, code=sformCode)
672

673
674
675
        affines, isneuro = Nifti.generateAffines(xform,
                                                 self.shape,
                                                 self.pixdim)
676
677
678

        self.__affines        = affines
        self.__isNeurological = isneuro
679

Paul McCarthy's avatar
Paul McCarthy committed
680
        log.debug('Affine changed:\npixdims: '
681
682
683
684
                  '%s\nsform: %s\nqform: %s',
                  header.get_zooms(),
                  header.get_sform(),
                  header.get_qform())
685
686
687
688

        self.notify(topic='transform')


689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
    @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')


710
711
    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
712
        underlying ``nibabel`` NIFTI image object.
713
714

        See the :meth:`__determineShape` method.
715

716
717
718
        :arg sliceobj: Something that can be used to slice a
                       multi-dimensional array, e.g. ``arr[sliceobj]``.
        """
719

720
721
722
        # How convenient - nibabel has a function
        # that does the dirty work for us.
        return fileslice.canonical_slicers(sliceobj, self.__origShape)
723
724


725
    def getXFormCode(self, code=None):
726
        """This method returns the code contained in the NIFTI header,
727
        indicating the space to which the (transformed) image is oriented.
728
729

        The ``code`` parameter may be either ``sform`` (the default) or
730
731
732
733
734
735
736
737
        ``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`
738
                    - :data:`~.constants.NIFTI_XFORM_ANALYZE`
739
740
        """

741
742
743
        if self.niftiVersion == 0:
            return constants.NIFTI_XFORM_ANALYZE

744
        if   code == 'sform' : code = 'sform_code'
745
        elif code == 'qform' : code = 'qform_code'
746
747
748
749
750
751
        elif code is not None:
            raise ValueError('code must be None, sform, or qform')

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

752
753
754
755
756
        # 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.
757
        else:
758

759
760
            sform_code = self.header['sform_code']
            qform_code = self.header['qform_code']
761

762
763
            if   sform_code != constants.NIFTI_XFORM_UNKNOWN: code = sform_code
            elif qform_code != constants.NIFTI_XFORM_UNKNOWN: code = qform_code
764

765
        # Invalid values
766
767
        if code not in range(5):
            code = constants.NIFTI_XFORM_UNKNOWN
768

769
        return int(code)
770

771

772
773
774
775
    # 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.
776
    @memoize.Instanceify(memoize.memoizeMD5)
777
778
779
780
781
782
783
784
785
786
787
788
    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)


789
    def isNeurological(self):
790
791
792
793
794
        """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
795
796
        radiological orientation.

Paul McCarthy's avatar
Paul McCarthy committed
797
798
799
800
801
802
        ..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
803
804
805
        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
806
        """
807
        return self.__isNeurological
808
809


810
811
812
813
814
    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
815
816
817
818
819
820
        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))
821
822


823
    def getOrientation(self, axis, xform):
824
        """Returns a code representing the orientation of the specified
825
        axis in the input coordinate system of the given transformation
826
827
828
        matrix.

        :arg xform: A transformation matrix which is assumed to transform
829
830
831
832
833
834
835
                    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.
836
837
838

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

840
841
842
843
844
845
846
          - :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
847
848

        The returned value is dictated by the XForm code contained in the
849
850
851
852
853
854
        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).
855
856
        """

857
        if self.getXFormCode() == constants.NIFTI_XFORM_UNKNOWN:
858
859
            return constants.ORIENT_UNKNOWN

860
        code = nib.orientations.aff2axcodes(
861
            xform,
862
863
864
            ((constants.ORIENT_R2L, constants.ORIENT_L2R),
             (constants.ORIENT_A2P, constants.ORIENT_P2A),
             (constants.ORIENT_S2I, constants.ORIENT_I2S)))[axis]
865

866
        return code
867
868


869
    def adjust(self, pixdim=None, shape=None, origin=None):
Paul McCarthy's avatar
Paul McCarthy committed
870
871
        """Return a new ``Nifti`` object with the specified ``pixdim`` or
        ``shape``. The affine of the new ``Nifti`` is adjusted accordingly.
872
873
874

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

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

Paul McCarthy's avatar
Paul McCarthy committed
877
878
879
880
        Only the spatial dimensions may be adjusted - use the functions in
        the :mod:`.image.resample` module if you need to adjust non-spatial
        dimensions.

881
882
883
884
885
886
887
888
889
890
891
892
893
        :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
894
895
896
897
898
899
900
901
902
903
        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])
904
905
906
907
        newShape  = shape
        newPixdim = pixdim

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

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


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

942

943
    In addition to the attributes added by the :meth:`Nifti.__init__` method,
944
    the following attributes/properties are present on an ``Image`` instance
945
    as properties (https://docs.python.org/2/library/functions.html#property):
946
947


948
    ============== ===========================================================
949
950
    ``name``       The name of this ``Image`` - defaults to the image
                   file name, sans-suffix.
951

952
953
954
    ``dataSource`` The data source of this ``Image`` - the name of the
                   file from where it was loaded, or some other string
                   describing its origin.
955

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

958
959
960
    ``saveState``  A boolean value which is ``True`` if this image is
                   saved to disk, ``False`` if it is in-memory, or has
                   been edited.
961

962
963
964
    ``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``
965
                   is created, or may be incrementally updated as more image
966
                   data is loaded from disk.
967
    ============== ===========================================================
968

969

970
971
972
973
    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):
974

975

976
977
    =============== ======================================================
    ``'data'``      This topic is notified whenever the image data changes
978
979
980
981
                    (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`).
982

983
    ``'saveState'`` This topic is notified whenever the saved state of the
984
985
                    image changes (i.e. data or ``voxToWorldMat`` is
                    edited, or the image saved to disk).
986

987
988
989
    ``'dataRange'`` This topic is notified whenever the image data range
                    is changed/adjusted.
    =============== ======================================================
990
991
992
    """


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

1006
1007
        :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
1008
                         object, or an ``Image`` object.
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
1043
1044
1045
1046
1047

        :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.
1048

1049
1050
1051
1052
1053
1054
        :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``.

1055
1056
        All other arguments are passed through to the ``nibabel.load`` function
        (if it is called).
1057
        """
1058

1059
        nibImage = None
1060
        saved    = False
1061

1062
1063
1064
        if loadData:
            threaded = False

1065
1066
1067
1068
1069
1070
1071
1072
1073
        # 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()

1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
        # 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)

1085
1086
        # The image parameter may be the name of an image file
        if isinstance(image, six.string_types):
1087
1088
            image      = op.abspath(addExt(image))
            nibImage   = nib.load(image, **kwargs)
1089
            dataSource = image
1090
            saved      = True
1091

1092
1093
1094
1095
1096
        # 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):

1097
1098
1099
            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
1100

1101
1102
            # default to NIFTI1 if FSLOUTPUTTYPE
            # is not set, just to be safe.
1103
            if header is None:
1104
1105
1106
                outputType = os.environ.get('FSLOUTPUTTYPE', 'NIFTI_GZ')
                if 'NIFTI2' in outputType: ctr = nib.Nifti2Image
                else:                      ctr = nib.Nifti1Image
1107

1108
1109
1110
1111
1112
1113
            # 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)

1114
1115
1116
1117
1118
            # 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):
1119
                ctr = nib.nifti1.Nifti1Image
1120
1121
1122
1123
            elif isinstance(header, nib.analyze.AnalyzeHeader):
                ctr = nib.analyze.AnalyzeImage

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

1125
1126
1127
1128
1129
1130
1131
        # 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
1132
        else:
1133
            nibImage = image
1134

1135
        # Figure out the name of this image, if
1136
1137
        # it has not beenbeen explicitly passed in
        if name is None:
1138

1139
1140
1141
1142
            # If this image was loaded
            # from disk, use the file name.
            if isinstance(image, six.string_types):
                name = removeExt(op.basename(image))
1143

1144
1145
1146
            # Or the image was created from a numpy array
            elif isinstance(image, np.ndarray):
                name = 'Numpy array'
1147

1148
1149
1150
            # Or image from a nibabel image
            else:
                name = 'Nibabel image'
1151

1152
        Nifti.__init__(self, nibImage.header)
1153

1154
1155
1156
1157
1158
        self.name           = name
        self.__lName        = '{}_{}'.format(id(self), self.name)
        self.__dataSource   = dataSource
        self.__threaded     = threaded
        self.__nibImage     = nibImage
1159
        self.__saveState    = saved
1160
1161
1162
1163
        self.__imageWrapper = imagewrapper.ImageWrapper(self.nibImage,
                                                        self.name,
                                                        loadData=loadData,
                                                        threaded=threaded)