Building a 100% Pure-Rust SIFT & Image Stitcher from Scratch (Zero OpenCV, Zero Unsafe)

Eddie — Thu, 12 Nov 2026 09:55:25 +0000

If you have ever needed to do classic feature matching or image alignment in a production Rust service, you know the drill: you either pull in heavy C++ OpenCV wrappers, fight system-level CMake configurations in Docker, or deal with non-trivial FFI safety hazards.

David Lowe’s landmark Scale-Invariant Feature Transform (SIFT) patent (US Patent 6,711,293 B1) expired in March 2020. Yet in the Rust ecosystem, almost every serious feature matching setup still leaned on external C/C++ libraries.

So I built pure-sift: an open-source, 100% pure-Rust library (#![forbid(unsafe_code)]) implementing the full SIFT pipeline, 2D Similarity & 4-Point DLT Homography RANSAC, and an automated image stitcher with zero OpenCV or C++ build dependencies.

GitHub: https://github.com/dannkkie/pure-sift
crates.io: https://crates.io/crates/pure-sift

Here are three concrete engineering lessons learned from reading the patent claims, paper math, and implementing it from first principles in safe Rust.

1. 3D Taylor Quadratic Sub-Pixel Localization Is Where the Battle Is Won

Extrema in the Difference-of-Gaussians (DoG) scale space are detected by comparing a sample point against its 26 neighbors (8 in the current scale, 9 in the scale above, and 9 in the scale below).

Finding a discrete integer coordinate (x, y, σ) is easy. But discrete coordinates give jagged, noisy keypoint locations that destroy subsequent homography calculations.

Claim 1 of US6711293B1 specifies expanding the scale-space function D(x, y, σ) via a 3D Taylor series truncated to quadratic terms:

D(x) = D + (∂D/∂x)ᵀ x + ½ xᵀ (∂²D/∂x²) x

Taking the derivative and setting it to zero yields the exact sub-pixel offset vector:

x̂ = - (∂²D/∂x²)⁻¹ (∂D/∂x)

The Rust Implementation Details:

3x3 Hessian Inversion: Because this is evaluated for thousands of candidate points, we solve the 3×3 linear system directly using an analytical cofactor matrix rather than bringing in heavy linear algebra libraries.
Convergence Limits: If any component of x̂ exceeds 0.5, the true extremum lies closer to an adjacent pixel. We iteratively shift the candidate center (up to 5 steps). In our unit tests against synthetic paraboloids, the sub-pixel error converges to < 0.12 pixels.
Contrast & Edge Thresholding: Extrema with low contrast (|D(x̂)| < 0.03) or high principal curvature ratio (governed by the Hessian trace and determinant ratio r > 10.0) are filtered out to reject poorly localized edge responses.

2. The >= 80% Secondary Orientation Peak Rule

A keypoint with only a single orientation breaks invariance under symmetric or ambiguous textures.

In Lowe’s formulation, we create a 36-bin orientation histogram weighted by gradient magnitude and a Gaussian window (σ_w = 1.5 × σ).

The subtle detail often missed in naive implementations:

Any local histogram peak within 80% of the dominant peak height splits off into a new, independent keypoint sharing the exact same (x, y, σ) coordinates but with the secondary orientation angle.
Each peak undergoes parabolic interpolation across its adjacent bins for smooth angular resolution:

θ = θ₀ + [ (h_left - h_right) / (2 × (h_left - 2×h_center + h_right)) ] × (2π / 36)

This orientation splitting accounts for roughly 15% of the total keypoint count, but it is responsible for nearly 40% of robust matches in difficult, repeating environments.

3. Why 4-Point DLT Homography Beats 2D Similarity on Perspective Slant

Early in the project, we tested a 2-point 2D Similarity model (estimating rotation θ, uniform scale s, and translation t_x, t_y via 4 degrees of freedom).

On planar aerial imagery with straight-down camera angles, 2D Similarity performed great:

Requires only 2 point correspondences per RANSAC sample.
Evaluates in under 10 ms.

However, the moment we tested heavy perspective slant (such as large out-of-plane camera tilt angles akin to the classic Oxford affine benchmarks), 2D Similarity failed completely: inlier count dropped to zero.

To solve this, we implemented the 4-Point Direct Linear Transform (DLT) with SVD:

Sets up an 8×9 algebraic matrix A from 4 point matches:

| -x  -y  -1   0   0   0   x·x'  y·x'  x' |
|  0   0   0  -x  -y  -1   x·y'  y·y'  y' |

Solves A h = 0 via pure-Rust Jacobi SVD.
Horizon & Degeneracy Guard: Checks that collinearity is avoided, points are convex, and enforces positive determinant (det(H) > 0) to discard unphysical "turned inside-out" backplane projections.
With 4-Point DLT Homography RANSAC, matching through heavy perspective angle shifts went from complete failure to hundreds of validated geometric inliers.

Architecture & Zero-Unsafe Guarantee

To keep the crate clean for production and embedded use:

#![forbid(unsafe_code)]: Zero raw pointer arithmetic or unchecked memory buffers.
Rayon Parallelism: Keypoint descriptor generation is inherently embarrassingly parallel; Rayon distributes the 128-D orientation binning across CPU threads.
Lightweight Optional Web Workbench: A built-in Axum/Tokio inspection UI is included under the default server feature for visual debugging. When using pure-sift as an embedded or WebAssembly (WASM) library, pass --no-default-features to compile with zero async runtime overhead.

[dependencies]
pure-sift = { version = "0.1.0", default-features = false }

Try It Out

cargo install pure-sift

# Extract SIFT keypoints and match two photos with RANSAC:
pure-sift match image1.jpg image2.jpg --model homography --output matches.jpg

# Automatically stitch into a panorama:
pure-sift stitch tile_left.jpg tile_right.jpg --output panorama.jpg

Code, tests, and documentation are up on GitHub. If you test it against your own datasets or have ideas for optimizations, issues and PRs are warmly welcome!