imagewrapper.py 33.9 KB
Newer Older
1
2
3
4
5
6
#!/usr/bin/env python
#
# imagewrapper.py - The ImageWrapper class.
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
7
8
"""This module provides the :class:`ImageWrapper` class, which can be used
to manage data access to ``nibabel`` NIFTI images.
9

10
11
12
13
14

Terminology
-----------


15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
There are some confusing terms used in this module, so it may be of use to
get their definitions straight:

  - *Coverage*:  The portion of an image that has been covered in the data
                 range calculation. The ``ImageWrapper`` keeps track of
                 the coverage for individual volumes within a 4D image (or
                 slices in a 3D image).

  - *Slice*:     Portion of the image data which is being accessed. A slice
                 comprises either a tuple of ``slice`` objects (or integers),
                 or a sequence of ``(low, high)`` tuples, specifying the
                 index range into each image dimension that is covered by
                 the slice.

  - *Expansion*: A sequence of ``(low, high)`` tuples, specifying an
                 index range into each image dimension, that is used to
                 *expand* the *coverage* of an image, based on a given set of
                 *slices*.
33
34
35
"""


36
37
import logging
import collections
38
import itertools as it
39

40
41
import numpy     as np
import nibabel   as nib
42
43

import fsl.utils.notifier as notifier
44
import fsl.utils.async    as async
45
46
47
48
49
50


log = logging.getLogger(__name__)


class ImageWrapper(notifier.Notifier):
51
52
53
    """The ``ImageWrapper`` class is a convenience class which manages data
    access to ``nibabel`` NIFTI images. The ``ImageWrapper`` class can be
    used to:
54

55
56
57
    
      - Control whether the image is loaded into memory, or kept on disk
    
58
      - Incrementally update the known image data range, as more image
59
60
61
        data is read in.


62
63
64
65
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
91
92
93
94
95
96
97
98
99
100
101
102
103
    *In memory or on disk?*

    The image data will be kept on disk, and accessed through the
    ``nibabel.Nifti1Image.dataobj`` array proxy, if:

     - The ``loadData`` parameter to :meth:`__init__` is ``False``.
     - The :meth:`loadData` method never gets called.
     - The image data is not modified (via :meth:`__setitem__`.

    If any of these conditions do not hold, the image data will be loaded into
    memory and accessed directly, via the ``nibabel.Nifti1Image.get_data``
    method.


    *Image dimensionality*

    
    The ``ImageWrapper`` abstracts away trailing image dimensions of length 1.
    This means that if the header for a NIFTI image specifies that the image
    has four dimensions, but the fourth dimension is of length 1, you do not
    need to worry about indexing that fourth dimension.


    *Data range*

    
    In order to avoid the computational overhead of calculating the image data
    range (its minimum/maximum values) when an image is first loaded in, an
    ``ImageWrapper`` incrementally updates the known image data range as data
    is accessed. The ``ImageWrapper`` keeps track of the image data _coverage_,
    the portion of the image which has already been considered in the data
    range calculation. When data from a region of the image not in the coverage
    is accessed, the coverage is expanded to include this region. The coverage
    is always expanded in a rectilinear manner, i.e. the coverage is always
    rectangular for a 2D image, or cuboid for a 3D image.

    
    For a 4D image, the ``ImageWrapper`` internally maintains a separate
    coverage and known data range for each 3D volume within the image. For a 3D
    image, separate coverages and data ranges are stored for each 2D slice.


104
105
106
107
    The ``ImageWrapper`` implements the :class:`.Notifier` interface.
    Listeners can register to be notified whenever the known image data range
    is updated. The data range can be accessed via the :attr:`dataRange`
    property.
108

109

110
111
112
113
114
115
116
117
118
119
    The ``ImageWrapper`` class uses the following functions (also defined in 
    this module) to keep track of the portion of the image that has currently
    been included in the data range calculation:

    .. autosummary::
       :nosignatures:

       sliceObjToSliceTuple
       sliceTupleToSliceObj
       sliceCovered
120
       calcExpansion
121
       adjustCoverage
122
123
    """

124
    
