image.py 30.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
9
"""This module provides the :class:`Nifti` and :class:`Image` classes, for
representing 3D/4D NIFI1 images. The ``nibabel`` package is used for file
10
I/O.
11

12
.. note:: Support for ANALYZE75 and NIFTI2 images has not been tested.
13
14
15
16
17
18
19
20


It is very easy to load a NIFTI image::

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


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

.. autosummary::
   :nosignatures:

27
   looksLikeImage
28
29
   removeExt
   addExt
30
   loadIndexedImageFile
31
"""
32

33

34
35
import               logging
import os.path    as op
36

Paul McCarthy's avatar
Paul McCarthy committed
37
import               six 
38
39
import numpy      as np

40

41
42
import fsl.utils.transform   as transform
import fsl.utils.notifier    as notifier
43
import fsl.utils.memoize     as memoize
44
45
46
import fsl.utils.path        as fslpath
import fsl.data.constants    as constants
import fsl.data.imagewrapper as imagewrapper
47

48

49
50
log = logging.getLogger(__name__)

51

52
53
54
55
56
class Nifti(object):
    """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
57
    a file, use the :class:`Image` class instead.
Paul McCarthy's avatar
Paul McCarthy committed
58

59

60
    When a ``Nifti`` instance is created, it adds the following attributes
61
    to itself:
62

63
    
64
    ================= ====================================================
65
    ``header``        The :mod:`nibabel.Nifti1Header` object.
66
    ``shape``         A list/tuple containing the number of voxels along
67
68
69
                      each image dimension.
    ``pixdim``        A list/tuple containing the length of one voxel 
                      along each image dimension. 
70
71
72
73
74
75
    ``voxToWorldMat`` A 4*4 array specifying the affine transformation
                      for transforming voxel coordinates into real world
                      coordinates.
    ``worldToVoxMat`` A 4*4 array specifying the affine transformation
                      for transforming real world coordinates into voxel
                      coordinates.
76
    ``intent``        The NIFTI intent code specified in the header.
77
    ================= ====================================================
78
79
80

    
    .. note:: The ``shape`` attribute may not precisely match the image shape
81
              as reported in the NIFTI header, because trailing dimensions of
82
83
              size 1 are squeezed out. See the :meth:`__determineShape` and
              :meth:`mapIndices` methods.
84
85
    """

86
    def __init__(self, header):
87
        """Create a ``Nifti`` object.
Paul McCarthy's avatar
Paul McCarthy committed
88

89
90
        :arg header:   A :class:`nibabel.nifti1.Nifti1Header` to be used as 
                       the image header. 
91
        """
92

93
94
        header                   = header.copy()
        origShape, shape, pixdim = self.__determineShape(header)
95

96
97
        if len(shape) < 3 or len(shape) > 4:
            raise RuntimeError('Only 3D or 4D images are supported')
98

99
100
101
102
103
        voxToWorldMat = self.__determineTransform(header)
        worldToVoxMat = transform.invert(voxToWorldMat)

        self.header        = header
        self.shape         = shape
104
105
        self.intent        = header.get('intent_code',
                                        constants.NIFTI_INTENT_NONE)
106
107
108
109
110
111
112
113
114
        self.__origShape   = origShape
        self.pixdim        = pixdim
        self.voxToWorldMat = voxToWorldMat
        self.worldToVoxMat = worldToVoxMat

        
    def __determineTransform(self, header):
        """Called by :meth:`__init__`. Figures out the voxel-to-world
        coordinate transformation matrix that is associated with this
115
        ``Nifti`` instance.
116
117
        """
        
118
119
120
        # We have to treat FSL/FNIRT images
        # specially, as FNIRT clobbers the
        # sform section of the NIFTI header
121
        # to store other data. 
122
        intent = header.get('intent_code', -1)
123
124
125
126
        if intent in (constants.FSL_FNIRT_DISPLACEMENT_FIELD,
                      constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
                      constants.FSL_DCT_COEFFICIENTS,
                      constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS):
127
            log.debug('FNIRT output image detected - using qform matrix')
128
            voxToWorldMat = np.array(header.get_qform())
129

130
131
132
133
134
135
136
137
        # 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.
        # 
        # n.b. For images like this, nibabel returns
        # a scaling matrix where the centre voxel
        # corresponds to world location (0, 0, 0).
138
        # This goes against the NIFTI spec - it
