test_label.py 36.8 KB
Newer Older
1
2
3
4
5
6
# -*- coding: utf-8 -*-

import unittest
import os
import tifffile
import numpy
7
8
import random

9

Edward Andò's avatar
Edward Andò committed
10
import spam.label
11
import spam.kalisphera
12
import spam.label.contacts as con
Gustavo Pinzon's avatar
Gustavo Pinzon committed
13
import spam.plotting
14
import skimage.morphology
15
VERBOSE = True
16

Edward Andò's avatar
Edward Andò committed
17
# small labelled volume is a 3x3x3 with a single 1 in the middle
18
19
20
21
22
threeCubedLabelVol = numpy.zeros((3, 3, 3), dtype=spam.label.labelType)
threeCubedLabelVol[1, 1, 1] = 1

longSingle = numpy.zeros((9, 6, 4), dtype=spam.label.labelType)
longSingle[1:9, 1:6, 1:3] = 1
Edward Andò's avatar
Edward Andò committed
23

24

Edward Andò's avatar
Edward Andò committed
25
26
class TestFunctionLabel(unittest.TestCase):

27
28
29
30
31
    # def tearDown(self):
    #     try:
    #         pass
    #     except OSError:
    #         pass
Edward Andò's avatar
Edward Andò committed
32
33

    def test_boundingBoxes(self):
34
        BB = spam.label.boundingBoxes(threeCubedLabelVol)
35
36
        self.assertEqual(numpy.zeros(6, dtype=spam.label.labelType).tolist(), BB[0].tolist())
        self.assertEqual(numpy.ones(6, dtype=spam.label.labelType).tolist(), BB[1].tolist())