125
126
127
128
129
130
    def __init__(self,
                 image,
                 name=None,
                 loadData=False,
                 dataRange=None,
                 threaded=False):
131
        """Create an ``ImageWrapper``.
132

133
        :arg image:     A ``nibabel.Nifti1Image``.
134

135
136
        :arg name:      A name for this ``ImageWrapper``, solely used for 
                        debug log messages.
137

138
139
140
141
142
143
        :arg loadData:  If ``True``, the image data is loaded into memory.
                        Otherwise it is kept on disk (and data access is
                        performed through the ``nibabel.Nifti1Image.dataobj``
                        array proxy).

        :arg dataRange: A tuple containing the initial ``(min, max)``  data
144
145
                        range to use. See the :meth:`reset` method for
                        important information about this parameter.
146
147
148
149

        :arg threaded:  If ``True``, the data range is updated on a
                        :class:`.TaskThread`. Otherwise (the default), the
                        data range is updated directly on reads/writes.
150
        """
151

152
153
154
        self.__image      = image
        self.__name       = name
        self.__taskThread = None
155

156
157
158
159
160
161
162
163
164
165
166
167
        # Save the number of 'real' dimensions,
        # that is the number of dimensions minus
        # any trailing dimensions of length 1
        self.__numRealDims = len(image.shape)
        for d in reversed(image.shape):
            if d == 1: self.__numRealDims -= 1
            else:      break

        # And save the number of
        # 'padding' dimensions too.
        self.__numPadDims = len(image.shape) - self.__numRealDims

168
169
170
171
172
173
174
175
176
177
178
179
180
        # The internal state is stored
        # in these attributes - they're
        # initialised in the reset method.
        self.__range     = None
        self.__coverage  = None
        self.__volRanges = None
        self.__covered   = False

        self.reset(dataRange)

        if loadData:
            self.loadData()

181
        if threaded:
182
            self.__taskThread = async.TaskThread()
183
            self.__taskThread.daemon = True
184
185
186
187
188
189
190
191
192
193
194
195
            self.__taskThread.start()


    def __del__(self):
        """If this ``ImageWrapper`` was created with ``threaded=True``,
        the :class:`.TaskThread` is stopped.
        """
        self.__image = None
        if self.__taskThread is not None:
            self.__taskThread.stop()
            self.__taskThraed = None

196
197
198
199
200
201
202
203
204
205
206
207

    def reset(self, dataRange=None):
        """Reset the internal state and known data range of this
        ``ImageWrapper``.

        
        :arg dataRange: A tuple containing the initial ``(min, max)``  data
                        range to use. See the :meth:`reset` method.


        .. note:: The ``dataRange`` parameter is intended for situations where
                  the image data range is known (e.g. it was calculated
208
                  earlier, and the image is being re-loaded). If a
209
210
211
212
213
214
215
216
                  ``dataRange`` is passed in, it will *not* be overwritten by
                  any range calculated from the data, unless the calculated
                  data range is wider than the provided ``dataRange``. 
        """
        
        if dataRange is None:
            dataRange = None, None

217
        image =             self.__image
218
219
220
        ndims =             self.__numRealDims - 1
        nvols = image.shape[self.__numRealDims - 1]

221
222
        # The current known image data range. This
        # gets updated as more image data gets read.
223
        self.__range = dataRange
224

225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
        # The coverage array is used to keep track of
        # the portions of the image which have been
        # considered in the data range calculation.
        # We use this coverage to avoid unnecessarily
        # re-calculating the data range on the same
        # part of the image.
        #
        # First of all, we're going to store a separate
        # 'coverage' for each 2D slice in the 3D image
        # (or 3D volume for 4D images). This effectively
        # means a seaprate coverage for each index in the
        # last 'real' image dimension (see above).
        # 
        # For each slice/volume, the the coverage is
        # stored as sequences of (low, high) indices, one
        # for each dimension in the slice/volume (e.g.
        # row/column for a slice, or row/column/depth
        # for a volume).
        #
244
        # All of these indices are stored in a numpy array:
245
246
        #   - first dimension:  low/high index
        #   - second dimension: image dimension
247
        #   - third dimension:  slice/volume index 
