Rebuilding SURF from the paper up

SURF is usually one line: cv2.xfeatures2d_SURF.create(). The last post already corrected a wrong claim in this repository’s README: neither FERN nor SURF comes from OpenCV here, both are written from the paper up, and only the AGAST detector used by FERN is a library call. This post is SURF’s turn: the box-filter approximation to a Gaussian second derivative, the mysterious 0.912 correction, a 3×3 Hessian solve for sub-pixel and sub-scale refinement that most from-scratch implementations quietly skip, and a 4×4×4 descriptor layout. Building it is where the library’s parameters (that hessianThreshold, that nOctaves) stop being magic numbers.

It also has a specific, teachable failure, and a fun one to catch years later. The report argues, correctly, that SURF chooses box filters over Gaussians because an integral image evaluates a box sum in four lookups regardless of size. Working back through the code behind that same report, I found it never actually implements that. Fixing it in Rust, and measuring the fix in your browser on the widget below, is this post’s payoff. Credit throughout: SURF is Bay, Tuytelaars and Van Gool, SURF: Speeded Up Robust Features (2006), reference [14]/[15] of the report’s own bibliography.

Detection: the determinant of a Hessian, at many blur levels

The report’s own §5.2, transcribed:

“The Speeded Up Robust Features algorithm (SURF), provides a method for both feature detection and description. This detector can be implemented as rotationally invariant (SURF) or not (USURF). Due to this, the upright variant was chosen as it is robust to variations in rotation of up to 15°, which is more than what is needed for the stereo correspondence problem. […] The feature detector is based on the Hessian-Laplace detector. Here the determinant of the Hessian matrix is used to extract key points. This is applied to the image with various levels of Gaussian blur $L(\mathbf{x}, \sigma)$ . The level of $\sigma$ can then be used to extract the feature’s scale.”

H(\mathbf{x}, \sigma) = \begin{pmatrix} L_{xx}(\mathbf{x},\sigma) & L_{xy}(\mathbf{x},\sigma) \\ L_{xy}(\mathbf{x},\sigma) & L_{yy}(\mathbf{x},\sigma) \end{pmatrix}

A point is a blob candidate where $\det(H)$ (which is large and positive only where the surface curves the same way in both directions, a peak or a pit, not a saddle) is large. Computing $L_{xx}$ , $L_{yy}$ , $L_{xy}$ exactly means convolving with the second derivatives of a Gaussian at the chosen $\sigma$ , once per scale, over the whole image. That’s the part SURF exists to avoid.

Why box filters: the report’s own argument

“In this paper, the use of the Gaussian blur as opposed to a box filter is bought into question. They argue that a Gaussian filter does not ensure that aliasing will not occur. Hence, the use of a Gaussian filter may not be important. They suggest the use of a Box filter as it can be easily evaluated at various levels with the use of an integral image. This allows the features to be detected at a much faster speed. To compensate for this change, the determinant must be weighted with a factor which is calculated to be approximately 0.912.”

\det(H_{\text{approx}}) = D_{xx} D_{yy} - (0.912\, D_{xy})^2

$D_{xx}$ , $D_{yy}$ , $D_{xy}$ are rectangle-shaped stand-ins for the true derivatives: cheap to build, coarse by construction, and the 0.912 factor exists purely to correct the resulting energy mismatch back toward the true Gaussian answer. SURF.py:62-77 builds exactly three such filters per scale:

for i, l in enumerate([3, 5, 7, 9, 13, 17, 25, 33, 49, 65]):
    KDyp = np.ones((3 * l, 2 * l + 1))
    KDyp[l:2 * l] = -2
    KDxp = KDyp.T
    KDxy = np.zeros((2 * l + 1, 2 * l + 1))
    KDxy[:l, :l] = 1
    KDxy[l + 1:, l + 1:] = 1
    KDxy[l + 1:, :l] = -1
    KDxy[:l, l + 1:] = -1
    Dy = cv2.filter2D(img_grey, cv2.CV_64F, KDyp, borderType=cv2.BORDER_REFLECT)
    Dx = cv2.filter2D(img_grey, cv2.CV_64F, KDxp, borderType=cv2.BORDER_REFLECT)
    Dxy = cv2.filter2D(img_grey, cv2.CV_64F, KDxy, borderType=cv2.BORDER_REFLECT)
    scale = 1 / (l * l * l * l)
    detH = scale * (Dx * Dy - (0.912 * Dxy) ** 2)

