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


8
from __future__ import print_function
9
10
11
12
13

import              collections
import              random
import itertools as it
import numpy     as np
14
import nibabel   as nib
15
import              pytest
16

Paul McCarthy's avatar
Paul McCarthy committed
17
import fsl.utils.naninfrange as nir
18
19
import fsl.data.imagewrapper as imagewrap

20
21
from . import random_voxels

22

23
24
25
26
27
real_print = print
def print(*args, **kwargs):
    pass


28
29
30
31
32
33
34
def setup_module():
    pass

def teardown_module():
    pass


35

36
def random_coverage(shape, vol_limit=None):
37

38
39
40
    ndims = len(shape) - 1
    nvols = shape[-1]

41
    # Generate a random coverage.
42
43
44
    # We use the same coverage for
    # each vector/slice/volume, so
    # are not fully testing the function.
45
    coverage = np.zeros((2, ndims, nvols), dtype=np.float32)
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
55
        # the coverage is complete.
56
57
58
59
60
        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
61
62
63

    if vol_limit is not None:
        coverage[:, :, vol_limit:] = np.nan
64

65
66
67
68
69
    return coverage


def random_slices(coverage, shape, mode):

70
71
72
73
    shape    = np.array(shape)
    ndims    = len(shape) - 1
    nvols    = shape[-1]
    slices   = np.zeros((2, len(shape)))
74
75
    origMode = mode

76
77
78
79
80
    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)

81
82
83
84
85
86
87
88
89
    # 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'

90
91
    for dim, size in enumerate(shape):

92
        # Volumes
93
94
        if dim == ndims:
            lowCover  = np.random.randint(0,            nvols)
95
            highCover = np.random.randint(lowCover + 1, nvols + 1)
96
97
98

            slices[:, dim] = lowCover, highCover
            continue
99

100
101
        if origMode == 'out':
            mode = dimModes[dim]
102
103
104
105

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

        if (np.isnan(lowCover) or np.isnan(highCover)) and mode in ('in', 'overlap'):
108
109
110
111
            if origMode == 'out':
                mode = 'out'
            else:
                raise RuntimeError('Can\'t generate in/overlapping slices on an empty coverage')
112

113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
        # 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':

132
            # No coverage, anything that
133
134
135
            # we generate will be outside
            if np.isnan(lowCover) or np.isnan(highCover):
                lowSlice  = np.random.randint(0,            size)
136
                highSlice = np.random.randint(lowSlice + 1, size + 1)
137

138
            # The coverage is full, so we can't
139
140
            # generate an outside range
            elif lowCover == 0 and highCover == size:
141
142
143
144
145
146
147
148
149
150
151
152
153
154
                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)
155

156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
            # 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


173
174
175
176
177
178
def rfloat(lo, hi):
    return lo + np.random.random() * (hi - lo)

def applyCoverage(wrapper, coverage):

    ndims = coverage.shape[1]
179
    nvols = coverage.shape[2]
180
181
182
183
184

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

185
186
187
    for vol in range(nvols):
        if np.any(np.isnan(coverage[..., vol])):
            continue
188

189
190
191
192
193
194
195
196
        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)]
197
        _ImageWraper_busy_wait(wrapper)
198
199
200

    # Check that the image wrapper
    # has covered what we just told
201
    # it to cover
202

203
204
205
    for vol in range(nvols):

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

207
208
209
210
        wcov = wrapper.coverage(vol)

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

212
        else:
213

214
215
216
217
218
219
            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):
220
221
222
223
224
225
226
227

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

    origcoverage = coverage

    if slices is not None:
228

229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
        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
246

247
248
249
        sliceobj = []

        for d in range(ndims):
Paul McCarthy's avatar
Paul McCarthy committed
250
            sliceobj.append(slice(int(cov[0, d]), int(cov[1, d]), 1))