248
        self.__coverage = np.zeros((2, ndims, nvols), dtype=np.float32)
249

250
251
252
253
254
255
        # Internally, we calculate and store the
        # data range for each volume/slice/vector
        self.__volRanges = np.zeros((nvols, 2), dtype=np.float32)

        self.__coverage[ :] = np.nan
        self.__volRanges[:] = np.nan
256

257
258
259
260
261
        # This flag is set to true if/when the
        # full image data range becomes known
        # (i.e. when all data has been loaded in).
        self.__covered = False

262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
        
    @property
    def dataRange(self):
        """Returns the currently known data range as a tuple of ``(min, max)``
        values.
        """
        # If no image data has been accessed, we
        # default to whatever is stored in the
        # header (which may or may not contain
        # useful values).
        low, high = self.__range
        hdr       = self.__image.get_header()

        if low  is None: low  = float(hdr['cal_min'])
        if high is None: high = float(hdr['cal_max'])

        return low, high

    
    @property
    def covered(self):
        """Returns ``True`` if this ``ImageWrapper`` has read the entire
        image data, ``False`` otherwise.
        """
        return self.__covered


    def coverage(self, vol):
        """Returns the current image data coverage for the specified volume
        (for a 4D image, slice for a 3D image, or vector for a 2D images).

        :arg vol: Index of the volume/slice/vector to return the coverage
                  for.

        :returns: The coverage for the specified volume, as a ``numpy``
                  array of shape ``(nd, 2)``, where ``nd`` is the number
                  of dimensions in the volume.
299
300
301

        .. note:: If the specified volume is not covered, the returned array
                  will contain ``np.nan`` values.
302
        """
303
        return np.array(self.__coverage[..., vol])
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320

    
    def loadData(self):
        """Forces all of the image data to be loaded into memory.

        .. note:: This method will be called by :meth:`__init__` if its
                  ``loadData`` parameter is ``True``.
        """

        # If the data is not already loaded, this will
        # cause nibabel to load it. By default, nibabel
        # will cache the numpy array that contains the
        # image data, so subsequent calls to this
        # method will not overwrite any changes that
        # have been made to the data array.
        self.__image.get_data()

321

322
    def __getData(self, sliceobj, isTuple=False):
323
        """Retrieves the image data at the location specified by ``sliceobj``.
324
325
326
327
328
329

        :arg sliceobj: Something which can be used to slice an array, or
                       a sequence of (low, high) index pairs.

        :arg isTuple:  Set to ``True`` if ``sliceobj`` is a sequence of
                       (low, high) index pairs.
330
331
        """

332
        if isTuple:
333
            sliceobj = sliceTupleToSliceObj(sliceobj)
334

335
        # If the image has not been loaded
336
        # into memory, we can use the nibabel
337
338
339
340
341
342
343
344
345
346
        # ArrayProxy. Otheriwse if it is in
        # memory, we can access it directly.
        #
        # Furthermore, if it is in memory and
        # has been modified, the ArrayProxy
        # will give us out-of-date values (as
        # the ArrayProxy reads from disk). So
        # we have to read from the in-memory
        # array to get changed values.
        if self.__image.in_memory: return self.__image.get_data()[sliceobj]
347
        else:                      return self.__image.dataobj[   sliceobj]
348
349


350
351
352
353
354
355
356
    def __imageIsCovered(self):
        """Returns ``True`` if all portions of the image have been covered
        in the data range calculation, ``False`` otherwise.
        """
        
        shape  = self.__image.shape
        slices = zip([0] * len(shape), shape)
357
        return sliceCovered(slices, self.__coverage)
358

359

360
361
362
    def __expandCoverage(self, slices):
        """Expands the current image data range and coverage to encompass the
        given ``slices``.
363
        """
364

365
366
        _, expansions = calcExpansion(slices, self.__coverage)
        expansions    = collapseExpansions(expansions, self.__numRealDims - 1)
367
368
369
370
371

        log.debug('Updating image {} data range [slice: {}] '
                  '(current range: [{}, {}]; '
                  'number of expansions: {}; '
                  'current coverage: {})'.format(
372
                      self.__name,
373
                      slices,
374
375
                      self.__range[0],
                      self.__range[1],
376
                      len(expansions),
377
                      self.__coverage))