37
38

        BB = spam.label.boundingBoxes(longSingle)
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
        self.assertEqual(numpy.array([1, 8, 1, 5, 1, 2], dtype=spam.label.labelType).tolist(), BB[1].tolist())

    def test_centresOfMass(self):
        # without pre-computation of BB
        COM = spam.label.centresOfMass(threeCubedLabelVol)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), COM[0].tolist())
        self.assertEqual(numpy.ones(3, dtype='<f4').tolist(), COM[1].tolist())

        # put in min vol filter
        COM = spam.label.centresOfMass(threeCubedLabelVol, minVol=2)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), COM[0].tolist())
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), COM[1].tolist())

        # with pre-computation of BB
        BB = spam.label.boundingBoxes(threeCubedLabelVol)
        COM = spam.label.centresOfMass(threeCubedLabelVol, boundingBoxes=BB)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), COM[0].tolist())
        self.assertEqual(numpy.ones(3, dtype='<f4').tolist(), COM[1].tolist())

    def test_volumes(self):
        # without pre-computation of BB
        volumes = spam.label.volumes(threeCubedLabelVol)
        self.assertEqual(0, volumes[0])
        self.assertEqual(1, volumes[1])

        # with pre-computation of BB
        BB = spam.label.boundingBoxes(threeCubedLabelVol)
        volumes = spam.label.volumes(threeCubedLabelVol, boundingBoxes=BB)
        self.assertEqual(0, volumes[0])
        self.assertEqual(1, volumes[1])

    def test_equivalentRadii(self):
        # without pre-computation of BB or volumes
        radii = spam.label.equivalentRadii(threeCubedLabelVol)
        self.assertEqual(0.0, radii[0])
        self.assertEqual(((3.0 * 1.0) / (4.0 * numpy.pi))**(1.0 / 3.0), radii[1])

        volumes = spam.label.volumes(threeCubedLabelVol)
        radii = spam.label.equivalentRadii(threeCubedLabelVol, volumes=volumes)
        self.assertEqual(0.0, radii[0])
        self.assertEqual(((3.0 * 1.0) / (4.0 * numpy.pi))**(1.0 / 3.0), radii[1])

        BB = spam.label.boundingBoxes(threeCubedLabelVol)
        radii = spam.label.equivalentRadii(threeCubedLabelVol, boundingBoxes=BB)
        self.assertEqual(0.0, radii[0])
        self.assertEqual(((3.0 * 1.0) / (4.0 * numpy.pi))**(1.0 / 3.0), radii[1])

    def test_momentOfInertia(self):
        MOIval, MOIvec = spam.label.momentOfInertia(longSingle)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), MOIval[0].tolist())
        self.assertEqual(numpy.zeros(9, dtype='<f4').tolist(), MOIvec[0].tolist())
        self.assertTrue(MOIval[1, 0] > MOIval[1, 1])
        self.assertTrue(MOIval[1, 1] > MOIval[1, 2])
        self.assertTrue(MOIval[1, 2] > 0.0)
        self.assertEqual(numpy.array([0, 0, 1], dtype='<f4').tolist(), MOIvec[1, 0:3].tolist())
        self.assertEqual(numpy.array([0, 1, 0], dtype='<f4').tolist(), MOIvec[1, 3:6].tolist())
        self.assertEqual(numpy.array([1, 0, 0], dtype='<f4').tolist(), MOIvec[1, 6:9].tolist())

        BB = spam.label.boundingBoxes(longSingle)
        MOIval, MOIvec = spam.label.momentOfInertia(longSingle, boundingBoxes=BB)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), MOIval[0].tolist())
        self.assertEqual(numpy.zeros(9, dtype='<f4').tolist(), MOIvec[0].tolist())
        self.assertTrue(MOIval[1, 0] > MOIval[1, 1])
        self.assertTrue(MOIval[1, 1] > MOIval[1, 2])
        self.assertTrue(MOIval[1, 2] > 0.0)
        self.assertEqual(numpy.array([0, 0, 1], dtype='<f4').tolist(), MOIvec[1, 0:3].tolist())
        self.assertEqual(numpy.array([0, 1, 0], dtype='<f4').tolist(), MOIvec[1, 3:6].tolist())
        self.assertEqual(numpy.array([1, 0, 0], dtype='<f4').tolist(), MOIvec[1, 6:9].tolist())

        COM = spam.label.centresOfMass(longSingle, boundingBoxes=BB)
        MOIval, MOIvec = spam.label.momentOfInertia(longSingle, centresOfMass=COM)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), MOIval[0].tolist())
        self.assertEqual(numpy.zeros(9, dtype='<f4').tolist(), MOIvec[0].tolist())
        self.assertTrue(MOIval[1, 0] > MOIval[1, 1])
        self.assertTrue(MOIval[1, 1] > MOIval[1, 2])
        self.assertTrue(MOIval[1, 2] > 0.0)
        self.assertEqual(numpy.array([0, 0, 1], dtype='<f4').tolist(), MOIvec[1, 0:3].tolist())
        self.assertEqual(numpy.array([0, 1, 0], dtype='<f4').tolist(), MOIvec[1, 3:6].tolist())
        self.assertEqual(numpy.array([1, 0, 0], dtype='<f4').tolist(), MOIvec[1, 6:9].tolist())

    def test_ellipseAxes(self):
        BB = spam.label.boundingBoxes(longSingle)
        EA = spam.label.ellipseAxes(longSingle)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), EA[0].tolist())
        self.assertTrue(EA[1, 0] > EA[1, 1])
        self.assertTrue(EA[1, 1] > EA[1, 2])
        self.assertTrue(EA[1, 2] > 0.0)
        self.assertTrue(EA[1, 0] < (BB[1, 1] - BB[1, 0]))
        self.assertTrue(EA[1, 0] > (BB[1, 1] - BB[1, 0]) / 2.0)
        self.assertTrue(EA[1, 1] < (BB[1, 3] - BB[1, 2]))
        self.assertTrue(EA[1, 1] > (BB[1, 3] - BB[1, 2]) / 2.0)
        # self.assertTrue( EA[1,2]<(BB[1,5]-BB[1,4]))
        self.assertTrue(EA[1, 2] > (BB[1, 5] - BB[1, 4]) / 2.0)

        EA = spam.label.ellipseAxes(longSingle, enforceVolume=False)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), EA[0].tolist())
        self.assertTrue(EA[1, 0] > EA[1, 1])
        self.assertTrue(EA[1, 1] > EA[1, 2])
        self.assertTrue(EA[1, 2] > 0.0)
        self.assertTrue(EA[1, 0] < (BB[1, 1] - BB[1, 0]))
        self.assertTrue(EA[1, 0] > (BB[1, 1] - BB[1, 0]) / 2.0)
        self.assertTrue(EA[1, 1] < (BB[1, 3] - BB[1, 2]))
        self.assertTrue(EA[1, 1] > (BB[1, 3] - BB[1, 2]) / 2.0)
        # self.assertTrue( EA[1,2]<(BB[1,5]-BB[1,4]))
        self.assertTrue(EA[1, 2] > (BB[1, 5] - BB[1, 4]) / 2.0)

        MOIval = spam.label.momentOfInertia(longSingle)[0]
        EA = spam.label.ellipseAxes(longSingle, MOIeigenValues=MOIval)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), EA[0].tolist())
        self.assertTrue(EA[1, 0] > EA[1, 1])
        self.assertTrue(EA[1, 1] > EA[1, 2])
        self.assertTrue(EA[1, 2] > 0.0)
        self.assertTrue(EA[1, 0] < (BB[1, 1] - BB[1, 0]))
        self.assertTrue(EA[1, 0] > (BB[1, 1] - BB[1, 0]) / 2.0)
        self.assertTrue(EA[1, 1] < (BB[1, 3] - BB[1, 2]))
        self.assertTrue(EA[1, 1] > (BB[1, 3] - BB[1, 2]) / 2.0)
        # self.assertTrue( EA[1,2]<(BB[1,5]-BB[1,4]))
        self.assertTrue(EA[1, 2] > (BB[1, 5] - BB[1, 4]) / 2.0)

        volumes = spam.label.volumes(longSingle)
        EA = spam.label.ellipseAxes(longSingle, volumes=volumes)
        self.assertEqual(numpy.zeros(3, dtype='<f4').tolist(), EA[0].tolist())
        self.assertTrue(EA[1, 0] > EA[1, 1])
        self.assertTrue(EA[1, 1] > EA[1, 2])
        self.assertTrue(EA[1, 2] > 0.0)
        self.assertTrue(EA[1, 0] < (BB[1, 1] - BB[1, 0]))
        self.assertTrue(EA[1, 0] > (BB[1, 1] - BB[1, 0]) / 2.0)
        self.assertTrue(EA[1, 1] < (BB[1, 3] - BB[1, 2]))
        self.assertTrue(EA[1, 1] > (BB[1, 3] - BB[1, 2]) / 2.0)
        # self.assertTrue( EA[1,2]<(BB[1,5]-BB[1,4]))
        self.assertTrue(EA[1, 2] > (BB[1, 5] - BB[1, 4]) / 2.0)

    def test_convertLabelToFloat(self):
        f = spam.label.convertLabelToFloat(threeCubedLabelVol, numpy.array([0.0, numpy.pi]))
        self.assertAlmostEqual(f[1, 1, 1], float(numpy.pi), places=5)

    def test_makeLabelsSequential(self):
        threeCubedLabelVolDouble = threeCubedLabelVol.copy() * 2
        f = spam.label.makeLabelsSequential(threeCubedLabelVol)
        self.assertEqual(f[1, 1, 1], 1)

    def test_getLabel(self):
Edward Andò's avatar
Edward Andò committed
181
        # No-option get label