(Confusingly, the code’s Dy is the filter that approximates $L_{yy}$ and Dx approximates $L_{xx}$ : the names describe which axis the three bands stack along, not which second derivative results.) Ten scales run, l = 3, 5, 7, 9, 13, 17, 25, 33, 49, 65; only the inner eight ever produce a keypoint, the two extremes exist purely so every usable scale has a scale-neighbour on both sides for refinement, below. These are literally three rectangles stacked (KDyp, KDxp) or four rectangles in a checkerboard (KDxy), and seeing them beats describing them, so here they are at the smallest scale in the bank, $l=5$ :

The middle row and column of KDxy are left at zero: a deliberate one-pixel gap between the four quadrants, visible above as the thin seam. A blob under KDyp/KDxp reads as a bright centre band flanked by two dark ones; a saddle under KDxy reads as two diagonal quadrants pulling opposite ways.

The argument the code doesn’t cash in

Here is the entire justification for box filters, one more time: “they can be easily evaluated at various levels with the use of an integral image. This allows the features to be detected at a much faster speed.” An integral image (a summed-area table) makes any rectangle sum four lookups, independent of its size. SURF.py:73-75 (quoted in full above) evaluates KDyp, KDxp and KDxy with cv2.filter2D, a general 2-D convolution, at every one of the ten scales. At $l=49$ that is a $147\times99$ and a $99\times147$ kernel plus a $99\times99$ one, applied by direct convolution: work proportional to kernel area, at every pixel. The integral image the report’s own paragraph describes is never built. The report explains, correctly, why SURF should be fast. The code it accompanies doesn’t implement the reason, which is a satisfying thing to catch years after writing it.

wasm/crates/stereo-wasm/src/surf.rs implements the Hessian response two ways behind one flag, compute_layer: direct rectangle summation (the honest equivalent of what filter2D does to a piecewise-constant kernel: touch every pixel in every rectangle) or a lookup against post 29’s summed-area table (image::integral_i64, integral_rect, already sitting in this crate specifically for this). Both paths read the same rectangles: response_at in surf.rs derives them once from SURF.py’s kernel geometry, so there is only one place either method could disagree with the other:

// Dyy ("Dy" in SURF.py): three l-tall, (2l+1)-wide bands stacked at y-l, y, y+l
let band_y = |cy: i32| s.sum(x - hw, cy - hh, x + hw, cy + hh);
let dyy = (band_y(y - l) - 2 * band_y(y) + band_y(y + l)) as f64;

Here, s.sum is either the direct or the integral rectangle sum, chosen once per call. cargo test -p stereo-wasm’s direct_and_integral_agree_on_every_rectangle checks they return bit-identical answers (everything here is an exact integer sum, so “identical” is a real claim, not an approximation converging).

The numbers, measured by the benchmark button below on the bundled sample image (300×220): at the widget’s default scale the direct pass runs at roughly the operation count 2·(3·(2l+1)·l) + 4·l² per pixel predicts, and the integral pass at roughly 10·4 + 9, independent of l. At the crate’s largest usable scale, $l=49$ , that is on the order of tens of thousands of operations per pixel for the direct pass against a few dozen for the integral one. Run it yourself; the exact multiplier depends on your machine, and reporting someone else’s number here would be exactly the kind of unverified claim this series tries not to make.

Non-maximum suppression: the honestly ugly part

Candidates above threshold across the eight inner scales are pooled, border pixels within 8px of the edge dropped, and then thinned. SURF.py:96-123:

len_before = len(all_pts) + 1
while len(all_pts) < len_before:
    i = 0
    len_before = len(all_pts)
    while i < len(all_pts):
        pt = all_pts[i]
        neighbors_i = np.where(np.logical_and(
            np.logical_and(all_pts[:, 0] <= pt[0] + non_max_d, all_pts[:, 0] >= pt[0] - non_max_d),
            np.logical_and(all_pts[:, 1] <= pt[1] + non_max_d, all_pts[:, 1] >= pt[1] - non_max_d)
        ))[0]
        best_i = np.argmax(all_scores[neighbors_i])
        keep_pt, refinedPt, refinedL = refine_pt(...)
        if keep_pt:
            all_pts[i] = refinedPt
            ...
            all_pts = np.delete(all_pts, items_to_delete, axis=0)
            i += 1
        else:
            all_pts = np.delete(all_pts, neighbors_i, axis=0)

An outer loop restarts from the top of the (shrinking) array every time anything changes, repeating work on points that already settled, until nothing moves. It is $O(n^2)$ in the number of candidates and it re-scans a lot of already-decided pixels. The Rust port keeps the same behaviour (best score in a box wins, gets refined, everything else in the box is dropped, and if refinement rejects the winner the whole box is dropped too), as one sort by score followed by a single greedy suppression pass, with no restart. surf.rs’s doc comment says exactly this rather than pretending the loop was never ugly.

A raw detection sits on whichever integer pixel and discrete scale happened to score highest. Brown & Lowe’s method (ref [16] of the bibliography) fits a 3×3 Hessian of $\det(H)$ in $(x, y, \text{scale})$ around that point and solves for the offset $\zeta$ that would zero its gradient:

\zeta = -H^{-1}\nabla,\qquad H=\begin{pmatrix}H_{xx}&H_{xy}&H_{xL}\\H_{xy}&H_{yy}&H_{yL}\\H_{xL}&H_{yL}&H_{LL}\end{pmatrix}

SURF.py:9-53 builds this from finite differences across the neighbouring pixels and neighbouring scales, rejects the point if the offset is larger than the sampling step (it overshot), and otherwise nudges the keypoint onto the sub-pixel location and sub-integer scale that the quadratic approximation predicts. This is the step that turns “the pixel where the loop happened to land” into “the peak of the surface underneath it.” Toggle it off in the widget below and watch keypoints snap back onto the detection grid.

Caveat

SURF.py:42-43 has a visible typo. The dy term of the gradient reads:

dy = s2 * (all_detH[j, pt[0], pt[1] + p] -
           all_detH[j, pt[0], pt[1]] - p)

It reuses the centre pixel (pt[1], offset zero) where the pattern set by every other term (and by dx two lines above) calls for the neighbour at pt[1] - p, then subtracts the raw integer step p (1, 2, 4 or 8) on top. I ported the corrected central difference, s2 * (above - below), and wrote a test rather than assume how much the typo would have mattered: the_refine_pt_typo_is_not_a_small_perturbation computes both at a real keypoint on the bundled synthetic blob. The finding surprised me: it is not a small, rounding-scale error. Because the wrong pixel is substituted rather than the right one merely mis-scaled, the buggy term behaves like a piece of the curvature itself, not like noise, and it is largest exactly where the correct term is smallest, near a real extremum, which is where keypoints actually sit. At the test’s keypoint the corrected dy is under 50 (near-flat, as expected at a symmetric peak); the buggy one differs from it by over 500. Ported faithfully, this typo would bias refinement’s $y$ -offset in a way that could plausibly flip the accept/reject check on a genuine detection.

The 3×3 solve needs no linear-algebra crate: a closed-form inverse is nine lines, surf.rs’s invert3.

The descriptor: a 20σ box, 4×4 grid, 5×5 samples

“Once the features with their associated scales are detected, the descriptor is formed by placing a square region over the feature of size 20s. This is then split up into 4×4 sub-regions. The first order derivatives $d_x$ and $d_y$ are computed in 5×5 locations across the sub-region. These are then weighted with a Gaussian function with $\sigma = 3.3s$ , which is centred at the key-point. The values are then summed per region. The absolute value of the derivatives are also summed.”

\mathbf{v} = \Big(\textstyle\sum_{i\in\mathcal R} w_i d_{x,i} \quad \sum_{i\in\mathcal R} w_i d_{y,i} \quad \sum_{i\in\mathcal R} |w_i d_{x,i}| \quad \sum_{i\in\mathcal R} |w_i d_{y,i}|\Big)

