scikit-image

1
import math
2

3
import numpy as np
4
import pytest
5
from numpy.testing import assert_almost_equal
6

7
from skimage import feature
8
from skimage.draw import disk
9
from skimage.draw.draw3d import ellipsoid
10
from skimage.feature import blob_dog, blob_doh, blob_log
11
from skimage.feature.blob import _blob_overlap
12

13

14
@pytest.mark.parametrize('dtype', [np.uint8, np.float16, np.float32, np.float64])
15
@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
16
def test_blob_dog(dtype, threshold_type):
17
    r2 = math.sqrt(2)
18
    img = np.ones((512, 512), dtype=dtype)
19

20
    xs, ys = disk((400, 130), 5)
21
    img[xs, ys] = 255
22

23
    xs, ys = disk((100, 300), 25)
24
    img[xs, ys] = 255
25

26
    xs, ys = disk((200, 350), 45)
27
    img[xs, ys] = 255
28

29
    if threshold_type == 'absolute':
30
        threshold = 2.0
31
        if img.dtype.kind != 'f':
32
            # account for internal scaling to [0, 1] by img_as_float
33
            threshold /= np.ptp(img)
34
        threshold_rel = None
35
    elif threshold_type == 'relative':
36
        threshold = None
37
        threshold_rel = 0.5
38

39
    blobs = blob_dog(
40
        img,
41
        min_sigma=4,
42
        max_sigma=50,
43
        threshold=threshold,
44
        threshold_rel=threshold_rel,
45
    )
46

47
    def radius(x):
48
        return r2 * x[2]
49

50
    s = sorted(blobs, key=radius)
51
    thresh = 5
52
    ratio_thresh = 0.25
53

54
    b = s[0]
55
    assert abs(b[0] - 400) <= thresh
56
    assert abs(b[1] - 130) <= thresh
57
    assert abs(radius(b) - 5) <= ratio_thresh * 5
58

59
    b = s[1]
60
    assert abs(b[0] - 100) <= thresh
61
    assert abs(b[1] - 300) <= thresh
62
    assert abs(radius(b) - 25) <= ratio_thresh * 25
63

64
    b = s[2]
65
    assert abs(b[0] - 200) <= thresh
66
    assert abs(b[1] - 350) <= thresh
67
    assert abs(radius(b) - 45) <= ratio_thresh * 45
68

69
    # Testing no peaks
70
    img_empty = np.zeros((100, 100), dtype=dtype)
71
    assert blob_dog(img_empty).size == 0
72

73

74
@pytest.mark.parametrize('dtype', [np.uint8, np.float16, np.float32, np.float64])
75
@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
76
def test_blob_dog_3d(dtype, threshold_type):
77
    # Testing 3D
78
    r = 10
79
    pad = 10
80
    im3 = ellipsoid(r, r, r)
81
    im3 = np.pad(im3, pad, mode='constant')
82

83
    if threshold_type == 'absolute':
84
        threshold = 0.001
85
        threshold_rel = 0
86
    elif threshold_type == 'relative':
87
        threshold = 0
88
        threshold_rel = 0.5
89

90
    blobs = blob_dog(
91
        im3,
92
        min_sigma=3,
93
        max_sigma=10,
94
        sigma_ratio=1.2,
95
        threshold=threshold,
96
        threshold_rel=threshold_rel,
97
    )
98
    b = blobs[0]
99

100
    assert b.shape == (4,)
101
    assert b[0] == r + pad + 1
102
    assert b[1] == r + pad + 1
103
    assert b[2] == r + pad + 1
104
    assert abs(math.sqrt(3) * b[3] - r) < 1.1
105

106

107
@pytest.mark.parametrize('dtype', [np.uint8, np.float16, np.float32, np.float64])
108
@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
109
def test_blob_dog_3d_anisotropic(dtype, threshold_type):
110
    # Testing 3D anisotropic
111
    r = 10
112
    pad = 10
113
    im3 = ellipsoid(r / 2, r, r)
114
    im3 = np.pad(im3, pad, mode='constant')
115

116
    if threshold_type == 'absolute':
117
        threshold = 0.001
118
        threshold_rel = None
119
    elif threshold_type == 'relative':
120
        threshold = None
121
        threshold_rel = 0.5
122

123
    blobs = blob_dog(
124
        im3.astype(dtype, copy=False),
125
        min_sigma=[1.5, 3, 3],
126
        max_sigma=[5, 10, 10],
127
        sigma_ratio=1.2,
128
        threshold=threshold,
129
        threshold_rel=threshold_rel,
130
    )
