scikit-image

1
import math
2

3
import numpy as np
4
import scipy.ndimage as ndi
5
from scipy import spatial
6

7
from .._shared.filters import gaussian
8
from .._shared.utils import _supported_float_type, check_nD
9
from ..transform import integral_image
10
from ..util import img_as_float
11
from ._hessian_det_appx import _hessian_matrix_det
12
from .peak import peak_local_max
13

14
# This basic blob detection algorithm is based on:
15
# http://www.cs.utah.edu/~jfishbau/advimproc/project1/ (04.04.2013)
16
# Theory behind: https://en.wikipedia.org/wiki/Blob_detection (04.04.2013)
17

18

19
def _compute_disk_overlap(d, r1, r2):
20
    """
21
    Compute fraction of surface overlap between two disks of radii
22
    ``r1`` and ``r2``, with centers separated by a distance ``d``.
23

24
    Parameters
25
    ----------
26
    d : float
27
        Distance between centers.
28
    r1 : float
29
        Radius of the first disk.
30
    r2 : float
31
        Radius of the second disk.
32

33
    Returns
34
    -------
35
    fraction: float
36
        Fraction of area of the overlap between the two disks.
37
    """
38

39
    ratio1 = (d**2 + r1**2 - r2**2) / (2 * d * r1)
40
    ratio1 = np.clip(ratio1, -1, 1)
41
    acos1 = math.acos(ratio1)
42

43
    ratio2 = (d**2 + r2**2 - r1**2) / (2 * d * r2)
44
    ratio2 = np.clip(ratio2, -1, 1)
45
    acos2 = math.acos(ratio2)
46

47
    a = -d + r2 + r1
48
    b = d - r2 + r1
49
    c = d + r2 - r1
50
    d = d + r2 + r1
51
    area = r1**2 * acos1 + r2**2 * acos2 - 0.5 * math.sqrt(abs(a * b * c * d))
52
    return area / (math.pi * (min(r1, r2) ** 2))
53

54

55
def _compute_sphere_overlap(d, r1, r2):
56
    """
57
    Compute volume overlap fraction between two spheres of radii
58
    ``r1`` and ``r2``, with centers separated by a distance ``d``.
59

60
    Parameters
61
    ----------
62
    d : float
63
        Distance between centers.
64
    r1 : float
65
        Radius of the first sphere.
66
    r2 : float
67
        Radius of the second sphere.
68

69
    Returns
70
    -------
71
    fraction: float
72
        Fraction of volume of the overlap between the two spheres.
73

74
    Notes
75
    -----
76
    See for example http://mathworld.wolfram.com/Sphere-SphereIntersection.html
77
    for more details.
78
    """
79
    vol = (
80
        math.pi
81
        / (12 * d)
82
        * (r1 + r2 - d) ** 2
83
        * (d**2 + 2 * d * (r1 + r2) - 3 * (r1**2 + r2**2) + 6 * r1 * r2)
84
    )
85
    return vol / (4.0 / 3 * math.pi * min(r1, r2) ** 3)
86

87

88
def _blob_overlap(blob1, blob2, *, sigma_dim=1):
89
    """Finds the overlapping area fraction between two blobs.
90

91
    Returns a float representing fraction of overlapped area. Note that 0.0
92
    is *always* returned for dimension greater than 3.
93

94
    Parameters
95
    ----------
96
    blob1 : sequence of arrays
97
        A sequence of ``(row, col, sigma)`` or ``(pln, row, col, sigma)``,
98
        where ``row, col`` (or ``(pln, row, col)``) are coordinates
99
        of blob and ``sigma`` is the standard deviation of the Gaussian kernel
100
        which detected the blob.
101
    blob2 : sequence of arrays
102
        A sequence of ``(row, col, sigma)`` or ``(pln, row, col, sigma)``,
103
        where ``row, col`` (or ``(pln, row, col)``) are coordinates
104
        of blob and ``sigma`` is the standard deviation of the Gaussian kernel
105
        which detected the blob.
106
    sigma_dim : int, optional
107
        The dimensionality of the sigma value. Can be 1 or the same as the
108
        dimensionality of the blob space (2 or 3).
109

110
    Returns
111
    -------
112
    f : float
113
        Fraction of overlapped area (or volume in 3D).
114
    """