Edward Andò's avatar
Edward Andò committed
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
        gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1)
        self.assertEqual(numpy.array([[[1]]], dtype=bool), gottenLabel['subvol'].tolist())
        self.assertEqual(gottenLabel['slice'][0].start, 1)
        self.assertEqual(gottenLabel['slice'][0].stop,  2)
        self.assertEqual(gottenLabel['slice'][1].start, 1)
        self.assertEqual(gottenLabel['slice'][1].stop,  2)
        self.assertEqual(gottenLabel['slice'][2].start, 1)
        self.assertEqual(gottenLabel['slice'][2].stop,  2)

        self.assertEqual(gottenLabel['boundingBox'][0], 1)
        self.assertEqual(gottenLabel['boundingBox'][1], 1)
        self.assertEqual(gottenLabel['boundingBox'][2], 1)
        self.assertEqual(gottenLabel['boundingBox'][3], 1)
        self.assertEqual(gottenLabel['boundingBox'][4], 1)
        self.assertEqual(gottenLabel['boundingBox'][5], 1)

198
        COM = spam.label.centresOfMass(threeCubedLabelVol)
Edward Andò's avatar
Edward Andò committed
199
200
201
        self.assertEqual(COM[1].tolist(), gottenLabel['centreOfMassABS'].tolist())
        self.assertEqual((COM[1] - numpy.array((1, 1, 1))).tolist(), gottenLabel['centreOfMassREL'].tolist())
        self.assertEqual(gottenLabel['volumeInitial'], 1)
202
203

        BB = spam.label.boundingBoxes(threeCubedLabelVol)
Edward Andò's avatar
Edward Andò committed
204
        # With extract cube
Edward Andò's avatar
Edward Andò committed
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
        gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, extractCube=True)
        self.assertEqual(gottenLabel['sliceCube'][0].start, 1)
        self.assertEqual(gottenLabel['sliceCube'][0].stop,  2)
        self.assertEqual(gottenLabel['sliceCube'][1].start, 1)
        self.assertEqual(gottenLabel['sliceCube'][1].stop,  2)
        self.assertEqual(gottenLabel['sliceCube'][2].start, 1)

        self.assertEqual(gottenLabel['boundingBoxCube'][0], 1)
        self.assertEqual(gottenLabel['boundingBoxCube'][1], 1)
        self.assertEqual(gottenLabel['boundingBoxCube'][2], 1)
        self.assertEqual(gottenLabel['boundingBoxCube'][3], 1)
        self.assertEqual(gottenLabel['boundingBoxCube'][4], 1)
        self.assertEqual(gottenLabel['boundingBoxCube'][5], 1)

        # print gottenLabel['slice']
220
        COM = spam.label.centresOfMass(threeCubedLabelVol)
Edward Andò's avatar
Edward Andò committed
221
222
223
        self.assertEqual(COM[1].tolist(), gottenLabel['centreOfMassABS'].tolist())
        self.assertEqual((COM[1] - numpy.array((1, 1, 1))).tolist(), gottenLabel['centreOfMassREL'].tolist())
        self.assertEqual(gottenLabel['volumeInitial'], 1)
224

Edward Andò's avatar
Edward Andò committed
225
        # Asking for a label that is not there
Edward Andò's avatar
Edward Andò committed
226
227
        gottenLabel = spam.label.getLabel(threeCubedLabelVol, 100, boundingBoxes=BB, centresOfMass=COM, extractCube=True)
        self.assertEqual(gottenLabel, None)
Edward Andò's avatar
Edward Andò committed
228
229

        # Just run the remaining cases
Edward Andò's avatar
Edward Andò committed
230
231
232
233
234
        gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, extractCube=True, extractCubeSize=10)
        gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, extractCube=True, extractCubeSize=0)
        gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, margin=1, labelDilate=-2)
        gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, margin=2, labelDilate=2)
        volRef = gottenLabel['volumeInitial']
235
236
237
238

        # New case for "labelDilateMaskOtherLabels"
        threeCubedLabelVol2 = threeCubedLabelVol.copy()
        threeCubedLabelVol2[0,1,1] = 2
Edward Andò's avatar
Edward Andò committed
239
240
241
        gottenLabel = spam.label.getLabel(threeCubedLabelVol2, 1, boundingBoxes=BB, centresOfMass=COM, margin=2, labelDilate=2, labelDilateMaskOtherLabels=True)
        #print(gottenLabel['volume'], volRef)
        self.assertEqual(gottenLabel['volumeDilated'] > volRef, True)
Edward Andò's avatar
Edward Andò committed
242

243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
    def test_labelsOnEdges(self):
        labelled = numpy.zeros((5, 5, 5), dtype="<u2")
        labelled[1, 2, 1] = 1
        labelled[1, 1, 1] = 2
        labelled[3, 3, 3] = 2
        labelled[3, 3, 2] = 4
        labelled[1, 1, 2] = 4
        labelled[2, 2, 2] = 5
        labelled[0, 0, 0] = 6
        labelled[0, 4, 4] = 7

        loe = spam.label.labelsOnEdges(labelled)
        self.assertEqual(loe.tolist(), [0, 6, 7])

    def test_removeLabels(self):
        f = spam.label.removeLabels(longSingle, numpy.array([1]))
        self.assertEqual(f.max(), 0)

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
    #def test_feretDiameters(self):
        #import scipy.ndimage
        ## Very simple test for radii tested on a sphere
        #radius = 20
        #im = numpy.zeros((50, 50, 50), dtype='<f8')
        #spam.kalisphera.makeSphere(im, (25,25,25), radius)
        #im = im > 0.5
        #diams = spam.label.feretDiameters(im)[0][1]
        #self.assertAlmostEqual(diams[0], radius*2, delta=1)
        #self.assertAlmostEqual(diams[1], radius*2, delta=1)

        #imSquish = scipy.ndimage.zoom(im, (0.5, 1, 1), order=0)
        #imSquish = imSquish > 0.5
        #diams = spam.label.feretDiameters(imSquish)[0][1]
        #self.assertAlmostEqual(diams[0], radius*2, delta=1)
        #self.assertAlmostEqual(diams[1], radius*2/2, delta=1)

    def test_trueSphericity(self):
        def kaliToTS(v,r):
            vol = numpy.zeros((v,v,v), dtype='f8')
            spam.kalisphera.makeSphere(vol, (v/2,v/2,v/2), r)
            vol = vol > 0.5
            return spam.label.trueSphericity(vol, gaussianFilterSigma=0.75)[1]

        # This is too small, should be ignored
        self.assertEqual(kaliToTS(50, 2), 0.0)
        self.assertAlmostEqual(kaliToTS(50,  5), 1.0, places=1)
        self.assertAlmostEqual(kaliToTS(50, 10), 1.0, places=1)
        self.assertAlmostEqual(kaliToTS(50, 15), 1.0, places=1)

        # Make sure with a box
        vol = numpy.zeros((50,50,50), dtype=spam.label.labelType)
        vol[0:40, 0:40, 0:40] = 1