Why keep both the signed sum and the absolute one? A signed sum alone cannot tell a smooth gradient from a fine alternating stripe pattern: both can sum to zero. The absolute sum can’t, because it never cancels. Four numbers per sub-region carries that distinction the whole way to the final vector, at the cost of doubling it.

The reshape bug, and why it’s worth its own section

Two implementations of the descriptor ship in this repository, and test_all_reconstruction.py runs both. Modules/SURF_DESCRIPTOR.py (v1) slices the 20σ×20σ patch and calls:

Gx_r = Gx[pt[0] - 10 * sigma_i: pt[0] + 10 * sigma_i,
          pt[1] - 10 * sigma_i: pt[1] + 10 * sigma_i].reshape(4, 4, -1)

For an $N\times N$ patch, .reshape(4, 4, -1) walks the array row-major: it flattens the whole patch into one long sequence, top row first, then the next, and slices that sequence into 16 equal chunks. The result is not sixteen spatial squares. surf.rs’s v1_group_of computes exactly which flattened chunk each pixel lands in, and the crate’s tests establish two things about it directly, not from description:

At a patch size not divisible by 16 (n=100, from σ=5), some rows genuinely straddle two different chunks: a row’s first half is one “sub-region,” its second half another. This is the interleaving the bug is usually described by.
At a patch size that is divisible by 16 (n=80, from σ=4), no row straddles a chunk boundary, but every pixel of a row lands in the same chunk regardless of its column. The “4×4 grid” degenerates into sixteen horizontal strips; column position never affects which sub-region a pixel belongs to, at all.

Neither case is a spatial 4×4 grid. Modules/SURF_DESCRIPTOR_2.py (v2) fixes it with explicit sampling grids built once, at import time:

ind = np.indices((5, 5))
indx = np.stack((ind[0]-10, ind[0]-10, ind[0]-10, ind[0]-10,
                 ind[0]-5,  ind[0]-5,  ind[0]-5,  ind[0]-5,
                 ind[0],    ind[0],    ind[0],    ind[0],
                 ind[0]+5,  ind[0]+5,  ind[0]+5,  ind[0]+5), axis=2).reshape(-1, 4, 4)

This, worked through, places a genuine $5\times5$ sample grid in each of the sixteen $5s\times5s$ sub-regions of the $20s$ box, offset from the keypoint by exact multiples of $5s$ in each axis, then weights every sample by a Gaussian centred on the keypoint (variance $10.89=3.3^2$ , matching $\sigma=3.3s$ ). v2 also hoists the gradient computation out of the per-keypoint loop entirely: v1 calls cv2.GaussianBlur on the whole image, at $\sigma_{\text{blur}}=2\sigma_i$ , inside the loop, once per keypoint; v2 computes one plain Sobel gradient for the whole image up front and never blurs.

The 64-D vector is finally normalised per sub-region: each of the sixteen 4-numbers cells divided by its own L2 norm, independently, + 1e-12 guarding a keypoint that landed on a perfectly flat patch:

holdD /= np.linalg.norm(holdD, axis=0)[None] + 1e-12
descriptors[i] = holdD.flat

Matching, and what the report found

Matching is a ratio test along the epipolar line, the same scanline restriction FERN uses, for the same reason:

if (match_distance_1 <= thresh * match_distance_2):
    was_matched.append(True)

thresh = 0.95 (SURF_DESCRIPTOR_2.py:33): accept the best match only if it beats the second-best by at least 5%, a much stricter bar than FERN’s classification-by-argmax. Table 6 is the result on Middlebury, bad-pixel percentages at 1/2/3/5 px tolerance against ground truth:

Dataset	Image	ERR1 (%)	ERR2 (%)	ERR3 (%)	ERR5 (%)
cones	2	31.7	21.8	15.9	10.5
cones	6	85.0	81.7	80.1	77.5
teddy	2	27.4	15.9	13.1	7.8
teddy	6	94.0	76.1	68.2	65.5
tsukuba	row 3 col 3	35.2	26.1	23.8	16.9

