Blob Detection

Blob detection finds compact image regions that are brighter or darker than their surroundings: spots, dots, nuclei, windows on a facade, the heads of flowers in a field. A blob detector returns each region’s centre and, in its scale-invariant form, its size. The standard detectors smooth the image with Gaussians of many widths, form a second-derivative measure at every scale (the scale-normalized Laplacian of Gaussian, its difference-of-Gaussians approximation, or the determinant of the Hessian), and keep the points where that measure is extremal in both position and scale. Blob detection is one family of feature detection, alongside corner and region detectors; its keypoints feed the descriptors of SIFT and SURF, and in microscopy and medical imaging the blobs are themselves the objects of interest.

The word has a second, older sense in industrial vision, sometimes called blob analysis: thresholding an image and treating each connected region as a blob. That sense is described under Variants below; most of this article is about the differential detectors.

Problem Definition

A blob has no single formal definition, and the differential detectors make the idea operational through second derivatives: at the centre of a bright blob the intensity surface curves downward in every direction, so the Hessian matrix of second derivatives,

\mathbf{H}L = \begin{bmatrix} L_{xx} & L_{xy} \\ L_{xy} & L_{yy} \end{bmatrix},

is negative definite, and at the centre of a dark blob it is positive definite. Here $L(x, y; t) = G_t * I$ is the image smoothed with a Gaussian of variance $t = \sigma^2$ , and subscripts are partial derivatives of $L$ , computed by convolving with Gaussian derivative kernels (see the Gaussian filter).

The size of a blob is defined by the scale at which a detector responds most strongly. A detector should therefore satisfy two requirements: the response should peak at the blob’s centre, and, when the image is rescaled by a factor $s$ , the selected scale should rescale by $s$ too. The second requirement is what scale space theory delivers with scale-normalized derivatives; this article takes that theory as given and concentrates on the detectors built on it.

Input and Output

The input is a grayscale image, or one channel of a colour or multichannel image, sometimes a 3D volume. The output is a list of detections, each with:

Position $(x, y)$ , at pixel precision or refined to subpixel precision.
Scale $\sigma$ , from which a radius follows. scikit-image’s documentation gives the radius of a blob_log or blob_dog detection as about $\sqrt{2}\,\sigma$ in 2D and of a blob_doh detection as about $\sigma$ (release 0.25.2).
Polarity: bright, dark or, for the determinant of the Hessian, saddle.
Strength, the value of the normalized measure at the extremum, used to threshold and rank detections.

Affine-adapted variants replace the circle with an ellipse.

Variants

Bright, dark and saddle points. With the scale-normalized Laplacian $\nabla^2_{\text{norm}} L = t\,(L_{xx} + L_{yy})$ , bright blobs are negative scale-space minima and dark blobs positive maxima. With the scale-normalized determinant of the Hessian $\det \mathbf{H}_{\text{norm}} L = t^2 (L_{xx}L_{yy} - L_{xy}^2)$ , both bright and dark blobs are positive maxima, told apart by the sign of the Laplacian, and negative minima are saddle-like points (Lindeberg, 2014; Lindeberg, 2015, §4.2). scikit-image’s blob_log and blob_dog detect one polarity only, bright blobs on a dark background, while blob_doh keeps positive maxima of the determinant and so finds either polarity.
Fixed-scale against scale-selecting detectors. Applying the Laplacian at one predetermined scale favours structures of matching size. Lindeberg (1998, §5.2, technical-report version) shows this bias with the 50 strongest responses at three scales and contrasts it with detection of scale-space maxima, which handles structures of all sizes between the inner and outer scale in the same way.
Affine-adapted blobs. Under a change of viewpoint a round blob becomes an ellipse. Affine shape adaptation iterates on the second-moment matrix until the normalized patch is isotropic, as the Harris-Affine and Hessian-Affine detectors described in feature detection do.
Threshold-based blob analysis. OpenCV’s SimpleBlobDetector (source, release 4.14.0) binarizes an 8-bit image at a series of thresholds (by default from 50 in steps of 10 while below 220, that is 50 to 210), extracts the contours of the connected regions at each, filters them by area, inertia ratio (a moment-based elongation measure), convexity and colour (filters on by default; a circularity filter is off), groups centres that recur at nearby positions across thresholds, and keeps groups seen at least twice (minRepeatability = 2). By default it looks for dark blobs (blobColor = 0). This approach suits high-contrast, well-separated objects; it has no scale selection, and touching objects merge, a failure that cell detection measures for Otsu thresholding with connected components.
Extremal regions. Maximally stable extremal regions, described in feature detection, take the components of thresholded images that change least over a range of thresholds, a region-based alternative to differential detectors.
3D blobs. The same detectors extend to volumes: blob_log and blob_dog accept 3D arrays, and the example in lesion detection uses scale-normalized 3D Laplacian maxima as nodule candidates.

