SURF

SURF (Speeded-Up Robust Features) is a keypoint detector and descriptor that follows the plan of SIFT, detection of blobs over position and scale, an orientation for each keypoint and a descriptor pooled over a grid of subregions, but builds every filter from rectangles. The detector looks for scale-space maxima of an approximation of the determinant of the Hessian in which the Gaussian second derivatives are replaced by box filters; the orientation and the descriptor use Haar wavelet responses, differences of two box sums. Every box sum is read from an integral image in four look-ups, so a filter costs the same at every size and the scale space is built by enlarging the filters rather than by smoothing and subsampling the image. SURF was presented by Herbert Bay, Tinne Tuytelaars and Luc Van Gool at the European Conference on Computer Vision in Graz in 2006, and the authors’ project page asks readers to cite the journal version by Bay, Andreas Ess, Tuytelaars and Van Gool in Computer Vision and Image Understanding (2008). Neither paper could be read for this article, so the algorithm is described here from the implementations in OpenCV and scikit-image, from Lindeberg’s analyses of scale selection, and from the published neighbours; nothing in it is attributed to the papers beyond their identity.

Problem

Matching images taken from different positions, at different zoom or with the camera rotated needs keypoints that are found again on the same scene points, and descriptors that stay similar for the same point and differ for different points; feature detection and feature descriptors set out the requirements. SIFT meets them with a pyramid of Gaussian-smoothed images and gradient histograms. SURF addresses the cost of that pipeline: every quantity it computes is a weighted sum of rectangles, so the cost per filter does not grow with scale. scikit-image’s documentation of its box-filter determinant makes the same point: the running time depends only on the size of the image, not on the scale.

Inputs and Outputs

The input is a grayscale image (OpenCV’s implementation takes 8-bit images and converts colour images to gray). The output is a set of keypoints, each with a position, a size (the side of the box filter that detected it), an orientation, the detector response and the sign of the Hessian’s trace, and for each keypoint a descriptor: 64 floating-point numbers of unit length, or 128 in the extended form. Descriptors are compared by Euclidean distance, as SIFT’s are, and matched as described in feature matching.

Intuition

The determinant of the Hessian, the product of the two principal curvatures of the intensity surface, is large at the centre of a bright or dark blob and small on edges and ridges, where only one curvature is large. Blob detection develops it as a blob detector with Gaussian derivatives. The Gaussian second derivative across a blob has a negative centre lobe and two positive side lobes; a box filter with three flat lobes weighted 1, −2 and 1 has the same shape in caricature, and the mixed derivative $L_{xy}$ , with its four signed quadrants, becomes four signed squares. Box filters are poor approximations in detail, but they are cheap: once the integral image is built, which takes one pass over the image (Haar-like features derives it), a box of any size costs four reads.

The same trick serves the descriptor. A Haar wavelet, the difference between the sums of two adjacent rectangles, measures the intensity change in $x$ or $y$ at the scale of the rectangles. SURF summarizes a keypoint’s neighbourhood by sums of such responses and of their absolute values, which record how much the intensity changes in each direction and with which sign, in place of SIFT’s gradient-orientation histograms.

Algorithm

The steps below follow OpenCV 4.14.0’s SURF (in opencv_contrib), whose source comments describe several changes, made by Ian Mahon to an earlier contribution by Liu Liu, as bringing the code in line with the description in the SURF paper. Scales are expressed through $s = 1.2\,L/9$ , the factor the code uses to relate a filter of side $L$ to the sizes of its orientation and descriptor windows; a $9 \times 9$ filter corresponds to $s = 1.2$ .

Integral image. Compute the integral image of the input once.
Determinant of the Hessian over scales. For filter sizes $L = (9 + 6j) \cdot 2^o$ for layer $j$ and octave $o$ , evaluate the box filters $D_{xx}$ , $D_{yy}$ and $D_{xy}$ at every sample position and form $D_{xx}D_{yy} - 0.81\,D_{xy}^2$ . The sample step is 1 pixel in the first octave and doubles in each further octave; the image is never subsampled. The sign of the trace $D_{xx} + D_{yy}$ is stored with each keypoint, telling bright blobs from dark ones.
Detection. Keep samples above the threshold that exceed all 26 neighbours in a $3 \times 3 \times 3$ neighbourhood over position and scale, the non-maximum suppression of scale-space detectors, and refine position and size by fitting a 3D quadratic, which the comments credit to Brown and Lowe.
Orientation. At sample points spaced $s$ apart within a radius of $6s$ , compute Haar wavelet responses of side $4s$ in $x$ and $y$ , weight them by a Gaussian of standard deviation $2.5s$ , and treat each pair as a vector. Slide a window of 60° around the circle in steps of 5°, sum the vectors whose angles fall inside it, and take the direction of the longest sum as the orientation.
Descriptor. Resample a square window of side about $20s$ , aligned with the orientation, into $20 \times 20$ cells of side $s$ ; compute Haar responses $d_x$ , $d_y$ of side $2s$ at each cell in the window’s own axes and weight them by a Gaussian of standard deviation $3.3$ cells. Split the window into $4 \times 4$ subregions of $5 \times 5$ cells and record, for each, $\left(\sum d_x, \sum d_y, \sum |d_x|, \sum |d_y|\right)$ , giving 64 numbers. Normalize to unit length, which the code’s comment calls essential for contrast invariance.
Upright and extended variants. With upright set, step 4 is skipped and every window is axis-aligned (U-SURF). With extended set, the sums of $d_x$ and $|d_x|$ are split by the sign of $d_y$ and those of $d_y$ and $|d_y|$ by the sign of $d_x$ , giving 128 numbers.

Mathematical Formulation

Box filters for the Hessian