115
    ndim = len(blob1) - sigma_dim
116
    if ndim > 3:
117
        return 0.0
118
    root_ndim = math.sqrt(ndim)
119

120
    # we divide coordinates by sigma * sqrt(ndim) to rescale space to isotropy,
121
    # giving spheres of radius = 1 or < 1.
122
    if blob1[-1] == blob2[-1] == 0:
123
        return 0.0
124
    elif blob1[-1] > blob2[-1]:
125
        max_sigma = blob1[-sigma_dim:]
126
        r1 = 1
127
        r2 = blob2[-1] / blob1[-1]
128
    else:
129
        max_sigma = blob2[-sigma_dim:]
130
        r2 = 1
131
        r1 = blob1[-1] / blob2[-1]
132
    pos1 = blob1[:ndim] / (max_sigma * root_ndim)
133
    pos2 = blob2[:ndim] / (max_sigma * root_ndim)
134

135
    d = np.sqrt(np.sum((pos2 - pos1) ** 2))
136
    if d > r1 + r2:  # centers farther than sum of radii, so no overlap
137
        return 0.0
138

139
    # one blob is inside the other
140
    if d <= abs(r1 - r2):
141
        return 1.0
142

143
    if ndim == 2:
144
        return _compute_disk_overlap(d, r1, r2)
145

146
    else:  # ndim=3 http://mathworld.wolfram.com/Sphere-SphereIntersection.html
147
        return _compute_sphere_overlap(d, r1, r2)
148

149

150
def _prune_blobs(blobs_array, overlap, *, sigma_dim=1):
151
    """Eliminated blobs with area overlap.
152

153
    Parameters
154
    ----------
155
    blobs_array : ndarray
156
        A 2d array with each row representing 3 (or 4) values,
157
        ``(row, col, sigma)`` or ``(pln, row, col, sigma)`` in 3D,
158
        where ``(row, col)`` (``(pln, row, col)``) are coordinates of the blob
159
        and ``sigma`` is the standard deviation of the Gaussian kernel which
160
        detected the blob.
161
        This array must not have a dimension of size 0.
162
    overlap : float
163
        A value between 0 and 1. If the fraction of area overlapping for 2
164
        blobs is greater than `overlap` the smaller blob is eliminated.
165
    sigma_dim : int, optional
166
        The number of columns in ``blobs_array`` corresponding to sigmas rather
167
        than positions.
168

169
    Returns
170
    -------
171
    A : ndarray
172
        `array` with overlapping blobs removed.
173
    """
174
    sigma = blobs_array[:, -sigma_dim:].max()
175
    distance = 2 * sigma * math.sqrt(blobs_array.shape[1] - sigma_dim)
176
    tree = spatial.cKDTree(blobs_array[:, :-sigma_dim])
177
    pairs = np.array(list(tree.query_pairs(distance)))
178
    if len(pairs) == 0:
179
        return blobs_array
180
    else:
181
        for i, j in pairs:
182
            blob1, blob2 = blobs_array[i], blobs_array[j]
183
            if _blob_overlap(blob1, blob2, sigma_dim=sigma_dim) > overlap:
184
                # note: this test works even in the anisotropic case because
185
                # all sigmas increase together.
186
                if blob1[-1] > blob2[-1]:
187
                    blob2[-1] = 0
188
                else:
189
                    blob1[-1] = 0
190

191
    return np.stack([b for b in blobs_array if b[-1] > 0])
192

193

194
def _format_exclude_border(img_ndim, exclude_border):
195
    """Format an ``exclude_border`` argument as a tuple of ints for calling
196
    ``peak_local_max``.
197
    """