378
379
380
381
382
383
384
385
386

        # As we access the data for each expansions,
        # we want it to have the same dimensionality
        # as the full image, so we can access data
        # for each volume in the image separately.
        # So we squeeze out the padding dimensions,
        # but not the volume dimension.
        squeezeDims = tuple(range(self.__numRealDims,
                                  self.__numRealDims + self.__numPadDims))
387
        
388
389
390
391
392
        # The calcExpansion function splits up the
        # expansions on volumes - here we calculate
        # the min/max per volume/expansion, and
        # iteratively update the stored per-volume
        # coverage and data range.
393
394
395
396
397
398
399
400
        for i, exp in enumerate(expansions):

            data = self.__getData(exp, isTuple=True)
            data = data.squeeze(squeezeDims)

            vlo, vhi = exp[self.__numRealDims - 1]

            for vi, vol in enumerate(range(vlo, vhi)):
401

402
403
                oldvlo, oldvhi = self.__volRanges[vol, :]
                voldata        = data[..., vi]
404

405
406
                newvlo = float(np.nanmin(voldata))
                newvhi = float(np.nanmax(voldata))
407

408
409
                if (not np.isnan(oldvlo)) and oldvlo < newvlo: newvlo = oldvlo
                if (not np.isnan(oldvhi)) and oldvhi > newvhi: newvhi = oldvhi
410

411
412
413
414
415
                # Update the stored range and
                # coverage for each volume 
                self.__volRanges[vol, :]  = newvlo, newvhi
                self.__coverage[..., vol] = adjustCoverage(
                    self.__coverage[..., vol], exp)
416

417
418
419
        # Calculate the new known data
        # range over the entire image
        # (i.e. over all volumes).
420
421
        newmin = float(np.nanmin(self.__volRanges[:, 0]))
        newmax = float(np.nanmax(self.__volRanges[:, 1]))
422

423
        oldmin, oldmax = self.__range
424
        self.__range   = (newmin, newmax)
425
426
        self.__covered = self.__imageIsCovered()

427
428
        if any((oldmin is None, oldmax is None)) or \
           not np.all(np.isclose([oldmin, oldmax], [newmin, newmax])):
429
430
431
432
433
434
            log.debug('Image {} range changed: [{}, {}] -> [{}, {}]'.format(
                self.__name,
                oldmin,
                oldmax,
                newmin,
                newmax))
435
436
            self.notify()

437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454

    def __updateDataRangeOnRead(self, slices, data):
        """Called by :meth:`__getitem__`. Calculates the minimum/maximum
        values of the given data (which has been extracted from the portion of
        the image specified by ``slices``), and updates the known data range
        of the image.

        :arg slices: A tuple of tuples, each tuple being a ``(low, high)``
                     index pair, one for each dimension in the image. 
        
        :arg data:   The image data at the given ``slices`` (as a ``numpy``
                     array).
        """

        # TODO You could do something with
        #      the provided data to avoid
        #      reading it in again.

455
456
457
458
459
        if self.__taskThread is None:
            self.__expandCoverage(slices)
        else:
            name = '{}_read_{}'.format(id(self), slices)
            if not self.__taskThread.isQueued(name):
460
                self.__taskThread.enqueue(
461
                    self.__expandCoverage, slices, taskName=name)
462

463

464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
    def __updateDataRangeOnWrite(self, slices, data):
        """Called by :meth:`__setitem__`. Assumes that the image data has
        been changed (the data at ``slices`` has been replaced with ``data``.
        Updates the image data coverage, and known data range accordingly.

        :arg slices: A tuple of tuples, each tuple being a ``(low, high)``
                     index pair, one for each dimension in the image. 
        
        :arg data:   The image data at the given ``slices`` (as a ``numpy``
                     array). 
        """

        overlap = sliceOverlap(slices, self.__coverage)

        # If there's no overlap between the written
        # area and the current coverage, then it's
        # easy - we just expand the coverage to
        # include the newly written area.
        if overlap in (OVERLAP_SOME, OVERLAP_ALL):

            # If there is overlap between the written
            # area and the current coverage, things are
            # more complicated, because the portion of
            # the image that has been written over may
            # have contained the currently known data
            # minimum/maximum. We have no way of knowing
            # this, so we have to reset the coverage (on
            # the affected volumes), and recalculate the
            # data range.

            # TODO Could you store the location of the
            #      data minimum/maximum (in each volume),
            #      so you know whether resetting the
            #      coverage is necessary?
            lowVol, highVol = slices[self.__numRealDims - 1]