251
252
253
254
255
256
257
258
259
        sliceobj.append(vol)

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

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


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
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
def test_expectedShape():

    tests = [
        ((slice(None), ), (10,),
         (1, (10, ))),

        ((slice(None), slice(None)),
         (10, 10), (2, (10, 10))),

        ((slice(None), slice(None), slice(None)),
         (10, 10, 10), (3, (10, 10, 10))),

        ((slice(None), slice(None), slice(None)),
         (10, 10, 10), (3, (10, 10, 10))),

        ((slice(None), slice(None), slice(None), slice(None)),
         (10, 10, 10, 10), (4, (10, 10, 10, 10))),

        ((1, slice(None), slice(None)),
         (10, 10, 10), (2, (10, 10))),

        ((slice(1, 3), slice(None), slice(None)),
         (10, 10, 10), (3, (2, 10, 10))),

        ((slice(None), 1, slice(None)),
         (10, 10, 10), (2, (10, 10))),

        ((slice(None), slice(1, 3), slice(None)),
         (10, 10, 10), (3, (10, 2, 10))),

        ((slice(None), slice(None), 1),
         (10, 10, 10), (2, (10, 10))),

        ((slice(None), slice(None), slice(1, 3), ),
         (10, 10, 10), (3, (10, 10, 2))),

        ((slice(None), slice(None), slice(1, 20), ),
         (10, 10, 10), (3, (10, 10, 9))),
    ]

    for slc, shape, exp in tests:

        explen, exp = exp
        gotlen, got = imagewrap.expectedShape(slc, shape)

        assert explen     == gotlen
        assert tuple(exp) == tuple(got)


309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
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


341

342
343
344
345
def test_adjustCoverage():

    # TODO Randomise

346
347
    n = np.nan

348
    # Each test is a tuple of (coverage, expansion, expectedResult)
349
350
351
352
353
    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]]),
354
355
356
357
358
359
360
361
362
363
    ]

    for coverage, expansion, result in tests:

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

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


364
def test_sliceOverlap(niters):
Paul McCarthy's avatar
Paul McCarthy committed
365
366

    # A bunch of random coverages
367
    for _ in range(niters):
Paul McCarthy's avatar
Paul McCarthy committed
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383

        # 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
384
        for _ in range(niters):
Paul McCarthy's avatar
Paul McCarthy committed
385
386
387
388
            slices = random_slices(coverage, shape, 'in')
            assert imagewrap.sliceOverlap(slices, coverage) == imagewrap.OVERLAP_ALL

        # Generate some slices that should
389
        # overlap with the coverage
390
        for _ in range(niters):
Paul McCarthy's avatar
Paul McCarthy committed
391
392
393
394
            slices = random_slices(coverage, shape, 'overlap')
            assert imagewrap.sliceOverlap(slices, coverage) == imagewrap.OVERLAP_SOME

        # Generate some slices that should
395
        # be outside of the coverage
396
        for _ in range(niters):
Paul McCarthy's avatar
Paul McCarthy committed
397
398
399
            slices = random_slices(coverage, shape, 'out')
            assert imagewrap.sliceOverlap(slices, coverage)  == imagewrap.OVERLAP_NONE

400

401
def test_sliceCovered(niters):
402
403

    # A bunch of random coverages
404
    for _ in range(niters):
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420

        # 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
421
        for _ in range(niters):
422
            slices = random_slices(coverage, shape, 'in')
Paul McCarthy's avatar
Paul McCarthy committed
423
            assert imagewrap.sliceCovered(slices, coverage)
424
425

        # Generate some slices that should
426
        # overlap with the coverage
427
        for _ in range(niters):
428
            slices = random_slices(coverage, shape, 'overlap')
Paul McCarthy's avatar
Paul McCarthy committed
429
            assert not imagewrap.sliceCovered(slices, coverage)
430
431

        # Generate some slices that should
432
        # be outside of the coverage
433
        for _ in range(niters):
434
            slices = random_slices(coverage, shape, 'out')
435
            assert not imagewrap.sliceCovered(slices, coverage)
436
437
438
439
440
441
442
443
444
445
446
447
448


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

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

452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
    # 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):
Paul McCarthy's avatar
Paul McCarthy committed
468
469
470
471
        dimranges.append(np.linspace(nc[0, d],
                                     nc[1, d],
                                     int(nc[1, d] / 5),
                                     dtype=np.uint32))
472
473
474

    points = it.product(*dimranges)

475
476
477
    # Bin the expansions by volume
    expsByVol = collections.defaultdict(list)
    for vol, exp in zip(volumes, expansions):
478
        print('  {:3d}:    "{}"'.format(
479
            int(vol),
480
            " ".join(["{:2d} {:2d}".format(int(l), int(h)) for l, h in exp])))
481
        expsByVol[vol].append(exp)
482

483
484
485
486
487
488
489
490
491
492
493
494
    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