198
    if isinstance(exclude_border, tuple):
199
        if len(exclude_border) != img_ndim:
200
            raise ValueError(
201
                "`exclude_border` should have the same length as the "
202
                "dimensionality of the image."
203
            )
204
        for exclude in exclude_border:
205
            if not isinstance(exclude, int):
206
                raise ValueError(
207
                    "exclude border, when expressed as a tuple, must only "
208
                    "contain ints."
209
                )
210
        return exclude_border + (0,)
211
    elif isinstance(exclude_border, int):
212
        return (exclude_border,) * img_ndim + (0,)
213
    elif exclude_border is True:
214
        raise ValueError("exclude_border cannot be True")
215
    elif exclude_border is False:
216
        return (0,) * (img_ndim + 1)
217
    else:
218
        raise ValueError(f'Unsupported value ({exclude_border}) for exclude_border')
219

220

221
def blob_dog(
222
    image,
223
    min_sigma=1,
224
    max_sigma=50,
225
    sigma_ratio=1.6,
226
    threshold=0.5,
227
    overlap=0.5,
228
    *,
229
    threshold_rel=None,
230
    exclude_border=False,
231
):
232
    r"""Finds blobs in the given grayscale image.
233

234
    Blobs are found using the Difference of Gaussian (DoG) method [1]_, [2]_.
235
    For each blob found, the method returns its coordinates and the standard
236
    deviation of the Gaussian kernel that detected the blob.
237

238
    Parameters
239
    ----------
240
    image : ndarray
241
        Input grayscale image, blobs are assumed to be light on dark
242
        background (white on black).
243
    min_sigma : scalar or sequence of scalars, optional
244
        The minimum standard deviation for Gaussian kernel. Keep this low to
245
        detect smaller blobs. The standard deviations of the Gaussian filter
246
        are given for each axis as a sequence, or as a single number, in
247
        which case it is equal for all axes.
248
    max_sigma : scalar or sequence of scalars, optional
249
        The maximum standard deviation for Gaussian kernel. Keep this high to
250
        detect larger blobs. The standard deviations of the Gaussian filter
251
        are given for each axis as a sequence, or as a single number, in
252
        which case it is equal for all axes.
253
    sigma_ratio : float, optional
254
        The ratio between the standard deviation of Gaussian Kernels used for
255
        computing the Difference of Gaussians
256
    threshold : float or None, optional
257
        The absolute lower bound for scale space maxima. Local maxima smaller
258
        than `threshold` are ignored. Reduce this to detect blobs with lower
259
        intensities. If `threshold_rel` is also specified, whichever threshold
260
        is larger will be used. If None, `threshold_rel` is used instead.
261
    overlap : float, optional
262
        A value between 0 and 1. If the area of two blobs overlaps by a
263
        fraction greater than `threshold`, the smaller blob is eliminated.
264
    threshold_rel : float or None, optional
265
        Minimum intensity of peaks, calculated as
266
        ``max(dog_space) * threshold_rel``, where ``dog_space`` refers to the
267
        stack of Difference-of-Gaussian (DoG) images computed internally. This
268
        should have a value between 0 and 1. If None, `threshold` is used
269
        instead.
270
    exclude_border : tuple of ints, int, or False, optional
271
        If tuple of ints, the length of the tuple must match the input array's
272
        dimensionality.  Each element of the tuple will exclude peaks from
273
        within `exclude_border`-pixels of the border of the image along that
274
        dimension.
275
        If nonzero int, `exclude_border` excludes peaks from within
276
        `exclude_border`-pixels of the border of the image.
277
        If zero or False, peaks are identified regardless of their
278
        distance from the border.
279

280
    Returns
281
    -------
282
    A : (n, image.ndim + sigma) ndarray
283
        A 2d array with each row representing 2 coordinate values for a 2D
284
        image, or 3 coordinate values for a 3D image, plus the sigma(s) used.
285
        When a single sigma is passed, outputs are:
286
        ``(r, c, sigma)`` or ``(p, r, c, sigma)`` where ``(r, c)`` or
287
        ``(p, r, c)`` are coordinates of the blob and ``sigma`` is the standard
288
        deviation of the Gaussian kernel which detected the blob. When an
289
        anisotropic gaussian is used (sigmas per dimension), the detected sigma
290
        is returned for each dimension.
291

292
    See also
293
    --------
294
    skimage.filters.difference_of_gaussians
295

296
    References
297
    ----------
298
    .. [1] https://en.wikipedia.org/wiki/Blob_detection#The_difference_of_Gaussians_approach
299
    .. [2] Lowe, D. G. "Distinctive Image Features from Scale-Invariant
300
        Keypoints." International Journal of Computer Vision 60, 91–110 (2004).
301
        https://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf
302
        :DOI:`10.1023/B:VISI.0000029664.99615.94`
303

304
    Examples
305
    --------
306
    >>> from skimage import data, feature
307
    >>> coins = data.coins()
308
    >>> feature.blob_dog(coins, threshold=.05, min_sigma=10, max_sigma=40)
309
    array([[128., 155.,  10.],
310
           [198., 155.,  10.],
311
           [124., 338.,  10.],
312
           [127., 102.,  10.],
313
           [193., 281.,  10.],
314
           [126., 208.,  10.],
315
           [267., 115.,  10.],
316
           [197., 102.,  10.],
317
           [198., 215.,  10.],
318
           [123., 279.,  10.],
319
           [126.,  46.,  10.],
320
           [259., 247.,  10.],
321
           [196.,  43.,  10.],
322
           [ 54., 276.,  10.],
323
           [267., 358.,  10.],
324
           [ 58., 100.,  10.],
325
           [259., 305.,  10.],
326
           [185., 347.,  16.],
327
           [261., 174.,  16.],
328
           [ 46., 336.,  16.],
329
           [ 54., 217.,  10.],
330
           [ 55., 157.,  10.],
331
           [ 57.,  41.,  10.],
332
           [260.,  47.,  16.]])
333

334
    Notes
335
    -----
336
    The radius of each blob is approximately :math:`\sqrt{2}\sigma` for
337
    a 2-D image and :math:`\sqrt{3}\sigma` for a 3-D image.
338
    """  # noqa: E501
