test_imagewrapper.py 33.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
#!/usr/bin/env python
#
# test_imagewrapper.py -
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>



import              collections
import              random
import itertools as it
import numpy     as np
13
import nibabel   as nib
14
15
16
17
18
19
20
21
22
23
24
25
26


import fsl.data.image        as image
import fsl.data.imagewrapper as imagewrap


def setup_module():
    pass

def teardown_module():
    pass


27
def random_coverage(shape, vol_limit=None):
28
29
30
31
    
    ndims = len(shape) - 1
    nvols = shape[-1]

32
33
    print '{}D (shape: {}, {} vectors/slices/volumes)'.format(
        ndims, shape, nvols)
34
35
36
37
38

    # Generate a random coverage. 
    # We use the same coverage for
    # each vector/slice/volume, so
    # are not fully testing the function.
39
    coverage = np.zeros((2, ndims, nvols), dtype=np.float32)
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

    for dim in range(ndims):
        dsize = shape[dim]

        # Random low/high indices for each dimension.
        low  = np.random.randint(0, dsize)

        # We have to make sure that the coverage is not
        # complete, as some of the tests will fail if
        # the coverage is complete. 
        if low == 0: high = np.random.randint(low + 1, dsize)
        else:        high = np.random.randint(low + 1, dsize + 1)

        coverage[0, dim, :] = low
        coverage[1, dim, :] = high
55
56
57

    if vol_limit is not None:
        coverage[:, :, vol_limit:] = np.nan
58
59
60
61
62
63
 
    return coverage


def random_slices(coverage, shape, mode):

64
65
66
67
    shape    = np.array(shape)
    ndims    = len(shape) - 1
    nvols    = shape[-1]
    slices   = np.zeros((2, len(shape)))
68
69
    origMode = mode

70
71
72
73
74
    if coverage is None:
        coverage          = np.zeros((2, ndims, nvols), dtype=np.float32)
        coverage[0, :, :] = 0
        coverage[1, :, :] = shape[:-1].reshape(ndims, 1).repeat(nvols, 1)

75
76
77
78
79
80
81
82
83
    # If we're generating an 'out' slice (i.e.
    # a slice which is not covered by the coverage),
    # then only one dimension needs to be out. The
    # other dimensions don't matter.
    if mode == 'out':
        dimModes = [random.choice(('in', 'out', 'overlap')) for i in range(ndims)]
        if not any([m == 'out' for m in dimModes]):
            dimModes[random.randint(0, ndims - 1)] = 'out'

84
85
86
87
88
89
90
91
92
    for dim, size in enumerate(shape):

        # Volumes 
        if dim == ndims:
            lowCover  = np.random.randint(0,            nvols)
            highCover = np.random.randint(lowCover + 1, nvols + 1) 

            slices[:, dim] = lowCover, highCover
            continue
93
94
95
 
        if origMode == 'out':
            mode = dimModes[dim]
96
97
98
99

        # Assuming that coverage is same for each volume
        lowCover  = coverage[0, dim, 0]
        highCover = coverage[1, dim, 0]
100
101

        if (np.isnan(lowCover) or np.isnan(highCover)) and mode in ('in', 'overlap'):
102
103
104
105
            if origMode == 'out':
                mode = 'out'
            else:
                raise RuntimeError('Can\'t generate in/overlapping slices on an empty coverage')
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
        
        # Generate some slices that will
        # be contained within the coverage
        if mode == 'in':
            lowSlice  = np.random.randint(lowCover,     highCover)
            highSlice = np.random.randint(lowSlice + 1, highCover + 1)

        # Generate some indices which will
        # randomly overlap with the coverage
        # (if it is possible to do so)
        elif mode == 'overlap':

            if highCover == size: lowSlice = np.random.randint(0, lowCover)
            else:                 lowSlice = np.random.randint(0, highCover)

            if lowSlice < lowCover: highSlice = np.random.randint(lowCover  + 1, size + 1)
            else:                   highSlice = np.random.randint(highCover + 1, size + 1)

        elif mode == 'out':

126
127
128
129
130
131
            # No coverage, anything that 
            # we generate will be outside
            if np.isnan(lowCover) or np.isnan(highCover):
                lowSlice  = np.random.randint(0,            size)
                highSlice = np.random.randint(lowSlice + 1, size + 1) 

132
            # The coverage is full, so we can't