Challenges

Edges and ridges. The Laplacian responds strongly beside straight edges and along thin lines, where the position along the structure is undefined. Lowe (2004, §4.1) notes that the difference of Gaussians has a strong response along edges, even where the location along the edge is poorly determined and unstable to small amounts of noise, and removes such points with a test on the ratio of principal curvatures (see SIFT).
Elongated blobs. A blob stretched by foreshortening has two widths, and different detectors assign it different sizes, as the example below measures.
Touching and nested blobs. Two touching blobs form a larger blob at coarse scales, and a structure can have several characteristic scales: Lindeberg (2014) shows a lamp whose scale-space signature has one maximum for the lamp and a finer one for the bulb. A detector must either keep several hypotheses or impose a prior on size, as the scale-range experiment in cell detection does.
Noise and discretization. Second derivatives amplify noise, and small scales suffer from the sampling of the Gaussian kernels; the smallest usable scale is limited by both.
Thresholds. Whether a weak response is a blob depends on contrast, and the determinant of the Hessian scales with the square of contrast while the Laplacian scales linearly, so thresholds are not interchangeable between them.

Approaches

Scale-normalized Laplacian of Gaussian

The Laplacian operator article derives the Laplacian of Gaussian, shows why its response must be multiplied by $\sigma^2$ before scales are compared, and works out the best scale for a disk. As a detector, the scale-normalized Laplacian is evaluated on a list of scales, the responses are stacked into a volume over $(x, y, \sigma)$ , and points that are extrema among their 26 neighbours in a $3 \times 3 \times 3$ neighbourhood and exceed a threshold are kept (Lindeberg, 2014). For a rotationally symmetric Gaussian blob of variance $t_0$ , the extremum lies at $t = t_0$ , so the selected scale equals the blob’s size (Lindeberg, 2014).

Lindeberg (2015, §2) traces Laplacian blob features as primitives for texture analysis to work by Voorhees and Poggio and by Blostein and Ahuja, and describes how his own earlier scale-space primal sketch (Lindeberg, 1993) linked grey-level blobs across scales and measured their significance by their life length and the volume they occupy in scale space. The 1998 paper generalized scale selection to extrema over scales of normalized derivatives and notes that in practice it suffices to detect maxima over space and scale, without tracking them along their paths (Lindeberg, 1998, §1 and §4.4, footnote 3, technical-report version).

Difference of Gaussians

The Laplacian operator article derives the approximation $G(x, y, k\sigma) - G(x, y, \sigma) \approx (k - 1)\sigma^2 \nabla^2 G$ (Lowe, 2004, §3): a difference of Gaussians (DoG) with a constant scale ratio $k$ already carries the $\sigma^2$ normalization, and the constant $k - 1$ does not move the extrema. Lowe reports that the approximation has almost no effect on the stability of extremum detection even for $k = \sqrt{2}$ , and Lindeberg (2015, Appendix A) writes it as $\mathrm{DoG} \approx \tfrac{k^2 - 1}{2} \nabla^2_{\text{norm}} L$ , whose factor $\tfrac{k^2-1}{2} = (k - 1)\tfrac{k + 1}{2}$ tends to Lowe’s $k - 1$ as $k \to 1$ . Which scale to report is a convention: scikit-image’s blob_dog reports the smaller of the two scales, while SIFT’s analysis of a Gaussian blob shows that the geometric mean $\sigma\sqrt{k}$ equals the blob’s standard deviation.

Determinant of the Hessian

The determinant of the Hessian is the product of its two eigenvalues, $\det \mathbf{H}L = \lambda_1 \lambda_2$ , while the Laplacian is their sum, $\lambda_1 + \lambda_2$ . Lindeberg (2015, §4.2) presents the two as the natural pair of rotation-invariant second-order measures. A large determinant requires strong curvature in two orthogonal directions, so it is small on structures that vary in one direction only. Lindeberg (2015, §4.1–4.2) reports that on a real image the Laplacian gives many responses along an oblique elongated ridge and outside the edges of objects, while the determinant of the Hessian gives none on such one-dimensional structures and localizes corners better. Lindeberg (2015, §2) attributes early uses of the trace and determinant of the Hessian as image operators to Beaudet (1978).