Edward Andò's avatar
Edward Andò committed
294
295
        self.assertAlmostEqual(spam.label.trueSphericity(vol, gaussianFilterSigma=0.5)[1],
                               (36.0*numpy.pi*(vol.sum()**2.0))**(1/3.0) / (6*vol.sum()**(2/3.0)), places=1)
296

Edward Andò's avatar
Edward Andò committed
297
298
299
    def test_loadDEMdata(self):
        import spam.datasets
        data = spam.datasets.loadDEMtouchingGrainsAndBranchVector()
300
301
302
303
304
305
306
307
308
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
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
#     def
#         print gl['subvol'].tolist()
#         print spam.label.ellipseAxes(longSingle, enforceVolume=False )
# class TestFonctionGet(unittest.TestCase):
#
#     def tearDown(self):
#         try:
#             pass
#         except OSError:
#             pass
#
#     def test_volumesCOM( self ):
#         self.assertEqual(0,0)
#
#         # Define test volume
#         testVolDims = ( 50, 50, 50 )
#         margin = 1
#         minRadius = 3
#
#         ### Kalisphera setup for random radius and centre
#         # Random parameters for sphere, making sure it fits in the box and is bigger than 3px:
#         radius = minRadius + ( numpy.random.rand() * ( ( numpy.min( testVolDims ) / 2.0 ) - minRadius ) * 0.8 )
#
#         #randomThree = numpy.array([0.5,0.5,0.5])
#         randomThree = numpy.random.rand(3)
#
#         # Generate random centre with margin
#         #      ---offset to corner---     --- three random numbers * difference between vol dims and 2*radius ---
#         centre = 3*[radius+margin]     +  ( randomThree * ( numpy.array( testVolDims ) - 2.0 * numpy.array( ( 3*[radius+margin] )  ) ) )
#         #print centre
#
#         # Create empy array
#         perfectSphere =  numpy.zeros( testVolDims, dtype="<f8" )
#
#         # Call kalisphera to fill in array
#         if VERBOSE: print( "\ttest_label: Creating greyscale sphere of radius {:0.2f} at {}".format( radius, centre ) )
#         kalisphera.run( perfectSphere, numpy.array( centre ).astype('<f8'), float(radius) )
#
#         ### Check kalisphera created greyscale volume
#         # Let's check the volume
#         theoreticalVolume = (4.0/3.0) * numpy.pi * radius**3.0
#         greyscaleVolume = perfectSphere.sum()
#         normGreyDifference = ( (theoreticalVolume-greyscaleVolume)/theoreticalVolume )
#
#         if VERBOSE: print( "\ttest_label: theoretical volume = {:0.2f}vx, created volume = {:0.2f}vx".format( theoreticalVolume, greyscaleVolume ) )
#         if VERBOSE: print( "\ttest_label: Difference = {:0.5f}%\n".format( 100.0*normGreyDifference ) )
#
#         ### Binary tests
#         # binarise this sphere at partial volume of 0.5s
#         binary = perfectSphere > 0.5
#         binaryVolume = binary.sum()
#         normBinDifference = ( (theoreticalVolume-binaryVolume)/theoreticalVolume )
#         if VERBOSE: print( "\ttest_label: binary volume = {:0.2f}vx".format( binaryVolume ) )
#         if VERBOSE: print( "\ttest_label: Difference = {:0.5f}%".format( 100.0*normBinDifference ) )
#         ltkVolume = ltk.getVolumes( binary )[1]
#         if VERBOSE: print( "\ttest_label: difference between binary sum and ltk.getVolumes(): {} (should be zero and integer)\n".format( binaryVolume - ltkVolume ) )
#
#         binaryCOM = ltk.getCentresOfMass( binary )[1]
#         if VERBOSE: print( "\ttest_label: difference between binary and imposed sphere centre Z,Y,X: {}".format( numpy.subtract( centre, binaryCOM  ) ) )
#
#     def test_boundingBoxes( self ):
#         ### new test: on a 3x3x3 array, test all single-pixel labels and check that label parameters
#         ###  i.e., boundingBoxes, COM, and volume are correctly labelled.
#         for z in [ 0,1,2 ]:
#             for y in [ 0,1,2 ]:
#                 for x in [ 0,1,2 ]:
#                     nugget = numpy.zeros( ( 3, 3, 3 ) )
#                     nugget[ z, y, x ] = 1
#
#                     bbo = ltk.getBoundingBoxes( nugget )[1]
#                     com = ltk.getCentresOfMass( nugget )[1]
#                     vol = ltk.getVolumes( nugget )[1]
#
#                     if ( bbo != numpy.array( [z,z,y,y,x,x] ) ).any():
#                         print "boundingBox mismatch", bbo
#
#                     if ( com != numpy.array( [z+0.5,y+0.5,x+0.5] )).any():
#                         print "com problem", com
#
#                     if ( vol != 1 ).any():
#                         print "vol != 1", vol
#
#     def test_missingLabel( self ):
#         nugget = numpy.zeros( ( 3, 3, 3 ) )
#         nugget[ 1, 1, 1 ] = 2
#         #print ltk.getBoundingBoxes( nugget )
#
#     def test_iterator( self ):
#         ### new test: on a 3x3x3 array, test all single-pixel labels and check that label parameters
#         ###  i.e., boundingBoxes, COM, and volume are correctly labelled.
#         labelled = numpy.zeros( (5,5,5), dtype="<u2" )
#         labelled[ 1,2,1 ] = 1
#         labelled[ 1,1,1 ] = 2
#         labelled[ 3,3,3 ] = 2
#         labelled[ 3,3,2 ] = 4
#         labelled[ 1,1,2 ] = 4
#         labelled[ 2,2,2 ] = 5
#
#         bbo = ltk.getBoundingBoxes( labelled )
#         com = ltk.getCentresOfMass( labelled )
#         vol = ltk.getVolumes( labelled )
#
#         #for label in range( 1, labelled.max()+3 ):
#         for label in [5]:
#             returns = ltk.getLabel( labelled, label, boundingBoxes = bbo, centresOfMass = com, extractCube=False, margin=1 )
#
#             if type(returns) == dict:
#                 for key in returns.keys():
#                     print key, returns[key]
#
#             #print returns['subvol'].shape
#             #print numpy.where( returns['subvol'] == label )
412