133
134
            # generate an outside range
            elif lowCover == 0 and highCover == size:
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
                lowSlice  = np.random.randint(lowCover,     highCover)
                highSlice = np.random.randint(lowSlice + 1, highCover + 1)

            # If low coverage is 0, the slice
            # must be above the coverage
            elif lowCover == 0:
                lowSlice  = np.random.randint(highCover,    size)
                highSlice = np.random.randint(lowSlice + 1, size + 1)

            # If high coverage is size, the
            # slice must be below the coverage
            elif highCover == size:
                lowSlice  = np.random.randint(0,            lowCover)
                highSlice = np.random.randint(lowSlice + 1, lowCover + 1)
                
            # Otherwise the slice could be
            # below or above the coverage
            else:
                lowSlice = random.choice((np.random.randint(0,         lowCover),
                                          np.random.randint(highCover, size)))

                if    lowSlice < lowCover: highSlice = np.random.randint(lowSlice + 1, lowCover + 1)
                else:                      highSlice = np.random.randint(lowSlice + 1, size     + 1)


        slices[0, dim] = lowSlice
        slices[1, dim] = highSlice

    slices = [tuple(map(int, pair)) for pair in slices.T]
    return slices


167
168
169
170
171
172
def rfloat(lo, hi):
    return lo + np.random.random() * (hi - lo)

def applyCoverage(wrapper, coverage):

    ndims = coverage.shape[1]
173
    nvols = coverage.shape[2]
174
175
176
177
178

    wrapper.reset()
    # 'Apply' that coverage to the image
    # wrapper by accessing the data in it

179
180
181
    for vol in range(nvols):
        if np.any(np.isnan(coverage[..., vol])):
            continue
182

183
184
185
186
187
188
189
190
        sliceobjs = []
        for dim in range(ndims):
            sliceobjs.append(
                slice(int(coverage[0, dim, vol]),
                      int(coverage[1, dim, vol]), 1))
        sliceobjs.append(vol)

        wrapper[tuple(sliceobjs)]
191
192
193

    # Check that the image wrapper
    # has covered what we just told
194
    # it to cover
195

196
197
198
199
200
201
202
203
    for vol in range(nvols):

        uncovered = np.any(np.isnan(coverage[..., vol]))
            
        wcov = wrapper.coverage(vol)

        if uncovered:
            assert np.any(np.isnan(wcov))
204

205
        else:
206

207
208
209
210
211
212
            for dim in range(ndims):
                assert coverage[0, dim, vol] == wcov[0, dim]
                assert coverage[1, dim, vol] == wcov[1, dim]


def coverageDataRange(data, coverage, slices=None):
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252

    # Assuming that adjustCoverage is working.
    ndims = coverage.shape[1]
    nvols = coverage.shape[2]

    origcoverage = coverage

    if slices is not None:
        
        coverage = np.copy(coverage)

        lowVol, highVol = slices[-1]

        for vol in range(lowVol, highVol):
            coverage[..., vol] = imagewrap.adjustCoverage(
                coverage[..., vol], slices[:ndims])

    volmin = []
    volmax = []

    for vol in range(nvols):

        cov = coverage[..., vol]

        if np.any(np.isnan(cov)):
            continue
        
        sliceobj = []

        for d in range(ndims):
            sliceobj.append(slice(cov[0, d], cov[1, d], 1))
        sliceobj.append(vol)

        voldata = data[tuple(sliceobj)]
        volmin.append(voldata.min())
        volmax.append(voldata.max())

    return np.min(volmin), np.max(volmax)


253
254
255
256
257
258
259
260
261
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
def test_sliceObjToSliceTuple():

    func  = imagewrap.sliceObjToSliceTuple
    shape = (10, 10, 10)


    assert func( 2,                                       shape) == ((2, 3),  (0, 10), (0, 10))
    assert func( slice(None),                             shape) == ((0, 10), (0, 10), (0, 10))
    assert func((slice(None), slice(None),  slice(None)), shape) == ((0, 10), (0, 10), (0, 10))
    assert func((9,           slice(None),  slice(None)), shape) == ((9, 10), (0, 10), (0, 10))
    assert func((slice(None), 5,            slice(None)), shape) == ((0, 10), (5, 6),  (0, 10))
    assert func((slice(None), slice(None),  3),           shape) == ((0, 10), (0, 10), (3, 4))
    assert func((slice(4, 6), slice(None),  slice(None)), shape) == ((4, 6),  (0, 10), (0, 10))
    assert func((8,           slice(1, 10), slice(None)), shape) == ((8, 9),  (1, 10), (0, 10))