500
501
502
503
504
505
506
            log.debug('Image {} data written - clearing known data '
                      'range on volumes {} - {} (write slice: {})'.format(
                          self.__name,
                          lowVol,
                          highVol,
                          slices))

507
            for vol in range(lowVol, highVol):
508
509
                self.__coverage[:, :, vol]    = np.nan
                self.__volRanges[     vol, :] = np.nan
510

511
512
513
514
515
516

        if self.__taskThread is None:
            self.__expandCoverage(slices)
        else:
            name = '{}_write_{}'.format(id(self), slices)
            if not self.__taskThread.isQueued(name):
517
                self.__taskThread.enqueue(
518
                    self.__expandCoverage, slices, taskName=name)
519

520
            
521
    def __getitem__(self, sliceobj):
522
523
524
525
        """Returns the image data for the given ``sliceobj``, and updates
        the known image data range if necessary.

        :arg sliceobj: Something which can slice the image data.
526
527
        """

528
529
        log.debug('Getting image data: {}'.format(sliceobj))

530
531
532
        sliceobj = nib.fileslice.canonical_slicers(
            sliceobj, self.__image.shape)

533
        # TODO Cache 3D images for large 4D volumes, 
534
        #      so you don't have to hit the disk?
535

536
537
        data = self.__getData(sliceobj)

538
        if not self.__covered:
539

540
            slices = sliceObjToSliceTuple(sliceobj, self.__image.shape)
541

542
            if not sliceCovered(slices, self.__coverage):
543
                self.__updateDataRangeOnRead(slices, data)
544
545

        return data
546
547


548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
    def __setitem__(self, sliceobj, values):
        """Writes the given ``values`` to the image at the given ``sliceobj``.

        
        :arg sliceobj: Something which can be used to slice the array.
        :arg values:   Data to write to the image.

        
        .. note:: Modifying image data will cause the entire image to be
                  loaded into memory. 
        """

        sliceobj = nib.fileslice.canonical_slicers(sliceobj,
                                                   self.__image.shape)
        slices   = sliceObjToSliceTuple(           sliceobj,
                                                   self.__image.shape)

        # The image data has to be in memory
        # for the data to be changed. If it's
        # already in memory, this call won't
        # have any effect.
        self.loadData()

        self.__image.get_data()[sliceobj] = values
        self.__updateDataRangeOnWrite(slices, values)


575
576
577
def sliceObjToSliceTuple(sliceobj, shape):
    """Turns an array slice object into a tuple of (low, high) index
    pairs, one pair for each dimension in the given shape
578
579
580
581
582

    :arg sliceobj: Something which can be used to slice an array of shape
                   ``shape``.

    :arg shape:    Shape of the array being sliced.
583
584
585
586
    """

    indices = []

587
588
    # The sliceobj could be a single sliceobj
    # or integer, instead of a tuple
589
590
591
    if not isinstance(sliceobj, collections.Sequence):
        sliceobj = [sliceobj]

592
    # Turn e.g. array[6] into array[6, :, :]
593
594
595
596
    if len(sliceobj) != len(shape):
        missing  = len(shape) - len(sliceobj)
        sliceobj = list(sliceobj) + [slice(None) for i in range(missing)]

597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
    for dim, s in enumerate(sliceobj):

        # each element in the slices tuple should 
        # be a slice object or an integer
        if isinstance(s, slice): i = [s.start, s.stop]
        else:                    i = [s,       s + 1]

        if i[0] is None: i[0] = 0
        if i[1] is None: i[1] = shape[dim]

        indices.append(tuple(i))

    return tuple(indices)