Two photographs of a bookshelf, lamp and mannequin head side by side, with lines connecting matched SURF keypoints between them, coloured green for correct and red for incorrect against ground truth — Figure 17(b): SURF’s own matches on the Tsukuba pair, green correct / red incorrect against ground truth. Compare FERN’s matches on the same pair: visibly fewer lines, and a higher fraction of them green.

The report’s §5.4 discussion, on why SURF trades count for accuracy:

“The SURF detector and descriptor achieved lower error percentages than the FERN descriptor. However, it was not able to detect as many features as FERN. This is probably due to its higher threshold. Furthermore, SURF required a large region in order to compute its descriptor. This often lead features close to the edge of the image to be removed as their features could not be extracted.”

That “large region” is the $20\sigma$ box, and get_surf_descriptors’s own border check makes the cost concrete: a keypoint at scale $l=49$ ( $\sigma=20$ ) needs a clear $400\times 400$ pixel neighbourhood, on images that are $384\times288$ to begin with:

could_not_process = ((pt[0] - 10 * sigma_i < 0) or (pt[1] - 10 * sigma_i < 0) or
                     (pt[0] + 10 * sigma_i >= img_grey.shape[0]) or
                     (pt[1] + 10 * sigma_i >= img_grey.shape[1]))

Try it: threshold a response surface, then watch two descriptors disagree

A synthetic test image with blobs and corners at several scales — InteractiveSURF detector and descriptor explorer

A few things worth doing with it.

Push the threshold down until the heatmap and the keypoints agree. The heatmap panel shows the raw det(H) surface for whichever scale is selected: a landscape of blob-shaped peaks. The threshold slider is a plane cutting through it; every peak above the plane is a candidate. Watch keypoints appear and vanish exactly where the plane crosses a peak, not approximately.

Toggle refinement off, then on, at a low NMS radius. Off, keypoints sit on whichever integer grid position scored highest in their neighbourhood, visibly blocky if you zoom the canvas. On, they snap onto the true local peak the quadratic fit predicts. This is Brown & Lowe’s correction made visible rather than asserted.

Run the benchmark at the largest scale on offer, then the smallest. The direct pass’s cost grows with $l^2$ -ish work per pixel; the integral pass barely moves. That growing gap, not a single number, is the report’s argument.

Switch to the Descriptor tab, click a keypoint, and flip v1/v2. No caption needed: v1’s overlay paints the patch by which reshape group each pixel actually fell into, and those groups visibly are not the four-by-four squares v2 draws next to them. The 64-D bar charts underneath usually disagree most in exactly the cells where the visual groups look most scrambled.

Where the Rust differs from the NumPy

No border reflection. cv2.filter2D(..., borderType=cv2.BORDER_REFLECT) produces a response everywhere by padding the image outward; surf.rs simply doesn’t evaluate a scale within its own margin of the border ((3l-1)/2 pixels, 73px at $l=49$ , far more than the code’s fixed 8px extra check). A large-scale keypoint near the edge of a small image is something this port cannot find at all, where the original might, imperfectly, via padding.
Non-maximum suppression keeps the original’s box-suppress-refine-or-drop behaviour but as one sorted, single-pass sweep rather than the restarting $O(n^2)$ loop. See above.
refine_pt’s dy typo is fixed, not ported. See the callout above for how much that measurably matters.
v1’s Gaussian blur is bounded, not unbounded. The original re-blurs the entire image per keypoint at $\sigma_{\text{blur}}=2\sigma_i$ ; this port blurs only a window around the patch, with the kernel radius capped at 80px, so the widget stays responsive. That’s a performance simplification only: the reshape bug itself, which is what this section is about, is ported exactly.
Detection always uses the integral-image path, deliberately, so the widget’s sliders stay live; only the benchmark button ever runs the direct path over a full response map.

What’s next

FERN covers this repository’s other matcher: the one that finds more points here at a worse accuracy, for the opposite reason SURF finds fewer at a better one. Both feed the same z = bf/d triangulation from the post before this one.