For a filter of side $L = 3k$ with $k$ odd, $D_{xx}$ has three lobes $k$ pixels wide and $2k - 1$ pixels tall, side by side along $x$ , with weights $1, -2, 1$ ; $D_{yy}$ is its transpose; $D_{xy}$ has four $k \times k$ squares in the quadrants, separated by a one-pixel gap through the centre, with weights $+1$ on one diagonal and $-1$ on the other. This is the geometry of OpenCV’s $9 \times 9$ pattern table ( $k = 3$ ) and of scikit-image’s code at every such size. The response at a pixel is

\det \mathcal{H}_{\text{approx}} = D_{xx} D_{yy} - w^2 D_{xy}^2, \qquad w^2 = 0.81,

the constant both libraries use. Lindeberg (2014) describes the SURF detector as an approximation of the scale-normalized determinant of the Hessian, $t^2 (L_{xx}L_{yy} - L_{xy}^2)$ with $t = \sigma^2$ , in which the filters are normalized to constant $\ell_1$ or Frobenius norm over scales, and argues that this normalization plays the role of the $\gamma$ -normalized derivatives of scale space theory; the scale invariance of the SURF descriptor then follows from that of the detection step, as for SIFT. Lindeberg (2015, §1) likewise describes SURF’s features as approximations of the determinant of the Hessian with Gaussian derivatives replaced by Haar wavelets.

The weight on the mixed term

Why $w = 0.9$ ? One way to balance the two terms (our calculation) is to scale $D_{xy}$ so that its energy relative to $D_{yy}$ matches that of the Gaussian derivatives it replaces:

w = \frac{\lVert G_{xy} \rVert_F \, \lVert D_{yy} \rVert_F}{\lVert G_{yy} \rVert_F \, \lVert D_{xy} \rVert_F}.

For a continuous Gaussian, $\lVert G_{xy} \rVert / \lVert G_{yy} \rVert = 1/\sqrt{3}$ at any scale. For the box filters, $\lVert D_{yy} \rVert_F = \sqrt{6k(2k-1)}$ (three lobes of $k(2k-1)$ pixels with squared weights $1 + 4 + 1$ ) and $\lVert D_{xy} \rVert_F = 2k$ , so

w = \frac{1}{\sqrt{3}} \cdot \frac{\sqrt{6k(2k-1)}}{2k} = \sqrt{\frac{2k - 1}{2k}},

which is $0.913$ for the $9 \times 9$ filter ( $k = 3$ ), close to 0.9, and rises toward 1 as the filters grow. The example below prints this weight with sampled kernels, together with a second criterion, the weight that makes the determinant exact on a quadratic image, which grows much faster with $L$ .

Normalization matters for the weight. scikit-image divides all weights by $L^2$ , which leaves the ratio of $D_{xx}$ to $D_{xy}$ that of the raw $\pm 1$ , $-2$ weights. OpenCV’s resizeHaarPattern divides each lobe’s weight by that lobe’s own area, so on the geometry above $D_{xx}$ is scaled by $1/(k(2k-1))$ and $D_{xy}$ by $1/k^2$ . Relative to the raw filters the fixed 0.81 then acts on $D_{xy}^2$ with an extra factor $\left((2k-1)/k\right)^2$ (our algebra), which is $2.78$ at $k = 3$ and tends to 4. OpenCV also obtains its larger filters by rescaling the $9 \times 9$ pattern and rounding the corners, which gives shorter $D_{xx}$ lobes and a wider $D_{xy}$ gap than the geometry above for every $L \ge 21$ of that form (for $L = 21$ , $D_{xx}$ lobes 11 pixels tall instead of 13 and a gap of 3 pixels; our reading of the code). The sizes of the second and later octaves, such as $18 = 3 \cdot 6$ and $30 = 3 \cdot 10$ , have even $k$ , which the geometry above does not cover.

Descriptor

With $d_x$ and $d_y$ the weighted Haar responses in subregion $i$ of the oriented window, the descriptor is

\mathbf{v} = \frac{\left(\ldots, \textstyle\sum_i d_x, \sum_i d_y, \sum_i |d_x|, \sum_i |d_y|, \ldots\right)}{\lVert \cdot \rVert_2} \in \mathbb{R}^{64}.

The signed sums record the dominant direction of change; the absolute sums distinguish, for example, a single edge from a periodic pattern with the same mean derivative. Because a constant added to the image cancels in every Haar response and a gain cancels in the normalization, the descriptor is invariant to affine changes of intensity.

Parameters

OpenCV 4.14.0’s SURF::create(hessianThreshold = 100, nOctaves = 4, nOctaveLayers = 3, extended = false, upright = false):

hessianThreshold. The minimum response of a keypoint. The class documentation suggests 300 to 500 as a good value depending on image contrast; the response grows with the square of contrast, so thresholds do not transfer between images of different contrast.
nOctaves. The number of octaves; 4 by default, which the code’s comments say matches the default of the original SURF binary.
nOctaveLayers. The number of layers per octave in which maxima are searched; create defaults to 3, while the class’s documentation comment and the internal constructor give 2.
extended and upright. The 128-dimensional descriptor and U-SURF. The documentation recommends upright for stereo pairs and image stitching, where matched features likely have similar angles.

Fixed in the code: the first filter size 9, the increment 6, the orientation radius $6s$ , its Gaussian $2.5s$ , the 60° window and 5° step, the descriptor window $20s$ and its Gaussian $3.3$ cells. scikit-image’s blob_doh(image, min_sigma, max_sigma, num_sigma, threshold, ...) uses the same box filters with a different scale convention: a “sigma” $\sigma$ there means a filter of side $\lfloor 3\sigma \rfloor$ , while OpenCV’s $s = 1.2L/9$ assigns $s = 1.2$ to the $9 \times 9$ filter. Its documentation states that blobs of radius below about 3 pixels cannot be detected because of the box filters.