def test_sliceTupleToSliceObj():

    func  = imagewrap.sliceTupleToSliceObj
    shape = (10, 10, 10)

    for x1, y1, z1 in it.product(*[range(d - 1) for d in shape]):

        for x2, y2, z2 in it.product(*[range(s + 1, d) for s, d in zip((x1, y1, z1), shape)]):

            slices   = [[x1, x2], [y1, y2], [z1, z2]]
            sliceobj = (slice(x1, x2, 1), slice(y1, y2, 1), slice(z1, z2, 1))

            assert func(slices) == sliceobj


def test_adjustCoverage():

    # TODO Randomise

289
290
    n = np.nan

291
    # Each test is a tuple of (coverage, expansion, expectedResult) 
292
293
294
295
296
    tests = [([[3, 5], [2, 6]], [[6, 7], [8,  10]],         [[3, 7], [2,  10]]),
             ([[0, 0], [0, 0]], [[1, 2], [3,  5]],          [[0, 2], [0,  5]]),
             ([[2, 3], [0, 6]], [[1, 5], [4,  10]],         [[1, 5], [0,  10]]),
             ([[0, 1], [0, 1]], [[0, 7], [19, 25], [0, 1]], [[0, 7], [0,  25]]),
             ([[n, n], [n, n]], [[0, 7], [19, 25], [0, 1]], [[0, 7], [19, 25]]),
297
298
299
300
301
302
303
304
305
306
    ]

    for coverage, expansion, result in tests:

        result   = np.array(result)  .T
        coverage = np.array(coverage).T

        assert np.all(imagewrap.adjustCoverage(coverage, expansion) == result)


307
def test_sliceOverlap(niters=150):
Paul McCarthy's avatar
Paul McCarthy committed
308
309

    # A bunch of random coverages
310
    for i in range(niters):
Paul McCarthy's avatar
Paul McCarthy committed
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326

        # 2D, 3D or 4D?
        # ndims is the number of dimensions
        # in one vector/slice/volume
        ndims = random.choice((2, 3, 4)) - 1

        # Shape of one vector[2D]/slice[3D]/volume[4D]
        shape = np.random.randint(5, 100, size=ndims + 1)

        # Number of vectors/slices/volumes
        nvols = shape[-1]

        coverage = random_coverage(shape)

        # Generate some slices that should
        # be contained within the coverage
327
        for j in range(niters):
Paul McCarthy's avatar
Paul McCarthy committed
328
329
330
331
332
            slices = random_slices(coverage, shape, 'in')
            assert imagewrap.sliceOverlap(slices, coverage) == imagewrap.OVERLAP_ALL

        # Generate some slices that should
        # overlap with the coverage 
333
        for j in range(niters):
Paul McCarthy's avatar
Paul McCarthy committed
334
335
336
337
338
            slices = random_slices(coverage, shape, 'overlap')
            assert imagewrap.sliceOverlap(slices, coverage) == imagewrap.OVERLAP_SOME

        # Generate some slices that should
        # be outside of the coverage 
339
        for j in range(niters):
Paul McCarthy's avatar
Paul McCarthy committed
340
341
342
343
            slices = random_slices(coverage, shape, 'out')
            assert imagewrap.sliceOverlap(slices, coverage)  == imagewrap.OVERLAP_NONE

        
344
def test_sliceCovered(niters=150):
345
346

    # A bunch of random coverages
347
    for i in range(niters):
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363

        # 2D, 3D or 4D?
        # ndims is the number of dimensions
        # in one vector/slice/volume
        ndims = random.choice((2, 3, 4)) - 1

        # Shape of one vector[2D]/slice[3D]/volume[4D]
        shape = np.random.randint(5, 100, size=ndims + 1)

        # Number of vectors/slices/volumes
        nvols = shape[-1]

        coverage = random_coverage(shape)

        # Generate some slices that should
        # be contained within the coverage
364
        for j in range(niters):
365
            slices = random_slices(coverage, shape, 'in')
Paul McCarthy's avatar
Paul McCarthy committed
366
            assert imagewrap.sliceCovered(slices, coverage)
367
368
369

        # Generate some slices that should
        # overlap with the coverage 