131
    b = blobs[0]
132

133
    assert b.shape == (6,)
134
    assert b[0] == r / 2 + pad + 1
135
    assert b[1] == r + pad + 1
136
    assert b[2] == r + pad + 1
137
    assert abs(math.sqrt(3) * b[3] - r / 2) < 1.1
138
    assert abs(math.sqrt(3) * b[4] - r) < 1.1
139
    assert abs(math.sqrt(3) * b[5] - r) < 1.1
140

141

142
@pytest.mark.parametrize("disc_center", [(5, 5), (5, 20)])
143
@pytest.mark.parametrize("exclude_border", [6, (6, 6), (4, 15)])
144
def test_blob_dog_exclude_border(disc_center, exclude_border):
145
    # Testing exclude border
146

147
    # image where blob is disc_center px from borders, radius 5
148
    img = np.ones((512, 512))
149
    xs, ys = disk(disc_center, 5)
150
    img[xs, ys] = 255
151
    blobs = blob_dog(
152
        img,
153
        min_sigma=1.5,
154
        max_sigma=5,
155
        sigma_ratio=1.2,
156
    )
157
    assert blobs.shape[0] == 1, "one blob should have been detected"
158
    b = blobs[0]
159
    assert b[0] == disc_center[0], f"blob should be {disc_center[0]} px from x border"
160
    assert b[1] == disc_center[1], f"blob should be {disc_center[1]} px from y border"
161

162
    blobs = blob_dog(
163
        img,
164
        min_sigma=1.5,
165
        max_sigma=5,
166
        sigma_ratio=1.2,
167
        exclude_border=exclude_border,
168
    )
169

170
    if disc_center == (5, 20) and exclude_border == (4, 15):
171
        assert blobs.shape[0] == 1, "one blob should have been detected"
172
        b = blobs[0]
173
        assert (
174
            b[0] == disc_center[0]
175
        ), f"blob should be {disc_center[0]} px from x border"
176
        assert (
177
            b[1] == disc_center[1]
178
        ), f"blob should be {disc_center[1]} px from y border"
179
    else:
180
        msg = "zero blobs should be detected, as only blob is 5 px from border"
181
        assert blobs.shape[0] == 0, msg
182

183

184
@pytest.mark.parametrize('anisotropic', [False, True])
185
@pytest.mark.parametrize('ndim', [1, 2, 3, 4])
186
@pytest.mark.parametrize('function_name', ['blob_dog', 'blob_log'])
187
def test_nd_blob_no_peaks_shape(function_name, ndim, anisotropic):
188
    # uniform image so no blobs will be found
189
    z = np.zeros((16,) * ndim, dtype=np.float32)
190
    if anisotropic:
191
        max_sigma = 8 + np.arange(ndim)
192
    else:
193
        max_sigma = 8
194
    blob_func = getattr(feature, function_name)
195
    blobs = blob_func(z, max_sigma=max_sigma)
196
    # z.ndim +  (z.ndim sigmas if anisotropic, only one sigma otherwise)
197
    expected_shape = 2 * z.ndim if anisotropic else z.ndim + 1
198
    assert blobs.shape == (0, expected_shape)
199

200

201
@pytest.mark.parametrize('dtype', [np.uint8, np.float16, np.float32, np.float64])
202
@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
203
def test_blob_log(dtype, threshold_type):
204
    r2 = math.sqrt(2)
205
    img = np.ones((256, 256), dtype=dtype)
206

207
    xs, ys = disk((200, 65), 5)
208
    img[xs, ys] = 255
209

210
    xs, ys = disk((80, 25), 15)
211
    img[xs, ys] = 255
212

213
    xs, ys = disk((50, 150), 25)
214
    img[xs, ys] = 255
215

216
    xs, ys = disk((100, 175), 30)
217
    img[xs, ys] = 255
218

219
    if threshold_type == 'absolute':
220
        threshold = 1
221
        if img.dtype.kind != 'f':
222
            # account for internal scaling to [0, 1] by img_as_float
223
            threshold /= np.ptp(img)
224
        threshold_rel = None
225
    elif threshold_type == 'relative':
226
        threshold = None
227
        threshold_rel = 0.5
228

229
    blobs = blob_log(
230
        img, min_sigma=5, max_sigma=20, threshold=threshold, threshold_rel=threshold_rel
231
    )
232

233
    def radius(x):
234
        return r2 * x[2]
235

236
    s = sorted(blobs, key=radius)
237
    thresh = 3
238

239
    b = s[0]