The determinant is normalized with $t^2$ (that is, $\sigma^4$ ), since it is a product of two second derivatives each normalized by $t$ . Lindeberg (1998, §5.1, Example 1, technical-report version) analyses the anisotropic Gaussian blob $f(x, y) = g(x; t_1)\, g(y; t_2)$ , with widths $\sqrt{t_1}$ and $\sqrt{t_2}$ along the axes. The normalized Laplacian at its centre has a unique maximum over scale, at $t = t_0$ when $t_1 = t_2 = t_0$ , and the normalized determinant has its maximum at

t_{\det \mathbf{H}} = \sqrt{t_1 t_2}.

Lindeberg (2013, §3.3) uses this model to study foreshortening. Writing $t_1 = s\,t_0$ and $t_2 = t_0 / s$ , a non-uniform scaling with determinant one, the determinant of the Hessian selects $t = t_0$ for every $s$ : its scale estimate is affine covariant for this model. For the Laplacian the selected scale at the centre is the root of a cubic equation (his eq. 85, here with $\gamma = 1$ ),

4t^3 + 2t^2 (t_1 + t_2) + t\,(t_1^2 - 6 t_1 t_2 + t_2^2) - 2 t_1 t_2 (t_1 + t_2) = 0,

whose Taylor expansion around $s = 1$ begins $t = \bigl(1 - \tfrac{1}{4}(s - 1)^2 + \dots\bigr)\, t_0$ (his eq. 92): the Laplacian’s estimate shrinks as the blob is stretched. Lindeberg (2014) summarizes that determinant-of-Hessian features often have better repeatability under affine deformations than Laplacian features.

SURF detects scale-space extrema of an approximation of the determinant of the Hessian, built from box filters on an integral image as feature detection describes, and Lindeberg (2014) regards its detection step as an approximation of the scale-normalized determinant; scikit-image’s blob_doh computes a box-filter approximation on an integral image, citing the SURF paper, and according to its documentation cannot detect blobs of radius below about 3 pixels for that reason. Haar-like features explains the integral image.

Complementary tests

A detector built on one measure can be filtered with another. Lowe’s edge test keeps a DoG extremum only if the ratio of the Hessian’s eigenvalues is below $r = 10$ , checked through $\operatorname{Tr}(\mathbf{H})^2 / \operatorname{Det}(\mathbf{H}) < (r + 1)^2 / r$ (Lowe, 2004, §4.1). Lindeberg (2015, §6.2), inspired by it, thresholds Laplacian features on a Hessian-based feature strength $L_{xx}L_{yy} - L_{xy}^2 - k(L_{xx} + L_{yy})^2 \ge 0$ , and reports that this removes Laplacian responses along elongated structures and many outside the edges of objects.

Which measure works best

The published comparisons disagree, and they measured different things. Lowe (2004, §3) cites a detailed comparison by Mikolajczyk (2002) in which extrema of $\sigma^2\nabla^2 G$ gave the most stable image features of the functions tested, including the gradient, the Hessian and the Harris function. Lindeberg (2015) matched features between photographs of 12 posters taken at relative scale factors from 1.25 to 6 and at slant angles of about 22.5°, 30° and 45°, with SIFT-like descriptors computed from Gaussian derivatives. Among detectors based on scale-space extrema, both with the complementary test above, the determinant of the Hessian had a higher efficiency (the fraction of interest points that led to accepted matches) than the Laplacian, averaging 0.7678 against 0.7498 over the scaling and foreshortening pairs, and a lower 1 − precision, 0.0339 against 0.0359 (Lindeberg, 2015, Table 2). In that evaluation the best results came from detectors that link features across scales rather than taking scale-space extrema, and Lindeberg argues that replacing SIFT’s DoG extrema by determinant-of-Hessian extrema could be expected to improve SIFT.

Implementations

scikit-image (release 0.25.2) provides the three detectors: blob_log and blob_dog, whose scale lists and normalization cell detection describes (blob_log can also space its scales logarithmically), and blob_doh, which evaluates the box-filter determinant. All three find maxima with peak_local_max over a $3 \times 3 \times 3$ neighbourhood and then remove the smaller of two blobs whose overlap exceeds the fraction overlap. A Gaussian-derivative Hessian is available from hessian_matrix(..., use_gaussian_derivatives=True). OpenCV provides SimpleBlobDetector for threshold-based blobs and the SIFT detector for DoG blobs.

Example: Laplacian, DoG and determinant of the Hessian side by side

