SIFT Scale-Invariant Features

SIFT: Scale-Invariant Feature Transform

A scale space is a family of images $L(x,y,\sigma)$ obtained by convolving the original image $I(x,y)$ with Gaussian kernels of increasing standard deviation $\sigma$:

Increasing $\sigma$ progressively suppresses fine detail, simulating what the image would look like if the same structures were captured at a coarser resolution or from farther away. Lindeberg's scale-space theory shows that the Gaussian is essentially the unique linear, shift-invariant kernel that does not create new structure (new extrema or zero-crossings) as $\sigma$ increases — a property called causality. This is why Gaussian blur, rather than an arbitrary low-pass filter, is the right tool for building a scale axis. Blob-like structures in an image are well detected by the Laplacian of Gaussian (LoG), $\sigma^2 \nabla^2 G$, whose scale-normalized extrema over $(x,y,\sigma)$ mark blobs together with the scale at which they are most salient — this is the theoretical basis of scale selection. Computing the LoG at many scales is expensive (it needs a second-derivative convolution at every scale), so SIFT uses the Difference-of-Gaussians,

which can be shown via the heat-diffusion equation to approximate $(k-1)\sigma^2\nabla^2 G$ up to a constant scale factor. Since $k$ is fixed, this constant does not change the location of extrema, so DoG extrema approximate scale-normalized LoG extrema at a fraction of the cost — DoG images fall out for free as byproducts of building the Gaussian pyramid, using only subtraction. A keypoint in SIFT is therefore a location $(x,y)$ and scale $\sigma$ at which $D$ is a local extremum with respect to all three coordinates. This joint space-scale search is exactly what makes SIFT scale-invariant: rescaling the input image simply shifts which scale index detects the same physical feature, but does not change whether it is detected.

Step 1 — Build the Gaussian pyramid in octaves. The scale axis is divided into octaves, each spanning a doubling of $\sigma$, and within each octave into $s$ intervals so that $k = 2^{1/s}$. Within an octave, $s+3$ Gaussian-blurred images are produced so that $s+2$ DoG images result, giving $s$ usable extrema-search layers (each needs a layer above and below). At the start of the next octave the image is downsampled by 2, and blurring resumes from an image already blurred to $2\sigma_0$ — this reuses computation and keeps the absolute blur cost per octave roughly constant, since a smaller image needs proportionally less work to reach the same relative blur.

Step 2 — Difference-of-Gaussians. Subtract adjacent Gaussian-blurred images in each octave to get the DoG stack $D(x,y,\sigma)$.

Step 3 — Scale-space extrema detection. For every pixel in every interior DoG layer, compare its value against its 8 spatial neighbors in the same layer, and 9 neighbors each in the layer above and below (26 total). A pixel is a candidate keypoint if it is strictly greater or strictly less than all 26.

Step 4 — Sub-pixel/sub-scale localization. Fit a 3-D quadratic (second-order Taylor expansion of $D$ around the candidate) in $(x,y,\sigma)$ and solve for the offset that maximizes/minimizes it, giving continuous-valued location and scale rather than a discrete grid point. Reject the point if the offset is large (re-center and retry, up to a small number of iterations) or if $|D|$ at the interpolated extremum is below a contrast threshold.

Step 5 — Eliminate edge responses. DoG, like the Laplacian, responds strongly along edges, where localization along the edge direction is poorly conditioned. Compute the $2\times2$ Hessian of $D$ in $x,y$ and reject points where the ratio of principal curvatures exceeds a threshold (equivalent to a Harris-like corner-vs-edge test).