370
        for j in range(niters):
371
            slices = random_slices(coverage, shape, 'overlap')
Paul McCarthy's avatar
Paul McCarthy committed
372
            assert not imagewrap.sliceCovered(slices, coverage)
373
374
375

        # Generate some slices that should
        # be outside of the coverage 
376
        for j in range(niters):
377
            slices = random_slices(coverage, shape, 'out')
Paul McCarthy's avatar
Paul McCarthy committed
378
            assert not imagewrap.sliceCovered(slices, coverage) 
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394


# The sum of the coverage ranges + the
# expansion ranges should be equal to
# the coverage, expanded to include the
# original slices (or the expansions
# - should be equivalent). Note that
# if imagewrapper.adjustCoverage is
# broken, this validation will also be
# broken.
def _test_expansion(coverage, slices, volumes, expansions):
    ndims = coverage.shape[1]

    print
    print 'Slice:    "{}"'.format(" ".join(["{:2d} {:2d}".format(l, h) for l, h in slices]))

395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
    # Figure out what the adjusted
    # coverage should look like (assumes
    # that adjustCoverage is working, and
    # the coverage is the same on all
    # volumes)
    oldCoverage  = coverage[..., 0]
    newCoverage  = imagewrap.adjustCoverage(oldCoverage, slices)

    nc = newCoverage

    # We're goint to test that every point
    # in the expected (expanded) coverage
    # is contained either in the original
    # coverage, or in one of the expansions.
    dimranges = []
    for d in range(ndims):
411
        dimranges.append(np.linspace(nc[0, d], nc[1, d], nc[1, d] / 5, dtype=np.uint32))
412
413
414

    points = it.product(*dimranges)

415
416
417
    # Bin the expansions by volume
    expsByVol = collections.defaultdict(list)
    for vol, exp in zip(volumes, expansions):
418
419
420
        print '  {:3d}:    "{}"'.format(
            int(vol),
            " ".join(["{:2d} {:2d}".format(int(l), int(h)) for l, h in exp]))
421
        expsByVol[vol].append(exp)
422
423
424
425
426
427
428
429
430
431
432
433
434
        
    for point in points:

        # Is this point in the old coverage?
        covered = True
        for dim in range(ndims):
            covLow, covHigh = oldCoverage[:, dim]

            if np.isnan(covLow)    or \
               np.isnan(covHigh)   or \
               point[dim] < covLow or \
               point[dim] > covHigh:
                covered = False
435
436
                break

437
438
439
440
441
        if covered:
            continue
        
        for vol, exps in expsByVol.items():

442
            # Is this point in any of the expansions
443
            coveredInExp = [False] * len(exps)
444
445
            for i, exp in enumerate(exps):
                
446
                coveredInExp[i] = True
447
448
449
450
451
                
                for dim in range(ndims):

                    expLow, expHigh = exp[dim]
                    if point[dim] < expLow or point[dim] > expHigh:
452
                        coveredInExp[i] = False
453
454
                        break
                    
455
456
        if not (covered or any(coveredInExp)):
            raise AssertionError(point)
457

458
            
459
def test_calcExpansionNoCoverage(niters=150):
460

461
    for i in range(niters):
462
463
464
465
466
        ndims       = random.choice((2, 3, 4)) - 1
        shape       = np.random.randint(5, 100, size=ndims + 1)
        shape[-1]   = np.random.randint(1, 8)
        coverage    = np.zeros((2, ndims, shape[-1]))
        coverage[:] = np.nan
467

468
469
        print
        print '-- Out --' 
470
        for j in range(niters):
471
472
473
474
475
            slices     = random_slices(coverage, shape, 'out')
            vols, exps = imagewrap.calcExpansion(slices, coverage)
            _test_expansion(coverage, slices, vols, exps)

                
476
def test_calcExpansion(niters=150):
477

478
    for i in range(niters):
479
480
481
482
483

        ndims     = random.choice((2, 3, 4)) - 1
        shape     = np.random.randint(5, 60, size=ndims + 1)
        shape[-1] = np.random.randint(1, 6)
        coverage  = random_coverage(shape)
484
485
486
487
488
489
490
491

        cov = [(lo, hi) for lo, hi in coverage[:, :, 0].T]

        print 'Shape:    {}'.format(shape)
        print 'Coverage: {}'.format(cov)

        print
        print '-- In --'