612
613
def sliceTupleToSliceObj(slices):
    """Turns a sequence of (low, high) index pairs into a tuple of array
614
    ``slice`` objects.
615
616

    :arg slices: A sequence of (low, high) index pairs.
617
618
619
620
621
622
623
624
625
626
    """

    sliceobj = []

    for lo, hi in slices:
        sliceobj.append(slice(lo, hi, 1))

    return tuple(sliceobj)


627
628
629
630
def adjustCoverage(oldCoverage, slices): 
    """Adjusts/expands the given ``oldCoverage`` so that it covers the
    given set of ``slices``.

631
    :arg oldCoverage: A ``numpy`` array of shape ``(2, n)`` containing
632
633
                      the (low, high) index pairs for ``n`` dimensions of
                      a single slice/volume in the image.
634
    
635
636
637
638
    :arg slices:      A sequence of (low, high) index pairs. If ``slices``
                      contains more dimensions than are specified in
                      ``oldCoverage``, the trailing dimensions are ignored.

639
    :return: A ``numpy`` array containing the adjusted/expanded coverage.
640
641
    """

642
    newCoverage = np.zeros(oldCoverage.shape, dtype=np.uint32)
643

644
    for dim in range(oldCoverage.shape[1]):
645

646
647
        low,      high      = slices[        dim]
        lowCover, highCover = oldCoverage[:, dim]
648

649
650
        if np.isnan(lowCover)  or low  < lowCover:  lowCover  = low
        if np.isnan(highCover) or high > highCover: highCover = high
651

652
        newCoverage[:, dim] = lowCover, highCover
653
654
655

    return newCoverage

656

657
658
659
660
661
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
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
OVERLAP_ALL = 0
"""Indicates that the slice is wholly contained within the coverage.  This is
a return code for the :func:`sliceOverlap` function.
"""


OVERLAP_SOME = 1
"""Indicates that the slice partially overlaps with the coverage. This is a
return code for the :func:`sliceOverlap` function.
"""


OVERLAP_NONE = 2
"""Indicates that the slice does not overlap with the coverage. This is a
return code for the :func:`sliceOverlap` function.
"""


def sliceOverlap(slices, coverage):
    """Determines whether the given ``slices`` overlap with the given
    ``coverage``.

    :arg slices:    A sequence of (low, high) index pairs, assumed to cover
                    all image dimensions.
    :arg coverage:  A ``numpy`` array of shape ``(2, nd, nv)`` (where ``nd``
                    is the number of dimensions being covered, and ``nv`` is
                    the number of volumes (or vectors/slices) in the image,
                    which contains the (low, high) index pairs describing
                    the current image coverage.

    :returns: One of the following codes:
              .. autosummary::

              OVERLAP_ALL
              OVERLAP_SOME
              OVERLAP_NONE
    """

    numDims         = coverage.shape[1]
    lowVol, highVol = slices[numDims]

    # Overlap state is calculated for each volume
    overlapStates = np.zeros(highVol - lowVol)

    for i, vol in enumerate(range(lowVol, highVol)):

        state = OVERLAP_ALL

        for dim in range(numDims):

            lowCover, highCover = coverage[:, dim, vol]
            lowSlice, highSlice = slices[     dim] 

            # No coverage
            if np.isnan(lowCover) or np.isnan(highCover):
                state = OVERLAP_NONE
                break

            # The slice is contained within the
            # coverage on this dimension - check
            # the other dimensions.
            if lowSlice >= lowCover and highSlice <= highCover:
                continue

            # The slice does not overlap at all
            # with the coverage on this dimension
            # (or at all). No overlap - no need
            # to check the other dimensions.
            if lowSlice >= highCover or highSlice <= lowCover:
                state = OVERLAP_NONE
                break

            # There is some overlap between the
            # slice and coverage on this dimension
            # - check the other dimensions.
            state = OVERLAP_SOME
            
        overlapStates[i] = state

    if   np.any(overlapStates == OVERLAP_SOME): return OVERLAP_SOME
    elif np.all(overlapStates == OVERLAP_NONE): return OVERLAP_NONE
    elif np.all(overlapStates == OVERLAP_ALL):  return OVERLAP_ALL