495
496
                break

497
498
        if covered:
            continue
499

500
501
        for vol, exps in expsByVol.items():

502
            # Is this point in any of the expansions
503
            coveredInExp = [False] * len(exps)
504
            for i, exp in enumerate(exps):
505

506
                coveredInExp[i] = True
507

508
509
510
511
                for dim in range(ndims):

                    expLow, expHigh = exp[dim]
                    if point[dim] < expLow or point[dim] > expHigh:
512
                        coveredInExp[i] = False
513
                        break
514

515
516
        if not (covered or any(coveredInExp)):
            raise AssertionError(point)
517

518

519
def test_calcExpansionNoCoverage(niters):
520

521
    for _ in range(niters):
522
523
524
525
526
        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
527

528
529
530
        print()
        print('-- Out --' )
        for _ in range(niters):
531
532
533
534
            slices     = random_slices(coverage, shape, 'out')
            vols, exps = imagewrap.calcExpansion(slices, coverage)
            _test_expansion(coverage, slices, vols, exps)

535

536
def test_calcExpansion(niters):
537

538
    for _ in range(niters):
539
540
541
542
543

        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)
544
545
546

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

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

550
551
552
        print()
        print('-- In --')
        for _ in range(niters):
553
554
            slices     = random_slices(coverage, shape, 'in')
            vols, exps = imagewrap.calcExpansion(slices, coverage)
555

556
            # There should be no expansions for a
557
558
            # slice that's inside the coverage
            assert len(vols) == 0 and len(exps) == 0
559

560
561
562
        print()
        print('-- Overlap --' )
        for _ in range(niters):
563
564
565
566
            slices     = random_slices(coverage, shape, 'overlap')
            vols, exps = imagewrap.calcExpansion(slices, coverage)
            _test_expansion(coverage, slices, vols, exps)

567
568
569
        print()
        print('-- Out --' )
        for _ in range(niters):
570
571
572
            slices     = random_slices(coverage, shape, 'out')
            vols, exps = imagewrap.calcExpansion(slices, coverage)
            _test_expansion(coverage, slices, vols, exps)
573

574

575
576
577
578
579
def _ImageWraper_busy_wait(wrapper, v=0):
    tt = wrapper.getTaskThread()
    if tt is not None:
        tt.waitUntilIdle()

580
@pytest.mark.longtest
581
582
def test_ImageWrapper_read_threaded(niters, seed):
    _test_ImageWrapper_read(niters, seed, True)
583
@pytest.mark.longtest
584
def test_ImageWrapper_read_unthreaded(niters, seed):
585
    _test_ImageWrapper_read(niters, seed, False)
586
@pytest.mark.longtest
587
588
def test_ImageWrapper_read_nans_threaded(niters, seed):
    _test_ImageWrapper_read(niters, seed, True, True)
589
@pytest.mark.longtest
590
591
592
593
def test_ImageWrapper_read_nans_unthreaded(niters, seed):
    _test_ImageWrapper_read(niters, seed, False, True)

def _test_ImageWrapper_read(niters, seed, threaded, addnans=False):
594

595
    for _ in range(niters):
596
597
598
599
600
601
602
603
604
605
606
607

        # 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])
608
609
        for vol in range(1, nvols):
            data[..., vol] = data[..., 0] * (vol + 1)
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629

        # Add 10-50% infs/nans to the image
        if addnans:

            infprop = 0.10 + 0.2 * np.random.random()
            nanprop = 0.10 + 0.2 * np.random.random()

            ninfs = int(infprop * np.prod(shape[:-1]))
            nnans = int(nanprop * np.prod(shape[:-1]))

            for vol in range(nvols):
                coords = random_voxels(shape[:-1], int(ninfs))
                coords = tuple([coords[..., i] for i in range(ndims)] + [vol])
                data[coords] = np.inf

                coords = random_voxels(shape[:-1], int(nnans))
                coords = tuple([coords[..., i] for i in range(ndims)] + [vol])
                data[coords] = np.nan


Paul McCarthy's avatar
Paul McCarthy committed
630
        volRanges = [nir.naninfrange(data[..., v]) for v in range(nvols)]
631
632
633
634
635

        rrs = []
        for vol in range(nvols):
            rrs.append('{:3d}: {: 3.0f} - {: 3.0f}'.format(
                vol, *volRanges[vol]))