Step 6 — Assign a dominant orientation. Around each surviving keypoint, at its own characteristic scale, build a histogram of gradient orientations (weighted by gradient magnitude and a Gaussian window). The histogram's peak (and any secondary peak within 80% of the maximum) defines the keypoint's canonical orientation, making the eventual descriptor rotation-invariant by always being computed relative to this direction. ```text function SIFT_DETECT(image I, octaves O, intervals s, sigma0, contrast_thresh, edge_thresh_r):

# ---- 1. Build Gaussian scale-space pyramid ---- k <- 2^(1/s) for o in 0 .. O-1: base <- (o == 0) ? I : DOWNSAMPLE(G[o-1][s], factor=2) G[o][0] <- GAUSSIAN_BLUR(base, sigma0) for i in 1 .. s+2: sigma_i <- sigma0 * k^i G[o][i] <- GAUSSIAN_BLUR(G[o][i-1], incremental_sigma(sigma_i))

# incremental blur: convolve previous image with the

# additional sigma needed, using sigma_new^2 = sigma_i^2 - sigma_{i-1}^2

# ---- 2. Difference-of-Gaussians ---- for o in 0 .. O-1: for i in 0 .. s+1: D[o][i] <- G[o][i+1] - G[o][i] keypoints <- []

# ---- 3. Scale-space extrema search ---- for o in 0 .. O-1: for i in 1 .. s: # skip top/bottom DoG layer (need neighbors) for each pixel (x, y) interior to D[o][i]: v <- D[o][i][x,y] neighbors <- 26 samples from D[o][i-1], D[o][i], D[o][i+1] in the 3x3x3 cube minus the center if v is strict max or strict min over neighbors: candidate <- (o, i, x, y)

# ---- 4. Sub-pixel / sub-scale localization ---- (offset, D_hat) <- QUADRATIC_REFINE(D, o, i, x, y) if |offset| too large in any dimension: retry at shifted integer position (bounded # of tries) continue if |D_hat| < contrast_thresh: discard candidate continue

# ---- 5. Reject edge-like responses ---- H <- HESSIAN_2x2(D[o][i], x, y) # only x,y components Tr, Det <- trace(H), det(H) if Det <= 0: discard candidate # curvatures have opposite sign: not an extremum continue r <- edge_thresh_r if (Tr*Tr) / Det >= ((r+1)^2)/r: discard candidate # edge-like, poorly localized continue

# ---- 6. Orientation assignment ---- scale_kp <- effective_sigma(o, i, offset) hist <- ORIENTATION_HISTOGRAM(G[o][i], x, y, scale_kp) for each peak within 80% of hist max: keypoints.append({x, y, scale_kp, orientation=peak, octave=o}) return keypoints ```

Time complexity. Building each Gaussian-blurred image is a separable convolution, costing $O(WH \cdot r)$ for an $r$-tap 1-D kernel applied in $x$ then $y$, i.e. $O(WHr)$ per layer rather than $O(WHr^2)$ for a naive 2-D convolution. Since $O(\log)$ octaves halve resolution each time, the total pixel count across all octaves is $WH(1 + 1/4 + 1/16 + \dots) = O(WH)$, a geometric series bounded by $\tfrac{4}{3}WH$. With $s+3$ blur levels per octave the whole pyramid costs $O(WH \cdot s \cdot r)$, and DoG subtraction is $O(WH \cdot s)$ — dominated by the blurring. Extrema search touches every interior DoG pixel a constant number of times (26 comparisons), so it is $O(WHs)$. Sub-pixel refinement is a constant number of $3\times3$ linear solves per candidate, and orientation histograms cost $O(k^2)$ per keypoint for a $k\times k$ window — negligible compared to pyramid construction since the number of surviving keypoints is far smaller than $WH$. Overall detection is linear in image size, $O(WH)$, with a modest constant factor from the number of scales sampled per octave — a critical practical property, since Harris-style multi-scale search done naively (recomputing a fixed-size operator at many independently chosen scales without pyramid reuse) can cost more per level.

Space complexity is $O(WHs)$ for storing the Gaussian and DoG stacks of one octave (subsequent octaves are geometrically smaller so total storage is again bounded by a constant multiple of a single full-resolution image), plus $O(K)$ for the $K$ detected keypoints.

Correctness argument. Three separate claims must hold. Claim 1 — DoG approximates scale-normalized LoG. The Gaussian satisfies the heat/diffusion equation $\partial L/\partial \sigma = \sigma \nabla^2 L$. Approximating the derivative by a finite difference between adjacent scales,

gives $D(x,y,\sigma) = L(x,y,k\sigma) - L(x,y,\sigma) \approx (k-1)\,\sigma^2 \nabla^2 L$. Since $k$ is a fixed constant chosen at pyramid-construction time, $D$ and the scale-normalized LoG differ only by a positive multiplicative constant that does not depend on $(x,y,\sigma)$ — hence they share the same extrema locations. Claim 2 — searching 26 neighbors correctly finds scale-space extrema. A point is a genuine local extremum of the sampled 3-D function $D$ iff it is more extreme than every point in its discrete neighborhood; on a regular 3-D grid with one sample per octave-fraction, the natural discrete neighborhood is exactly the 26 grid-adjacent points (a 3×3×3 cube minus the center). This is the direct 3-D analogue of the 1-D condition "$f(x)$ is a local max iff $f(x) > f(x-1)$ and $f(x) > f(x+1)$," so no weaker test would correctly capture extrema jointly in space and scale. Claim 3 — the quadratic (Taylor) refinement step is a locally valid Newton step. Near a true extremum, $D$ is well approximated by its second-order Taylor expansion