413
    def test_meanOrientation(self):
414
        # Generate random main direction
415
416
        theta = numpy.radians(random.randrange(0, 360, 1))
        phi = numpy.radians(random.randrange(0, 90, 1))
417
        # Generate n random vector near the main direction
418
419
420
421
        n = 1000 #number of vectors
        R = 10 #Radius of vectors
        vect = numpy.zeros((n, 3))
        for x in range(n):
422
423
            thetaI = numpy.random.normal(theta, numpy.radians(5))
            phiI = numpy.random.normal(phi, numpy.radians(5))
424
            RI = numpy.random.normal(R, 2)
425
426
427
428
429
430
            vect[x,0] = RI*numpy.cos(phiI)
            vect[x,1] = RI*numpy.sin(phiI)*numpy.sin(thetaI)
            vect[x,2] = RI*numpy.sin(phiI)*numpy.cos(thetaI)
        projectedVectors,mainAxis,_,_ = spam.label.meanOrientation(vect)
        phiPCA = numpy.degrees(numpy.arctan(numpy.sqrt(mainAxis[1]**2 + mainAxis[2]**2) / mainAxis[0]))
        thetaPCA = numpy.degrees(numpy.arctan2(mainAxis[1], mainAxis[2]))
Edward Andò's avatar
Edward Andò committed
431
        if thetaPCA<0: thetaPCA = 360 + thetaPCA 
432
        # Check if the difference between the angles is less than 5 degree
433
434
        self.assertLess(numpy.abs(phiPCA - numpy.degrees(phi)), 5)
        self.assertLess(numpy.abs(thetaPCA - numpy.degrees(theta)), 5)
435
        # Check if the vector are normalized
436
        normVectors = numpy.sum(numpy.linalg.norm(projectedVectors, axis = 1)) / len(projectedVectors)
437
        self.assertAlmostEqual(normVectors, 1, places = 3)
438

439
    def test_Spheroid(self):
440
        # Test for Oblate (Lentil-shaped)
Edward Andò's avatar
Edward Andò committed
441

442
        # Generate random direction
443
444
445
446
        theta = numpy.radians(random.randrange(0,360,1))
        phi = numpy.radians(random.randrange(0,90,1))
        vect = numpy.array([numpy.cos(phi), numpy.sin(phi)*numpy.sin(theta), numpy.sin(phi)*numpy.cos(theta)])
        if vect[0]<0: vect[:]=-1*vect[:] 
447
        # Generate two random semi-axis values
448
449
        a = random.randrange(20,40,1)
        c = random.randrange(10,20,1)
450
451
        #Generate the spheroid
        spheroid = spam.label.label.Spheroid(a, c, numpy.asarray(vect)).digitize()
452
        # Compute its semi-axis
453
        semiAxis = spam.label.ellipseAxes(spheroid)
454
        # Compare the semi-axis
455
456
        self.assertLess(numpy.abs(numpy.max(semiAxis[1])-numpy.maximum(a,c)),2) 
        self.assertLess(numpy.abs(numpy.min(semiAxis[1])-numpy.minimum(a,c)),2) 
457
        # Compute the orientation
458
459
        eigenVal, eigenVect = spam.label.momentOfInertia(spheroid)
        eigenVal = eigenVal / numpy.max(eigenVal)
460
        # Get main orientation
461
462
        mainVect = eigenVect[1][0:3]
        if mainVect[0]<0: mainVect[:]=-1*mainVect[:] 
463
        # Compute the angle between them
464
465
        c = numpy.dot(vect,mainVect)/numpy.linalg.norm(vect)/numpy.linalg.norm(mainVect)
        angle = numpy.degrees(numpy.arccos(numpy.clip(c, -1, 1)))
466
        # Check angle less than 2 degree
467
        self.assertLess(angle,2)
Edward Andò's avatar
Edward Andò committed
468

469
        # Test for Prolate (Rice-shaped)
Edward Andò's avatar
Edward Andò committed
470

471
        # Generate random direction
472
473
474
475
        theta = numpy.radians(random.randrange(0,360,1))
        phi = numpy.radians(random.randrange(0,90,1))
        vect = numpy.array([numpy.cos(phi), numpy.sin(phi)*numpy.sin(theta), numpy.sin(phi)*numpy.cos(theta)])
        if vect[0]<0: vect[:]=-1*vect[:] 
476
        # Generate two random semi-axis values
477
478
        a = random.randrange(10,20,1)
        b = random.randrange(20,40,1)
479
        # Generate the spheroid
480
        spheroid = spam.label.label.Spheroid(a, b, numpy.asarray(vect)).digitize()
481
        # Compute its semi-axis
482
        semiAxis = spam.label.ellipseAxes(spheroid)