The script compares the three detectors on synthetic images with known structure. It uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2, runs in about 10 seconds, and gives identical output on repeated runs. Images are floating-point arrays with background 0 and structures of peak value 1; rows are $y$ and columns $x$ , and all blobs and lines are axis-aligned. The measures are the negated scale-normalized Laplacian $-\sigma^2 \nabla^2 L$ (positive at bright blobs, as in blob_log), the DoG $(L(\sigma) - L(k\sigma))/(k - 1)$ with $k = 2^{1/3}$ (as in blob_dog), and $\sigma^4 \det \mathbf{H}L$ from Gaussian derivative filters (both polarities, as in blob_doh). The Laplacian and determinant use 19 logarithmically spaced scales from 2 to 16 pixels; the DoG uses the series $2 k^i$ that blob_dog builds for the same range and reports the smaller scale of each pair. Detections are scale-space maxima over a $3 \times 3 \times 3$ neighbourhood above a threshold, without any pruning of overlapping detections.

The script has four parts:

Responses. Isolated noise-free images of a Gaussian blob ( $\sigma = 6$ ), a straight step edge, a straight ridge (a line with Gaussian profile, $\sigma = 2$ ) and a corner (one quadrant bright). Each measure’s largest value over all pixels and scales is divided by its value on the blob. It also prints how much SciPy’s second-derivative filter responds to a constant image.
Detections. A $384 \times 800$ image with three blobs of equal “size” $\sigma_0 = \sqrt{\sigma_x \sigma_y} = 6$ and aspect ratios 1, 2 and 4 (B, C and D), a vertical ridge, and a bright quadrant whose corner and two edges run to the image border, with Gaussian noise of standard deviation 0.02, for noise seeds 0, 1 and 2. Thresholds are $f = 0.1$ and $0.3$ times each measure’s noise-free peak on blob B. Because the determinant is quadratic in contrast, it is also run at $f^2$ , the threshold matching $f$ in contrast. Detections are counted by region: within 6 pixels of a blob centre, within 80 pixels of the corner, within 30 pixels of an edge (and farther from the corner), within 30 pixels of the ridge, or elsewhere.
Scale under stretching. Anisotropic Gaussian blobs with $t_1 = s\,t_0$ , $t_2 = t_0/s$ , $\sigma_0 = 6$ and $s = 1, 2, 4$ , noise-free; the response at the centre is computed on 241 logarithmically spaced scales from 2.5 to 12, and the selected scale is divided by $\sigma_0$ . The predictions in brackets are the root of the cubic above, $(t_1 t_2)^{1/4}$ for the determinant, and the continuous DoG on the same grid.
Library checks. The hand-written Laplacian and DoG detections are compared with blob_log and blob_dog run with overlap=1.0 (no pruning), and the hand-written determinant with scikit-image’s Gaussian-derivative Hessian.

import os
os.environ["OMP_NUM_THREADS"] = "1"          # identical sums on every run
os.environ["OPENBLAS_NUM_THREADS"] = "1"
import numpy as np
from scipy import ndimage as ndi
from scipy.special import erf
from skimage.feature import blob_log, blob_dog, hessian_matrix

Phi = lambda z: 0.5 * (1 + erf(z / np.sqrt(2)))       # soft step, blur 1 px
SIG = np.logspace(np.log10(2), np.log10(16), 19)       # LoG and DoH scales
K = 2 ** (1 / 3)                                       # DoG scale ratio
nd = int(np.mean(np.log(16 / 2) / np.log(K) + 1))      # blob_dog's rule
SIGD = np.array([2 * K ** i for i in range(nd + 1)])

def log_stack(im):        # -sigma^2 Laplacian of Gaussian (bright blobs positive)
    return np.stack([-ndi.gaussian_laplace(im, s) * s * s for s in SIG])

def dog_stack(im):        # (L(sigma) - L(k sigma)) / (k - 1), as blob_dog
    g = [ndi.gaussian_filter(im, s) for s in SIGD]
    return np.stack([(g[i] - g[i + 1]) / (K - 1) for i in range(nd)])

def doh_stack(im, sigmas=SIG):   # sigma^4 (Lxx Lyy - Lxy^2), Gaussian derivatives
    out = []
    for s in sigmas:
        a = ndi.gaussian_filter(im, s, order=(2, 0))
        b = ndi.gaussian_filter(im, s, order=(0, 2))
        c = ndi.gaussian_filter(im, s, order=(1, 1))
        out.append(s ** 4 * (a * b - c * c))
    return np.stack(out)