whose stationary point is $\hat{\mathbf{x}} = -H^{-1}\nabla D$, the standard Newton step; since $D$ is locally quadratic near a genuine extremum this converges in one or two iterations, and iterating with re-centering when $\hat{\mathbf{x}}$ leaves the sampled cell simply re-applies Newton's method around the nearest grid point, which is guaranteed to terminate because the search domain is finite. Rejecting edges is sound because the ratio test $\mathrm{Tr}(H)^2/\det(H) \geq (r+1)^2/r$ is exactly the Harris corner-vs-edge criterion applied to the Hessian of $D$: if the two eigenvalues (principal curvatures) $\alpha,\beta$ satisfy $\alpha = r\beta$, then $\mathrm{Tr}^2/\det = (r+1)^2/r$, so the threshold directly bounds the eigenvalue ratio without needing to compute eigenvalues explicitly — the same trick used in Harris Corner Detection to avoid an eigendecomposition.

Tradeoffs. SIFT trades computation (many blurred images per octave, floating-point convolutions) for invariance and repeatability; it is far more expensive than a single-scale Harris or FAST detector, and its scale invariance is only approximate since DoG is a discretely-sampled approximation to a continuous scale axis with a finite number of intervals $s$ per octave (Lowe's original experiments found $s=3$ near-optimal for repeatability versus cost). It is also not invariant to full affine or perspective distortion — only similarity transforms (scale + rotation + translation) are handled exactly, with affine-covariant extensions such as MSER or Hessian-Affine addressing the rest.

Consider a $512 \times 512$ grayscale photo of a checkerboard pattern with squares at two different physical densities: a coarse pattern near the camera (squares ≈ 40 px) and a fine pattern receding into the distance (squares ≈ 8 px), plus one small circular logo of radius ≈ 6 px pasted in a corner. Using $\sigma_0 = 1.6$, $s = 3$ intervals per octave (so $k = 2^{1/3} \approx 1.26$), and $O = 4$ octaves:

Trace a single candidate: at octave 1, layer $i=2$ ($\sigma \approx 5.0$), pixel $(x,y) = (130, 74)$ has DoG value $D = -0.087$ (a minimum, since it is a dark checkerboard corner against a light background — dark blobs give negative DoG extrema). All 8 same-layer neighbors are $> -0.087$; all 9 neighbors in layer $i=1$ and all 9 in layer $i=3$ are also $> -0.087$. It passes the 26-neighbor test. Quadratic refinement solves the $3\times3$ linear system for $\hat{\mathbf{x}} = (\Delta x, \Delta y, \Delta\sigma\text{-index})$, yielding $\hat{\mathbf{x}} = (0.18, -0.06, 0.31)$ — small, so no re-centering iteration is needed. The interpolated extremum value is $\hat{D} = -0.091$, comfortably past the contrast threshold (commonly 0.03 on an intensity range $[0,1]$), so the point survives. The Hessian of $D$ at this pixel gives $\mathrm{Tr}(H) = -0.31$, $\det(H) = 0.021$, so $\mathrm{Tr}^2/\det = 4.58$. With edge threshold $r = 10$, the bound $(r+1)^2/r = 12.1$ is not exceeded, so the point is corner-like enough (this is a checkerboard corner, as expected) and is kept. Finally, an orientation histogram is built in a Gaussian-weighted window of radius $\approx 3\sigma_{kp}$ around the refined location, using gradients of the octave-1 image at the appropriate blur level. Suppose the histogram (36 bins of 10°) peaks sharply at bin 14 (140°) with no secondary peak above 80% of the max — the keypoint is assigned orientation $140°$, one keypoint emitted: $(x=130.18,\, y=73.94,\, \sigma \approx 5.6,\, \theta = 140°,\, \text{octave}=1)$. This single trace illustrates the key invariance property: if the same photo were retaken twice as close (doubling all feature sizes), the same checkerboard corner would instead be detected in octave 2 at a proportionally doubled $\sigma$, with the same orientation $140°$ (orientation is invariant to scale, only to rotation) and an essentially identical descriptor once built — which is precisely why matching across zoom levels works.

SIFT-style scale-space keypoint detection underlies classical panorama stitching software, where images taken at different zooms and angles must be aligned by matching stable keypoints before estimating a homography. Structure-from-motion and early visual SLAM systems used SIFT (or its faster relatives, SURF and ORB) to establish correspondences across frames captured at varying distances, feeding into bundle adjustment for 3-D reconstruction. Content-based image retrieval and "bag of visual words" systems built large indexes of SIFT descriptors so that a query photo could be matched against a million-image database purely from local scale-invariant features. Robotics and augmented-reality tracking systems used SIFT keypoints to relocalize a camera against a known map. Even after deep-learned descriptors and detectors (e.g., learned keypoint networks, SuperPoint-style detectors) became competitive, the scale-space extremum idea remains the conceptual backbone that motivates learned scale-covariant detectors, and SIFT is still shipped as a baseline in OpenCV and used in benchmarking new Feature Descriptors and Matching methods.

Ask the tutor: How would SIFT's detection stage need to change to handle a fisheye or wide-angle lens where scale itself varies non-uniformly across the image?