SIFT

SIFT, the Scale-Invariant Feature Transform, is a method for finding distinctive points in an image and describing them so that they can be matched across images taken at different scales, orientations, and viewpoints. Introduced by David Lowe in 1999 and described in full in his 2004 journal paper, it combines a scale-space keypoint detector with a gradient-histogram descriptor. SIFT made wide-baseline feature matching practical and became the reference local feature for object recognition, panorama stitching, and structure from motion.

Problem

Corner detectors such as the Harris detector work at one fixed scale, which Lowe noted makes them a poor basis for matching images of different sizes. SIFT addresses two problems together: detecting points whose location and characteristic size can be recovered in any image of the scene, and describing their neighborhoods in a way that is unaffected by scale and rotation and robust to moderate changes in viewpoint and illumination.

Inputs and Outputs

The input is a grayscale image. The output is a set of keypoints, each with a subpixel position $(x, y)$, a scale $\sigma$, an orientation $\theta$, and a 128-dimensional descriptor vector of unit length. Lowe reported that a typical $500 \times 500$ image gives about 2,000 stable keypoints, depending on content and parameters.

Intuition

A blob, such as a dark window on a bright wall, has a natural size. A center-surround filter responds most strongly when its size matches the blob’s, so searching for the strongest responses over position and blur finds the blob and its size. A dominant gradient direction then gives the point an orientation. Measuring the neighborhood in a frame centered at the point, scaled by its size, and rotated to its orientation makes the measurements independent of all three, and a grid of small orientation histograms summarizes it coarsely enough to tolerate a few pixels of misalignment.

Algorithm

The steps form a cascade: the expensive ones run only where the cheap ones pass.

1. Scale-space construction. Double the input image with linear interpolation, blur it to an initial scale of $\sigma = 1.6$, and build a Gaussian scale space organized in octaves. Each octave contains $s + 3$ images blurred by successive factors $k = 2^{1/s}$; adjacent images are subtracted to give $s + 2$ difference-of-Gaussian (DoG) images. The next octave starts from the Gaussian image with twice the initial blur, subsampled by taking every second pixel.
2. Extrema detection. Compare each DoG sample with its 8 neighbors at the same scale and the 9 neighbors in each of the scales above and below. Keep it if it is larger or smaller than all 26.
3. Keypoint localization. Fit a quadratic to the DoG values around each candidate to estimate the extremum’s subpixel position and scale. If the offset exceeds half a sample in any dimension, move to the neighboring sample and repeat.
4. Low-contrast rejection. Discard extrema whose interpolated DoG value satisfies $|D(\hat{\mathbf{x}})| < 0.03$, for pixel values in $[0, 1]$.
5. Edge rejection. Discard extrema whose ratio of principal curvatures exceeds $r = 10$, measured with the $2 \times 2$ Hessian of the DoG.
6. Orientation assignment. In the Gaussian image closest to the keypoint’s scale, build a 36-bin histogram of gradient orientations (10° per bin) around the keypoint, weighted by gradient magnitude and a Gaussian window with $\sigma$ equal to 1.5 times the keypoint scale. Assign the highest peak as the orientation, interpolated with a parabola through the three nearest bins. Every other local peak within 80% of the highest creates an additional keypoint at the same location and scale with that orientation; Lowe reported that only about 15% of points receive more than one.
7. Descriptor. Rotate the sampling grid and gradient orientations to the keypoint orientation. Take a $16 \times 16$ array of gradient samples, weight them by a Gaussian with $\sigma$ equal to half the window width, and accumulate them into a $4 \times 4$ grid of 8-bin orientation histograms, distributing each sample with trilinear interpolation over neighboring cells and orientation bins. Concatenate the histograms into a 128-dimensional vector, normalize it to unit length, clip every entry at 0.2, and normalize again. Normalization removes affine illumination changes; the clipping, with a value Lowe chose experimentally, limits large gradient magnitudes caused by non-linear effects such as saturation.

Mathematical Formulation

Scale space and the difference of Gaussians

The scale space of an image $I$ is the family of blurred images

$$L(x, y, \sigma) = G(x, y, \sigma) * I(x, y), \qquad G(x, y, \sigma) = \frac{1}{2\pi\sigma^2} e^{-(x^2 + y^2)/2\sigma^2},$$

and the DoG function is the difference of two nearby scales,