DET = {"LoG": (log_stack, SIG), "DoG": (dog_stack, SIGD[:-1]), "DoH": (doh_stack, SIG)}

def maxima(st, sig, thr):  # 3x3x3 scale-space maxima above thr -> rows (r, c, sigma)
    m = (st == ndi.maximum_filter(st, size=3, mode="nearest")) & (st > thr)
    s, r, c = np.nonzero(m)
    return np.column_stack([r, c, sig[s]]).astype(float)

# ---- 1. Peak responses to single structures (noise-free 256 x 256 canvases)
n = 256
r, c = np.mgrid[:n, :n].astype(float)
canvases = {
    "blob": np.exp(-((r - 128) ** 2 + (c - 128) ** 2) / (2 * 6.0 ** 2)),
    "edge": Phi(c - 127.5),
    "ridge": np.exp(-((c - 128) ** 2) / (2 * 2.0 ** 2)),
    "corner": Phi(r - 127.5) * Phi(c - 127.5),
}
print("1. peak response / peak at the blob (sigma at the peak)")
print("       " + "".join(f"{k:>16s}" for k in canvases))
for name, (fn, sig) in DET.items():
    peaks = {}
    for k, im in canvases.items():
        st = fn(im)
        i, y, x = np.unravel_index(st.argmax(), st.shape)
        peaks[k] = (st.max(), sig[i], (y, x))
    ref = peaks["blob"][0]
    print(f"{name}   " + "".join(f"{v / ref:9.4f} ({s:5.2f})" for v, s, _ in peaks.values()),
          f" blob peak at ({int(peaks['blob'][2][0])}, {int(peaks['blob'][2][1])})")
for s in (2.83, 10.08):
    dc = ndi.gaussian_filter1d(np.ones(400), s, order=2)[200]
    print(f"sigma^2 x response of the 2nd-derivative filter to a constant, sigma={s}: {dc * s * s:.1e}")

# ---- 2. Detections on a composite image with noise (seeds 0, 1, 2)
H, W = 384, 800
r, c = np.mgrid[:H, :W].astype(float)
def gblob(y0, x0, sx, sy):
    return np.exp(-((c - x0) ** 2) / (2 * sx ** 2) - ((r - y0) ** 2) / (2 * sy ** 2))
centres = {"B": (64, 80), "C": (64, 200), "D": (64, 330)}
clean = (gblob(64, 80, 6, 6) + gblob(64, 200, 6 * 2 ** 0.5, 6 / 2 ** 0.5)
         + gblob(64, 330, 12, 3) + np.exp(-((c - 440) ** 2) / (2 * 2.0 ** 2))
         + Phi(r - 159.5) * Phi(c - 559.5))

def region(y, x):
    for k, (y0, x0) in centres.items():
        if np.hypot(y - y0, x - x0) <= 6:
            return k
    dcor = np.hypot(y - 159.5, x - 559.5)
    if dcor <= 80:
        return "corner"
    if (abs(y - 159.5) <= 30 and x >= 559.5) or (abs(x - 559.5) <= 30 and y >= 159.5):
        return "edge"
    if abs(x - 440) <= 30:
        return "ridge"
    return "other"

cols = ["B", "C", "D", "corner", "edge", "ridge", "other"]
ref = {k: fn(clean).max(axis=0)[58:71, 74:87].max() for k, (fn, _) in DET.items()}
rows = [(f"{d} f={f}", d, f * ref[d]) for f in (0.1, 0.3) for d in ("LoG", "DoG", "DoH")]
rows += [(f"DoH f^2={f * f:.2f}", "DoH", f * f * ref["DoH"]) for f in (0.1, 0.3)]
counts = {lab: dict.fromkeys(cols, 0) for lab, _, _ in rows}
hitsB = {lab: [] for lab, _, _ in rows}
stacks0 = None
for seed in (0, 1, 2):
    img = clean + 0.02 * np.random.default_rng(seed).standard_normal(clean.shape)
    st = {d: fn(img) for d, (fn, _) in DET.items()}
    if seed == 0:
        img0, stacks0 = img, st
        print(f"data check: noise std {np.std(img - clean):.4f}, ridge peak {clean[300, 440]:.3f},"
              f" quadrant {clean[300, 700]:.3f}, background {clean[300, 100]:.4f}")
    for lab, d, thr in rows:
        det = maxima(st[d], DET[d][1], thr)
        b = [p for p in det if region(p[0], p[1]) == "B"]
        hitsB[lab].append(f"{len(b)}@{max([np.hypot(p[0]-64, p[1]-80) for p in b], default=-1):.0f}px")
        for p in det:
            counts[lab][region(p[0], p[1])] += 1