139
140
141
142
143
        # should just be a straight scaling matrix.
        elif header['qform_code'] == 0 or header['sform_code'] == 0:
            pixdims       = header.get_zooms()
            voxToWorldMat = transform.scaleOffsetXform(pixdims, 0)

144
145
146
        # Otherwise we let nibabel decide
        # which transform to use.
        else:
147
            voxToWorldMat = np.array(header.get_best_affine())
148

149
        return voxToWorldMat
150

151

152
153
154
155
156
    def __determineShape(self, header):
        """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:
157

158
159
160
161
         - 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.
162
163
        """

164
165
166
        origShape = list(header.get_data_shape())
        shape     = list(origShape)
        pixdims   = list(header.get_zooms())
167
168
169
170
171
172
173
174

        # Squeeze out empty dimensions, as
        # 3D image can sometimes be listed
        # as having 4 or more dimensions 
        for i in reversed(range(len(shape))):
            if shape[i] == 1: shape = shape[:i]
            else:             break

175
176
        # But make sure the shape 
        # has at 3 least dimensions
177
178
179
180
181
182
183
        if len(shape) < 3:
            shape = shape + [1] * (3 - len(shape))

        # The same goes for the pixdim - if get_zooms()
        # doesn't return at least 3 values, we'll fall
        # back to the pixdim field in the header.
        if len(pixdims) < 3:
184
            pixdims = header['pixdim'][1:]
185
186
187

        pixdims = pixdims[:len(shape)]
        
188
        return origShape, shape, pixdims
189
190


191
192
    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
193
        underlying ``nibabel.Nifti1Image`` object.
194
195

        See the :meth:`__determineShape` method.
196

197
198
199
        :arg sliceobj: Something that can be used to slice a
                       multi-dimensional array, e.g. ``arr[sliceobj]``.
        """
200

201
202
203
204
205
        # How convenient - nibabel has a function
        # that does the dirty work for us.
        import nibabel.fileslice as fileslice
        return fileslice.canonical_slicers(sliceobj, self.__origShape)
 
206
        
207
    # TODO: Remove this method, and use the shape attribute directly
208
    def is4DImage(self):
209
        """Returns ``True`` if this image is 4D, ``False`` otherwise. """
210
        return len(self.shape) > 3 and self.shape[3] > 1 
211

212
    
213
    def getXFormCode(self, code=None):
214
        """This method returns the code contained in the NIFTI header,
215
        indicating the space to which the (transformed) image is oriented.
216
217

        The ``code`` parameter may be either ``sform`` (the default) or