339
    image = img_as_float(image)
340
    float_dtype = _supported_float_type(image.dtype)
341
    image = image.astype(float_dtype, copy=False)
342

343
    # if both min and max sigma are scalar, function returns only one sigma
344
    scalar_sigma = np.isscalar(max_sigma) and np.isscalar(min_sigma)
345

346
    # Gaussian filter requires that sequence-type sigmas have same
347
    # dimensionality as image. This broadcasts scalar kernels
348
    if np.isscalar(max_sigma):
349
        max_sigma = np.full(image.ndim, max_sigma, dtype=float_dtype)
350
    if np.isscalar(min_sigma):
351
        min_sigma = np.full(image.ndim, min_sigma, dtype=float_dtype)
352

353
    # Convert sequence types to array
354
    min_sigma = np.asarray(min_sigma, dtype=float_dtype)
355
    max_sigma = np.asarray(max_sigma, dtype=float_dtype)
356

357
    if sigma_ratio <= 1.0:
358
        raise ValueError('sigma_ratio must be > 1.0')
359

360
    # k such that min_sigma*(sigma_ratio**k) > max_sigma
361
    k = int(np.mean(np.log(max_sigma / min_sigma) / np.log(sigma_ratio) + 1))
362

363
    # a geometric progression of standard deviations for gaussian kernels
364
    sigma_list = np.array([min_sigma * (sigma_ratio**i) for i in range(k + 1)])
365

366
    # computing difference between two successive Gaussian blurred images
367
    # to obtain an approximation of the scale invariant Laplacian of the
368
    # Gaussian operator