print("\n2. detections summed over noise seeds 0-2 (blob B per seed: count@max distance)")
print(f"{'':16s}" + "".join(f"{k:>7s}" for k in cols) + "   blob B")
for lab, _, _ in rows:
    print(f"{lab:16s}" + "".join(f"{counts[lab][k]:7d}" for k in cols) + "   " + " ".join(hitsB[lab]))

# ---- 3. Selected scale at the centre of anisotropic Gaussian blobs, sigma0 = 6
sf = np.logspace(np.log10(2.5), np.log10(12), 241)
r, c = np.mgrid[:201, :201].astype(float)
print("\n3. selected sigma / sigma0 at the centre (prediction in brackets)")
print("  s     LoG              DoG sigma   DoG sigma*sqrt(k)        DoH")
for s in (1, 2, 4):
    t0 = 36.0; t1, t2 = s * t0, t0 / s
    im = np.exp(-((c - 100) ** 2) / (2 * t1) - ((r - 100) ** 2) / (2 * t2))
    mx, my = np.sqrt((im * (c - 100) ** 2).sum() / im.sum()), np.sqrt((im * (r - 100) ** 2).sum() / im.sum())
    lap = [-ndi.gaussian_laplace(im, q)[100, 100] * q * q for q in sf]
    dog = [(ndi.gaussian_filter(im, q)[100, 100] - ndi.gaussian_filter(im, K * q)[100, 100]) / (K - 1)
           for q in sf]
    doh = [doh_stack(im, [q])[0, 100, 100] for q in sf]
    # exact model: Lindeberg (2013) eq. (85) for the Laplacian, eq. (96) for det H,
    # and the continuous DoG on the same grid
    roots = np.roots([4, 2 * (t1 + t2), t1 ** 2 - 6 * t1 * t2 + t2 ** 2, -2 * t1 * t2 * (t1 + t2)])
    p_lap = np.sqrt(roots[(abs(roots.imag) < 1e-9) & (roots.real > 0)].real[0]) / 6
    L0 = lambda t: np.sqrt(t1 * t2 / ((t1 + t) * (t2 + t)))
    p_dog = sf[np.argmax(L0(sf ** 2) - L0((K * sf) ** 2))] * np.sqrt(K) / 6
    p_doh = (t1 * t2) ** 0.25 / 6
    g = lambda v: sf[int(np.argmax(v))] / 6
    print(f"  {s}   {g(lap):.4f} ({p_lap:.4f})   {g(dog):.4f}      {g(dog) * np.sqrt(K):.4f} ({p_dog:.4f})"
          f"   {g(doh):.4f} ({p_doh:.4f})   data: sx {mx:.2f} sy {my:.2f}")

# ---- 4. Hand implementations against scikit-image (noise seed 0, f = 0.1)
print("\n4. checks against scikit-image 0.25.2")
tl, td = 0.1 * ref["LoG"], 0.1 * ref["DoG"]
lib = blob_log(img0, 2, 16, num_sigma=19, threshold=tl, overlap=1.0, log_scale=True)
mine = maxima(stacks0["LoG"], SIG, tl)
srt = lambda a: a[np.lexsort(a.T[::-1])]
same = lambda a, b: a.shape == b.shape and np.array_equal(srt(a), srt(b))
print(f"blob_log : {len(lib)} blobs, identical to hand LoG: {same(lib, mine)}")
lib = blob_dog(img0, 2, 16, sigma_ratio=K, threshold=td, overlap=1.0)
mine = maxima(stacks0["DoG"], SIGD[:-1], td)
print(f"blob_dog : {len(lib)} blobs, identical to hand DoG: {same(lib, mine)}")
worst, mine = 0.0, doh_stack(clean)      # noise-free composite
for i, s in enumerate(SIG):
    hrr, hrc, hcc = hessian_matrix(clean, s, mode="reflect", use_gaussian_derivatives=True)
    worst = max(worst, np.abs(s ** 4 * (hrr * hcc - hrc ** 2) - mine[i]).max())
print(f"hessian_matrix det (noise-free image): max |difference| / blob DoH peak = {worst / ref['DoH']:.1e}")

Output:

1. peak response / peak at the blob (sigma at the peak)
                   blob            edge           ridge          corner