$$D(x, y, \sigma) = \big(G(x, y, k\sigma) - G(x, y, \sigma)\big) * I(x, y) = L(x, y, k\sigma) - L(x, y, \sigma).$$

From the diffusion equation $\partial G / \partial \sigma = \sigma \nabla^2 G$, Lowe showed that

$$G(x, y, k\sigma) - G(x, y, \sigma) \approx (k - 1)\, \sigma^2 \nabla^2 G ,$$

so the DoG approximates the scale-normalized Laplacian of Gaussian studied by Lindeberg, with the factor $\sigma^2$ built in (see the Laplacian operator); the constant $k - 1$ does not move the extrema, and the blurred images are needed anyway.

For a Gaussian blob of standard deviation $\sigma_b$, the magnitude of the DoG response at the blob’s center is proportional to $1/(\sigma_b^2 + \sigma^2) - 1/(\sigma_b^2 + k^2\sigma^2)$, which is extremal at $\sigma = \sigma_b / \sqrt{k}$. The geometric mean of the two blur levels, $\sigma\sqrt{k}$, then equals $\sigma_b$: the detected scale measures the blob’s size, which the example below confirms.

Subpixel localization

With $\mathbf{x} = (x, y, \sigma)^\top$ the offset from a sample point, the second-order Taylor expansion of $D$ is

$$D(\mathbf{x}) = D + \frac{\partial D}{\partial \mathbf{x}}^{\top} \mathbf{x} + \frac{1}{2} \mathbf{x}^{\top} \frac{\partial^2 D}{\partial \mathbf{x}^2} \mathbf{x},$$

with the derivatives approximated by differences of neighboring samples. Setting its derivative to zero gives the extremum and its value,

$$\hat{\mathbf{x}} = -\left(\frac{\partial^2 D}{\partial \mathbf{x}^2}\right)^{-1} \frac{\partial D}{\partial \mathbf{x}}, \qquad D(\hat{\mathbf{x}}) = D + \frac{1}{2} \frac{\partial D}{\partial \mathbf{x}}^{\top} \hat{\mathbf{x}} .$$

Edge rejection

The DoG also responds along edges, where position along the edge is poorly determined. With the spatial Hessian

$$\mathbf{H} = \begin{bmatrix} D_{xx} & D_{xy} \\ D_{xy} & D_{yy} \end{bmatrix}$$

and eigenvalues $\alpha = r\beta$, the ratio $\operatorname{Tr}(\mathbf{H})^2 / \operatorname{Det}(\mathbf{H}) = (r + 1)^2 / r$ depends only on $r$ and grows with it. A keypoint is kept only if

$$\frac{\operatorname{Tr}(\mathbf{H})^2}{\operatorname{Det}(\mathbf{H})} < \frac{(r + 1)^2}{r},$$

and discarded if the determinant is negative. Like the Harris corner measure, from which Lowe borrowed the idea, the test avoids computing eigenvalues.

Matching and Recognition

Lowe matched each descriptor to its Euclidean nearest neighbor and kept the match only if it was closer than 0.8 times the second-nearest, the ratio test described in feature matching. Because each keypoint carries a location, scale, and orientation, a single match predicts an object’s pose. For recognition, matches vote for pose in a coarse Hough transform, which Lowe preferred because the correct matches for one object can be under 1% of all matches, far too few for RANSAC; clusters of at least three votes are then verified by a least-squares affine fit and a probabilistic acceptance test. Modern pipelines usually replace these steps with RANSAC on a homography or epipolar geometry.

Parameters

Parameter	Lowe (2004)	Effect of the choice
Intervals per octave $s$	3	Best repeatability in Lowe’s tests; more intervals give more but less stable keypoints
Initial blur $\sigma$	1.6	Repeatability rises with $\sigma$, at growing cost
Input doubling	Yes	Almost four times as many stable keypoints
Contrast threshold	0.03	Removes weak, noise-sensitive extrema
Edge ratio $r$	10	Removes extrema elongated along edges
Descriptor grid	$4 \times 4 \times 8$	Smaller grids discriminate worse; larger ones are more sensitive to distortion
Ratio threshold	0.8	Trades false matches against missed ones

The contrast threshold depends on the intensity range and on $s$, since the DoG magnitude shrinks as the scale step $k - 1$ shrinks; an implementation that changes $s$ must adjust the threshold to match.

Complexity