369
    dog_image_cube = np.empty(image.shape + (k,), dtype=float_dtype)
370
    gaussian_previous = gaussian(image, sigma=sigma_list[0], mode='reflect')
371
    for i, s in enumerate(sigma_list[1:]):
372
        gaussian_current = gaussian(image, sigma=s, mode='reflect')
373
        dog_image_cube[..., i] = gaussian_previous - gaussian_current
374
        gaussian_previous = gaussian_current
375

376
    # normalization factor for consistency in DoG magnitude
377
    sf = 1 / (sigma_ratio - 1)
378
    dog_image_cube *= sf
379

380
    exclude_border = _format_exclude_border(image.ndim, exclude_border)
381
    local_maxima = peak_local_max(
382
        dog_image_cube,
383
        threshold_abs=threshold,
384
        threshold_rel=threshold_rel,
385
        exclude_border=exclude_border,
386
        footprint=np.ones((3,) * (image.ndim + 1)),
387
    )
388

389
    # Catch no peaks
390
    if local_maxima.size == 0:
391
        return np.empty((0, image.ndim + (1 if scalar_sigma else image.ndim)))
392

393
    # Convert local_maxima to float64
394
    lm = local_maxima.astype(float_dtype)
395

396
    # translate final column of lm, which contains the index of the
397
    # sigma that produced the maximum intensity value, into the sigma
398
    sigmas_of_peaks = sigma_list[local_maxima[:, -1]]
399

400
    if scalar_sigma:
401
        # select one sigma column, keeping dimension
402
        sigmas_of_peaks = sigmas_of_peaks[:, 0:1]
403

404
    # Remove sigma index and replace with sigmas
405
    lm = np.hstack([lm[:, :-1], sigmas_of_peaks])
406

407
    sigma_dim = sigmas_of_peaks.shape[1]
408

409
    return _prune_blobs(lm, overlap, sigma_dim=sigma_dim)
410

411