def sliceCovered(slices, coverage):
742
743
    """Returns ``True`` if the portion of the image data calculated by
    the given ``slices` has already been calculated, ``False`` otherwise.
744

745
746
747
748
749
750
751
    :arg slices:    A sequence of (low, high) index pairs, assumed to cover
                    all image dimensions.
    :arg coverage:  A ``numpy`` array of shape ``(2, nd, nv)`` (where ``nd``
                    is the number of dimensions being covered, and ``nv`` is
                    the number of volumes (or vectors/slices) in the image,
                    which contains the (low, high) index pairs describing
                    the current image coverage.
752
753
    """

754
755
    numDims         = coverage.shape[1]
    lowVol, highVol = slices[numDims]
756
757
758

    for vol in range(lowVol, highVol):

759
        for dim in range(numDims):
760

761
762
            lowCover, highCover = coverage[:, dim, vol]
            lowSlice, highSlice = slices[     dim] 
763

764
            if np.isnan(lowCover) or np.isnan(highCover):
765
766
767
768
769
770
771
772
                return False

            if lowSlice  < lowCover:  return False
            if highSlice > highCover: return False

    return True


773
def calcExpansion(slices, coverage):
774
775
776
777
778
779
780
781
    """Calculates a series of *expansion* slices, which can be used to expand
    the given ``coverage`` so that it includes the given ``slices``.

    :arg slices:   Slices that the coverage needs to be expanded to cover.
    :arg coverage: Current image coverage.

    :returns: A list of volume indices, and a corresponding list of
              expansions.
782
783
    """

784
785
786
    numDims         = coverage.shape[1]
    padDims         = len(slices) - numDims - 1
    lowVol, highVol = slices[numDims] 
787
788

    expansions = []
789
    volumes    = []
790

791
792
793
794
795
796
797
798
799
800
    # Finish off an expansion by
    # adding indices for the vector/
    # slice/volume dimension, and for
    # 'padding' dimensions of size 1.
    def finishExpansion(exp, vol):
        exp.append((vol, vol + 1))
        for i in range(padDims):
            exp.append((0, 1))
        return exp
    
801
    for vol in range(lowVol, highVol):
802

803
804
805
        # No coverage of this volume - 
        # we need the whole slice.
        if np.any(np.isnan(coverage[:, :, vol])):
806
807
            exp = [(s[0], s[1]) for s in slices[:numDims]]
            exp = finishExpansion(exp, vol)
808
            volumes   .append(vol)
809
            expansions.append(exp)
810
811
            continue

812
813
814
815
816
817
818
819
        # First we'll figure out the index
        # range for each dimension that
        # needs to be added to the coverage.
        # We build a list of required ranges,
        # where each entry is a tuple
        # containing:
        #   (dimension, lowIndex, highIndex)
        reqRanges = []
820

821
        for dim in range(numDims):
822

823
            lowCover, highCover = coverage[:, dim, vol]
824
            lowSlice, highSlice = slices[     dim]
825

826
827
828
            # The slice covers a region
            # below the current coverage
            if lowCover - lowSlice > 0:
829
                reqRanges.append((dim, int(lowSlice), int(lowCover)))
830
831
832
833
                
            # The slice covers a region
            # above the current coverage
            if highCover - highSlice < 0:
834
                reqRanges.append((dim, int(highCover), int(highSlice)))
835
836

        # Now we generate an expansion for
837
838
        # each of those ranges.
        volExpansions = []
839
840
        for dimx, xlo, xhi in reqRanges:

841
            expansion = [[np.nan, np.nan] for d in range(numDims)]
842
843
844
845
846
847
848
849

            # The expansion for each
            # dimension will span the range
            # for that dimension...
            expansion[dimx][0] = xlo
            expansion[dimx][1] = xhi
                
            # And will span the union of
850
851
            # the coverage, and calculated
            # range for every other dimension.
852
853
854
855
856
            for dimy, ylo, yhi in reqRanges:
                if dimy == dimx:
                    continue

                yLowCover, yHighCover = coverage[:, dimy, vol]
857
                expLow,    expHigh    = expansion[  dimy]
858