LoG      1.0000 ( 6.35)   0.4857 (16.00)   0.7731 ( 2.83)   0.8297 (16.00)  blob peak at (128, 128)
DoG      1.0000 ( 5.04)   0.4831 (12.70)   0.7736 ( 2.52)   0.8259 (12.70)  blob peak at (128, 128)
DoH      1.0000 ( 6.35)   0.0033 (10.08)   0.0028 ( 2.83)   0.6603 (16.00)  blob peak at (128, 128)
sigma^2 x response of the 2nd-derivative filter to a constant, sigma=2.83: -7.8e-04
sigma^2 x response of the 2nd-derivative filter to a constant, sigma=10.08: -1.0e-03
data check: noise std 0.0200, ridge peak 1.000, quadrant 1.000, background 0.0000

2. detections summed over noise seeds 0-2 (blob B per seed: count@max distance)
                      B      C      D corner   edge  ridge  other   blob B
LoG f=0.1             3      3      3     27     88    111      0   1@0px 1@0px 1@0px
DoG f=0.1             3      3      3     23     76    114      0   1@0px 1@0px 1@0px
DoH f=0.1             3      3      3     12      0      0      0   1@0px 1@0px 1@0px
LoG f=0.3             3      3      3     27     88    111      0   1@0px 1@0px 1@0px
DoG f=0.3             3      3      3     23     76    114      0   1@0px 1@0px 1@0px
DoH f=0.3             3      3      3     12      0      0      0   1@0px 1@0px 1@0px
DoH f^2=0.01          3      3      3     70    133    246      0   1@0px 1@0px 1@0px
DoH f^2=0.09          3      3      3     12      0      0      0   1@0px 1@0px 1@0px

3. selected sigma / sigma0 at the centre (prediction in brackets)
  s     LoG              DoG sigma   DoG sigma*sqrt(k)        DoH
  1   1.0003 (1.0000)   0.8951      1.0048 (0.9982)   1.0003 (1.0000)   data: sx 6.00 sy 6.00
  2   0.9371 (0.9362)   0.8440      0.9474 (0.9351)   1.0003 (1.0000)   data: sx 8.49 sy 4.24
  4   0.7262 (0.7336)   0.6456      0.7247 (0.7342)   1.0003 (1.0000)   data: sx 12.00 sy 3.00

4. checks against scikit-image 0.25.2
blob_log : 72 blobs, identical to hand LoG: True
blob_dog : 68 blobs, identical to hand DoG: True
hessian_matrix det (noise-free image): max |difference| / blob DoH peak = 3.3e-03

Responses. On the isolated structures, the Laplacian and the DoG respond to the ridge at about 0.77 of their blob response and to the edge at about 0.48, close to the continuous values for these profiles (our calculation: $2/3^{3/2} \div \tfrac{1}{2} = 0.770$ for a Gaussian ridge at $\sigma = \sqrt{2}\,\sigma_r$ , and about 0.48 for the edge at $\sigma = 16$ ). Their edge peaks sit at the largest scale of each series (12.70 is the largest scale the DoG reports); in the continuous model (our calculation) the normalized response to this slightly blurred edge grows slowly with scale. The determinant of the Hessian responds to the edge and the ridge at about 0.3% of its blob response. In the continuous theory the value is exactly zero, because these images vary along one axis only. In our reading, the residual comes from the discrete filters: the second-derivative kernel along the line should give zero on a constant, but the printed line shows that SciPy’s kernel, truncated at four standard deviations, gives about $-10^{-3}$ after normalization, and this small error multiplies the strong curvature across the line. All three measures respond to the corner, the Laplacian and DoG more strongly than the determinant, relative to their blob responses.

Detections. Every detector finds blob B exactly once in every seed, on its centre, and makes three detections near each of blobs C and D over the three seeds. The Laplacian and the DoG also fire along the edges (88 and 76 detections over the three seeds) and the ridge (111 and 114), at both thresholds, because the edge and ridge responses are 0.48 to 0.77 of the blob response, above either threshold. The determinant of the Hessian fires only at the corner (12 detections) and the blobs at $f = 0.1$ , $f = 0.3$ and $f^2 = 0.09$ . At the lowest threshold, $f^2 = 0.01$ , the planned hypothesis that the determinant makes no detections on the edges and the ridge failed: it made 133 on the edges, 246 on the ridge and 70 at the corner. In our reading, noise gives the second derivative along the line a small random value, which the large curvature across the line multiplies into a determinant above that threshold. Nothing was detected outside the labelled regions at any threshold. The many detections per edge come from counting raw maxima: on a noise-free straight edge the response is constant along the edge, and in our reading noise breaks that plateau into many small maxima.