Complexity

Building the integral image costs one pass over the $N$ pixels. Each box filter response then costs a constant number of reads: 3 rectangles for $D_{xx}$ and $D_{yy}$ , 4 for $D_{xy}$ . A layer sampled with step $p$ costs $O(N/p^2)$ regardless of its filter size, and because the step doubles each octave the whole detector is dominated by the first octave, $O(N)$ per layer. Orientation costs about 113 pairs of Haar responses per keypoint and the descriptor 400, plus the resampling of the window. Matching $N_1$ against $N_2$ descriptors by brute force costs $N_1 N_2$ distances of 64 dimensions; the ORB paper’s timings of ORB, SURF and SIFT are summarized in ORB.

Implementation

The example compares SURF’s box filters with the Gaussian derivatives they replace, and its Haar-wavelet orientation with a Gaussian-derivative one. It uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2, runs in about 3 seconds, and printed identical output on two runs. Rows are $y$ (downward) and columns $x$ ; angles are counter-clockwise as seen on screen. All images are rendered analytically with peak values of order 1. Three detectors are compared: the Gaussian determinant $\sigma^4 \det \mathcal{H}$ from sampled Gaussian derivative kernels, truncated at $4\sigma$ ; the box determinant with weights divided by $L^2$ (“box raw”, scikit-image’s normalization and geometry); and the same box geometry with each lobe divided by its own area (“box per-lobe”: per-lobe normalization as in OpenCV, on scikit-image’s box geometry; OpenCV’s own geometry differs for $L \ge 21$ , as described above, and the variant is not a reproduction of OpenCV’s detector). Both box detectors use $w^2 = 0.81$ , and a box of side $L$ is assigned the Gaussian scale $\sigma_L = 1.2L/9$ .

Kernels. For $L = 9, 15, 27, 51$ : the cosine similarity of $D_{xx}$ with $G_{xx}$ and of $D_{xy}$ with $G_{xy}$ at $\sigma_L$ ; the $\sigma$ in $[0.7, 1.4]\,\sigma_L$ (141 values) at which $D_{xx}$ is most similar to $G_{xx}$ ; the Frobenius weight above, with sampled kernels; and the quadratic weight $D_{xx}(x^2/2) / D_{xy}(xy)$ , which makes the determinant exact on quadratic images, for both normalizations.
Blobs. (a) Isotropic Gaussian blobs of standard deviation $\sigma_b = 2$ to 8 on a $257 \times 257$ image: the response at the centre for 16 box sizes $L = 9, 15, \ldots, 99$ , the selected scale by a parabola through the maximum and its neighbours, as a fraction of $\sigma_b$ , and the peak value. The exact model is the continuous $t^2 t_0^2 / (t_0 + t)^4$ at the centre, with $t = \sigma_L^2$ and $t_0 = \sigma_b^2$ , on the same grid. (b) An elongated Gaussian blob with axis variances $2\sigma_0^2$ and $\sigma_0^2/2$ , rotated by 15°, 30° and 45°: the peak over scale at the centre relative to its value unrotated, and the spread max/min − 1. A library check compares scikit-image’s hessian_matrix_det with the hand-written box determinant on the $\sigma_b = 4$ image.
Rotation. A texture of 800 Gaussian bumps (standard deviations 1.5 to 4 pixels, signed amplitudes), $256 \times 256$ , rendered exactly at rotations of 0° to 90° about the image centre, with 100 keypoints at fixed scale $s = 2$ ( $L = 15$ ). Orientation follows OpenCV’s recipe (Haar wavelets of side $4s$ at rounded positions), and is repeated with Gaussian derivatives of standard deviation $4s/\sqrt{12}$ , the standard deviation of a box $4s$ wide. The descriptor follows OpenCV’s recipe (a $42 \times 42$ window sampled bilinearly along the orientation and averaged to $21 \times 21$ cells). The script prints the median orientation error and the median descriptor distance to the unrotated image for U-SURF, for SURF with its estimated orientation, and for SURF given the true orientation (the unrotated estimate plus the rotation).

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 skimage.feature import hessian_matrix_det

LS = np.arange(9, 100, 6)            # box filter sizes 9, 15, ..., 99
SIG = 1.2 * LS / 9                   # Gaussian scale of each size, as OpenCV's s = 1.2 L / 9

