Scale Space

A scale space is a family of images derived from one image by smoothing it with Gaussians of every width. The width becomes a continuous third coordinate: structures much smaller than the current width have been smoothed away, while larger ones survive. Scale-space theory explains why the smoothing kernel should be a Gaussian, how derivatives are measured at a given scale, and how a characteristic scale can be selected for each image structure, without knowing in advance how large the objects in the image are. The last idea, automatic scale selection, is what makes detectors such as SIFT return a scale for every keypoint along with its position.

Definition

The scale-space representation of an image $f$ is the one-parameter family of images $L(\cdot\,; t)$ obtained by convolving $f$ with Gaussian kernels of variance $t$, for all $t \ge 0$. Derivatives taken from this family are called scale-space derivatives or Gaussian derivatives. Scale selection picks, at each point, the values of $t$ at which a suitably normalized combination of such derivatives has a local maximum over scale; the chosen scale is then treated as a characteristic size of the structure there.

Intuition

Two facts motivate the construction (Lindeberg, 2008). Objects are meaningful only over a range of scales: a branch makes sense from centimetres to metres, not at the scale of molecules or of forests. And a camera adds a second scale problem through perspective, since the same object looks larger when it is nearer. An operator applied at a fixed size sees only what matches that size, and in an unknown scene there is no way to know the right size beforehand. Lindeberg’s conclusion is that an uncommitted system should represent the image at all scales at once.

Simply blurring by different amounts is not enough. Coarse levels should be simplifications of fine levels, not new structures invented by the smoothing; Lindeberg’s reviews present this requirement, in several formal versions, as the reason the Gaussian is chosen. Once the family exists, a second question follows: given all scales, which ones matter at a particular point? A blob of radius 10 and a blob of radius 40 should each be described at their own scale, and if the image is enlarged twofold, the chosen scales should double. Automatic scale selection answers this with derivatives normalized so that their responses can be compared across scales.

A Gaussian image pyramid samples the same family at octave steps and subsamples each level. Lindeberg (2008) describes pyramids as a predecessor of scale space that is computationally efficient but makes it harder to relate structures across scales, because the levels are discrete; pyramids are still used for fast approximations of scale-space algorithms.

Formal Definition

Gaussian scale space

For an image $f:\mathbb{R}^2 \to \mathbb{R}$, the scale-space representation $L:\mathbb{R}^2 \times \mathbb{R}_+ \to \mathbb{R}$ is

$$L(x, y; t) = \iint f(x - \xi, y - \eta)\, g(\xi, \eta; t)\, d\xi\, d\eta, \qquad g(x, y; t) = \frac{1}{2\pi t}\, e^{-(x^2 + y^2)/2t},$$

with $L(\cdot\,; 0) = f$. The scale parameter $t = \sigma^2$ is the variance of the Gaussian filter; $\sigma = \sqrt{t}$ has the dimension of length and is the natural unit for comparing scales with sizes. Equivalently, $L$ solves the diffusion (heat) equation

$$\partial_t L = \tfrac{1}{2} \nabla^2 L = \tfrac{1}{2}(L_{xx} + L_{yy})$$

with initial condition $f$ (Lindeberg, 2008, eqs. 1–3). This article uses $t$ as the variance and the factor $\tfrac{1}{2}$ throughout, as Lindeberg does; some texts write $\partial_t L = \nabla^2 L$ with $t = \sigma^2/2$.

Scale-space derivatives

Derivatives at scale $t$ are obtained either by differentiating $L$ or by convolving $f$ with derivatives of the Gaussian:

$$L_{x^\alpha y^\beta}(\cdot\,; t) = \partial_{x^\alpha y^\beta} L(\cdot\,; t) = \big(\partial_{x^\alpha y^\beta}\, g(\cdot\,; t)\big) * f .$$

Defining derivatives this way makes the ill-posed problem of differentiating sampled, noisy data well-posed (Lindeberg, 2008). The set of derivatives up to order $N$ at a point and scale is the $N$-jet, a truncated Taylor expansion of the smoothed image. From the 2-jet $(L_x, L_y, L_{xx}, L_{xy}, L_{yy})$ one can form rotationally invariant quantities: the gradient magnitude $|\nabla L|$, the Laplacian $\nabla^2 L = L_{xx} + L_{yy}$, the determinant of the Hessian $\det \mathcal{H}L = L_{xx}L_{yy} - L_{xy}^2$, and the rescaled level-curve curvature $\tilde\kappa(L) = L_x^2 L_{yy} + L_y^2 L_{xx} - 2 L_x L_y L_{xy}$. Lindeberg (2008) builds the classical detectors from them: edges where the second directional derivative along the gradient vanishes and the third is negative, blobs at extrema of the Laplacian or of the determinant of the Hessian, corners at extrema of $\tilde\kappa$, and ridges where the first derivative along the principal curvature direction vanishes.