412
def blob_log(
413
    image,
414
    min_sigma=1,
415
    max_sigma=50,
416
    num_sigma=10,
417
    threshold=0.2,
418
    overlap=0.5,
419
    log_scale=False,
420
    *,
421
    threshold_rel=None,
422
    exclude_border=False,
423
):
424
    r"""Finds blobs in the given grayscale image.
425

426
    Blobs are found using the Laplacian of Gaussian (LoG) method [1]_.
427
    For each blob found, the method returns its coordinates and the standard
428
    deviation of the Gaussian kernel that detected the blob.
429

430
    Parameters
431
    ----------
432
    image : ndarray
433
        Input grayscale image, blobs are assumed to be light on dark
434
        background (white on black).
435
    min_sigma : scalar or sequence of scalars, optional
436
        the minimum standard deviation for Gaussian kernel. Keep this low to
437
        detect smaller blobs. The standard deviations of the Gaussian filter
438
        are given for each axis as a sequence, or as a single number, in
439
        which case it is equal for all axes.
440
    max_sigma : scalar or sequence of scalars, optional
441
        The maximum standard deviation for Gaussian kernel. Keep this high to
442
        detect larger blobs. The standard deviations of the Gaussian filter
443
        are given for each axis as a sequence, or as a single number, in
444
        which case it is equal for all axes.
445
    num_sigma : int, optional
446
        The number of intermediate values of standard deviations to consider
447
        between `min_sigma` and `max_sigma`.
448
    threshold : float or None, optional
449
        The absolute lower bound for scale space maxima. Local maxima smaller
450
        than `threshold` are ignored. Reduce this to detect blobs with lower
451
        intensities. If `threshold_rel` is also specified, whichever threshold
452
        is larger will be used. If None, `threshold_rel` is used instead.
453
    overlap : float, optional
454
        A value between 0 and 1. If the area of two blobs overlaps by a
455
        fraction greater than `threshold`, the smaller blob is eliminated.
456
    log_scale : bool, optional
457
        If set intermediate values of standard deviations are interpolated
458
        using a logarithmic scale to the base `10`. If not, linear
459
        interpolation is used.
460
    threshold_rel : float or None, optional
461
        Minimum intensity of peaks, calculated as
462
        ``max(log_space) * threshold_rel``, where ``log_space`` refers to the
463
        stack of Laplacian-of-Gaussian (LoG) images computed internally. This
464
        should have a value between 0 and 1. If None, `threshold` is used
465
        instead.
466
    exclude_border : tuple of ints, int, or False, optional
467
        If tuple of ints, the length of the tuple must match the input array's
468
        dimensionality.  Each element of the tuple will exclude peaks from
469
        within `exclude_border`-pixels of the border of the image along that
470
        dimension.
471
        If nonzero int, `exclude_border` excludes peaks from within
472
        `exclude_border`-pixels of the border of the image.
473
        If zero or False, peaks are identified regardless of their
474
        distance from the border.
475

476
    Returns
477
    -------
478
    A : (n, image.ndim + sigma) ndarray
479
        A 2d array with each row representing 2 coordinate values for a 2D
480
        image, or 3 coordinate values for a 3D image, plus the sigma(s) used.
481
        When a single sigma is passed, outputs are:
482
        ``(r, c, sigma)`` or ``(p, r, c, sigma)`` where ``(r, c)`` or
483
        ``(p, r, c)`` are coordinates of the blob and ``sigma`` is the standard
484
        deviation of the Gaussian kernel which detected the blob. When an
485
        anisotropic gaussian is used (sigmas per dimension), the detected sigma
486
        is returned for each dimension.
487

488
    References
489
    ----------
490
    .. [1] https://en.wikipedia.org/wiki/Blob_detection#The_Laplacian_of_Gaussian
491

492
    Examples
493
    --------
494
    >>> from skimage import data, feature, exposure
495
    >>> img = data.coins()
496
    >>> img = exposure.equalize_hist(img)  # improves detection
497
    >>> feature.blob_log(img, threshold = .3)
498
    array([[124.        , 336.        ,  11.88888889],
499
           [198.        , 155.        ,  11.88888889],
500
           [194.        , 213.        ,  17.33333333],
501
           [121.        , 272.        ,  17.33333333],
502
           [263.        , 244.        ,  17.33333333],
503
           [194.        , 276.        ,  17.33333333],
504
           [266.        , 115.        ,  11.88888889],
505
           [128.        , 154.        ,  11.88888889],
506
           [260.        , 174.        ,  17.33333333],
507
           [198.        , 103.        ,  11.88888889],
508
           [126.        , 208.        ,  11.88888889],
509
           [127.        , 102.        ,  11.88888889],
510
           [263.        , 302.        ,  17.33333333],
511
           [197.        ,  44.        ,  11.88888889],
512
           [185.        , 344.        ,  17.33333333],
513
           [126.        ,  46.        ,  11.88888889],
514
           [113.        , 323.        ,   1.        ]])
515

516
    Notes
517
    -----
518
    The radius of each blob is approximately :math:`\sqrt{2}\sigma` for
519
    a 2-D image and :math:`\sqrt{3}\sigma` for a 3-D image.
520
    """  # noqa: E501
521
    image = img_as_float(image)
522
    float_dtype = _supported_float_type(image.dtype)
523
    image = image.astype(float_dtype, copy=False)
524

525
    # if both min and max sigma are scalar, function returns only one sigma
526
    scalar_sigma = True if np.isscalar(max_sigma) and np.isscalar(min_sigma) else False
527

528
    # Gaussian filter requires that sequence-type sigmas have same
529
    # dimensionality as image. This broadcasts scalar kernels
530
    if np.isscalar(max_sigma):
531
        max_sigma = np.full(image.ndim, max_sigma, dtype=float_dtype)
532
    if np.isscalar(min_sigma):
533
        min_sigma = np.full(image.ndim, min_sigma, dtype=float_dtype)
534

535
    # Convert sequence types to array