636
637

        img     = nib.Nifti1Image(data, np.eye(4))
638
639
640
        wrapper = imagewrap.ImageWrapper(img,
                                         loadData=False,
                                         threaded=threaded)
641

642
        # We're going to access data volumes
643
644
        # through the image wrapper with a
        # bunch of random volume orderings.
645
        for _ in range(niters):
646
647
648
649
650
651
652
653
654
655
656
657
658

            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
Paul McCarthy's avatar
Paul McCarthy committed
659
660
                if len(data.shape) >= 3: wrapper[..., vol]
                else:                    wrapper[:, vol, 0]
661

662
663
                _ImageWraper_busy_wait(wrapper, vol)

664
665
666
667
668
669
670
671
672
673
                # 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
674

675

676
@pytest.mark.longtest
677
678
def test_ImageWrapper_write_out_threaded(niters, seed):
    _test_ImageWrapper_write_out(niters, seed, True)
679
@pytest.mark.longtest
680
681
682
def test_ImageWrapper_write_out_unthreaded(niters, seed):
    _test_ImageWrapper_write_out(niters, seed, False)
def _test_ImageWrapper_write_out(niters, seed, threaded):
683
    # This is HORRIBLE
684
685
686

    loop = 0

687

688
    # Generate a bunch of random coverages
689
    for _ in range(niters):
690

691
692
693
694
695
696
697
698
699
700
701
702
        # 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])
703
704
        for vol in range(1, nvols):
            data[..., vol] = data[..., 0] * (vol + 1)
705
706
707
708
709

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

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

714
715
        loop += 1

716
717
        # Now, we'll simulate some writes
        # outside of the coverage area.
718
        for _ in range(niters):
719
720
721
722

            # Generate some slices outside
            # of the coverage area, making
            # sure that the slice covers
723
            # at least one element
724
725
726
727
728
            while True:
                slices     = random_slices(cov, shape, 'out')
                slices[-1] = [0, 1]
                sliceshape = [hi - lo for lo, hi in slices]

729
                if np.prod(sliceshape) == 0:
730
731
732
733
734
735
736
                    continue

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

737
738
            # print('---------------')
            # print('Slice {}'.format(slices))
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755

            # 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)
756
757
758
759
760
761
762
763
764
765
766
767
            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),
768
                        rfloat(loRanges[5], elo)]
769

770
771
772
773
774
775
            for rlo, rhi in zip(loRanges, hiRanges):

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

776
777
778
                # print('ndims', ndims)
                # print('sliceshape', sliceshape)
                # print('sliceobjs', sliceobjs)
779
780
781
782
783
784
785

                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

786
                # print('Writing data (shape: {})'.format(newData.shape))
787
788
789
790

                oldCov = wrapper.coverage(0)

                wrapper[tuple(sliceobjs)] = newData
791
                _ImageWraper_busy_wait(wrapper)
792

Paul McCarthy's avatar
Paul McCarthy committed
793
                expLo, expHi = coverageDataRange(np.asanyarray(img.dataobj), cov, slices)
794
795
                newLo, newHi = wrapper.dataRange

796
797
798
799
800
801
802
                # 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))
803
804

                newCov = wrapper.coverage(0)
805
806
807
808
809
                # print('Old coverage:      {}'.format(oldCov))
                # print('New coverage:      {}'.format(newCov))
                # print('Expected coverage: {}'.format(expCov))
                # print()
                # print()
810
811
812
813
814

                assert np.all(newCov == expCov)

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

817

818
@pytest.mark.longtest
819
820
def test_ImageWrapper_write_in_overlap_threaded(niters, seed):
    _test_ImageWrapper_write_in_overlap(niters, seed, True)
821
@pytest.mark.longtest
822
def test_ImageWrapper_write_in_overlap_unthreaded(niters, seed):
823
    _test_ImageWrapper_write_in_overlap(niters, seed, False)
824
def _test_ImageWrapper_write_in_overlap(niters, seed, threaded):
825
826

    # Generate a bunch of random coverages
827
    for _ in range(niters):
828

829
830
831
832
833
834
835
836
837
838
839
840
        # 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])
841
842
        for vol in range(1, nvols):
            data[..., vol] = data[..., 0] * (vol + 1)
843

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