492
        for j in range(niters):
493
494
            slices     = random_slices(coverage, shape, 'in')
            vols, exps = imagewrap.calcExpansion(slices, coverage)
495
496
497
498

            # There should be no expansions for a 
            # slice that's inside the coverage
            assert len(vols) == 0 and len(exps) == 0
499
500
501
            
        print
        print '-- Overlap --' 
502
        for j in range(niters):
503
504
505
506
507
508
            slices     = random_slices(coverage, shape, 'overlap')
            vols, exps = imagewrap.calcExpansion(slices, coverage)
            _test_expansion(coverage, slices, vols, exps)

        print
        print '-- Out --' 
509
        for j in range(niters):
510
511
512
            slices     = random_slices(coverage, shape, 'out')
            vols, exps = imagewrap.calcExpansion(slices, coverage)
            _test_expansion(coverage, slices, vols, exps)
513

514

515
def test_ImageWrapper_read(niters=150):
516

517
    for i in range(niters):
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544

        # Generate an image with a number of volumes
        ndims     = random.choice((2, 3, 4)) - 1
        shape     = np.random.randint(5, 60, size=ndims + 1)
        shape[-1] = np.random.randint(5, 15)
        nvols     = shape[-1]

        data = np.zeros(shape)

        # The data range of each volume
        # increases sequentially
        data[..., 0] = np.random.randint(-5, 6, shape[:-1])
        volRanges    = [(np.min(data[..., 0]),
                         np.max(data[..., 0]))]

        for i in range(1, nvols):
            data[..., i] = data[..., 0] * (i + 1)
            volRanges.append((np.min(data[..., i]),
                              np.max(data[..., i])))

        img     = nib.Nifti1Image(data, np.eye(4))
        wrapper = imagewrap.ImageWrapper(img, loadData=False)

        # We're going to access data volumes 
        # through the image wrapper with a
        # bunch of random volume orderings.

545
        for i in range(niters):
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

            ordering = list(range(nvols))
            random.shuffle(ordering)

            ranges = [volRanges[o] for o in ordering]

            wrapper.reset()

            assert wrapper.dataRange == (0.0, 0.0)

            for j, (vol, r) in enumerate(zip(ordering, ranges)):

                # Access the volume
                wrapper[..., vol]

                # The current known data range
                # should be the min/max of
                # all acccessed volumes so far
                expMin = min([r[0] for r in ranges[:j + 1]])
                expMax = max([r[1] for r in ranges[:j + 1]])

                assert wrapper.dataRange == (expMin, expMax)

                if j < nvols - 1: assert not wrapper.covered
                else:             assert     wrapper.covered