536
    min_sigma = np.asarray(min_sigma, dtype=float_dtype)
537
    max_sigma = np.asarray(max_sigma, dtype=float_dtype)
538

539
    if log_scale:
540
        start = np.log10(min_sigma)
541
        stop = np.log10(max_sigma)
542
        sigma_list = np.logspace(start, stop, num_sigma)
543
    else:
544
        sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)
545

546
    # computing gaussian laplace
547
    image_cube = np.empty(image.shape + (len(sigma_list),), dtype=float_dtype)
548
    for i, s in enumerate(sigma_list):
549
        # average s**2 provides scale invariance
550
        image_cube[..., i] = -ndi.gaussian_laplace(image, s) * np.mean(s) ** 2
551

552
    exclude_border = _format_exclude_border(image.ndim, exclude_border)
553
    local_maxima = peak_local_max(
554
        image_cube,
555
        threshold_abs=threshold,
556
        threshold_rel=threshold_rel,
557
        exclude_border=exclude_border,
558
        footprint=np.ones((3,) * (image.ndim + 1)),
559
    )
560

561
    # Catch no peaks
562
    if local_maxima.size == 0:
563
        return np.empty((0, image.ndim + (1 if scalar_sigma else image.ndim)))
564

565
    # Convert local_maxima to float64
566
    lm = local_maxima.astype(float_dtype)
567

568
    # translate final column of lm, which contains the index of the
569
    # sigma that produced the maximum intensity value, into the sigma
570
    sigmas_of_peaks = sigma_list[local_maxima[:, -1]]
571

572
    if scalar_sigma:
573
        # select one sigma column, keeping dimension
574
        sigmas_of_peaks = sigmas_of_peaks[:, 0:1]
575

576
    # Remove sigma index and replace with sigmas
577
    lm = np.hstack([lm[:, :-1], sigmas_of_peaks])
578

579
    sigma_dim = sigmas_of_peaks.shape[1]
580

581
    return _prune_blobs(lm, overlap, sigma_dim=sigma_dim)
582

583