240
    assert abs(b[0] - 200) <= thresh
241
    assert abs(b[1] - 65) <= thresh
242
    assert abs(radius(b) - 5) <= thresh
243

244
    b = s[1]
245
    assert abs(b[0] - 80) <= thresh
246
    assert abs(b[1] - 25) <= thresh
247
    assert abs(radius(b) - 15) <= thresh
248

249
    b = s[2]
250
    assert abs(b[0] - 50) <= thresh
251
    assert abs(b[1] - 150) <= thresh
252
    assert abs(radius(b) - 25) <= thresh
253

254
    b = s[3]
255
    assert abs(b[0] - 100) <= thresh
256
    assert abs(b[1] - 175) <= thresh
257
    assert abs(radius(b) - 30) <= thresh
258

259
    # Testing log scale
260
    blobs = blob_log(
261
        img,
262
        min_sigma=5,
263
        max_sigma=20,
264
        threshold=threshold,
265
        threshold_rel=threshold_rel,
266
        log_scale=True,
267
    )
268

269
    b = s[0]
270
    assert abs(b[0] - 200) <= thresh
271
    assert abs(b[1] - 65) <= thresh
272
    assert abs(radius(b) - 5) <= thresh
273

274
    b = s[1]
275
    assert abs(b[0] - 80) <= thresh
276
    assert abs(b[1] - 25) <= thresh
277
    assert abs(radius(b) - 15) <= thresh
278

279
    b = s[2]
280
    assert abs(b[0] - 50) <= thresh
281
    assert abs(b[1] - 150) <= thresh
282
    assert abs(radius(b) - 25) <= thresh
283

284
    b = s[3]
285
    assert abs(b[0] - 100) <= thresh
286
    assert abs(b[1] - 175) <= thresh
287
    assert abs(radius(b) - 30) <= thresh
288

289
    # Testing no peaks
290
    img_empty = np.zeros((100, 100))
291
    assert blob_log(img_empty).size == 0
292

293

294
def test_blob_log_no_warnings():
295
    img = np.ones((11, 11))
296

297
    xs, ys = disk((5, 5), 2)
298
    img[xs, ys] = 255
299

300
    xs, ys = disk((7, 6), 2)
301
    img[xs, ys] = 255
302

303
    blob_log(img, max_sigma=20, num_sigma=10, threshold=0.1)
304

305

306
def test_blob_log_3d():
307
    # Testing 3D
308
    r = 6
309
    pad = 10
310
    im3 = ellipsoid(r, r, r)
311
    im3 = np.pad(im3, pad, mode='constant')
312

313
    blobs = blob_log(im3, min_sigma=3, max_sigma=10)
314
    b = blobs[0]
315

316
    assert b.shape == (4,)
317
    assert b[0] == r + pad + 1
318
    assert b[1] == r + pad + 1
319
    assert b[2] == r + pad + 1
320
    assert abs(math.sqrt(3) * b[3] - r) < 1
321

322

323
def test_blob_log_3d_anisotropic():
324
    # Testing 3D anisotropic
325
    r = 6
326
    pad = 10
327
    im3 = ellipsoid(r / 2, r, r)
328
    im3 = np.pad(im3, pad, mode='constant')
329

330
    blobs = blob_log(
331
        im3,
332
        min_sigma=[1, 2, 2],
333
        max_sigma=[5, 10, 10],
334
    )
335

336
    b = blobs[0]
337
    assert b.shape == (6,)
338
    assert b[0] == r / 2 + pad + 1
339
    assert b[1] == r + pad + 1
340
    assert b[2] == r + pad + 1
341
    assert abs(math.sqrt(3) * b[3] - r / 2) < 1
342
    assert abs(math.sqrt(3) * b[4] - r) < 1
343
    assert abs(math.sqrt(3) * b[5] - r) < 1
344

345

346
@pytest.mark.parametrize("disc_center", [(5, 5), (5, 20)])
347
@pytest.mark.parametrize("exclude_border", [6, (6, 6), (4, 15)])
348
def test_blob_log_exclude_border(disc_center, exclude_border):
349
    # image where blob is disc_center px from borders, radius 5
350
    img = np.ones((512, 512))
351
    xs, ys = disk(disc_center, 5)
352
    img[xs, ys] = 255
353

354
    blobs = blob_log(
355
        img,
356
        min_sigma=1.5,
357
        max_sigma=5,
358
    )
359
    assert blobs.shape[0] == 1
360
    b = blobs[0]
361
    assert b[0] == disc_center[0], f"blob should be {disc_center[0]} px from x border"