Paul McCarthy's avatar
Paul McCarthy committed
847
848
849
        print('Shape:    {}'.format(shape))
        print('Coverage: {}'.format(cov))
        print('Data:     {}'.format(data))
850

851
        # Now, we'll simulate some writes
852
853
        # which are contained within, or
        # overlap with, the initial coverage
854
        for _ in range(niters):
855

856
857
858
            # Generate some slices outside
            # of the coverage area, making
            # sure that the slice covers
859
            # at least one element
860
            while True:
861
                slices     = random_slices(cov, shape, random.choice(('in', 'overlap')))
862
863
864
                slices[-1] = [0, 1]
                sliceshape = [hi - lo for lo, hi in slices]

865
                if np.prod(sliceshape) == 0:
866
867
868
869
870
871
872
                    continue

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

873
874
            # Expected wrapper coverage after the
            # write is the union of the original
Paul McCarthy's avatar
Paul McCarthy committed
875
876
            # coverage and the write slice.
            expCov = imagewrap.adjustCoverage(cov[..., 0], slices)
877

878
            for _ in range(10):
879

880
881
                rlo = rfloat(data.min() - 100, data.max() + 100)
                rhi = rfloat(rlo,              data.max() + 100)
882

Paul McCarthy's avatar
Paul McCarthy committed
883
884
885
                if np.prod(sliceshape) == 1:
                    rhi = rlo

886
                img     = nib.Nifti1Image(np.copy(data), np.eye(4))
887
                wrapper = imagewrap.ImageWrapper(img, threaded=threaded)
888
889
890
891
892
                applyCoverage(wrapper, cov)

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

Paul McCarthy's avatar
Paul McCarthy committed
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
                print('Old coverage:      {}'.format(cov[..., 0]))
                print('Slice:             {}'.format(sliceobjs[:-1]))
                print('Expected coverage: {}'.format(expCov))
                print('Old range:         {} - {}'.format(*wrapper.dataRange))
                print('New data range:    {} - {}'.format(newData.min(), newData.max()))

                # We figure out the expected data
                # range by creating a copy of the
                # data, and doing the same write
                expData = np.copy(data[..., 0])
                expData[sliceobjs[:-1]] = newData

                # Then calcultaing the min/max
                # on this copy
                expCovSlice = [slice(int(lo), int(hi)) for lo, hi in expCov.T]

909
910
                expLo = expData[tuple(expCovSlice)].min()
                expHi = expData[tuple(expCovSlice)].max()
911
912

                wrapper[tuple(sliceobjs)] = newData
913
                _ImageWraper_busy_wait(wrapper)
914

915
                newCov       = wrapper.coverage(0)
916
917
                newLo, newHi = wrapper.dataRange

Paul McCarthy's avatar
Paul McCarthy committed
918
919
                print('Expected range:    {} - {}'.format(expLo, expHi))
                print('New range:         {} - {}'.format(newLo, newHi))
Paul McCarthy's avatar
Paul McCarthy committed
920
921
922
923
                print('Slice min/max:     {} - {}'.format(np.asanyarray(img.dataobj)[tuple(sliceobjs)].min(),
                                                          np.asanyarray(img.dataobj)[tuple(sliceobjs)].max()))
                print('Data min/max:      {} - {}'.format(np.asanyarray(img.dataobj).min(),
                                                          np.asanyarray(img.dataobj).max()))
924
925
926

                assert np.all(newCov == expCov)

Paul McCarthy's avatar
Paul McCarthy committed
927
928
                assert np.isclose(newLo, expLo)
                assert np.isclose(newHi, expHi)
929

930

931
@pytest.mark.longtest
932
933
def test_ImageWrapper_write_different_volume_threaded(niters, seed):
    _test_ImageWrapper_write_different_volume(niters, seed, True)
934
@pytest.mark.longtest
935
936
937
def test_ImageWrapper_write_different_volume_unthreaded(niters, seed):
    _test_ImageWrapper_write_different_volume(niters, seed, False)
def _test_ImageWrapper_write_different_volume(niters, seed, threaded):
938

939
    for _ in range(niters):
940

941
        # Generate an image with several volumes.
942
943
944
945
946
947
948
949
950
951
        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])
952
953
        for vol in range(1, nvols):
            data[..., vol] = data[..., 0] * (vol + 1)
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976


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

977
        print('Coverage: {} [{} - {}]'.format(
978
            [(lo, hi) for lo, hi in cov[:, :, covvlo].T],
979
            covvlo, covvhi))