584
def blob_doh(
585
    image,
586
    min_sigma=1,
587
    max_sigma=30,
588
    num_sigma=10,
589
    threshold=0.01,
590
    overlap=0.5,
591
    log_scale=False,
592
    *,
593
    threshold_rel=None,
594
):
595
    """Finds blobs in the given grayscale image.
596

597
    Blobs are found using the Determinant of Hessian method [1]_. For each blob
598
    found, the method returns its coordinates and the standard deviation
599
    of the Gaussian Kernel used for the Hessian matrix whose determinant
600
    detected the blob. Determinant of Hessians is approximated using [2]_.
601

602
    Parameters
603
    ----------
604
    image : 2D ndarray
605
        Input grayscale image.Blobs can either be light on dark or vice versa.
606
    min_sigma : float, optional
607
        The minimum standard deviation for Gaussian Kernel used to compute
608
        Hessian matrix. Keep this low to detect smaller blobs.
609
    max_sigma : float, optional
610
        The maximum standard deviation for Gaussian Kernel used to compute
611
        Hessian matrix. Keep this high to detect larger blobs.
612
    num_sigma : int, optional
613
        The number of intermediate values of standard deviations to consider
614
        between `min_sigma` and `max_sigma`.
615
    threshold : float or None, optional
616
        The absolute lower bound for scale space maxima. Local maxima smaller
617
        than `threshold` are ignored. Reduce this to detect blobs with lower
618
        intensities. If `threshold_rel` is also specified, whichever threshold
619
        is larger will be used. If None, `threshold_rel` is used instead.
620
    overlap : float, optional
621
        A value between 0 and 1. If the area of two blobs overlaps by a
622
        fraction greater than `threshold`, the smaller blob is eliminated.
623
    log_scale : bool, optional
624
        If set intermediate values of standard deviations are interpolated
625
        using a logarithmic scale to the base `10`. If not, linear
626
        interpolation is used.
627
    threshold_rel : float or None, optional
628
        Minimum intensity of peaks, calculated as
629
        ``max(doh_space) * threshold_rel``, where ``doh_space`` refers to the
630
        stack of Determinant-of-Hessian (DoH) images computed internally. This
631
        should have a value between 0 and 1. If None, `threshold` is used
632
        instead.
633

634
    Returns
635
    -------
636
    A : (n, 3) ndarray
637
        A 2d array with each row representing 3 values, ``(y,x,sigma)``
638
        where ``(y,x)`` are coordinates of the blob and ``sigma`` is the
639
        standard deviation of the Gaussian kernel of the Hessian Matrix whose
640
        determinant detected the blob.
641

642
    References
643
    ----------
644
    .. [1] https://en.wikipedia.org/wiki/Blob_detection#The_determinant_of_the_Hessian
645
    .. [2] Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool,
646
           "SURF: Speeded Up Robust Features"
647
           ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf
648

649
    Examples
650
    --------
651
    >>> from skimage import data, feature
652
    >>> img = data.coins()
653
    >>> feature.blob_doh(img)
654
    array([[197.        , 153.        ,  20.33333333],
655
           [124.        , 336.        ,  20.33333333],
656
           [126.        , 153.        ,  20.33333333],
657
           [195.        , 100.        ,  23.55555556],
658
           [192.        , 212.        ,  23.55555556],
659
           [121.        , 271.        ,  30.        ],
660
           [126.        , 101.        ,  20.33333333],
661
           [193.        , 275.        ,  23.55555556],
662
           [123.        , 205.        ,  20.33333333],
663
           [270.        , 363.        ,  30.        ],
664
           [265.        , 113.        ,  23.55555556],
665
           [262.        , 243.        ,  23.55555556],
666
           [185.        , 348.        ,  30.        ],
667
           [156.        , 302.        ,  30.        ],
668
           [123.        ,  44.        ,  23.55555556],
669
           [260.        , 173.        ,  30.        ],
670
           [197.        ,  44.        ,  20.33333333]])
671

672
    Notes
673
    -----
674
    The radius of each blob is approximately `sigma`.
675
    Computation of Determinant of Hessians is independent of the standard
676
    deviation. Therefore detecting larger blobs won't take more time. In
677
    methods line :py:meth:`blob_dog` and :py:meth:`blob_log` the computation
678
    of Gaussians for larger `sigma` takes more time. The downside is that
679
    this method can't be used for detecting blobs of radius less than `3px`
680
    due to the box filters used in the approximation of Hessian Determinant.
681
    """  # noqa: E501
682
    check_nD(image, 2)
683

684
    image = img_as_float(image)
685
    float_dtype = _supported_float_type(image.dtype)
686
    image = image.astype(float_dtype, copy=False)
687

688
    image = integral_image(image)
689

690
    if log_scale:
691
        start, stop = math.log(min_sigma, 10), math.log(max_sigma, 10)
692
        sigma_list = np.logspace(start, stop, num_sigma)
693
    else:
694
        sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)
695

696
    image_cube = np.empty(shape=image.shape + (len(sigma_list),), dtype=float_dtype)
697
    for j, s in enumerate(sigma_list):
698
        image_cube[..., j] = _hessian_matrix_det(image, s)
699

700
    local_maxima = peak_local_max(
701
        image_cube,
702
        threshold_abs=threshold,
703
        threshold_rel=threshold_rel,
704
        exclude_border=False,
705
        footprint=np.ones((3,) * image_cube.ndim),
706
    )
707

708
    # Catch no peaks
709
    if local_maxima.size == 0:
710
        return np.empty((0, 3))
711
    # Convert local_maxima to float64
712
    lm = local_maxima.astype(np.float64)
713
    # Convert the last index to its corresponding scale value
714
    lm[:, -1] = sigma_list[local_maxima[:, -1]]
715
    return _prune_blobs(lm, overlap)
716