362
    assert b[1] == disc_center[1], f"blob should be {disc_center[1]} px from y border"
363

364
    blobs = blob_log(
365
        img,
366
        min_sigma=1.5,
367
        max_sigma=5,
368
        exclude_border=exclude_border,
369
    )
370

371
    if disc_center == (5, 20) and exclude_border == (4, 15):
372
        assert blobs.shape[0] == 1, "one blob should have been detected"
373
        b = blobs[0]
374
        assert (
375
            b[0] == disc_center[0]
376
        ), f"blob should be {disc_center[0]} px from x border"
377
        assert (
378
            b[1] == disc_center[1]
379
        ), f"blob should be {disc_center[1]} px from y border"
380
    else:
381
        msg = "zero blobs should be detected, as only blob is 5 px from border"
382
        assert blobs.shape[0] == 0, msg
383

384

385
@pytest.mark.parametrize("dtype", [np.uint8, np.float16, np.float32])
386
@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
387
def test_blob_doh(dtype, threshold_type):
388
    img = np.ones((512, 512), dtype=dtype)
389

390
    xs, ys = disk((400, 130), 20)
391
    img[xs, ys] = 255
392

393
    xs, ys = disk((460, 50), 30)
394
    img[xs, ys] = 255
395

396
    xs, ys = disk((100, 300), 40)
397
    img[xs, ys] = 255
398

399
    xs, ys = disk((200, 350), 50)
400
    img[xs, ys] = 255
401

402
    if threshold_type == 'absolute':
403
        # Note: have to either scale up threshold or rescale the image to the
404
        #       range [0, 1] internally.
405
        threshold = 0.05
406
        if img.dtype.kind == 'f':
407
            # account for lack of internal scaling to [0, 1] by img_as_float
408
            ptp = np.ptp(img)
409
            threshold *= ptp**2
410
        threshold_rel = None
411
    elif threshold_type == 'relative':
412
        threshold = None
413
        threshold_rel = 0.5
414

415
    blobs = blob_doh(
416
        img,
417
        min_sigma=1,
418
        max_sigma=60,
419
        num_sigma=10,
420
        threshold=threshold,
421
        threshold_rel=threshold_rel,
422
    )
423

424
    def radius(x):
425
        return x[2]
426

427
    s = sorted(blobs, key=radius)
428
    thresh = 4
429

430
    b = s[0]
431
    assert abs(b[0] - 400) <= thresh
432
    assert abs(b[1] - 130) <= thresh
433
    assert abs(radius(b) - 20) <= thresh
434

435
    b = s[1]
436
    assert abs(b[0] - 460) <= thresh
437
    assert abs(b[1] - 50) <= thresh
438
    assert abs(radius(b) - 30) <= thresh
439

440
    b = s[2]
441
    assert abs(b[0] - 100) <= thresh
442
    assert abs(b[1] - 300) <= thresh
443
    assert abs(radius(b) - 40) <= thresh
444

445
    b = s[3]
446
    assert abs(b[0] - 200) <= thresh
447
    assert abs(b[1] - 350) <= thresh
448
    assert abs(radius(b) - 50) <= thresh
449

450

451
def test_blob_doh_log_scale():
452
    img = np.ones((512, 512), dtype=np.uint8)
453

454
    xs, ys = disk((400, 130), 20)
455
    img[xs, ys] = 255
456

457
    xs, ys = disk((460, 50), 30)
458
    img[xs, ys] = 255
459

460
    xs, ys = disk((100, 300), 40)
461
    img[xs, ys] = 255
462

463
    xs, ys = disk((200, 350), 50)
464
    img[xs, ys] = 255
465

466
    blobs = blob_doh(
467
        img, min_sigma=1, max_sigma=60, num_sigma=10, log_scale=True, threshold=0.05
468
    )
469

470
    def radius(x):
471
        return x[2]
472

473
    s = sorted(blobs, key=radius)
474
    thresh = 10
475

476
    b = s[0]
477
    assert abs(b[0] - 400) <= thresh
478
    assert abs(b[1] - 130) <= thresh
479
    assert abs(radius(b) - 20) <= thresh
480

481
    b = s[2]
482
    assert abs(b[0] - 460) <= thresh
483
    assert abs(b[1] - 50) <= thresh
484
    assert abs(radius(b) - 30) <= thresh
485

486
    b = s[1]
487
    assert abs(b[0] - 100) <= thresh
488
    assert abs(b[1] - 300) <= thresh
489
    assert abs(radius(b) - 40) <= thresh