980
981
982
983

        # Now, we'll simulate some writes
        # on volumes that are not in the
        # coverage
984
        for _ in range(niters):
985
986
987

            # Generate a slice, making
            # sure that the slice covers
Paul McCarthy's avatar
Paul McCarthy committed
988
            # at least one element
989
990
991
992
993
            while True:
                slices = random_slices(fullcov,
                                       shape,
                                       random.choice(('out', 'in', 'overlap')))

994
                # print(slices)
995
996
997
998
999
1000
1001

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

1003
1004
                    if vhi < covvlo or vlo > covvhi:
                        break
1005

1006
                slices[-1] = vlo, vhi
1007

1008
1009
                sliceshape = [hi - lo for lo, hi in slices]

1010
                if np.prod(sliceshape) == 0:
1011
1012
1013
1014
1015
                    continue

                sliceobjs = imagewrap.sliceTupleToSliceObj(slices)
                break

1016
            # Calculate what we expect the
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
            # 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)
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042

            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),
1043
1044
                        rfloat(loRanges[5], covlo)]

1045
            # What we expect the range to
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
            # 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))
1057
                wrapper = imagewrap.ImageWrapper(img, threaded=threaded)
1058
1059
                applyCoverage(wrapper, cov)

Paul McCarthy's avatar
Paul McCarthy committed
1060
1061
1062
                oldLo, oldHi = wrapper.dataRange


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

Paul McCarthy's avatar
Paul McCarthy committed
1066
                if np.prod(sliceshape) == 1:
1067
                    ehi = max(newData.max(), oldHi)
Paul McCarthy's avatar
Paul McCarthy committed
1068

1069
                wrapper[tuple(sliceobjs)] = newData
1070
                _ImageWraper_busy_wait(wrapper)
1071
1072
1073
1074
1075

                newLo, newHi = wrapper.dataRange

                for vol in range(nvols):
                    np.all(wrapper.coverage(vol) == expCov[..., vol])
Paul McCarthy's avatar
Paul McCarthy committed
1076
1077
1078
1079
1080

                print('Old range:      {} - {}'.format(oldLo, oldHi))
                print('Newdata range:  {} - {}'.format(newData.min(), newData.max()))
                print('Expected range: {} - {}'.format(elo,   ehi))
                print('New range:      {} - {}'.format(newLo, newHi))
1081

1082
1083
                assert np.isclose(newLo, elo)
                assert np.isclose(newHi, ehi)
1084
1085


1086
def test_collapseExpansions(niters):
1087
1088

    def expEq(exp1, exp2):
1089

1090
1091
        if len(exp1) != len(exp2):
            return False
1092

1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
        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))

1106
    for _ in range(niters):
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119

        # 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 = []

1120
        for _ in range(10):
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140

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

1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
                    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)
Paul McCarthy's avatar
Paul McCarthy committed
1170
1171
1172
1173
1174
1175


def test_3D_indexing(shape=None, img=None):

    # Test that a 3D image looks like a 3D image

1176
1177
1178
1179
    if   shape is None:
        shape = (21, 22, 23)
    elif len(shape) < 3:
        shape = tuple(list(shape) + [1] * (3 - len(shape)))
1180

Paul McCarthy's avatar
Paul McCarthy committed
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
    if img is None:
        data   = np.random.random(shape)
        nibImg = nib.Nifti1Image(data, np.eye(4))
        img    = imagewrap.ImageWrapper(nibImg, loadData=True)

    assert tuple(img[:]      .shape) == tuple(shape)
    assert tuple(img[:, :]   .shape) == tuple(shape)
    assert tuple(img[:, :, :].shape) == tuple(shape)
    assert tuple(img[:, 0, 0].shape) == (shape[0], )
    assert tuple(img[0, :, 0].shape) == (shape[1], )
    assert tuple(img[0, 0, :].shape) == (shape[2], )
    assert tuple(img[0, :, :].shape) == (shape[1], shape[2])
    assert tuple(img[:, 0, :].shape) == (shape[0], shape[2])
    assert tuple(img[:, :, 0].shape) == (shape[0], shape[1])

    assert type(img[0, 0, 0]) == np.float64

1198
    mask1 = np.zeros(shape, dtype=bool)
Paul McCarthy's avatar
Paul McCarthy committed
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210