SIFT Architecture & Implementation Details

For researchers and advanced users aiming to master computer vision, understanding how SIFT works under the hood is crucial. This page meticulously breaks down the mathematical architecture of the SIFT algorithm and maps it directly to the NumPy operations implemented in the pysift.py codebase.

SIFT is a multi-stage, mathematically dense pipeline. PythonSIFT organizes these stages logically to mirror the steps outlined in David Lowe's 2004 paper.

Stage 1: Scale-Space Extrema Detection

To find features that are invariant to the distance of the camera (scale invariance), SIFT searches for stable features across a continuous "scale-space". It simulates looking at the image from increasingly further distances.

  1. Image Pyramid: The input image is repeatedly blurred using a Gaussian kernel and then downsampled (halved in resolution) to form an "octave". This process is repeated across multiple octaves.
  2. Difference-of-Gaussians (DoG): Finding sharp edges and corners across multiple scales normally requires calculating the Laplacian of Gaussian (LoG), which is computationally devastating. Lowe proved that subtracting adjacent blurred images in an octave creates a highly accurate approximation of the LoG.

Code Mapping in pysift.py:

  • generateGaussianKernels(): Calculates the exact delta standard deviation required to move from one scale to the next. The math sqrt(sigma_total**2 - sigma_previous**2) ensures the convolutions accumulate correctly.
  • generateGaussianImages(): Loops through the octaves, applying cv2.GaussianBlur, and uses cv2.resize(..., interpolation=INTER_NEAREST) to step down the resolution for the next octave.
  • generateDoGImages(): Executes the subtraction (cv2.subtract) to create the 3D volume of responses.
  • findScaleSpaceExtrema(): Scans the DoG pyramids. For every pixel, it compares its intensity against its 8 neighbors in the current scale, 9 neighbors in the scale above, and 9 in the scale below (26 neighbors total). Only strict local maxima or minima proceed.

Stage 2: Sub-Pixel Keypoint Localization

Once candidate extrema are found on the discrete pixel grid, they must be refined. Pixel-level accuracy is insufficient for stable tracking. Furthermore, the DoG operator naturally fires on long, straight edges, which are terrible for tracking because they suffer from the "aperture problem" (sliding along the edge looks identical).

Code Mapping in pysift.py:

  • Taylor Expansion: localizeExtremumViaQuadraticFit() fits a 3D quadratic function to the local sample points to determine the interpolated, continuous location of the extremum (sub-pixel x, y and sub-scale s).
  • Solving the Linear System: It calculates the 3D Gradient vector and the 3x3 Hessian matrix (computeHessianAtCenterPixel) using central difference approximations. It calculates the offset using numpy.linalg.lstsq(hessian, gradient, rcond=None). If the offset is larger than 0.5, the extremum has shifted to a different pixel, and the algorithm iterates.
  • Edge Suppression: To eliminate flat regions and straight edges, the code evaluates the principal curvatures of the 2D Hessian matrix (just the X and Y components). It calculates the trace and determinant. If (trace**2) / det is greater than a specific threshold (derived from eigenvalue_ratio=10), the feature is deemed an unstable edge and discarded.

Stage 3: Orientation Assignment

To achieve true rotation invariance, a local dominant orientation is assigned to every surviving keypoint. All future operations will be calculated relative to this angle.

Code Mapping in pysift.py:

  • computeKeypointsWithOrientations() computes the gradient magnitude and direction for all pixels within a circular neighborhood around the keypoint. The size of this neighborhood scales with the keypoint's detected size.
  • Histogram Creation: It builds a 36-bin array, representing 360 degrees (10 degrees per bin). Each pixel votes for its orientation bin. The weight of the vote is determined by the pixel's gradient magnitude, multiplied by a Gaussian weight (closer pixels matter more).
  • Peak Interpolation: The highest peak in the histogram becomes the primary angle. A parabola is fitted to the 3 bins closest to the peak to find the precise, sub-degree angle.
  • Duplicate Spawning: If a secondary peak exists that is at least 80% as high as the primary peak (peak_ratio=0.8), a duplicate keypoint is spawned at the exact same location but with the secondary orientation. This drastically improves stability in complex textures.

Stage 4: Descriptor Generation

Finally, SIFT computes a highly distinctive signature (descriptor) for the local image region. This vector must remain stable even if the image is subjected to affine distortion, noise, or illumination changes.

Code Mapping in pysift.py:

  • generateDescriptors() establishes a 16x16 pixel window around the keypoint. Crucially, the axes of this window are rotated to align with the keypoint's calculated angle.
  • Spatial Bins: The 16x16 window is divided into sixteen 4x4 sub-regions. Within each sub-region, an 8-bin orientation histogram is calculated. (16 regions * 8 bins = 128 dimensions).
  • Inverse Trilinear Interpolation: When the image slightly rotates or translates, a pixel near the boundary of a bin might "snap" into an adjacent bin, causing a sudden spike in the descriptor. To prevent this, the gradient magnitude of every pixel is distributed among its 8 adjacent bins in 3D space (row, col, orientation) using weights proportional to the distance to the bin centers (c000 through c111 in the code).
  • Illumination Normalization: The final 128-dimensional vector is normalized to unit length. This cancels out uniform illumination changes (e.g., the sun comes out). Finally, to counter non-linear illumination changes (e.g., camera glare), all values above 0.2 are capped, and the vector is normalized to unit length one last time.