490

491
    b = s[3]
492
    assert abs(b[0] - 200) <= thresh
493
    assert abs(b[1] - 350) <= thresh
494
    assert abs(radius(b) - 50) <= thresh
495

496

497
def test_blob_doh_no_peaks():
498
    # Testing no peaks
499
    img_empty = np.zeros((100, 100))
500
    assert blob_doh(img_empty).size == 0
501

502

503
def test_blob_doh_overlap():
504
    img = np.ones((256, 256), dtype=np.uint8)
505

506
    xs, ys = disk((100, 100), 20)
507
    img[xs, ys] = 255
508

509
    xs, ys = disk((120, 100), 30)
510
    img[xs, ys] = 255
511

512
    blobs = blob_doh(img, min_sigma=1, max_sigma=60, num_sigma=10, threshold=0.05)
513

514
    assert len(blobs) == 1
515

516

517
def test_blob_log_overlap_3d():
518
    r1, r2 = 7, 6
519
    pad1, pad2 = 11, 12
520
    blob1 = ellipsoid(r1, r1, r1)
521
    blob1 = np.pad(blob1, pad1, mode='constant')
522
    blob2 = ellipsoid(r2, r2, r2)
523
    blob2 = np.pad(
524
        blob2, [(pad2, pad2), (pad2 - 9, pad2 + 9), (pad2, pad2)], mode='constant'
525
    )
526
    im3 = np.logical_or(blob1, blob2)
527

528
    blobs = blob_log(im3, min_sigma=2, max_sigma=10, overlap=0.1)
529
    assert len(blobs) == 1
530

531

532
def test_blob_overlap_3d_anisotropic():
533
    # Two spheres with distance between centers equal to radius
534
    # One sphere is much smaller than the other so about half of it is within
535
    # the bigger sphere.
536
    s3 = math.sqrt(3)
537
    overlap = _blob_overlap(
538
        np.array([0, 0, 0, 2 / s3, 10 / s3, 10 / s3]),
539
        np.array([0, 0, 10, 0.2 / s3, 1 / s3, 1 / s3]),
540
        sigma_dim=3,
541
    )
542
    assert_almost_equal(overlap, 0.48125)
543
    overlap = _blob_overlap(
544
        np.array([0, 0, 0, 2 / s3, 10 / s3, 10 / s3]),
545
        np.array([2, 0, 0, 0.2 / s3, 1 / s3, 1 / s3]),
546
        sigma_dim=3,
547
    )
548
    assert_almost_equal(overlap, 0.48125)
549

550

551
def test_blob_log_anisotropic():
552
    image = np.zeros((50, 50))
553
    image[20, 10:20] = 1
554
    isotropic_blobs = blob_log(image, min_sigma=0.5, max_sigma=2, num_sigma=3)
555
    assert len(isotropic_blobs) > 1  # many small blobs found in line
556
    ani_blobs = blob_log(
557
        image, min_sigma=[0.5, 5], max_sigma=[2, 20], num_sigma=3
558
    )  # 10x anisotropy, line is 1x10
559
    assert len(ani_blobs) == 1  # single anisotropic blob found
560

561

562
def test_blob_log_overlap_3d_anisotropic():
563
    r1, r2 = 7, 6
564
    pad1, pad2 = 11, 12
565
    blob1 = ellipsoid(r1, r1, r1)
566
    blob1 = np.pad(blob1, pad1, mode='constant')
567
    blob2 = ellipsoid(r2, r2, r2)
568
    blob2 = np.pad(
569
        blob2, [(pad2, pad2), (pad2 - 9, pad2 + 9), (pad2, pad2)], mode='constant'
570
    )
571
    im3 = np.logical_or(blob1, blob2)
572

573
    blobs = blob_log(im3, min_sigma=[2, 2.01, 2.005], max_sigma=10, overlap=0.1)
574
    assert len(blobs) == 1
575

576
    # Two circles with distance between centers equal to radius
577
    overlap = _blob_overlap(
578
        np.array([0, 0, 10 / math.sqrt(2)]), np.array([0, 10, 10 / math.sqrt(2)])
579
    )
580
    assert_almost_equal(
581
        overlap, 1.0 / math.pi * (2 * math.acos(1.0 / 2) - math.sqrt(3) / 2.0)
582
    )
583

584

585
def test_no_blob():
586
    im = np.zeros((10, 10))
587
    blobs = blob_log(im, min_sigma=2, max_sigma=5, num_sigma=4)
588
    assert len(blobs) == 0
589