571
572
573
574
575
576
577

                
def test_ImageWrapper_write_out(niters=150):
    # This is HORRIBLE
 
    # Generate a bunch of random coverages
    for i in range(niters):
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
        
        # Generate an image with just two volumes. We're only
        # testing within-volume modifications here.
        ndims     = random.choice((2, 3, 4)) - 1
        shape     = np.random.randint(5, 60, size=ndims + 1)
        shape[-1] = np.random.randint(2, 3)
        nvols     = shape[-1]

        data = np.zeros(shape)

        # The data range of each volume
        # increases sequentially
        data[..., 0] = np.random.randint(-5, 6, shape[:-1])
        for i in range(1, nvols):
            data[..., i] = data[..., 0] * (i + 1)
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
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

        # Generate a random coverage
        cov = random_coverage(shape, vol_limit=1)

        print 'This is the coverage: {}'.format(cov)

        img     = nib.Nifti1Image(data, np.eye(4))
        wrapper = imagewrap.ImageWrapper(img)
        applyCoverage(wrapper, cov)
        clo, chi = wrapper.dataRange

        # Now, we'll simulate some writes
        # outside of the coverage area.
        for i in range(niters):

            # Generate some slices outside
            # of the coverage area, making
            # sure that the slice covers
            # at least two elements
            while True:
                slices     = random_slices(cov, shape, 'out')
                slices[-1] = [0, 1]
                sliceshape = [hi - lo for lo, hi in slices]

                if np.prod(sliceshape) == 1:
                    continue

                sliceobjs = imagewrap.sliceTupleToSliceObj(slices)
                sliceobjs  = tuple(list(sliceobjs[:-1]) + [0])
                sliceshape = sliceshape[:-1]
                break

            print '---------------'
            print 'Slice {}'.format(slices)

            # Expected wrapper coverage after the write
            expCov = imagewrap.adjustCoverage(cov[..., 0], slices)

            # Figure out the data range of the
            # expanded coverage (the original
            # coverage expanded to include this
            # slice).
            elo, ehi = coverageDataRange(data, cov, slices)

            # Test all data range possibilities:
            #  - inside the existing range       (clo < rlo < rhi < chi)
            #  - encompassing the existing range (rlo < clo < chi < rhi)
            #  - Overlapping the existing range  (rlo < clo < rhi < chi)
            #                                    (clo < rlo < chi < rhi)
            #  - Outside of the existing range   (clo < chi < rlo < rhi)
            #                                    (rlo < rhi < clo < chi)
            loRanges = [rfloat(clo,         chi), 
                        rfloat(elo - 100,   elo), 
                        rfloat(elo - 100,   elo), 
                        rfloat(clo,         chi), 
                        rfloat(ehi,         ehi + 100), 
                        rfloat(elo - 100,   elo)] 

            hiRanges = [rfloat(loRanges[0], chi), 
                        rfloat(ehi,         ehi + 100), 
                        rfloat(clo,         chi), 
                        rfloat(ehi,         ehi + 100), 
                        rfloat(loRanges[4], ehi + 100), 
                        rfloat(loRanges[5], elo)]
            
            for rlo, rhi in zip(loRanges, hiRanges):

                img     = nib.Nifti1Image(np.copy(data), np.eye(4))
                wrapper = imagewrap.ImageWrapper(img)
                applyCoverage(wrapper, cov)

                print 'ndims', ndims
                print 'sliceshape', sliceshape
                print 'sliceobjs', sliceobjs

                newData = np.linspace(rlo, rhi, np.prod(sliceshape))
                newData = newData.reshape(sliceshape)

                # Make sure that the expected low/high values
                # are present in the data being written

                print 'Writing data (shape: {})'.format(newData.shape)

                oldCov = wrapper.coverage(0)

                wrapper[tuple(sliceobjs)] = newData

                expLo, expHi = coverageDataRange(img.get_data(), cov, slices)
                newLo, newHi = wrapper.dataRange

                print 'Old    range: {} - {}'.format(clo,   chi)
                print 'Sim    range: {} - {}'.format(rlo,   rhi)
                print 'Exp    range: {} - {}'.format(expLo, expHi)
                print 'NewDat range: {} - {}'.format(newData.min(), newData.max())
                print 'Data   range: {} - {}'.format(data.min(),   data.max())
                print 'Expand range: {} - {}'.format(elo, ehi)
                print 'New    range: {} - {}'.format(newLo, newHi)

                newCov = wrapper.coverage(0)
                print 'Old coverage:      {}'.format(oldCov)
                print 'New coverage:      {}'.format(newCov)
                print 'Expected coverage: {}'.format(expCov)
                print
                print

                assert np.all(newCov == expCov)

                assert np.isclose(newLo, expLo)
                assert np.isclose(newHi, expHi)
            print '--------------'

704
705
            
def test_ImageWrapper_write_in_overlap(niters=150):
706
707

    # Generate a bunch of random coverages
708
    for i in range(niters):
709

710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
        # Generate an image with just two volumes. We're only
        # testing within-volume modifications here.
        ndims     = random.choice((2, 3, 4)) - 1
        shape     = np.random.randint(5, 60, size=ndims + 1)
        shape[-1] = np.random.randint(2, 3)
        nvols     = shape[-1]

        data = np.zeros(shape)

        # The data range of each volume
        # increases sequentially
        data[..., 0] = np.random.randint(-5, 6, shape[:-1])
        for i in range(1, nvols):
            data[..., i] = data[..., 0] * (i + 1)        

725
726
727
728
        # Generate a random coverage
        cov = random_coverage(shape, vol_limit=1)

        # Now, we'll simulate some writes
729
730
731
        # which are contained within, or
        # overlap with, the initial coverage
        for i in range(niters):
732

733
734
735
736
            # Generate some slices outside
            # of the coverage area, making
            # sure that the slice covers
            # at least two elements
737
            while True:
738
                slices     = random_slices(cov, shape, random.choice(('in', 'overlap')))