Building the scale space dominates. Each octave has a quarter of the pixels of the one before, so the total is a constant factor times the first octave, which has four times the original pixels. Extremum detection is a linear scan, and the later steps cost a constant amount per keypoint. Lowe reported recognizing all objects in his example scenes in under 0.3 seconds on a 2 GHz Pentium 4.

Implementation

The example implements SIFT’s detection stages (steps 2–5) on a synthetic image with four Gaussian blobs of standard deviation 2, 4, 8, and 16 pixels, an elongated ridge, and mild noise. For clarity it uses a single-resolution scale space ($s = 3$, $\sigma = 1.6$) without image doubling or per-octave subsampling, and omits orientation and the descriptor. The reported scale is $\sigma\sqrt{k}$, which should equal the blob’s standard deviation.

import numpy as np
from scipy import ndimage

# Synthetic image in [0, 1]: Gaussian blobs of known size plus one elongated ridge.
H, W = 256, 640
y, x = np.mgrid[0:H, 0:W]
blobs = [(128, 60, 2.0), (128, 140, 4.0), (128, 260, 8.0), (128, 440, 16.0)]
image = np.zeros((H, W))
for cy, cx, s in blobs:
    image += np.exp(-((x - cx) ** 2 + (y - cy) ** 2) / (2 * s**2))
image += np.exp(-((x - 580) ** 2 / (2 * 3.0**2) + (y - 128) ** 2 / (2 * 30.0**2)))
image += np.random.default_rng(0).normal(0, 0.02, image.shape)

# Gaussian scale space, s = 3 intervals per octave starting at sigma = 1.6.
# (Single resolution for clarity; SIFT halves the image after each octave.)
s_per_octave, sigma0, levels = 3, 1.6, 16
k = 2 ** (1 / s_per_octave)
sigmas = sigma0 * k ** np.arange(levels)
L = np.stack([ndimage.gaussian_filter(image, s) for s in sigmas])
D = L[1:] - L[:-1]  # difference of Gaussians; D[i] pairs sigmas[i] and sigmas[i + 1]

# 1. Candidates: larger or smaller than all 26 neighbours in space and scale.
is_max = D == ndimage.maximum_filter(D, size=3)
is_min = D == ndimage.minimum_filter(D, size=3)
cand = np.argwhere((is_max | is_min) & (np.abs(D) > 0.5 * 0.03))  # pre-filter (shortcut)
inner = (cand > 0) & (cand < np.array(D.shape) - 1)  # skip the outermost samples
cand = cand[inner.all(axis=1)]

def derivatives(i, r, c):
    """Gradient and Hessian of D at a sample, by finite differences (s, y, x)."""
    d = D[i - 1:i + 2, r - 1:r + 2, c - 1:c + 2]
    g = 0.5 * np.array([d[2, 1, 1] - d[0, 1, 1], d[1, 2, 1] - d[1, 0, 1], d[1, 1, 2] - d[1, 1, 0]])
    Hm = np.empty((3, 3))
    Hm[0, 0] = d[2, 1, 1] - 2 * d[1, 1, 1] + d[0, 1, 1]
    Hm[1, 1] = d[1, 2, 1] - 2 * d[1, 1, 1] + d[1, 0, 1]
    Hm[2, 2] = d[1, 1, 2] - 2 * d[1, 1, 1] + d[1, 1, 0]
    Hm[0, 1] = Hm[1, 0] = 0.25 * (d[2, 2, 1] - d[2, 0, 1] - d[0, 2, 1] + d[0, 0, 1])
    Hm[0, 2] = Hm[2, 0] = 0.25 * (d[2, 1, 2] - d[2, 1, 0] - d[0, 1, 2] + d[0, 1, 0])
    Hm[1, 2] = Hm[2, 1] = 0.25 * (d[1, 2, 2] - d[1, 2, 0] - d[1, 0, 2] + d[1, 0, 0])
    return g, Hm