483
        # Compare the semi-axis
484
485
        self.assertLess(numpy.abs(numpy.max(semiAxis[1])-numpy.maximum(a,b)),2) 
        self.assertLess(numpy.abs(numpy.min(semiAxis[1])-numpy.minimum(a,b)),2)
486
        # Compute the orientation
487
488
        eigenVal, eigenVect = spam.label.momentOfInertia(spheroid)
        eigenVal = eigenVal / numpy.max(eigenVal)
489
        # Get main orientation
490
491
        mainVect = eigenVect[1][6:9]
        if mainVect[0]<0: mainVect[:]=-1*mainVect[:] 
492
        # Compute the angle between them
493
494
        c = numpy.dot(vect,mainVect)/numpy.linalg.norm(vect)/numpy.linalg.norm(mainVect)
        angle = numpy.degrees(numpy.arccos(numpy.clip(c, -1, 1)))
495
        # Check angle less than 2 degree
496
        self.assertLess(angle,2)
Gustavo Pinzon's avatar
Gustavo Pinzon committed
497
        
498
        # Check that raises an error when the vector and dim are passed along
Gustavo Pinzon's avatar
Gustavo Pinzon committed
499
        with self.assertRaises(ValueError): spam.label.label.Spheroid(10, 20, numpy.asarray([0,1,0]), dim = 1).digitize()
500
        # Check that it runs even without a vector
Gustavo Pinzon's avatar
Gustavo Pinzon committed
501
502
        res = spam.label.label.Spheroid(10, 20, dim = 3).digitize()
        self.assertIsNotNone(res)
Edward Andò's avatar
Edward Andò committed
503

504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
    def _test_FixUnderSegmentation(self):
        # Generate two prolate grains (rice-like)
        grain1 = spam.label.label.Spheroid(10, 20, numpy.asarray([0,1,0])).digitize()
        grain2 = spam.label.label.Spheroid(10, 20, numpy.asarray([0,1,0])).digitize()
        # Create the bigger  labelled image
        grainIm = numpy.concatenate((grain1,grain2))
        grainIm = numpy.zeros(grainIm.shape)
        # Add the grains to the bigger image
        grainIm[:grain1.shape[0]-1,:,:] = grain1[:grain1.shape[0]-1,:,:]
        grainIm[grain2.shape[0]-5:-5,:,:] =  grainIm[grain2.shape[0]-5:-5,:,:] + grain1[:,:,:]
        # Set all the labels to 1
        grainIm = numpy.where(grainIm >= 1, 3, grainIm)
        # Pad a border
        grainIm = numpy.pad(grainIm, pad_width=10, mode='constant', constant_values = 0)
        # Create the 'greyScale' image
        greyIm = numpy.where(grainIm == 3, 30000, 10000)

        # Check that the greyscale image is normalized
        res1 = spam.label.label.fixUndersegmentation(grainIm, greyIm, [3], 10*0.5, 20*0.5)
        self.assertEqual(res1, None)
        # Check that a or c is valid
        res2 = spam.label.label.fixUndersegmentation(grainIm, greyIm, [3], numpy.nan, numpy.nan)
        self.assertEqual(res2, None)
        # Check that a or c are positive
        res3 = spam.label.label.fixUndersegmentation(grainIm, greyIm, [3], -1, 10)
        self.assertEqual(res3, None)

        # Run fixUnderSegmentation
        greyIm = numpy.where(grainIm == 3, 0.75, 0.25)
        res4 = spam.label.label.fixUndersegmentation(grainIm, greyIm, [3], 10*0.8, 20*0.8, vect = [[0,1,0]])
        # Check that there are two grains
        self.assertEqual(numpy.max(numpy.unique(res4)), 2)
        # Check that it runs even if the label does not exist
        res5 = spam.label.label.fixUndersegmentation(grainIm, greyIm, [3], 10*0.8, 20*0.8, vect = [[0,1,0]])
        self.assertIsNotNone(res5)
        # Check for a vect that is not a list
        res6 = spam.label.label.fixUndersegmentation(grainIm, greyIm, [3], 10*0.8, 20*0.8, vect = (0,1,0))
        self.assertIsNone(res6)
        # Check that it works without the input vect
        res7 = spam.label.label.fixUndersegmentation(grainIm, greyIm, [3], 10*0.8, 20*0.8, numVect = 1)
        self.assertIsNotNone(res7)
        # Check that it works with verbose = True
        res8 = spam.label.label.fixUndersegmentation(grainIm, greyIm, [3], 10*0.8, 20*0.8, numVect = 1, verbose = True)
        self.assertIsNotNone(res8)
        # Check that it works even for a non-existing label
        res9 = spam.label.label.fixUndersegmentation(grainIm, greyIm, [5], 10*0.8, 20*0.8, numVect = 1)
        self.assertIsNotNone(res9)
Gustavo Pinzon's avatar
Gustavo Pinzon committed
551
        
Edward Andò's avatar
Edward Andò committed
552

553
    def test_contactOrientations(self):
554
        # Generate two oblate grains
555
556
557
558
559
560
561
        grain1 = spam.label.label.Spheroid(7, 15, numpy.asarray([0,1,0])).digitize()
        grain2 = spam.label.label.Spheroid(7, 15, numpy.asarray([0,1,0])).digitize()
        grainIm = numpy.concatenate((grain1,grain2))
        grainIm = numpy.zeros(grainIm.shape)
        grainIm[:grain1.shape[0]-1,:,:] = grain1[:grain1.shape[0]-1,:,:]
        grainIm[grain2.shape[0]-5:-5,:,:] =  grainIm[grain2.shape[0]-5:-5,:,:] + grain2[:,:,:]
        grainIm = numpy.where(grainIm >1,1,grainIm)
562
        # Generate a new grain and create contact