739
740
741
742
743
744
745
746
747
748
749
750
                slices[-1] = [0, 1]
                sliceshape = [hi - lo for lo, hi in slices]

                if np.prod(sliceshape) == 1:
                    continue

                sliceobjs = imagewrap.sliceTupleToSliceObj(slices)
                sliceobjs  = tuple(list(sliceobjs[:-1]) + [0])
                sliceshape = sliceshape[:-1]
                break


751
752
753
754
755
756
757
758
            # Expected wrapper coverage after the write 
            # The write will invalidate the current
            # known data range and coverage, so the
            # expected coverage after the write will
            # be the area covered by the write slice
            nullCov    = np.zeros(cov.shape)
            nullCov[:] = np.nan
            expCov     = imagewrap.adjustCoverage(nullCov[..., 0], slices)
759

760
            for i in range(10):
761

762
763
                rlo = rfloat(data.min() - 100, data.max() + 100)
                rhi = rfloat(rlo,              data.max() + 100)
764
765
766
767
768
769
770
771
772

                img     = nib.Nifti1Image(np.copy(data), np.eye(4))
                wrapper = imagewrap.ImageWrapper(img)
                applyCoverage(wrapper, cov)

                newData = np.linspace(rlo, rhi, np.prod(sliceshape))
                newData = newData.reshape(sliceshape)


773
                print 'Old range:      {} - {}'.format(*wrapper.dataRange)
774
775
776

                wrapper[tuple(sliceobjs)] = newData

777
                newCov       = wrapper.coverage(0)
778
779
                newLo, newHi = wrapper.dataRange

780
781
                print 'Expected range: {} - {}'.format(rlo,   rhi)
                print 'New range:      {} - {}'.format(newLo, newHi)
782
783
784

                assert np.all(newCov == expCov)

785
786
                assert np.isclose(newLo, rlo)
                assert np.isclose(newHi, rhi)
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
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

            
def test_ImageWrapper_write_different_volume(niters=150):

    for i in range(niters):

        # Generate an image with several volumes. 
        ndims     = random.choice((2, 3, 4)) - 1
        nvols     = np.random.randint(10, 40)
        shape     = np.random.randint(5, 60, size=ndims + 1)
        shape[-1] = nvols

        data = np.zeros(shape)

        # The data range of each volume
        # increases sequentially
        data[..., 0] = np.random.randint(-5, 6, shape[:-1])
        for i in range(1, nvols):
            data[..., i] = data[..., 0] * (i + 1)


        # Generate a random coverage
        fullcov = random_coverage(shape)
        cov     = np.copy(fullcov)

        # Choose some consecutive volumes
        # to limit that coverage to.
        while True:
            covvlo = np.random.randint(0,          nvols - 1)
            covvhi = np.random.randint(covvlo + 1, nvols + 1)

            # Only include up to 4
            # volumes in the coverage
            if covvhi - covvlo <= 3:
                break

        for v in range(nvols):
            if v < covvlo or v >= covvhi:
                cov[..., v] = np.nan

        covlo, covhi = coverageDataRange(data, cov)

        print 'Coverage: {} [{} - {}]'.format(
            [(lo, hi) for lo, hi in cov[:, :, covvlo].T],
            covvlo, covvhi)

        # Now, we'll simulate some writes
        # on volumes that are not in the
        # coverage
        for i in range(niters):

            # Generate a slice, making
            # sure that the slice covers
            # at least two elements
            while True:
                slices = random_slices(fullcov,
                                       shape,
                                       random.choice(('out', 'in', 'overlap')))

                # print slices

                # Clobber the slice volume range
                # so it does not overlap with the
                # coverage volume range
                while True:
                    vlo = np.random.randint(0,       nvols)
                    vhi = np.random.randint(vlo + 1, nvols + 1)
                    
                    if vhi < covvlo or vlo > covvhi:
                        break
                    
                slices[-1] = vlo, vhi
                
                sliceshape = [hi - lo for lo, hi in slices]

                if np.prod(sliceshape) == 1:
                    continue

                sliceobjs = imagewrap.sliceTupleToSliceObj(slices)
                break

            # Calculate what we expect the 
            # coverage to be after the write
            expCov = np.copy(cov)
            for vol in range(slices[-1][0], slices[-1][1]):
                expCov[..., vol] = imagewrap.adjustCoverage(
                    expCov[..., vol], slices)

            # Test all data range possibilities:
            #  - inside the existing range       (clo < rlo < rhi < chi)
            #  - encompassing the existing range (rlo < clo < chi < rhi)
            #  - Overlapping the existing range  (rlo < clo < rhi < chi)
            #                                    (clo < rlo < chi < rhi)
            #  - Outside of the existing range   (clo < chi < rlo < rhi)
            #                                    (rlo < rhi < clo < chi)
            
            loRanges = [rfloat(covlo,         covhi), 
                        rfloat(covlo - 100,   covlo), 
                        rfloat(covlo - 100,   covlo), 
                        rfloat(covlo,         covhi), 
                        rfloat(covhi,         covhi + 100), 
                        rfloat(covlo - 100,   covlo)] 

            hiRanges = [rfloat(loRanges[0], covhi), 
                        rfloat(covhi,       covhi + 100), 
                        rfloat(covlo,       covhi), 
                        rfloat(covhi,       covhi + 100), 
                        rfloat(loRanges[4], covhi + 100), 
                        rfloat(loRanges[5], covlo)]

            # What we expect the range to 
            # be after the data write
            expected = [(covlo,       covhi),
                        (loRanges[1], hiRanges[1]),
                        (loRanges[2], covhi),
                        (covlo,       hiRanges[3]),
                        (covlo,       hiRanges[4]),
                        (loRanges[5], covhi)]

            for rlo, rhi, (elo, ehi) in zip(loRanges, hiRanges, expected):

                img     = nib.Nifti1Image(np.copy(data), np.eye(4))
                wrapper = imagewrap.ImageWrapper(img)
                applyCoverage(wrapper, cov)

                newData = np.linspace(rlo, rhi, np.prod(sliceshape))
                newData = newData.reshape(sliceshape)

                wrapper[tuple(sliceobjs)] = newData

                newLo, newHi = wrapper.dataRange

                for vol in range(nvols):
                    np.all(wrapper.coverage(vol) == expCov[..., vol])
                    
                assert np.isclose(newLo, elo)
                assert np.isclose(newHi, ehi)
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
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
1001
1002
1003
1004
1005
1006
1007
1008
1009