kept, n_unstable, n_low, n_edge = [], 0, 0, 0
for i, r, c in cand:
    # 2. Quadratic (Taylor) interpolation, moving to a neighbour if the offset exceeds 0.5.
    for _ in range(5):
        g, Hm = derivatives(i, r, c)
        offset = -np.linalg.solve(Hm, g)
        if np.all(np.abs(offset) <= 0.5):
            break
        i, r, c = np.array([i, r, c]) + np.round(offset).astype(int)
        if not (0 < i < len(D) - 1 and 0 < r < H - 1 and 0 < c < W - 1):
            break
    if np.any(np.abs(offset) > 0.5):  # did not settle on a sample: discard
        n_unstable += 1
        continue
    value = D[i, r, c] + 0.5 * g @ offset
    # 3. Contrast test |D(x_hat)| >= 0.03, then 4. edge test with r = 10.
    if abs(value) < 0.03:
        n_low += 1
        continue
    Dxx, Dyy, Dxy = Hm[2, 2], Hm[1, 1], Hm[1, 2]
    tr, det = Dxx + Dyy, Dxx * Dyy - Dxy**2
    if det <= 0 or tr**2 / det >= (10 + 1) ** 2 / 10:
        n_edge += 1
        continue
    sigma = sigma0 * k ** (i + offset[0])  # interpolated scale of the lower Gaussian
    kept.append((r + offset[1], c + offset[2], sigma * np.sqrt(k), value))

print(f"candidates {len(cand)}, unstable {n_unstable}, low contrast {n_low}, edge-like {n_edge}, kept {len(kept)}")
for row, col, scale, value in sorted(kept, key=lambda p: p[1]):
    print(f"x = {col:6.1f}  y = {row:6.1f}  scale = {scale:5.2f}  D = {value:+.3f}")

Output:

candidates 142, unstable 68, low contrast 66, edge-like 4, kept 4
x =   60.0  y =  128.0  scale =  2.03  D = -0.113
x =  140.0  y =  128.0  scale =  3.99  D = -0.115
x =  260.0  y =  128.0  scale =  8.02  D = -0.115
x =  440.0  y =  128.0  scale = 15.99  D = -0.115

Each blob yields one keypoint at its center, with a scale within 1.5% of its size. The DoG value is negative because a bright blob is a minimum of $L(k\sigma) - L(\sigma)$. Most candidates are noise extrema that fail to settle or have too little contrast, and the ridge, despite strong contrast, fails the edge test. OpenCV (cv2.SIFT_create) and VLFeat provide complete, tested implementations.

Properties and Behavior

• Scale and rotation invariance, from detection over all scales and description in the keypoint’s normalized frame.
• Partial viewpoint robustness. Lowe reported reliable recognition of textured planar surfaces up to about 50° of rotation in depth, and of 3D objects up to about 30°.
• Noise robustness. In Lowe’s tests orientation assignment stayed accurate 95% of the time even with 10% pixel noise; most matching failures came from location and scale detection.
• Distinctiveness. A single descriptor was often matched correctly against a database of 40,000 keypoints, which made recognition from as few as three matches possible.

Limitations

• Cost. The scale space and 128-dimensional floating-point descriptors are expensive compared with FAST detection and binary descriptors such as ORB, designed as faster alternatives.
• Limited affine invariance. The frame normalizes only position, scale, and rotation, so strong out-of-plane rotation distorts the patch beyond what pooling absorbs.
• Textureless and repetitive scenes. Uniform regions give few keypoints, and repeated structure produces descriptors that the ratio test rejects as ambiguous.
• Orientation dependence. An unstable dominant orientation puts the descriptor in the wrong frame; the extra keypoints at secondary peaks only partly compensate.

Patent and Availability

The University of British Columbia held U.S. patent 6,711,293 on the method (filed March 2000, granted March 2004, Lowe as inventor), which expired in March 2020. OpenCV had kept SIFT in its non-free contributed modules; release 4.4.0 (July 2020) moved it into the main repository.

Variants

• RootSIFT (Arandjelović and Zisserman, 2012) normalizes the SIFT vector to unit L1 norm and takes the element-wise square root, so that Euclidean distance between the results corresponds to the Hellinger kernel on the original histograms. The authors found it markedly improved object retrieval at negligible cost and with no extra storage.
• PCA-SIFT (Ke and Sukthankar, 2004) keeps SIFT’s detector but projects the normalized gradient patch onto a low-dimensional basis learned with principal component analysis.
• DSP-SIFT (Dong and Soatto, 2015) pools gradient orientations across several domain sizes as well as locations, keeping SIFT’s dimension without training.
• ASIFT (Morel and Yu, 2009) applies SIFT to simulated views over the two angles of the camera axis, extending invariance to all six parameters of an affine transformation at higher computational cost.
• GLOH, SURF, and binary alternatives such as ORB are covered in feature descriptors.