Scale under stretching. For the round blob ( $s = 1$ ) all three measures recover $\sigma_0$ to within half a percent. As the blob is stretched, the determinant of the Hessian keeps selecting $\sigma_0$ (1.0003 is the grid point nearest 1), as Lindeberg’s analysis predicts, while the Laplacian selects 0.94 $\sigma_0$ at $s = 2$ and 0.73 $\sigma_0$ at $s = 4$ , within 0.01 of the predicted root. The DoG follows the Laplacian when its scale is read as $\sigma\sqrt{k}$ , within 0.02 of the continuous DoG prediction; read as the smaller scale, as blob_dog reports it, every value is a further factor $\sqrt{k} \approx 1.12$ smaller. The scale the Laplacian assigns to a foreshortened blob therefore depends on the foreshortening, which matters when scales are used to normalize patches for description.

Library checks. The hand-written Laplacian and DoG detectors give exactly the same detections as blob_log and blob_dog, and the hand-written determinant agrees with scikit-image’s Gaussian-derivative Hessian to within 0.33% of the blob response, about the size of the edge and ridge residuals above (scikit-image applies two first-derivative filters of width $\sigma/\sqrt{2}$ with a wider truncation instead of one second-derivative filter). The library check of the determinant, as first planned, failed: it compared the two on the noisy image, while its tolerance had been computed on the noise-free image, and it was repaired by editor decision before the first run, and piloted once more, to use the noise-free image, as shown.

Datasets and Benchmarks

Blob detectors used as keypoint detectors are evaluated on the same image pairs as other detectors, with a known homography between images: the Oxford affine sequences and HPatches, described in feature detection, and purpose-built sets such as Lindeberg’s (2015) poster photographs, taken at known scale factors and slant angles. Where blobs are the objects of interest, benchmarks come from the application: cell and nucleus detection with point annotations, described in cell detection, and lung nodule detection in CT, described in lesion detection. Synthetic images with blobs of known size and shape, as in the example above and in the analyses of Lindeberg (1998, 2013), test scale estimates directly, since real images rarely come with ground-truth scales.

Evaluation Metrics

As keypoint detectors, blob detectors are scored by repeatability, the fraction of detections found again after a known transformation, with a tolerance on position and on the overlap of the detected regions (feature detection gives the protocol), and downstream by matching performance: Lindeberg (2015) reports efficiency, the fraction of interest points that lead to accepted matches, and 1 − precision, the fraction of rejected matches. As object detectors, they are scored by matching detected and annotated centres within a radius and computing precision, recall and F1, as in cell detection. Scale accuracy is measured against known sizes on synthetic data.

Applications

Keypoints for matching. SIFT’s DoG extrema and SURF’s box-filter Hessian extrema are blob detectors; their keypoints are matched in feature matching for image stitching, structure from motion and image registration, and the detected scale sets the size of the patch the descriptor summarizes.
Microscopy. Nuclei, spots and particles in fluorescence microscopy images are bright blobs on a dark background; cell detection compares LoG and DoG detectors with thresholding on such images, one of the tasks that bioimage analysis surveys.
Medical imaging. In the example of lesion detection, scale-normalized Laplacian maxima generate nodule candidates in synthetic CT-like volumes, and a classifier that uses a ratio of Hessian eigenvalues rejects vessel-like tubes.
Industrial inspection and counting. Threshold-based blob analysis measures well-separated, high-contrast parts and defects, with area, circularity and moment filters as in OpenCV’s SimpleBlobDetector.
Texture and scale cues. Lindeberg (1998, §5) shows blobs detected in an image of a sunflower field, whose detected sizes reflect the size variations in the image, and Lindeberg (2015, §2) lists uses of scale-invariant blobs for feature tracking and gesture recognition.

Open Problems

Choosing the measure. The Laplacian, its DoG approximation and the determinant of the Hessian behave differently on edges, ridges and foreshortened blobs, and the published comparisons point in different directions depending on the data and protocol.
Elongated and touching structures. Real blobs are rarely Gaussian; stretched, irregular and touching objects shift the selected scale or merge detections, and no single scale range suits populations of very different sizes.
Discretization at small scales. At scales of a few pixels, sampled Gaussian kernels and box-filter approximations depart from the continuous theory, which limits how small a blob can be detected and sized reliably.
Learned detectors. Networks trained for keypoint or object detection replace hand-designed measures; in our reading, their scale behaviour has no closed-form analysis of the kind that makes the differential detectors predictable.