563
564
565
566
567
568
        grain3 = 2*spam.label.label.Spheroid(7, 15, numpy.asarray([0,1,0])).digitize()
        imLab = numpy.concatenate((grain3,grainIm))
        imLab = numpy.zeros(imLab.shape)
        imLab[:grain3.shape[0]-1,:,:] = grain3[:grain3.shape[0]-1,:,:]
        imLab[grain3.shape[0]-3:-3,:,:] = imLab[grain3.shape[0]-3:-3,:,:] + grainIm[:,:,:]
        imLab = numpy.where(imLab >2, 2, imLab)
569
        # Create bin image
570
        imBin = numpy.where(imLab >0, 1, 0)
571
        # Run ITK and check angle of contact
572
573
574
575
576
577
        contactNormal, intervox, NotTreatedContact = con.contactOrientations(imBin, imLab, watershed="ITK")
        if contactNormal[0]<0: contactNormal = contactNormal*-1
        c = numpy.dot(contactNormal, [1, 0, 0])/numpy.linalg.norm(contactNormal)/numpy.linalg.norm([1, 0, 0])
        angle = numpy.degrees(numpy.arccos(numpy.clip(c, -1, 1)))
        self.assertFalse(NotTreatedContact)
        self.assertLess(angle, 1.)
578
        # Run RW and check angle of contact
579
580
581
582
583
584
585
        contactNormal, intervox, NotTreatedContact = con.contactOrientations(imBin, imLab, watershed="RW")
        if contactNormal[0]<0: contactNormal = contactNormal*-1
        c = numpy.dot(contactNormal, [1, 0, 0])/numpy.linalg.norm(contactNormal)/numpy.linalg.norm([1, 0, 0])
        angle = numpy.degrees(numpy.arccos(numpy.clip(c, -1, 1)))
        self.assertFalse(NotTreatedContact)
        self.assertLess(angle, 1.)

Gustavo Pinzon's avatar
Gustavo Pinzon committed
586
    def test_fabricTensor(self):
587
        # Define number of vector
Gustavo Pinzon's avatar
Gustavo Pinzon committed
588
        n = 1000
589
        # Create equally-spaced vector
Gustavo Pinzon's avatar
Gustavo Pinzon committed
590
        vectEq = spam.plotting.orientationPlotter.SaffAndKuijlaarsSpiral(n)
591
        # Compute fabric
Gustavo Pinzon's avatar
Gustavo Pinzon committed
592
        NEq,FEq,aEq = spam.label.label.fabricTensor(vectEq)
593
        # Create a set of vectors using a normal distribution centered on a random orientation
Gustavo Pinzon's avatar
Gustavo Pinzon committed
594
595
596
597
598
599
600
601
602
        theta = numpy.radians(random.randrange(0, 360, 1))
        phi = numpy.radians(random.randrange(0, 90, 1))
        vectRand = numpy.zeros((n, 3))
        for x in range(n):
            thetaI = numpy.random.normal(theta, numpy.radians(5))
            phiI = numpy.random.normal(phi, numpy.radians(5))
            vectRand[x,0] = numpy.cos(phiI)
            vectRand[x,1] = numpy.sin(phiI)*numpy.sin(thetaI)
            vectRand[x,2] = numpy.sin(phiI)*numpy.cos(thetaI)
603
        # Compute fabric
Gustavo Pinzon's avatar
Gustavo Pinzon committed
604
        NRand, FRand, aRand = spam.label.label.fabricTensor(vectRand)
605
        # Test the results
Gustavo Pinzon's avatar
Gustavo Pinzon committed
606
607
        self.assertGreater(aRand, aEq)

Edward Andò's avatar
Edward Andò committed
608
609
610
611
    def test_setVoronoi(self):
        box = numpy.zeros((100, 100, 100), dtype='<f8')
        spam.kalisphera.makeSphere(box, [[50, 50, 50-20], [50, 50, 50+20]], [20, 20])
        box = box > 0.5
612
        lab = spam.label.watershed(box, verbose=True)
Edward Andò's avatar
Edward Andò committed
613
614
615
616
617
618
        labVolumes = spam.label.volumes(lab)

        setVoronoi = spam.label.setVoronoi(lab)
        setVoronoiVolumes = spam.label.volumes(setVoronoi)
        self.assertEqual(setVoronoiVolumes[1]>labVolumes[1], True)
        self.assertEqual(setVoronoiVolumes[2]>labVolumes[2], True)
619

620
    def test_convexVolume(self):
621
        # Get a random separation between 1 and 10
622
        dist = int(10*numpy.random.rand())
623
        # Create a rectangular image
624
        im = numpy.zeros((15,15,30))
625
        # Create a first block of 5x5x5
626
        im[5:10,5:10,5:10] = 1
627
        # Create second block of 5x5x5 located at dist
628
        im[5:10,5:10,10+dist:15+dist] = 1
629
        # Compute convex volume
630
631
        im = im.astype('<u4')
        convexVol = spam.label.convexVolume(im)
632
        # Compute theoretical volume
633
634
        volOr = numpy.sum(im) + 5*5*dist
        self.assertEqual(volOr, convexVol[-1])
635
        # Check that it shows the error for particles smaller than 3 voxels
636
637
        im = numpy.zeros((15,15,30))
        im[10,10,10] = 1
638
        # Compute convex volume