def test_collapseExpansions(niters=500):

    def expEq(exp1, exp2):
        
        if len(exp1) != len(exp2):
            return False
        
        for (e1lo, e1hi), (e2lo, e2hi) in zip(exp1, exp2):
            if e1lo != e2lo: return False
            if e1hi != e2hi: return False
        return True

    def rangesOverlapOrAdjacent(r1, r2):

        r1lo, r1hi = r1
        r2lo, r2hi = r2

        return not ((r1lo > r2hi) or (r1hi < r2lo) or \
                    (r2lo > r1hi) or (r2hi < r1lo))

    for i in range(niters):

        # Generate a random coverage shape
        ndims     = random.choice((2, 3, 4)) - 1
        nvols     = np.random.randint(10, 40)
        shape     = np.random.randint(5, 60, size=ndims + 1)
        shape[-1] = nvols

        # Generate a bunch of random slices, and
        # split each of them up by volume to
        # turn them each into a set of expansions
        exps     = []
        expected = []

        for i in range(10):

            # Generate a random slice with a volume
            # range that doesn't overlap with, and
            # is not adjacent to, any of the other
            # generated random slices.
            #
            # The only reason I'm doing this is to
            # overocme a chicken-and-egg problem -
            # if I want to test overlapping/adjacent
            # region, I basically have to
            # re-implement the collapseExpansions
            # function.
            while True:
                slices = random_slices(None, shape, 'in')
                vlo, vhi = slices[-1]

                for exp in expected:

                    if not expEq(exp[:-1], slices[:-1]):
                        continue
                    
                    evlo, evhi = exp[-1]

                    # Overlap
                    if rangesOverlapOrAdjacent((vlo, vhi), (evlo, evhi)):
                        break
                else:
                    break

            expected.append(slices)

            for vol in range(slices[-1][0], slices[-1][1]):

                exp = tuple(list(slices[:-1]) + [(vol, vol + 1)])
                exps.append(exp)

        # Now shuffle them up randomly
        np.random.shuffle(exps)

        # And attempt to collapse them
        collapsed = imagewrap.collapseExpansions(exps, ndims)

        collapsed = sorted(collapsed)
        expected  = sorted(expected)

        assert len(expected) == len(collapsed)

        for exp, col in zip(expected, collapsed):
            assert expEq(exp, col)