γ-normalized derivatives

Under smoothing, derivative amplitudes fall with scale: for a sine wave $\sin(\omega_0 x)$, $L(x; t) = e^{-\omega_0^2 t/2} \sin(\omega_0 x)$, and every derivative decays exponentially in $t$. Raw responses at different scales therefore cannot be compared. Lindeberg treated scale selection in two papers of 1998, “Feature detection with automatic scale selection” (1998a here) and “Edge detection and ridge detection with automatic scale selection” (1998b); section and equation numbers below are those of the authors’ technical-report versions. Lindeberg (1998a, eq. 9) introduced the γ-normalized derivative

$$\partial_\xi = t^{\gamma/2}\, \partial_x ,$$

so an $m$-th order derivative is multiplied by $t^{m\gamma/2}$. For the sine wave, the normalized $m$-th derivative $t^{m\gamma/2}\omega_0^m e^{-\omega_0^2 t/2}$ first rises and then falls, with its maximum at

$$\sigma_{\max} = \frac{\sqrt{\gamma m}}{2\pi}\, \lambda_0 ,$$

proportional to the wavelength $\lambda_0 = 2\pi/\omega_0$ (Lindeberg, 1998a, eq. 12). A maximum over scale of a normalized derivative thus measures a length at a single point, without a window or a ruler.

Automatic scale selection

Lindeberg (1998a, §4) stated a principle for scale selection, whose usefulness he says must be verified empirically: in the absence of other evidence, a scale at which some combination of normalized derivatives has a local maximum over scales can be treated as reflecting a characteristic length of a structure in the data. A scale-space extremum is a point $(x_0; t_0)$ that is a local extremum of a normalized differential expression $\mathcal{D}_{\text{norm}}L$ with respect to both position and scale:

$$\nabla\big(\mathcal{D}_{\text{norm}}L\big)(x_0; t_0) = 0, \qquad \partial_t\big(\mathcal{D}_{\text{norm}}L\big)(x_0; t_0) = 0 .$$

Two expressions with $\gamma = 1$ serve as blob measures (Lindeberg, 2014): the scale-normalized Laplacian $\nabla^2_{\text{norm}}L = t(L_{xx} + L_{yy})$ and the scale-normalized determinant of the Hessian $\det\mathcal{H}_{\text{norm}}L = t^2(L_{xx}L_{yy} - L_{xy}^2)$.

Properties

The scale-space axioms

The Gaussian is singled out by several different sets of requirements; Lindeberg (1996) reviews them and notes that similar conclusions follow from different starting points. All combine linearity and shift invariance with some formalization of “no new structure at coarser scales”:

• Non-creation of local extrema in one dimension. Witkin (1983), who introduced the term scale-space, observed for one-dimensional signals that smoothing creates no new local extrema, as Lindeberg (1996) describes; because differentiation commutes with convolution, the same holds for zero-crossings of every derivative, so zero-crossing paths in the $(x, t)$ plane are never closed from below. As Lindeberg (1996) points out, this property does not extend to two or more dimensions.
• Causality. Koenderink (1984), as Lindeberg (1996) describes it, required that new level surfaces must not be created in the scale-space volume as $t$ increases, and with isotropy and homogeneity (all positions and scales treated alike) showed that the family must satisfy the diffusion equation. The Gaussian filter article gives this result in brief.
• Semigroup and non-enhancement. Lindeberg combined the semigroup property $g(\cdot\,; t_1) * g(\cdot\,; t_2) = g(\cdot\,; t_1 + t_2)$ with non-creation of local extrema in one dimension, and in higher dimensions with non-enhancement of local extrema: the value at a spatial local maximum must not increase, and at a local minimum must not decrease, as $t$ grows (Lindeberg, 1996). The semigroup property means any coarse level can be computed from any finer one, $L(\cdot\,; t_2) = g(\cdot\,; t_2 - t_1) * L(\cdot\,; t_1)$.
• Scale invariance. Formulations by Florack and by Pauwels and colleagues, as Lindeberg (1996) summarizes them, restrict the kernel by requiring scale invariance together with linearity, shift and rotation invariance and the semigroup. This leaves a one-parameter family of kernels; adding separability (Florack, in two or more dimensions), or a local infinitesimal generator and a non-negative kernel (Pauwels and colleagues), leaves the Gaussian.