639
640
641
        im = im.astype('<u4')
        convexVol = spam.label.convexVolume(im)
        self.assertEqual(0, convexVol[-1])
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
704
705
706
707
708
709
710
        
    def test_moveLabels(self):
        # Create the sphere
        imLab = skimage.morphology.ball(20)
        # Pad with zero at each boundaries
        imLab = numpy.pad(imLab, (10), 'constant', constant_values=(0))
        # Compute initial volume
        iniVol = spam.label.volumes(imLab)[-1]
        # Compute initial COM
        iniCOM = spam.label.centresOfMass(imLab)[-1]
        # Compute boundingBoxes
        boundingBoxes = spam.label.boundingBoxes(imLab)
        # Compute centre of Mass
        centresOfMass = spam.label.centresOfMass(imLab)
        # Create empty PhiField
        PhiField = numpy.zeros((2, 4, 4))

        # Test #1 -> COM=iniCOM, vol = iniVol for no dilate and Phi = I
        transformation = {'t': [0, 0, 0]}
        PhiField[1] = spam.deformation.computePhi(transformation)
        imLab2 = spam.label.moveLabels(imLab, PhiField, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
        # Compute Volume and COM
        newVol = spam.label.volumes(imLab2)[-1]
        newCOM = spam.label.centresOfMass(imLab2)[-1]
        self.assertEqual(0, numpy.sum(iniVol - newVol))
        self.assertEqual(0, numpy.sum(iniCOM - newCOM))
        # Test #2 -> COM=iniCOM, vol = +10% volIni for a dilation of 10%
        transformation = {'z': [1.1, 1, 1]}
        PhiField[1] = spam.deformation.computePhi(transformation)
        imLab2 = spam.label.moveLabels(imLab, PhiField, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
        # Compute Volume and COM
        newVol = spam.label.volumes(imLab2)[-1]
        newCOM = spam.label.centresOfMass(imLab2)[-1]
        self.assertEqual(0, numpy.sum(iniCOM - newCOM))
        self.assertAlmostEqual(1.1, newVol/iniVol, places = 1)
        # Test #3 -> COM=iniCOM and vol/ iniVol > 1 for dilate = 1 and Phi = I
        transformation = {'t': [0, 0, 0]}
        PhiField[1] = spam.deformation.computePhi(transformation)
        imLab2 = spam.label.moveLabels(imLab, PhiField, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass, labelDilate=1)
        # Compute Volume and COM
        newVol = spam.label.volumes(imLab2)[-1]
        newCOM = spam.label.centresOfMass(imLab2)[-1]
        self.assertEqual(0, numpy.sum(iniCOM - newCOM))
        self.assertGreater(newVol/iniVol, 1.0)
        # Test#4 -> vol=iniVol and COM follows the 5vox displacement
        transformation = {'t': [5, 5, 5]}
        PhiField[1] = spam.deformation.computePhi(transformation)
        imLab2 = spam.label.moveLabels(imLab, PhiField, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
        # Compute Volume and COM
        newVol = spam.label.volumes(imLab2)[-1]
        newCOM = spam.label.centresOfMass(imLab2)[-1]
        self.assertEqual(0, numpy.sum(iniVol - newVol))
        self.assertTrue(((newCOM - iniCOM) == 5).all())
        # Test 5 -> move grain half out and check volume is near 0.5
        transformation = {'t': [30, 0, 0,]}
        PhiField[1] = spam.deformation.computePhi(transformation)
        imLab2 = spam.label.moveLabels(imLab, PhiField, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
        # Compute Volume and COM
        newVol = spam.label.volumes(imLab2)[-1]
        newCOM = spam.label.centresOfMass(imLab2)[-1]
        self.assertAlmostEqual(0.5, newVol/iniVol, places = 1)
        # Test 6A -> Check that it runs when the PhiField has less labels than the labelled image
        PhiField = numpy.zeros((1, 4, 4))
        imLab2 = spam.label.moveLabels(imLab, PhiField, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
        self.assertIsNotNone(imLab2)
        # Test 6B -> Check that it runs when the PhiField has more labels than the labelled image
        PhiField = numpy.zeros((3, 4, 4))
        imLab2 = spam.label.moveLabels(imLab, PhiField, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
        self.assertIsNotNone(imLab2)
711
        
712
713
714
715
716
717
718
719
720
721
722
723
724
        # Test 7A -> Check that it works if the return status is equal to 2
        transformation = {'t': [5, 5, 5]}
        PhiField[1] = spam.deformation.computePhi(transformation)
        imLab2 = spam.label.moveLabels(imLab, PhiField, returnStatus = numpy.array([0, 2]),boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
        # Compute Volume and COM
        newVol = spam.label.volumes(imLab2)[-1]
        newCOM = spam.label.centresOfMass(imLab2)[-1]
        self.assertEqual(0, numpy.sum(iniVol - newVol))
        self.assertTrue(((newCOM - iniCOM) == 5).all())
        # Test 7B -> Check that it works if the return status is different to 2
        imLab2 = spam.label.moveLabels(imLab, PhiField, returnStatus = numpy.array([0, 0]),boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
        self.assertEqual(1, len(numpy.unique(imLab2)))
        
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
    def test_erodeLabels(self):
        # Create the sphere
        imLab = skimage.morphology.ball(20)
        # Pad with zero at each boundaries
        imLab = numpy.pad(imLab, (10), 'constant', constant_values=(0))
        # Compute initial volume
        iniVol = spam.label.volumes(imLab)[-1]
        # Compute initial COM
        iniCOM = spam.label.centresOfMass(imLab)[-1]
        # Test #1 -> COM=iniCOM, vol = iniVol for no dilate and Phi = I
        imLab2 = spam.label.erodeLabels(imLab, erosion=0)
        # Compute Volume and COM
        newVol = spam.label.volumes(imLab2)[-1]
        newCOM = spam.label.centresOfMass(imLab2)[-1]
        self.assertEqual(0, numpy.sum(iniVol - newVol))
        self.assertEqual(0, numpy.sum(iniCOM - newCOM))
        # Test #2 -> COM=iniCOM and vol/ iniVol < 1 for dilate = 1 and Phi = I
        imLab2 = spam.label.erodeLabels(imLab, erosion=1)
        # Compute Volume and COM
        newVol = spam.label.volumes(imLab2)[-1]
        newCOM = spam.label.centresOfMass(imLab2)[-1]
        self.assertEqual(0, numpy.sum(iniCOM - newCOM))
        self.assertGreater(1.0, newVol/iniVol)
748
if __name__ == '__main__':
749
    unittest.main()