859
860
861
862
863
864
865
866
                if np.isnan(expLow):  expLow  = yLowCover
                if np.isnan(expHigh): expHigh = yHighCover

                expLow  = min((ylo, yLowCover,  expLow))
                expHigh = max((yhi, yHighCover, expHigh))

                expansion[dimy][0] = expLow
                expansion[dimy][1] = expHigh
867

868
869
870
871
872
873
874
875
876
877
878
879

            # If no range exists for any of the
            # other dimensions, the range for
            # all expansions will be the current
            # coverage
            for dimy in range(numDims):
                if dimy == dimx:
                    continue

                if np.any(np.isnan(expansion[dimy])):
                    expansion[dimy] = coverage[:, dimy, vol]

880
881
            # Finish off this expansion
            expansion = finishExpansion(expansion, vol)
882

883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
            volumes.      append(vol)
            volExpansions.append(expansion)

        # We do a final run through all pairs
        # of expansions, and adjust their
        # range if they overlap with each other.
        for exp1, exp2 in it.product(volExpansions, volExpansions):

            # Check each dimension
            for dimx in range(numDims):

                xlo1, xhi1 = exp1[dimx]
                xlo2, xhi2 = exp2[dimx]

                # These expansions do not
                # overlap with each other
                # on this dimension (or at
                # all). No need to check
                # the other dimensions.
                if xhi1 <= xlo2: break
                if xlo1 >= xhi2: break

                # These expansions overlap on
                # this dimension - check to see
                # if exp1 is wholly contained
                # within exp2 in all other
                # dimensions.
                adjustable = True

                for dimy in range(numDims):

                    if dimy == dimx:
                        continue

                    ylo1, yhi1 = exp1[dimy]
                    ylo2, yhi2 = exp2[dimy]

                    # Exp1 is not contained within
                    # exp2 on another dimension -
                    # we can't reduce the overlap.
                    if ylo1 < ylo2 or yhi1 > yhi2:
                        adjustable = False
                        break

                # The x dimension range of exp1
                # can be reduced, as it is covered
                # by exp2.
                if adjustable:
                    if   xlo1 <  xlo2 and xhi1 <= xhi2 and xhi1 > xlo2:
                        xhi1 = xlo2

                    elif xlo1 >= xlo2 and xhi1 >  xhi2 and xlo1 < xhi2:
                        xlo1 = xhi2

                    exp1[dimx] = xlo1, xhi1

        expansions.extend(volExpansions)
940
941

    return volumes, expansions
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000


def collapseExpansions(expansions, numDims):
    """Scans through the given list of expansions (each assumed to pertain 
    to a single 3D image), and combines any which cover the same
    image area, and cover adjacent volumes.

    :args expansions: A list of expansion slices - see :func:`calcExpansions`.

    :args numDims:    Number of dimensions covered by each expansion,
                      not including the volume dimension (i.e. 3 for a 4D
                      image).

    :returns: A list of expansions, with equivalent expansions that cover
              adjacent images collapsed down.

    .. note:: For one expansion ``exp`` in the ``expansions`` list, this
              function assumes that the range at ``exp[numDims]`` contains
              the image to which ``exp`` pertains (i.e.
              ``exp[numDims] == (vol, vol + 1)``).
    """
    if len(expansions) == 0:
        return []

    commonExpansions = collections.OrderedDict()
    expansions       = sorted(expansions)

    for exp in expansions:

        vol        = exp[numDims][0]
        exp        = tuple(exp[:numDims])
        commonExps = commonExpansions.get(exp, None)
        
        if commonExps is None:
            commonExps            = []
            commonExpansions[exp] = commonExps

        for i, (vlo, vhi) in enumerate(commonExps):

            if vol >= vlo and vol < vhi:
                break
            
            elif vol == vlo - 1:
                commonExps[i] = vol, vhi
                break
            elif vol == vhi:
                commonExps[i] = vlo, vol + 1
                break
            
        else:
            commonExps.append((vol, vol + 1))
            
    collapsed = []

    for exp, volRanges in commonExpansions.items():
        for vlo, vhi in volRanges:
            newExp = list(exp) + [(vlo, vhi)]
            collapsed.append(newExp)

For faster browsing, not all history is shown. View entire blame