Non-enhancement can be read off the diffusion equation (our reasoning): at a spatial maximum the Hessian is negative semidefinite, so $\nabla^2 L \le 0$ and $\partial_t L \le 0$ there. It follows, as Lindeberg (1998a) notes, that derivative amplitudes generally decrease with scale, which is why normalization is needed.

Covariance under rescaling

The central property of normalized derivatives is how they behave when the image is rescaled. Let $f'(sx) = f(x)$, so $f'$ is $f$ enlarged by a factor $s$, and match points and scales by $x' = sx$, $t' = s^2 t$. Then $L(x; t) = L'(x'; t')$ and (Lindeberg, 1998a, eq. 23)

$$\partial_{\xi^m} L(x; t) = s^{m(1 - \gamma)}\, \partial_{\xi'^m} L'(x'; t') .$$

Three consequences follow (Lindeberg, 1998a, §4.1):

• With $\gamma = 1$ the normalized derivatives are equal at corresponding points: perfect scale invariance, for any expression built from the $N$-jet.
• With $\gamma \ne 1$ the values differ by a power of $s$, but for an expression that is a homogeneous polynomial of total order $M$ the factor $s^{M(1-\gamma)}$ is constant, so a maximum over scale at $(x_0; t_0)$ in $f$ is still a maximum at $(s x_0; s^2 t_0)$ in $f'$. The selected scale is covariant: it moves with the image.
• The magnitude at the maximum depends on the structure’s size unless $\gamma = 1$; multiplying by $t^{M(1-\gamma)/2}$ compensates.

Lindeberg (1998a, appendix A.1) also shows that γ-normalization is necessary: if each derivative is multiplied by a smooth function of scale, detecting maxima over scale commutes with rescaling only when that function is a power of $t$.

What γ means

The free parameter $\gamma$ fixes which scale is selected for a given kind of structure. Lindeberg (1998a, §§9.1–9.2) gives two readings. Normalizing by $t^{m\gamma/2}$ keeps the $L_p$ norm of the $m$-th order Gaussian derivative kernel constant over scale, with $p = 1/(1 + \tfrac{m}{D}(1 - \gamma))$ in $D$ dimensions, so $\gamma = 1$ is $L_1$ normalization for every order. And γ-normalized derivatives are neutral to images whose power spectrum falls off as a particular power of frequency (Lindeberg, 2014), a behaviour natural images show qualitatively.

For idealized structures the selected scale can be computed in closed form:

Structure	Measure	Selected scale	Choice making it equal the structure’s own scale
Gaussian blob, variance $t_0$	$t\,\nabla^2 L$ at the centre ($\gamma = 1$)	$t = t_0$	$\gamma = 1$
Anisotropic Gaussian blob, $t_1, t_2$	$t^2 \det\mathcal{H}L$ ($\gamma = 1$)	$t = \sqrt{t_1 t_2}$	—
Diffuse step edge, $\Phi(x; t_0)$	$t^{\gamma}(L_x^2 + L_y^2)$	$t = \frac{\gamma}{1 - \gamma}\, t_0$	$\gamma = 1/2$
Gaussian ridge, $g(x; t_0)$	$t^{\gamma} \lvert L_{pp} \rvert$	$t = \frac{2\gamma}{3 - 2\gamma}\, t_0$	$\gamma = 3/4$

The blob rows are from Lindeberg (1998a, §5.1, eqs. 37 and 39), the edge row from Lindeberg (1998b, §4.5.1, eq. 22) and the ridge row from Lindeberg (1998b, §5.6.1, eq. 62). For the edge with $\gamma = 1$, the normalized gradient grows monotonically with scale toward the height of the step, so the selected scale would be infinite: ordinary normalization responds to the infinite extent of a step rather than to its diffuseness (Lindeberg, 1998b, footnote 7). For the ridge, $\gamma = 1$ selects $t = 2t_0$, a scale too coarse by $\sqrt{2}$ in $\sigma$. For corners from the rescaled level-curve curvature, Lindeberg (2014) reports $\gamma = 7/8$ as a good choice. These values are why a single normalization does not serve every feature: $\gamma = 1$ for blobs, $\gamma = 1/2$ for edges and $\gamma = 3/4$ for ridges (Lindeberg, 2008).

Deep structure and multiple scales

A single point can have several maxima over scale. Lindeberg (2014) shows a lamp whose signature has two: a coarse one for the lamp and a finer one for the bulb inside, so a system must keep several scale hypotheses. The relations between structures at different scales, such as small blobs that merge into one larger blob as $t$ grows, are called the deep structure of scale space (Lindeberg, 1996); in a continuous scale space the implicit function theorem defines feature trajectories, and singularity theory describes the bifurcations, which is harder to do across the discrete levels of a pyramid (Lindeberg, 2008).

Affine scale space

Rotationally symmetric Gaussians make the representation closed under translations, rotations and rescalings, not under the affine distortions caused by viewing a surface obliquely. An affine scale space smooths with Gaussians of covariance matrix $\Sigma$; if $f'(Ax) = f(x)$, the two representations correspond when $\Sigma' = A\Sigma A^\top$ (Lindeberg, 2008, eq. 14–15). Estimating $\Sigma$ from the second-moment matrix by iterative shape adaptation gives affine-covariant interest points up to a scale factor and a rotation (Lindeberg, 2008); feature detection surveys such affine-covariant regions.

Examples

The example tests the scaling properties and the role of $\gamma$ on synthetic images, then the one-dimensional causality property. Three structures are drawn in $513 \times 513$ images, centred on pixel $(256, 256)$: a bright square of side $10s$, rendered with exact pixel-area coverage; a vertical diffuse step edge $\tfrac12(1 + \operatorname{erf}(x/\sqrt{2}\,3s))$, the Gaussian-blurred step of standard deviation $3s$; and a vertical Gaussian ridge of standard deviation $3s$. The factor $s \in \{1, 1.5, 2, 3, 4\}$ enlarges each structure; each image is drawn afresh from the continuous model, not resampled from a smaller image, so the rescaling is exact. Axis 0 of the arrays is $y$.

For each image the code computes the derivatives at the centre pixel over 241 scales from $\sigma = 0.5$ to 48, forms the γ-normalized measure of the table above for $\gamma = 1, 3/4, 1/2$ (for the square, $-t^\gamma \nabla^2 L$, positive at a bright blob), and refines the scale of the largest value with a parabola in $\log\sigma$. A maximum at either end of the grid counts as none. The printed table gives the selected $\sigma$ divided by $s$, which covariance predicts to be the same for every $s$; the value predicted at $s = 1$ by the continuous models (the closed forms for edge and ridge, and for the square the exact centre value $\operatorname{erf}(u)^2$, $u = A/(2\sqrt{2t})$, of the smoothed square); and the ratio of the maximum values at $s = 4$ and $s = 1$ against the law $s^{m(\gamma-1)}$, with derivative order $m = 2$ for square and ridge and $m = 1$ for the edge.

The discrete scale space is Lindeberg’s: the discrete analogue of the Gaussian, $T(n; t) = e^{-t} I_n(t)$ with $I_n$ the modified Bessel function, followed by the central differences $[-\tfrac12, 0, \tfrac12]$ and $[1, -2, 1]$. Lindeberg (2014) describes this kernel as a particularly convenient implementation, because the semigroup property holds exactly for it. The first printed lines show a further reason that matters here: at $\sigma = 0.5$ the sampled second-derivative Gaussian kernel of scipy.ndimage.gaussian_filter1d (truncated at $4\sigma$) does not sum to zero, so it responds to a constant region; in our reading, the normalized measures at fine scales would pick up that response. A check on a quadratic image confirms that the difference weights return its exact derivatives.

The second part smooths a synthetic scanline, 40 bright Gaussian spots with random positions, widths and heights, for five seeds. Because variances add under Gaussian smoothing, its scale space is exact: each spot of variance $v$ becomes one of variance $v + t$. Local extrema are counted as sign changes of the analytic $L_x$ on a grid of step 0.01; since all spots are positive, every extremum lies between the outermost centres. The count is printed at seven scales, with the number of increases over 161 scales and, as a check of the grid, the count at $t = 0$ with half the step.

import os
os.environ["OMP_NUM_THREADS"] = "1"          # one BLAS thread, set before NumPy loads
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
from scipy import ndimage, special

N, C = 513, 256                               # image size, centre pixel
coord = np.arange(N) - C
SIGMAS = np.geomspace(0.5, 48.0, 241)         # scale grid, sigma = sqrt(t)
FACTORS = [1.0, 1.5, 2.0, 3.0, 4.0]
GAMMAS = [1.0, 0.75, 0.5]

# 0. Second-derivative kernels: response to a constant (the kernel sum).
for s in (0.5, 1.0, 2.0):
    e = np.zeros(201); e[100] = 1
    k_scipy = ndimage.gaussian_filter1d(e, s, order=2)
    k_disc = np.convolve(special.ive(np.abs(np.arange(-100, 101)), s**2), [1, -2, 1], "same")
    print(f"sigma {s}: sum of sampled 2nd-derivative kernel {k_scipy.sum():+.1e}, "
          f"discrete analogue {k_disc.sum():+.1e}")


def weights(sigma, order):
    """Weights on f(C + k) of the order-th central difference of T(.; t) * f at C,
    with T(n; t) = exp(-t) I_n(t), Lindeberg's discrete analogue of the Gaussian."""
    T = special.ive(np.abs(coord), sigma**2)
    Tp, Tm = np.roll(T, -1), np.roll(T, 1)
    return [T, (Tm - Tp) / 2, Tm - 2 * T + Tp][order]


def derivatives(F, sigma):
    """Lx, Ly, Lxx, Lyy, Lxy at the centre pixel (axis 0 is y)."""
    w = [weights(sigma, k) for k in range(3)]
    return {n: w[oy] @ F @ w[ox] for n, (oy, ox) in
            {"x": (0, 1), "y": (1, 0), "xx": (0, 2), "yy": (2, 0), "xy": (1, 1)}.items()}


yy, xx = np.meshgrid(coord.astype(float), coord.astype(float), indexing="ij")
quad = xx + 2 * yy + 0.5 * xx**2 + 1.5 * yy**2 + xx * yy
for sg in (0.5, 5.0, 30.0):
    d = derivatives(quad, sg)
    print(f"sigma {sg}: derivatives of a quadratic (exact 1, 2, 1, 3, 1):",
          " ".join(f"{d[k]:.6f}" for k in ("x", "y", "xx", "yy", "xy")))


def render(kind, s):
    """The continuous structure drawn at scale factor s (re-rendered, not resampled)."""
    if kind == "square":                      # side 10 s, exact pixel-area coverage
        h = 5.0 * s
        c = np.clip(coord + 0.5, -h, h) - np.clip(coord - 0.5, -h, h)
        return np.outer(c, c)
    if kind == "edge":                        # diffuse step, std 3 s
        return 0.5 * (1 + special.erf(xx / (np.sqrt(2) * 3 * s)))
    return np.exp(-xx**2 / (2 * (3 * s) ** 2))   # ridge, std 3 s


def raw_measure(kind, F):
    """Unnormalized measure at the centre for every sigma, and its derivative order m."""
    out = []
    for sg in SIGMAS:
        d = derivatives(F, sg)
        lap = d["xx"] + d["yy"]
        if kind == "square":
            out.append(-lap)
        elif kind == "edge":
            out.append(np.hypot(d["x"], d["y"]))
        else:                                 # bright ridge: -L_pp, L_pp the smaller eigenvalue
            out.append(max(0.0, -(0.5 * lap - np.hypot(0.5 * (d["xx"] - d["yy"]), d["xy"]))))
    return np.array(out), {"square": 2, "edge": 1, "ridge": 2}[kind]


def select(raw, m, gamma):
    """Maximum over scale of t^(m gamma / 2) * raw, refined by a parabola in log sigma."""
    r = SIGMAS ** (m * gamma) * raw
    i = int(np.argmax(r))
    if i in (0, len(r) - 1):
        return None, None                     # no interior maximum on the grid
    a, b, c = np.polyfit(np.log(SIGMAS[i - 1:i + 2] / SIGMAS[i]), r[i - 1:i + 2], 2)
    return SIGMAS[i] * np.exp(-b / (2 * a)), c - b**2 / (4 * a)


def predicted(kind, gamma):
    """Selected sigma at s = 1 from the continuous models."""
    if kind == "edge":
        return None if gamma == 1 else 3 * np.sqrt(gamma / (1 - gamma))
    if kind == "ridge":
        return 3 * np.sqrt(2 * gamma / (3 - 2 * gamma))
    u = np.linspace(0.05, 5, 200001)          # square: centre value erf(u)^2,
    t = (10 / (2 * np.sqrt(2) * u)) ** 2      # u = A / (2 sqrt(2t)), A = 10
    resp = t ** (gamma - 1) * u * np.exp(-u**2) * special.erf(u)
    return np.sqrt(t[np.argmax(resp)])


print("\n                 selected sigma / s at s =            ratio of maxima, s=4 / s=1")
print("structure gamma   1.0    1.5    2.0    3.0    4.0  predicted  measured  s^(m(g-1))")
cov, err, mag, edge1 = [], [], [], []
for kind in ("square", "edge", "ridge"):
    raws = {s: raw_measure(kind, render(kind, s)) for s in FACTORS}
    for g in GAMMAS:
        sel = [select(*raws[s], g) for s in FACTORS]
        if sel[0][0] is None:
            edge1.append(all(v[0] is None for v in sel))
            print(f"{kind:9s} {g:4.2f}  none at any s (measure still rising at sigma = 48)")
            continue
        m, p = raws[1.0][1], predicted(kind, g)
        ratios = [v[0] / s for s, v in zip(FACTORS, sel)]
        cov += [abs(r / ratios[0] - 1) for r in ratios]
        err += [abs(r / p - 1) for r in ratios]
        mag += [abs(v[1] / sel[0][1] / s ** (m * (g - 1)) - 1) for s, v in zip(FACTORS, sel)]
        print(f"{kind:9s} {g:4.2f} " + "".join(f"{r:7.3f}" for r in ratios)
              + f"  {p:8.3f}  {sel[-1][1] / sel[0][1]:8.3f}  {4.0 ** (m * (g - 1)):9.3f}")
print(f"largest deviation from covariance {max(cov):.4f}, from the prediction {max(err):.4f}, "
      f"of the magnitude ratio {max(mag):.4f}; edge with gamma = 1 has no maximum: {all(edge1)}")


# 1D causality: a scanline of 40 bright Gaussian spots, smoothed exactly
# (variances add), local extrema counted as sign changes of L_x.
def count_extrema(mu, var, amp, t, dx=0.01):
    x = np.arange(mu.min(), mu.max() + dx, dx)
    v = var[:, None] + t
    lx = (amp[:, None] * -(x - mu[:, None]) / v
          * np.exp(-(x - mu[:, None]) ** 2 / (2 * v)) / np.sqrt(2 * np.pi * v)).sum(0)
    return int(np.sum(np.sign(lx[:-1]) * np.sign(lx[1:]) < 0))


ts = np.concatenate([[0.0], np.geomspace(0.01, 1600, 160)])
shown = [0, 1, 2, 4, 8, 16, 32]
print("\nlocal extrema at sigma =", shown, "| increases over 161 scales | count at t=0 with dx/2")
for seed in range(5):
    rng = np.random.default_rng(seed)
    mu, var, amp = rng.uniform(0, 400, 40), rng.uniform(1, 4, 40) ** 2, rng.uniform(0.2, 1, 40)
    counts = [count_extrema(mu, var, amp, t) for t in ts]
    at = [count_extrema(mu, var, amp, s**2) for s in shown]
    print(f"seed {seed}: {at} | {int(np.sum(np.diff(counts) > 0))} | "
          f"{count_extrema(mu, var, amp, 0.0, 0.005)}")

Output:

sigma 0.5: sum of sampled 2nd-derivative kernel -5.6e-01, discrete analogue +2.1e-16
sigma 1.0: sum of sampled 2nd-derivative kernel -7.2e-05, discrete analogue +4.5e-17
sigma 2.0: sum of sampled 2nd-derivative kernel -8.7e-05, discrete analogue +2.4e-17
sigma 0.5: derivatives of a quadratic (exact 1, 2, 1, 3, 1): 1.000000 2.000000 1.000000 3.000000 1.000000
sigma 5.0: derivatives of a quadratic (exact 1, 2, 1, 3, 1): 1.000000 2.000000 1.000000 3.000000 1.000000
sigma 30.0: derivatives of a quadratic (exact 1, 2, 1, 3, 1): 1.000000 2.000000 1.000000 3.000000 1.000000

                 selected sigma / s at s =            ratio of maxima, s=4 / s=1
structure gamma   1.0    1.5    2.0    3.0    4.0  predicted  measured  s^(m(g-1))
square    1.00   4.036  4.006  3.996  3.989  3.986     3.983     1.011      1.000
square    0.75   3.591  3.570  3.554  3.547  3.545     3.541     0.509      0.500
square    0.50   3.256  3.245  3.225  3.219  3.217     3.214     0.256      0.250
edge      1.00  none at any s (measure still rising at sigma = 48)
edge      0.75   5.257  5.223  5.211  5.203  5.200     5.196     0.708      0.707
edge      0.50   3.056  3.025  3.014  3.006  3.003     3.000     0.503      0.500
ridge     1.00   4.249  4.246  4.244  4.244  4.243     4.243     0.994      1.000
ridge     0.75   3.028  3.012  3.007  3.003  3.002     3.000     0.498      0.500
ridge     0.50   2.158  2.137  2.130  2.125  2.123     2.121     0.251      0.250
largest deviation from covariance 0.0171, from the prediction 0.0186, of the magnitude ratio 0.0238; edge with gamma = 1 has no maximum: True

local extrema at sigma = [0, 1, 2, 4, 8, 16, 32] | increases over 161 scales | count at t=0 with dx/2
seed 0: [47, 45, 39, 31, 19, 13, 5] | 0 | 47
seed 1: [43, 41, 39, 27, 13, 9, 3] | 0 | 43
seed 2: [51, 47, 45, 37, 15, 9, 7] | 0 | 51
seed 3: [49, 45, 39, 29, 17, 11, 3] | 0 | 49
seed 4: [45, 45, 39, 27, 15, 7, 5] | 0 | 45

The selected scale moves with the image for every structure and every $\gamma$ that has a maximum: $\sigma/s$ stays within 1.8% of its value at $s = 1$ and within 1.9% of the continuous prediction. The deviations are largest at $s = 1$ and shrink as the structure grows, which in our reading is the pixel grid’s effect at the smallest scales rather than a failure of covariance. Covariance does not depend on $\gamma$: with $\gamma = 1/2$ the square is still assigned a scale proportional to its size, only a smaller multiple of it.

What $\gamma$ changes is which multiple. With $\gamma = 1/2$ the edge’s selected $\sigma$ equals its blur of $3s$, and with $\gamma = 3/4$ the ridge’s equals its width $3s$, as the closed forms say; $\gamma = 1$ assigns the ridge $\sqrt{2}$ times its width (4.24 for width 3) and gives the edge no scale at all, since its normalized gradient is still rising at $\sigma = 48$. The square, which has no single “true” scale, is assigned $0.40$, $0.35$ or $0.32$ times its side depending on $\gamma$, and the measured values match the exact model of the smoothed square.

The last columns show the magnitude law. With $\gamma = 1$ the largest normalized response is the same for the small and the four-times-larger structure (ratios 1.011 and 0.994), so the strength of a feature can be compared across sizes. With $\gamma = 3/4$ the second-order measures fall by $4^{-1/2} = 0.5$ and the first-order edge measure by $4^{-1/4} = 0.707$, as predicted; a detector that ranks features by strength with $\gamma \ne 1$ has to apply the compensation factor.

In one dimension the number of local extrema falls steadily with scale for all five scanlines, from 43–51 at $\sigma = 0$ to 3–7 at $\sigma = 32$, and never increases over the 161 scales, as non-creation requires. The halved grid step finds the same initial counts, which in our reading means the grid separates even the closest pairs of extrema. The script runs in about 10 seconds with NumPy 2.3.5 and SciPy 1.16.3.

Common Misconceptions

• “Scale space is the same as an image pyramid.” A pyramid samples the family at octave steps and subsamples it; scale space is the continuous family at full resolution, in which features can be followed continuously across scales (see image pyramid).
• “Gaussian smoothing never creates new extrema.” That holds for one-dimensional signals. In two or more dimensions new local extrema can appear as scale increases; what the Gaussian guarantees there is non-enhancement, that maxima do not grow and minima do not deepen (Lindeberg, 1996, 2008).
• “Normalizing by $\sigma^2$ is always right.” Multiplying second derivatives by $t = \sigma^2$ is the $\gamma = 1$ choice that suits blobs. Edges need $\gamma = 1/2$ to select a finite scale at all, and ridges $\gamma = 3/4$ to select their width; the example shows $\gamma = 1$ failing for both.
• “Only $\gamma = 1$ gives scale-invariant features.” Any $\gamma$ gives covariant scale selection for homogeneous differential expressions; $\gamma = 1$ additionally makes the response magnitudes comparable across sizes, and other values need a compensation factor.
• “The selected scale is the object’s size.” It is proportional to a characteristic length, with a constant set by the measure, $\gamma$ and the shape: the square in the example is assigned 0.32 to 0.40 times its side. In our reading, a detector that reports a radius has to convert the scale with a constant chosen for the shapes it expects.
• “Unnormalized responses can be compared across scales.” Raw derivative responses decrease with scale, so their maximum is biased toward fine scales; the Laplacian operator article shows this on a disk.
• “A small sampled Gaussian kernel is a faithful scale space.” At $\sigma$ below about one pixel, a sampled Gaussian approximates the continuous kernel poorly, as the Gaussian filter article notes (the example’s sampled second-derivative kernel responds to a constant at $\sigma = 0.5$). The discrete analogue $e^{-t}I_n(t)$ keeps the semigroup property exactly (Lindeberg, 2014).

Implementation

Lindeberg (2014) describes the practical steps. Scale-space extrema of a 2D image are found by comparing each sample with its 26 neighbours in a $3 \times 3 \times 3$ neighbourhood over space and scale, and a parabola fitted around an extremum refines position and scale. Discrete normalized derivatives can multiply difference approximations by $t^{m\gamma/2}$ or, alternatively, scale them so that their discrete $l_p$ norms equal the $L_p$ norms of the continuous normalized kernels; experiments that Lindeberg (2014) reports show that the $l_p$ normalization gives more accurate scale estimates where discretization matters. The discrete analogue of the Gaussian allows a set of incremental convolutions with small kernels instead of one large convolution per level, and scale selection can be transferred to a pyramid for real-time use.

The Gaussian filter article covers sampling, truncation and cascades of the kernel itself. Scale-space signatures are cheap at a single point, as in the example; full detection needs every derivative at every scale, which in our reading is why practical detectors subsample the coarser octaves.

Where It Is Used

• Interest points with scale. Scale-space extrema of $\nabla^2_{\text{norm}}L$ give dark and bright blobs, and those of $\det\mathcal{H}_{\text{norm}}L$ give blobs and saddles; Lindeberg (2014) notes that the determinant of the Hessian often has better repeatability under affine deformations. SIFT’s difference of Gaussians approximates the scale-normalized Laplacian (the Laplacian operator article gives the approximation), and SURF’s detector, which approximates the determinant of the Hessian with Haar wavelets, can be seen as an approximation of the scale-normalized determinant of the Hessian (Lindeberg, 2014). The Harris–Laplace detector, described in feature detection, localizes points with the Harris corner measure and selects their scale with the normalized Laplacian.
• Scale-invariant descriptors. Computing a descriptor in a window proportional to the selected $\sigma$, or resampling the patch by it, makes the descriptor scale invariant; Lindeberg (2014) attributes the scale invariance of SIFT and SURF descriptors to this step.
• Edges and ridges with automatic scale. Lindeberg’s scale-space edges select a scale along each edge, fine for sharp edges and coarse for diffuse ones (edge detection), which Lindeberg (2014) relates to focus blur and shadow edges; ridge detection with $\gamma = 3/4$ tunes the scale to the width of elongated structures (Lindeberg, 1998b). Multi-scale vesselness filters used in ophthalmic image analysis take a maximum over scales too, but that article notes its responses are not scale-normalized.
• Blob detection in microscopy and medical imaging. Scale-normalized Laplacian maxima generate nucleus and nodule candidates in cell detection and lesion detection.
• Feature tracking and frequency estimation. Lindeberg (2014) describes tracking with features whose scales adapt over time, robust enough to estimate time-to-collision from scale changes; Lindeberg (1998a, §8) estimates local frequency densely from maxima over scale.
• A model of early vision. Lindeberg (1998a) cites neurophysiological measurements by Young of receptive fields in the mammalian retina and visual cortex that are well modelled by Gaussian derivatives up to order four, and presents the scale-space operators as an idealized visual front-end. The scale of a learned network’s features is the subject of receptive fields, and multi-scale feature maps in detectors are the subject of feature pyramid networks.
• Other scale selection methods. Lindeberg (2014) lists alternatives, including peaks over scale of entropy measures, minimum reliable scales from a noise model, and stopping times for non-linear diffusion. He contrasts them with maxima of normalized derivatives: a noise-driven method chooses a scale set by the noise level, and on noise-free data no smoothing at all, while maxima over scale always reflect a characteristic length of the image structure.