# ---------------------------------------------------------------- kernels
def box_kernels(L, perlobe=False):
    """Dxx, Dyy, Dxy box filters of size L (rows y, cols x). Lobes are L/3 wide; the
    Dxx lobes are 2L/3 - 1 tall, the Dxy lobes are L/3 squares with a one-pixel gap.
    perlobe=False: weights 1, -2, 1 and +-1 divided by L^2 (scikit-image);
    perlobe=True: each lobe divided by its own area (OpenCV 4.14.0)."""
    k, c = L // 3, L // 2
    h = 2 * k - 1
    xx = [(-(h // 2), h // 2 + 1, -c, -c + k, 1), (-(h // 2), h // 2 + 1, -c + k, -c + 2 * k, -2),
          (-(h // 2), h // 2 + 1, -c + 2 * k, c + 1, 1)]
    yy = [(c0, c1, r0, r1, w) for (r0, r1, c0, c1, w) in xx]
    xy = [(-k, 0, -k, 0, 1), (-k, 0, 1, k + 1, -1), (1, k + 1, -k, 0, -1), (1, k + 1, 1, k + 1, 1)]
    out = []
    for lobes in (xx, yy, xy):
        K = np.zeros((L, L))
        for r0, r1, c0, c1, w in lobes:
            K[c + r0:c + r1, c + c0:c + c1] += w / ((r1 - r0) * (c1 - c0) if perlobe else L * L)
        out.append(K)
    return out

def gauss_kernels(sig, R=None):
    """Sampled Gxx, Gyy, Gxy of a unit-mass Gaussian, truncated at radius R (default 4 sigma)."""
    R = int(np.ceil(4 * sig)) if R is None else R
    y, x = np.mgrid[-R:R + 1, -R:R + 1].astype(float)
    g = np.exp(-(x ** 2 + y ** 2) / (2 * sig ** 2)) / (2 * np.pi * sig ** 2)
    return [(x ** 2 / sig ** 4 - 1 / sig ** 2) * g, (y ** 2 / sig ** 4 - 1 / sig ** 2) * g, x * y / sig ** 4 * g]

def pad(K, R):
    p = R - K.shape[0] // 2
    return np.pad(K, p)

cos = lambda a, b: (a * b).sum() / np.sqrt((a * a).sum() * (b * b).sum())

# ---- 1. Box filters against Gaussian second derivatives
print("1. box filters against Gaussian second derivatives at sigma = 1.2 L / 9")
print("    L  sigma  corr xx  corr xy  best sigma xx  Frobenius w  quadratic w (raw, per-lobe)")
for L in (9, 15, 27, 51):
    s = 1.2 * L / 9
    R = max(int(np.ceil(4 * 1.4 * s)), L // 2)        # room for the widest sigma searched
    gxx, gyy, gxy = gauss_kernels(s, R)
    dxx, dyy, dxy = (pad(K, R) for K in box_kernels(L))
    grid = s * np.linspace(0.7, 1.4, 141)
    best = grid[np.argmax([cos(dxx, gauss_kernels(q, R)[0]) for q in grid])]
    wF = np.linalg.norm(gxy) * np.linalg.norm(dyy) / (np.linalg.norm(gyy) * np.linalg.norm(dxy))
    y, x = np.mgrid[-R:R + 1, -R:R + 1].astype(float)
    wq = [(K[0] * x ** 2 / 2).sum() / (K[2] * x * y).sum()
          for K in ([pad(B, R) for B in box_kernels(L, p)] for p in (False, True))]
    print(f"  {L:3d}  {s:5.2f}   {cos(dxx, gxx):.4f}   {cos(dxy, gxy):.4f}   {best / s:.3f} x sigma"
          f"      {wF:.4f}      {wq[0]:.4f}, {wq[1]:.4f}")

# ---- 2. Responses at the centre of blobs, over scale
def centre_responses(img, c=128):
    """sigma^4 det H (Gaussian), and box det with w^2 = 0.81 (raw and per-lobe), for every L."""
    g, braw, blobe = [], [], []
    for L, s in zip(LS, SIG):
        kxx, kyy, kxy = gauss_kernels(s)
        R = kxx.shape[0] // 2
        p = img[c - R:c + R + 1, c - R:c + R + 1]
        a, b, d = (p * kxx).sum(), (p * kyy).sum(), (p * kxy).sum()
        g.append(s ** 4 * (a * b - d * d))
        p = img[c - L // 2:c + L // 2 + 1, c - L // 2:c + L // 2 + 1]
        for out, pl in ((braw, False), (blobe, True)):
            a, b, d = ((p * K).sum() for K in box_kernels(L, pl))
            out.append(a * b - 0.81 * d * d)
    return np.array(g), np.array(braw), np.array(blobe)

def interp_peak(v):
    """Peak position in index units and value, by a parabola through the max and its neighbours."""
    i = int(np.argmax(v))
    if 0 < i < len(v) - 1:
        a, b, c = v[i - 1], v[i], v[i + 1]
        d = 0.5 * (a - c) / (a - 2 * b + c)
        return i + d, b - 0.25 * (a - c) * d
    return float(i), v[i]

n = 257
yy, xx = np.mgrid[:n, :n].astype(float) - 128      # centre pixel (128, 128)
print("\n2a. isotropic Gaussian blobs (peak 1): selected sigma / sigma_b, and peak response")
print("  sigma_b   Gaussian   box      exact model    peak: Gaussian  box raw  box per-lobe")
for sb in (2, 3, 4, 6, 8):
    img = np.exp(-(xx ** 2 + yy ** 2) / (2 * sb ** 2))
    g, braw, blobe = centre_responses(img)
    t0, t = sb ** 2, SIG ** 2
    sel = lambda v: 1.2 * (9 + 6 * interp_peak(v)[0]) / 9 / sb
    print(f"  {sb:5d}     {sel(g):.4f}   {sel(braw):.4f}   {sel(t ** 2 * t0 ** 2 / (t0 + t) ** 4):.4f}"
          f"          {interp_peak(g)[1]:.5f}  {interp_peak(braw)[1]:.5f}  {interp_peak(blobe)[1]:.4f}")

print("\n2b. elongated blob (axis ratio 2) rotated by theta: peak response / its value at theta = 0")
print("  sigma0  detector       15 deg   30 deg   45 deg   max/min - 1")
for s0 in (2, 4, 8):
    rows = {"Gaussian": [], "box raw": [], "box per-lobe": []}
    for th in np.deg2rad([0, 15, 30, 45]):
        u = np.cos(th) * xx + np.sin(th) * yy
        v = -np.sin(th) * xx + np.cos(th) * yy
        img = np.exp(-u ** 2 / (2 * 2 * s0 ** 2) - v ** 2 / (2 * s0 ** 2 / 2))
        for key, resp in zip(rows, centre_responses(img)):
            rows[key].append(interp_peak(resp)[1])
    for key, vals in rows.items():
        r = np.array(vals) / vals[0]
        print(f"  {s0:5d}   {key:13s} " + "".join(f"{x:8.4f} " for x in r[1:]) + f"  {r.max() / r.min() - 1:.2e}")

# ---- library check: hand box determinant against scikit-image's hessian_matrix_det
img = np.exp(-(xx ** 2 + yy ** 2) / (2 * 4.0 ** 2))
print("\nlibrary check: scikit-image hessian_matrix_det(image, sigma=L/3) against the hand raw box det")
for L in (9, 15, 27):
    kxx, kyy, kxy = box_kernels(L)
    f = lambda K: ndi.correlate(img, K, mode="constant")
    mine = f(kxx) * f(kyy) - 0.81 * f(kxy) ** 2
    sk = hessian_matrix_det(img, sigma=L / 3)
    m = 60
    d0 = np.abs(sk - mine)[m:-m, m:-m].max()
    d1 = np.abs(sk[m:-m, m:-m] - mine[m + 1:n - m + 1, m + 1:n - m + 1]).max()
    print(f"  L={L}: max |difference| / peak, same pixel {d0 / mine.max():.1e}, "
          f"scikit-image at (r, c) vs hand at (r+1, c+1) {d1 / mine.max():.1e}")

# ---- 3. Orientation and descriptor under rotation (s = 2, filter size 15)
s, c0 = 2, 127.5                                   # keypoint scale; rotation centre (image centre)
rng = np.random.default_rng(0)
nb = 800
rad, ang = 200 * np.sqrt(rng.random(nb)), 2 * np.pi * rng.random(nb)
bx, by = c0 + rad * np.cos(ang), c0 + rad * np.sin(ang)
bs, ba = rng.uniform(1.5, 4.0, nb), rng.uniform(-1, 1, nb)
kr, ka = 60 * np.sqrt(rng.random(100)), 2 * np.pi * rng.random(100)
kx, ky = np.round(c0 + kr * np.cos(ka)), np.round(c0 + kr * np.sin(ka))   # integer keypoints

def rot(x, y, th):
    """Rotate points about the image centre by th, counter-clockwise on screen (y points down)."""
    dx, dy = x - c0, y - c0
    return c0 + np.cos(th) * dx + np.sin(th) * dy, c0 - np.sin(th) * dx + np.cos(th) * dy

def render(th):
    X, Y = rot(bx, by, th)
    yy3, xx3 = np.mgrid[:256, :256].astype(float)
    im = np.zeros((256, 256))
    for x, y, sg, a in zip(X, Y, bs, ba):
        im += a * np.exp(-((xx3 - x) ** 2 + (yy3 - y) ** 2) / (2 * sg ** 2))
    return im

# orientation samples: grid points (i, j) with i^2 + j^2 <= 36, step s, Gaussian weight sigma 2.5 s
oi, oj = np.mgrid[-6:7, -6:7]
inside = oi ** 2 + oj ** 2 <= 36
oi, oj = oi[inside], oj[inside]
g1 = np.exp(-np.arange(-6, 7) ** 2 / (2 * 2.5 ** 2)); g1 /= g1.sum()
ow = g1[oj + 6] * g1[oi + 6]
d1 = np.exp(-(np.arange(20) - 9.5) ** 2 / (2 * 3.3 ** 2)); d1 /= d1.sum()
DW = np.outer(d1, d1)                              # descriptor weights, sigma 3.3 samples

def haar(S, x0, y0, size):
    """Haar responses of a size x size wavelet with top-left (x0, y0), from the integral image S
    (zero row and column prepended). Returns (vx, vy) as differences of means, vy positive upwards."""
    h = size // 2
    box = lambda r0, c0_, r1, c1: S[r1, c1] - S[r0, c1] - S[r1, c0_] + S[r0, c0_]
    right = box(y0, x0 + h, y0 + size, x0 + size); left = box(y0, x0, y0 + size, x0 + h)
    top = box(y0, x0, y0 + h, x0 + size); bottom = box(y0 + h, x0, y0 + size, x0 + size)
    a = h * size
    return (right - left) / a, (top - bottom) / a

def window_angle(X, Y):
    """Sliding 60-degree window over sample angles in 5-degree steps, as OpenCV 4.14.0."""
    angle = np.degrees(np.arctan2(Y, X)) % 360
    best, bx_, by_ = -1.0, 0.0, 0.0
    for i in range(0, 360, 5):
        d = np.abs(np.round(angle) - i)
        m = (d < 30) | (d > 330)
        sx, sy = X[m].sum(), Y[m].sum()
        if sx * sx + sy * sy > best:
            best, bx_, by_ = sx * sx + sy * sy, sx, sy
    return np.arctan2(by_, bx_)                      # counter-clockwise on screen

def orient_haar(S, x, y):
    size = 2 * int(np.round(2 * s))                  # wavelet size 4s, even
    xs = np.round(x + oi * s - (size - 1) / 2).astype(int)
    ys = np.round(y + oj * s - (size - 1) / 2).astype(int)
    vx, vy = haar(S, xs, ys, size)
    return window_angle(vx * ow, vy * ow)

def orient_gauss(Gx, Gy, x, y):
    pts = [y + oj * s, x + oi * s]
    vx = ndi.map_coordinates(Gx, pts, order=1)
    vy = -ndi.map_coordinates(Gy, pts, order=1)       # positive upwards
    return window_angle(vx * ow, vy * ow)

def describe(im, x, y, phi):
    """64-D SURF descriptor: a 42 x 42 window along the orientation phi, averaged to 21 x 21 cells of
    size s, Haar responses of size 2s, Gaussian weights, 4 x 4 subregions of (dx, dy, |dx|, |dy|)."""
    a = np.arange(42) - 20.5
    A, B = np.meshgrid(a, a)                         # A along the orientation, B perpendicular (down)
    px = x + A * np.cos(phi) + B * np.sin(phi)
    py = y - A * np.sin(phi) + B * np.cos(phi)
    w = ndi.map_coordinates(im, [py, px], order=1, mode="nearest")
    P = w.reshape(21, 2, 21, 2).mean(axis=(1, 3))
    DX = (P[:-1, 1:] - P[:-1, :-1] + P[1:, 1:] - P[1:, :-1]) * DW
    DY = (P[1:, :-1] - P[:-1, :-1] + P[1:, 1:] - P[:-1, 1:]) * DW
    f = lambda M: M.reshape(4, 5, 4, 5).sum(axis=(1, 3)).ravel()
    v = np.stack([f(DX), f(DY), f(np.abs(DX)), f(np.abs(DY))], axis=1).ravel()
    return v / np.linalg.norm(v)

def features(im, xs, ys, phis_oracle=None):
    S = np.pad(im.cumsum(0).cumsum(1), ((1, 0), (1, 0)))
    sg = 4 * s / np.sqrt(12)                         # std of a 4s-wide box
    Gx, Gy = ndi.gaussian_filter(im, sg, order=(0, 1)), ndi.gaussian_filter(im, sg, order=(1, 0))
    ph = np.array([orient_haar(S, x, y) for x, y in zip(xs, ys)])
    pg = np.array([orient_gauss(Gx, Gy, x, y) for x, y in zip(xs, ys)])
    up = np.array([describe(im, x, y, 0.0) for x, y in zip(xs, ys)])
    est = np.array([describe(im, x, y, p) for x, y, p in zip(xs, ys, ph)])
    orc = None if phis_oracle is None else np.array([describe(im, x, y, p) for x, y, p in zip(xs, ys, phis_oracle)])
    return ph, pg, up, est, orc

wrap = lambda a: np.degrees(np.angle(np.exp(1j * a)))
im0 = render(0.0)
ph0, pg0, up0, est0, _ = features(im0, kx, ky)
inframe = ((bx > -10) & (bx < 266) & (by > -10) & (by < 266)).sum()
reach = (np.hypot(kx - c0, ky - c0) + 21 * np.sqrt(2) + 1).max()
print(f"\n3. data check: {inframe} bumps within 10 px of the frame, image std {im0.std():.3f}, "
      f"farthest window corner {reach:.1f} px from the centre (limit 127.5)")
nm = np.median([np.linalg.norm(est0[i] - est0[j]) for i in range(100) for j in range(i + 1, 100)])
print(f"   median descriptor distance between different keypoints (SURF, theta = 0): {nm:.3f}")
print("   theta  |orientation error| median, deg   descriptor distance to theta = 0, median")
print("            Haar     Gaussian derivative     U-SURF   SURF   SURF, true orientation")
for deg in (0, 15, 30, 45, 60, 75, 90):
    th = np.deg2rad(deg)
    X, Y = rot(kx, ky, th)
    im = render(th)
    ph, pg, up, est, orc = features(im, X, Y, ph0 + th)
    eh, eg = np.abs(wrap(ph - ph0 - th)), np.abs(wrap(pg - pg0 - th))
    dist = lambda a, b: np.median(np.linalg.norm(a - b, axis=1))
    print(f"   {deg:3d}     {np.median(eh):6.3f}   {np.median(eg):6.3f}                  "
          f"{dist(up, up0):.3f}   {dist(est, est0):.3f}  {dist(orc, est0):.3f}")

Output:

1. box filters against Gaussian second derivatives at sigma = 1.2 L / 9
    L  sigma  corr xx  corr xy  best sigma xx  Frobenius w  quadratic w (raw, per-lobe)
    9   1.20   0.7147   0.7784   1.400 x sigma      0.9127      0.9375, 0.5625
   15   2.00   0.6757   0.8278   1.400 x sigma      0.9487      1.2500, 0.6944
   27   3.60   0.6562   0.8378   1.400 x sigma      0.9718      1.5300, 0.8100
   51   6.80   0.6469   0.8353   1.400 x sigma      0.9852      1.7315, 0.8920

2a. isotropic Gaussian blobs (peak 1): selected sigma / sigma_b, and peak response
  sigma_b   Gaussian   box      exact model    peak: Gaussian  box raw  box per-lobe
      2     1.0656   0.6000   1.0656          0.06300  0.02783  0.8116
      3     1.0211   0.7331   1.0211          0.06299  0.03151  0.7781
      4     1.0070   0.7214   1.0070          0.06214  0.03152  0.7386
      6     1.0087   0.6987   1.0087          0.06251  0.03174  0.7094
      8     1.0015   0.6882   1.0015          0.06245  0.03220  0.7025

2b. elongated blob (axis ratio 2) rotated by theta: peak response / its value at theta = 0
  sigma0  detector       15 deg   30 deg   45 deg   max/min - 1
      2   Gaussian        1.0000   1.0000   1.0000   4.94e-06
      2   box raw         0.9972   1.0010   1.0076   1.04e-02
      2   box per-lobe    0.9103   0.7318   0.6423   5.57e-01
      4   Gaussian        1.0000   1.0000   1.0000   9.53e-06
      4   box raw         0.9906   0.9821   0.9838   1.82e-02
      4   box per-lobe    0.9126   0.7458   0.6671   4.99e-01
      8   Gaussian        1.0000   1.0000   1.0000   2.13e-05
      8   box raw         0.9918   0.9863   0.9900   1.39e-02
      8   box per-lobe    0.9111   0.7388   0.6558   5.25e-01

library check: scikit-image hessian_matrix_det(image, sigma=L/3) against the hand raw box det
  L=9: max |difference| / peak, same pixel 4.3e-01, scikit-image at (r, c) vs hand at (r+1, c+1) 8.5e-15
  L=15: max |difference| / peak, same pixel 3.7e-01, scikit-image at (r, c) vs hand at (r+1, c+1) 3.0e-15
  L=27: max |difference| / peak, same pixel 2.9e-01, scikit-image at (r, c) vs hand at (r+1, c+1) 2.0e-15

3. data check: 468 bumps within 10 px of the frame, image std 0.214, farthest window corner 91.0 px from the centre (limit 127.5)
   median descriptor distance between different keypoints (SURF, theta = 0): 0.998
   theta  |orientation error| median, deg   descriptor distance to theta = 0, median
            Haar     Gaussian derivative     U-SURF   SURF   SURF, true orientation
     0      0.000    0.000                  0.000   0.000  0.000
    15      8.632    4.095                  0.396   0.225  0.001
    30     13.042    4.069                  0.700   0.333  0.002
    45     17.576    4.469                  0.917   0.420  0.002
    60     13.482    4.099                  1.065   0.323  0.002
    75      6.427    5.446                  1.111   0.171  0.001
    90      0.000    0.000                  1.135   0.000  0.000

The example was planned with hypotheses fixed before its pilot; the outcomes are given here under every printed condition.

Box filters against Gaussian derivatives. The box kernels are crude: at $\sigma_L$ their cosine similarity with the Gaussian derivatives is 0.65 to 0.71 for $D_{xx}$ and 0.78 to 0.84 for $D_{xy}$ , and the $D_{xx}$ similarity falls slightly as $L$ grows. The planned hypothesis that both similarities exceed 0.9 at every $L$ failed, and so did the hypothesis that the best-matching $\sigma$ lies within 10% of $\sigma_L$ : at every $L$ the best $\sigma$ is the top of the searched range, so it is at least $1.4\,\sigma_L$ within the range searched. The Frobenius weight with sampled kernels, 0.9127 to 0.9852, agrees with the continuous formula above to within 0.0002 and rises toward 1, as predicted; the constant 0.9 is close to the value for the $9 \times 9$ filter. The quadratic weights, as predicted, depend strongly on the normalization and on $L$ (0.9375 to 1.7315 raw, 0.5625 to 0.8920 per-lobe): a single constant cannot make the determinant exact on quadratic images at every size under either normalization.

Scale and response on blobs. The Gaussian determinant selects exactly the scale the continuous model predicts on the same grid (the two columns agree to four decimals, including the grid-interpolation offset of 1.0656 at $\sigma_b = 2$ ), and its peaks lie within 1% of the model’s maximum, $1/16$ at $t = t_0$ . The box detector, read with $\sigma_L = 1.2L/9$ , selects 0.69 to 0.73 of $\sigma_b$ for $\sigma_b \ge 3$ and 0.60 at $\sigma_b = 2$ ; the planned hypothesis that it selects within 10% of $\sigma_b$ for $\sigma_b \ge 3$ failed. In our reading this agrees with part 1: a box of side $L$ behaves like a Gaussian wider than $1.2L/9$ , so a blob is selected at a smaller nominal scale. The ratio is not constant either: it drifts from 0.733 to 0.688 between $\sigma_b = 3$ and 8. The other clause held: the box raw peaks vary by 2.2% for $\sigma_b = 3$ to 8 (0.03151 to 0.03220), so the $1/L^2$ normalization keeps responses comparable across scales; the smallest blob gives 0.02783. The per-lobe values are on a different scale, and fall from 0.81 to 0.70 as $\sigma_b$ grows.

Rotation of an elongated blob. The Gaussian determinant is rotation invariant in theory, and its printed spread is at most $2.1 \times 10^{-5}$ . Both box detectors vary by more than ten times as much at every $\sigma_0$ , as predicted: the raw detector by 1.0% to 1.8%, the per-lobe variant by 50% to 56%, its response at 45° falling to 0.64 to 0.67 of its unrotated value. The planned hypothesis that this per-lobe variant varies less than the raw detector at $\sigma_0 = 4$ and 8 failed at both. The derivation above shows that, relative to the raw filters, the per-lobe normalization multiplies the weight on $D_{xy}^2$ by $((2k-1)/k)^2$ , between 2.78 and 4; in our reading this over-weighted mixed term, largest for a blob at 45°, is what lowers the response. These numbers describe the variant in this script, which has scikit-image’s lobe geometry; they are not a measurement of OpenCV’s detector, whose larger filters have different lobes.

Library check. scikit-image’s hessian_matrix_det agrees with the hand-written raw determinant to rounding error, below $10^{-14}$ of the peak, but only when its value at pixel $(r, c)$ is compared with the hand value at $(r + 1, c + 1)$ ; at the same pixel the difference is 29% to 43% of the peak. The cause is in its code: _integ(img, r, c, rl, cl) reads the integral image at $(r, c)$ and $(r + rl, c + cl)$ without a zero row and column, so the box covers rows $r + 1$ to $r + rl$ , and every lobe ends up centred one pixel down and to the right of the pixel it is stored at. Positions from blob_doh in release 0.25.2 carry the same offset.

Orientation and descriptor under rotation. At 0° every error and distance is exactly zero, and at 90° the Gaussian-derivative orientation error is zero, as required. The Haar-wavelet orientation error is larger than the Gaussian-derivative one at every angle from 15° to 75°, as predicted: it grows from 8.6° at 15° to 17.6° at 45° and falls back to 6.4° at 75°, while the Gaussian estimator stays between 4.1° and 5.4°. The prediction that the Haar error stays above zero at 90° failed: it is 0.000. The descriptor computed with the true orientation stays within 0.002 of the unrotated one at every angle, so in this example the descriptor’s own construction, axis-aligned Haar wavelets on a resampled window, is almost exactly rotation invariant, and the change in the SURF descriptor (0.17 to 0.42, largest at 45°) comes from the orientation error. As predicted, SURF with its estimated orientation stays below half the median distance between different keypoints (0.998) at every angle, and the true orientation is never worse than the estimated one. The upright descriptor’s distance rises with every step of rotation, as predicted, reaching 0.40 already at 15°; the prediction that it exceeds the different-keypoint median by 45° failed (0.917 at 45°), and it does so from 60° (1.065). Feature descriptors compares a gradient-histogram descriptor with and without orientation assignment, and ORB steered and unsteered BRIEF; here the upright column serves only as a reference, and what this example adds is the true-orientation column, which separates the orientation estimate from the descriptor.

The plan first described part 3 as “the descriptor’s sums of Haar responses under rotation with and without orientation assignment”; before the pilot it was changed to the angle sweep with Haar and Gaussian orientation and the true-orientation column, to avoid repeating those two examples. Nothing was changed after the pilot or the first run.

Properties and Behavior

Scale invariance. Detection at scale-space maxima of a scale-normalized measure makes the selected filter size proportional to the size of the structure, and computing the orientation and descriptor in windows proportional to it makes them scale invariant (Lindeberg, 2014). The example shows the box detector keeping its responses comparable across blob sizes, with a selected nominal scale that is a nearly, but not exactly, fixed fraction of the blob’s.
Rotation. The Hessian determinant is rotation invariant only for the Gaussian filters; the box filters, which are aligned with the pixel grid, make the response depend on the orientation of elongated structures, slightly with the raw weights in the example and strongly in the per-lobe variant. The descriptor relies on the orientation estimate, whose Haar-wavelet version was the larger error source in the example. Rublee et al. (2011), as described in the ORB paper article, found in their synthetic rotation test that SIFT outperformed SURF, which showed quantization effects at multiples of 45° that they attribute to its Haar-wavelet construction.
Edges and ridges. As a determinant-of-Hessian detector, SURF responds weakly to straight edges and ridges, where one principal curvature is near zero; blob detection measures this for the Gaussian determinant.
Polarity. The determinant is positive at bright and dark blobs alike, and the sign of the trace tells them apart; saddle points give negative determinants and are not detected as maxima.
Contrast. The descriptor is invariant to affine intensity changes. The detector’s response scales with the square of contrast.
Descriptor quality. Lindeberg (2015) replaced SIFT’s and SURF’s image measurements by Gaussian derivatives (“Gauss-SIFT” and “Gauss-SURF”, the latter with the sums of $L_x$ , $|L_x|$ , $L_y$ , $|L_y|$ over $4 \times 4$ subwindows and SIFT-like orientation normalization) to compare the descriptors on equal interest points. On photographs of 12 posters with known homographies, Gauss-SIFT gave higher efficiency (the fraction of interest points that led to accepted matches) and lower 1 − precision than Gauss-SURF for every interest point detector. For his analogue of SURF’s detector, scale-space extrema of the determinant of the Hessian with a complementary threshold on his $\mathcal{D}_1 L > 0$ , the average efficiency was 0.7529 with Gauss-SURF against 0.7678 with Gauss-SIFT, and the average 1 − precision 0.0391 against 0.0339 (his Tables 3 and 2, averages over scaling and foreshortening).

Limitations

Approximation error. Box filters are coarse approximations of Gaussian derivatives: the example measures cosine similarities of 0.65 to 0.84 and a mapping from filter size to Gaussian scale that the conventional $1.2L/9$ does not capture for blobs. scikit-image warns that its box determinant is not accurate for small scales and cannot detect blobs of radius below about 3 pixels.
Implementation differences. The determinant’s weighting depends on how the lobes are normalized, and the two common implementations differ in normalization, lobe geometry at larger sizes, scale convention and (in scikit-image 0.25.2) a one-pixel position offset, so their responses and scales are not interchangeable.
Orientation errors. The single dominant orientation fails on keypoints without a clear gradient direction; a wrong orientation spoils the descriptor regardless of its own invariance.
Viewpoint. Nothing in SURF models affine distortion beyond what the determinant’s scale estimate tolerates; affine-adapted detectors are described in feature detection.
Licensing. OpenCV builds SURF only with the OPENCV_ENABLE_NONFREE option; otherwise SURF::create raises an error stating that the algorithm is patented and excluded in that configuration (as of November 2026, opencv_contrib 4.14.0). As the ORB paper article describes, Rublee et al. (2011) presented ORB as free of the licensing restrictions of SIFT and SURF.

Variants

U-SURF. The upright descriptor skips orientation assignment; it suits cameras that do not rotate about their optical axis, a trade-off feature descriptors discusses.
Extended SURF. The 128-dimensional descriptor splits each sum by the sign of the other derivative.
Gauss-SURF. Lindeberg’s (2015) analogue computes the descriptor from Gaussian derivatives at the detection scale, and pairs it with generalized scale-space interest points, some linked over scales; in his evaluation scale linking generally raised efficiency over scale-space extrema.
scikit-image’s blob_doh and hessian_matrix_det. The box-filter determinant as a blob detector, without orientation or descriptor, with the scale convention and offset described above.
GPU versions. OpenCV’s opencv_contrib also contains CUDA and OpenCL implementations of SURF (surf.cuda.cpp, surf.ocl.cpp).
Binary successors. ORB and the other binary descriptors replaced floating-point descriptors where speed and memory matter, and match by Hamming distance.

Use

Several published articles describe uses of SURF. Wang and Schmid’s improved dense trajectories match SURF descriptors between consecutive frames, together with optical flow, to estimate camera motion as a homography; FAB-MAP builds its bag of words from SIFT or SURF features, and Gálvez-López and Tardós replaced SIFT and SURF by FAST and BRIEF to cut the cost of loop closure; and ORB-SLAM’s authors ruled out SIFT and SURF as too slow for their per-frame budget, as described in ORB.