218
219
220
221
222
223
224
225
        ``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`
226
227
        """

228
229
230
231
        if   code is None:     code = 'sform_code'
        elif code == 'sform' : code = 'sform_code'
        elif code == 'qform' : code = 'qform_code'
        else: raise ValueError('code must be None, sform, or qform')
232

233
        code = self.header[code]
234

235
236
237
238
        # Invalid values
        if   code > 4: code = constants.NIFTI_XFORM_UNKNOWN
        elif code < 0: code = constants.NIFTI_XFORM_UNKNOWN
        
239
        return int(code)
240

241
242
243
244
    # 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.
245
    @memoize.Instanceify(memoize.memoizeMD5)
246
247
248
249
250
251
252
253
254
255
256
257
258
259
    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.
        """

        import nibabel as nib

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

        return nib.orientations.aff2axcodes(xform, inaxes)


260
    @memoize.Instanceify(memoize.memoize)
261
    def isNeurological(self):
262
        """Returns ``True`` if it looks like this ``Nifti`` object is in
263
        neurological orientation, ``False`` otherwise. This test is purely
264
        based on the determinant of the voxel-to-mm transformation matrix -
265
266
267
268
        if it has a positive determinant, the image is assumed to be in
        neurological orientation, otherwise it is assumed to be in
        radiological orientation.

Paul McCarthy's avatar
Paul McCarthy committed
269
270
271
        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
272
273
274
275
276
        """
        import numpy.linalg as npla
        return npla.det(self.voxToWorldMat) > 0


277
278
279
280
281
282
283
284
    def getOrientation(self, axis, xform):
        """Returns a code representing the orientation of the specified data
        axis in the coordinate system defined by the given transformation
        matrix.

        :arg xform: A transformation matrix which is assumed to transform
                    coordinates from the image world coordinate system to
                    some other coordinate system.
285
286
287

        This method returns one of the following values, indicating the
        direction in which coordinates along the specified axis increase:
288
289
290
291
292
293
294
295
        
          - :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
296
297

        The returned value is dictated by the XForm code contained in the
298
299
300
301
302
303
        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).
304
305
        """

306
307
        if self.getXFormCode() == constants.NIFTI_XFORM_UNKNOWN:
            return constants.ORIENT_UNKNOWN 
308
        
309
        import nibabel as nib
310
        code = nib.orientations.aff2axcodes(
311
            xform,
312
313
314
            ((constants.ORIENT_R2L, constants.ORIENT_L2R),
             (constants.ORIENT_A2P, constants.ORIENT_P2A),
             (constants.ORIENT_S2I, constants.ORIENT_I2S)))[axis]
315

316
317
318
        return code 


319
320
class Image(Nifti, notifier.Notifier):
    """Class which represents a 3D/4D NIFTI image. Internally, the image
321
322
    is loaded/stored using a :mod:`nibabel.nifti1.Nifti1Image`, and data
    access managed by a :class:`.ImageWrapper`.
323
324

    
325
    In addition to the attributes added by the :meth:`Nifti.__init__` method,
326
327
    the following attributes/properties are present on an ``Image`` instance 
    as properties (https://docs.python.org/2/library/functions.html#property):
328
329


330
    ============== ===========================================================
331
332
    ``name``       The name of this ``Image`` - defaults to the image
                   file name, sans-suffix.
333

334
335
336
    ``dataSource`` The data source of this ``Image`` - the name of the
                   file from where it was loaded, or some other string
                   describing its origin.
337

338
339
340
341
342
343
    ``nibImage``   A reference to the ``nibabel.Nifti1Image`` object.
    
    ``saveState``  A boolean value which is ``True`` if this image is
                   saved to disk, ``False`` if it is in-memory, or has
                   been edited.
    
344
345
346
347
348
    ``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``
                   is created, or may be incrementally updated as more image 
                   data is loaded from disk.
349
    ============== ===========================================================
350

351
    
352
353
354
355
    The ``Image`` class implements the :class:`.Notifier` interface -
    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):
356

357
    
358
359
    =============== ======================================================
    ``'data'``      This topic is notified whenever the image data changes
360
361
                    (via the :meth:`__setitem__` method).
    
362
363
    ``'saveState'`` This topic is notified whenever the saved state of the
                    image changes (i.e. it is edited, or saved to disk).
364
    
365
366
367
    ``'dataRange'`` This topic is notified whenever the image data range
                    is changed/adjusted.
    =============== ======================================================
368
369
370
    """


371
372
373
374
375
    def __init__(self,
                 image,
                 name=None,
                 header=None,
                 xform=None,
376
                 loadData=True,
377
                 calcRange=True,
378
379
                 indexed=False,
                 threaded=False):
380
381
        """Create an ``Image`` object with the given image data or file name.

382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
        :arg image:     A string containing the name of an image file to load, 
                        or a :mod:`numpy` array, or a :mod:`nibabel` image
                        object.

        :arg name:      A name for the image.

        :arg header:    If not ``None``, assumed to be a
                        :class:`nibabel.nifti1.Nifti1Header` 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. Only used if ``image`` is a ``numpy``
                        array, and ``header`` is ``None``.

        :arg loadData:  If ``True`` (the default) the image data is loaded
399
400
401
402
403
404
405
                        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. 

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

        :arg indexed:   If ``True``, and the file is gzipped, it is opened 
                        using the :mod:`indexed_gzip` package. Otherwise the
414
415
                        file is opened by ``nibabel``. Ignored if ``loadData``
                        is ``True``.
416
417
418

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

        import nibabel as nib

424
425
        nibImage   = None
        dataSource = None
426
        fileobj    = None
427

428
429
430
431
        if loadData:
            indexed  = False
            threaded = False

432
433
434
        # The image parameter may be the name of an image file
        if isinstance(image, six.string_types):

435
436
437
438
439
440
441
442
443
444
445
            image = op.abspath(addExt(image))

            # Use indexed_gzip to open gzip files
            # if requested - this provides fast
            # on-disk access to the compressed
            # data.
            #
            # If using indexed_gzip, we store a
            # ref to the file object - we'll close
            # it when we are destroyed.
            if indexed and image.endswith('.gz'):
446
                nibImage, fileobj = loadIndexedImageFile(image)
447
448
449
450
451
452

            # Otherwise we let nibabel
            # manage the file reference(s)
            else:
                nibImage  = nib.load(image)
                
453
            dataSource = image
454
455
456
457
458
459
460
461
 
        # 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):

            if   header is not None: xform = header.get_best_affine()
            elif xform  is     None: xform = np.identity(4)
462
                    
463
464
465
            nibImage = nib.nifti1.Nifti1Image(image,
                                              xform,
                                              header=header)
466
467
468
            
        # otherwise, we assume that it is a nibabel image
        else:
469
            nibImage = image
470

471
472
473
        # Figure out the name of this image, if 
        # it has not beenbeen explicitly passed in
        if name is None:
474
            
475
476
477
478
            # If this image was loaded
            # from disk, use the file name.
            if isinstance(image, six.string_types):
                name = removeExt(op.basename(image))
479
            
480
481
482
483
484
485
486
            # Or the image was created from a numpy array
            elif isinstance(image, np.ndarray):
                name = 'Numpy array'
            
            # Or image from a nibabel image
            else:
                name = 'Nibabel image'
487
 
488
        Nifti.__init__(self, nibImage.get_header())
489

490
491
        self.name                = name
        self.__dataSource        = dataSource
492
        self.__fileobj           = fileobj
493
494
495
496
497
        self.__nibImage          = nibImage
        self.__saveState         = dataSource is not None
        self.__suppressDataRange = False
        self.__imageWrapper      = imagewrapper.ImageWrapper(self.nibImage,
                                                             self.name,
498
499
                                                             loadData=loadData,
                                                             threaded=threaded)
500

501
502
503
504
505
506
        if calcRange:
            self.calcRange()

        self.__imageWrapper.register(
            '{}_{}'.format(id(self), self.name),
            self.__dataRangeChanged) 
507

508
        
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
    def __hash__(self):
        """Returns a number which uniquely idenfities this ``Image`` instance
        (the result of ``id(self)``).
        """
        return id(self)


    def __str__(self):
        """Return a string representation of this ``Image`` instance."""
        return '{}({}, {})'.format(self.__class__.__name__,
                                   self.name,
                                   self.dataSource)

        
    def __repr__(self):
        """See the :meth:`__str__` method."""
        return self.__str__()

527
528

    def __del__(self):
529
530
531
532
533
        """Closes any open file handles, and clears some references. """
        
        self.__nibImage     = None
        self.__imageWrapper = None
        
534
535
536
        if self.__fileobj is not None:
            self.__fileobj.close()
        
537
538
539
    
    @property
    def dataSource(self):
540
541
        """Returns the data source (e.g. file name) that this ``Image`` was
        loaded from (``None`` if this image only exists in memory).
542
543
        """
        return self.__dataSource
544

545
546
547
    
    @property
    def nibImage(self):
548
        """Returns a reference to the ``nibabel.nifti1.Nifti1Image`` instance.
549
550
        """
        return self.__nibImage
551

552
553
554
555
556
557
558
    
    @property
    def saveState(self):
        """Returns ``True`` if this ``Image`` has been saved to disk, ``False``
        otherwise.
        """
        return self.__saveState
559

560
561
562
    
    @property
    def dataRange(self):
563
564
565
566
567
568
569
        """Returns the image data range as a  ``(min, max)`` tuple. If the
        ``calcRange`` parameter to :meth:`__init__` was ``False``, these
        values may not be accurate, and may change as more image data is
        accessed.

        If the data range has not been no data has been accessed,
        ``(None, None)`` is returned.
570
        """
571
572
573
574
        if self.__imageWrapper is None: drange = (None, None)
        else:                           drange = self.__imageWrapper.dataRange

        # Fall back to the cal_min/max
575
        # fields in the NIFTI header
576
577
578
        # if we don't yet know anything
        # about the image data range.
        if drange[0] is None or drange[1] is None:
579
580
            drange = (float(self.header['cal_min']),
                      float(self.header['cal_max']))
581
582
583
584
585
586
587
588
589
590
591
592

        return drange

    
    @property
    def dtype(self):
        """Returns the ``numpy`` data type of the image data. """
        
        # Get the data type from the
        # first voxel in the image
        coords = [0] * len(self.__nibImage.shape)
        return self.__nibImage.dataobj[tuple(coords)].dtype
593

594

595
    def __dataRangeChanged(self, *args, **kwargs):
596
597
598
599
        """Called when the :class:`.ImageWrapper` data range changes.
        Notifies any listeners of this ``Image`` (registered through the
        :class:`.Notifier` interface) on the ``'dataRange'`` topic.
        """
600
601
        if not self.__suppressDataRange:
            self.notify(notifier_topic='dataRange')
602

603

604
    def calcRange(self, sizethres=None):
605
606
607
        """Forces calculation of the image data range.

        :arg sizethres: If not ``None``, specifies an image size threshold
608
609
610
611
                        (total number of bytes). If the number of bytes in
                        the image is greater than this threshold, the range 
                        is calculated on a sample (the first volume for a
                        4D image, or slice for a 3D image).
612
        """
613

614
615
616
617
        # The ImageWrapper automatically calculates
        # the range of the specified slice, whenever
        # it gets indexed. All we have to do is
        # access a portion of the data to trigger the
618
619
        # range calculation.
        nbytes = np.prod(self.shape) * self.dtype.itemsize
620
621
622

        # If an image size threshold has not been specified,
        # then we'll calculate the full data range right now.
623
624
625
        if sizethres is None or nbytes < sizethres:
            log.debug('{}: Forcing calculation of full '
                      'data range'.format(self.name))
626
            self.__imageWrapper[:]
627
628
629
630
            
        else:
            log.debug('{}: Calculating data range '
                      'from sample'.format(self.name))
631

632
633
634
            # Otherwise if the number of values in the
            # image is bigger than the size threshold, 
            # we'll calculate the range from a sample:
635
636
            if len(self.shape) == 3: self.__imageWrapper[:, :, 0]
            else:                    self.__imageWrapper[:, :, :, 0]
637
638


639
640
641
642
643
644
645
    def loadData(self):
        """Makes sure that the image data is loaded into memory.
        See :meth:`.ImageWrapper.loadData`.
        """
        self.__imageWrapper.loadData()


646
    def save(self, filename=None):
Paul McCarthy's avatar
Paul McCarthy committed
647
648
649
650
        """Saves this ``Image`` to the specifed file, or the :attr:`dataSource`
        if ``filename`` is ``None``.
        """

651
652
        import nibabel as nib
        
Paul McCarthy's avatar
Paul McCarthy committed
653
654
655
656
657
658
        if self.__dataSource is None and filename is None:
            raise ValueError('A file name must be specified')

        if filename is None:
            filename = self.__dataSource

659
660
        filename = op.abspath(filename)

Paul McCarthy's avatar
Paul McCarthy committed
661
        log.debug('Saving {} to {}'.format(self.name, filename))
662

663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
        # If this Image is not managing its
        # own file object, nibabel does all
        # of the hard work.
        if self.__fileobj is None:
            nib.save(self.__nibImage, filename)

        # Otherwise we've got our own file
        # handle to an IndexedGzipFile
        else:
            # Currently indexed_gzip does not support
            # writing. So we're going to use nibabel
            # to save the image, then close and re-open
            # the file.
            #
            # Unfortunately this means that we'll
            # lose the file index (and fast random
            # access) - I'll fix this when I get a
            # chance to work on indexed_gzip a bit
            # more.
            #
            # Hopefully I should be able to add write
            # support to indexed_gzip, such that it
            # re-builds the index while writing the
            # compressed data. And then be able to
            # transfer the index generated from the
            # write to a new read-only file handle.
            nib.save(self.__nibImage, filename)
            self.__fileobj.close()
            self.__nibImage, self.__fileobj = loadIndexedImageFile(filename)

        self.__dataSource = filename
        self.__saveState  = True
        
        self.notify(notifier_topic='saveState')
Paul McCarthy's avatar
Paul McCarthy committed
697

698

699
    def __getitem__(self, sliceobj):
700
701
702
        """Access the image data with the specified ``sliceobj``.

        :arg sliceobj: Something which can slice the image data.
703
        """
704
705
706
707


        log.debug('{}: __getitem__ [{}]'.format(self.name, sliceobj))
        
708
709
710
711
712
713
714
715
716
        data = self.__imageWrapper.__getitem__(sliceobj)

        if len(data.shape) > len(self.shape):

            shape = data.shape[:len(self.shape)]
            data  = np.reshape(data, shape)

        return data

717

718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
    def __setitem__(self, sliceobj, values):
        """Set the image data at ``sliceobj`` to ``values``.

        :arg sliceobj: Something which can slice the image data.
        :arg values:   New image data.

        .. note:: Modifying image data will force the entire image to be 
                  loaded into memory if it has not already been loaded.
        """

        log.debug('{}: __setitem__ [{} = {}]'.format(self.name,
                                                     sliceobj,
                                                     values.shape))

        self.__suppressDataRange = True
        oldRange = self.__imageWrapper.dataRange

        self.__imageWrapper.__setitem__(sliceobj, values)

        newRange = self.__imageWrapper.dataRange
        self.__suppressDataRange = False

        if values.size > 0:
741

742
            self.notify(notifier_topic='data')
743
744
745
746

            if self.__saveState:
                self.__saveState = False
                self.notify(notifier_topic='saveState')
747
748
749
750
751

            if not np.all(np.isclose(oldRange, newRange)):
                self.notify(notifier_topic='dataRange') 


752
ALLOWED_EXTENSIONS = ['.nii.gz', '.nii', '.img', '.hdr', '.img.gz', '.hdr.gz']
753
"""The file extensions which we understand. This list is used as the default
754
if the ``allowedExts`` parameter is not passed to any of the functions
755
below.
756
757
"""

758

759
760
EXTENSION_DESCRIPTIONS = ['Compressed NIFTI images',
                          'NIFTI images',
761
                          'ANALYZE75 images',
762
763
                          'NIFTI/ANALYZE75 headers',
                          'Compressed NIFTI/ANALYZE75 images',
764
765
766
767
                          'Compressed images']
"""Descriptions for each of the extensions in :data:`ALLOWED_EXTENSIONS`. """


768
REPLACEMENTS = {'.hdr' : ['.img'], '.hdr.gz' : ['.img.gz']}
769
770
771
772
773
"""Suffix replacements used by :func:`addExt` to resolve file path
ambiguities - see :func:`fsl.utils.path.addExt`.
"""


774
775
776
777
DEFAULT_EXTENSION  = '.nii.gz'
"""The default file extension (TODO read this from ``$FSLOUTPUTTYPE``)."""


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


784
785
def looksLikeImage(filename, allowedExts=None):
    """Returns ``True`` if the given file looks like an image, ``False``
786
787
    otherwise.

788
789
790
791
    .. note:: The ``filename`` cannot just be a file prefix - it must
              include the file suffix (e.g. ``myfile.nii.gz``, not
              ``myfile``).

792
793
794
    :arg filename:    The file name to test.
    
    :arg allowedExts: A list of strings containing the allowed file
795
                      extensions - defaults to :attr:`ALLOWED_EXTENSIONS`.
796
797
798
799
    """

    if allowedExts is None: allowedExts = ALLOWED_EXTENSIONS

800
801
802
    # TODO A much more robust approach would be
    #      to try loading the file using nibabel.

Paul McCarthy's avatar
Paul McCarthy committed
803
    return any([filename.endswith(ext) for ext in allowedExts])
804
805


806
def removeExt(filename):
807
    """Removes the extension from the given file name. Returns the filename
808
809
    unmodified if it does not have a supported extension.

810
    See :func:`~fsl.utils.path.removeExt`.
811

812
813
814
    :arg filename: The file name to strip.
    """
    return fslpath.removeExt(filename, ALLOWED_EXTENSIONS)
815
816


817
def addExt(prefix, mustExist=True):
818
819
    """Adds a file extension to the given file ``prefix``.

820
    See :func:`~fsl.utils.path.addExt`.
821
    """
822
823
824
    return fslpath.addExt(prefix,
                          ALLOWED_EXTENSIONS,
                          mustExist,
825
826
                          DEFAULT_EXTENSION,
                          replace=REPLACEMENTS)
827
828
829
830
831


def loadIndexedImageFile(filename):
    """Loads the given image file using ``nibabel`` and ``indexed_gzip``.

832
    Returns a tuple containing the ``nibabel.Nifti1Image``, and the open
833
834
835
836
837
838
839
840
    ``IndexedGzipFile`` handle.
    """
    
    import nibabel      as nib
    import indexed_gzip as igzip

    log.debug('Loading {} using indexed gzip'.format(filename))

841
    fobj = igzip.SafeIndexedGzipFile(
842
843
844
845
846
847
848
849
850
        filename=filename,
        spacing=4194304,
        readbuf_size=131072)

    fmap = nib.Nifti1Image.make_file_map()
    fmap['image'].fileobj = fobj
    image = nib.Nifti1Image.from_file_map(fmap)

    return image, fobj