PySIFT

The 5-Stage GPU Pipeline

PySIFT implements Lowe's SIFT algorithm (IJCV 2004) entirely on the GPU. Every stage — from image preprocessing to descriptor output — executes as CUDA kernels via CuPy and Numba. The image never leaves VRAM.

Stage 1: CLAHE Preprocessing

Adaptive histogram equalization (CLAHE) normalizes local contrast before building the scale-space pyramid. This yields 40-80% more usable keypoints on night, haze, and low-contrast images vs. raw input — critical for drone imagery at dawn/dusk and medical scans with uneven illumination.

Stage 2: Gaussian Scale-Space Pyramid

The image is progressively blurred and downsampled across octaves. PySIFT supports fp16_pyramid=True for half-precision storage (2x memory savings on VRAM-constrained GPUs like Jetson or RTX 3050), with full fp32 precision during DoG subtraction to avoid catastrophic cancellation.

Stage 3: DoG Extrema Detection

Difference-of-Gaussians (DoG) volumes are computed by subtracting adjacent Gaussian levels. Local extrema in 3D (x, y, scale) are candidate keypoints. Sub-pixel refinement via Taylor expansion and edge-response rejection (Hessian ratio test) filter the candidates.

Stage 4: Orientation Assignment

A gradient orientation histogram (36 bins) around each keypoint determines its dominant orientation, making descriptors rotation-invariant. PySIFT also supports learned orientation via OriNet (orientation='orinet').

Stage 5: Descriptor Computation

A 128-dimensional histogram of oriented gradients over a 16x16 patch (4x4 spatial bins, 8 orientation bins). PySIFT applies RootSIFT normalization by default (L1 + sqrt), which improves matching over raw L2-normalized SIFT. With dsp=True, multi-scale descriptor pooling further boosts robustness.

GPU-Residency: Why It Matters

In a typical vision pipeline (detect → match → estimate → reconstruct), CPU-based SIFT forces a GPU→CPU→GPU round-trip at every stage boundary. PySIFT eliminates this overhead entirely:

Keypoint Detection

PySIFT detects scale-invariant keypoints across a Gaussian pyramid. The default configuration (4 octaves, 3 scales per octave, contrast threshold 0.04) closely matches OpenCV SIFT's behavior while running entirely on the GPU.

PySIFT feature matching on landmark images. Green = PySIFT inliers, Red = OpenCV inliers.

Basic Usage

from pysift import PySIFT
import cv2

sift = PySIFT()
gray = cv2.imread("image.jpg", cv2.IMREAD_GRAYSCALE)

# Returns OpenCV-compatible KeyPoints + NumPy descriptors
keypoints, descriptors = sift.detectAndCompute(gray, None)
print(f"Detected {len(keypoints)} keypoints")
print(f"Descriptor shape: {descriptors.shape}")  # (N, 128)

GPU-Resident Output

# Keep everything on GPU as CuPy arrays
kp_gpu, desc_gpu = sift.detectAndCompute(gray, None, gpu_output=True)
# kp_gpu: CuPy (N, 4) -- [x, y, size, angle]
# desc_gpu: CuPy (N, 128) -- ready for GPU matching

Tunable Parameters

Feature Matching

PySIFT descriptors are L2-compatible (or RootSIFT by default), so any standard matcher works. For best results, combine Lowe's ratio test with mutual nearest-neighbor (MNN) filtering.

CPU Matching (OpenCV BFMatcher)

bf = cv2.BFMatcher(cv2.NORM_L2)
raw = bf.knnMatch(desc_a, desc_b, k=2)

# Lowe's ratio test
good = [m for m, n in raw if m.distance < 0.85 * n.distance]

GPU Matching (Zero PCIe Transfer)

When using gpu_output=True, descriptors are CuPy arrays in VRAM. Match entirely on GPU using matmul-based brute-force:

import cupy as cp

def gpu_bf_match_mutual(desc_a, desc_b, ratio=0.85):
    """Mutual nearest-neighbor matching on GPU."""
    ratio_sq = ratio * ratio

    # L2 distance matrix via matmul
    a_sq = cp.sum(desc_a * desc_a, axis=1, keepdims=True)
    b_sq = cp.sum(desc_b * desc_b, axis=1, keepdims=True).T
    dist_sq = a_sq + b_sq - 2.0 * cp.matmul(desc_a, desc_b.T)
    cp.maximum(dist_sq, 0.0, out=dist_sq)

    # Forward: Lowe's ratio test (A -> B)
    top2 = cp.argpartition(dist_sq, 1, axis=1)[:, :2]
    rows = cp.arange(dist_sq.shape[0])
    d0 = dist_sq[rows, top2[:, 0]]
    d1 = dist_sq[rows, top2[:, 1]]
    swap = d0 > d1
    best_f = cp.where(swap, top2[:, 1], top2[:, 0])
    d_best = cp.where(swap, d1, d0)
    d_second = cp.where(swap, d0, d1)
    pass_ratio = (d_best < ratio_sq * d_second) & (d_second > 0)

    # Backward: nearest neighbor (B -> A)
    best_b = cp.argmin(dist_sq, axis=0)

    # Mutual check
    mutual = pass_ratio & (best_b[best_f] == rows)
    q = cp.asnumpy(cp.where(mutual)[0])
    t = cp.asnumpy(best_f[mutual])
    return q, t

Why Mutual Nearest-Neighbor?

RANSAC: From Matches to Geometry

# MAGSAC++ with aggressive settings
F, mask = cv2.findFundamentalMat(
    pts_a, pts_b,
    cv2.USAC_MAGSAC,
    ransacReprojThreshold=2.0,
    confidence=0.9999,
    maxIters=10000
)

DSP-SIFT: Multi-Scale Descriptor Pooling

DSP-SIFT (Dong & Soatto, CVPR 2015) computes descriptors at multiple Gaussian scale radii around each keypoint and averages the histograms. This makes each descriptor an implicit mixture over scales, dramatically improving robustness to scale estimation errors and viewpoint changes — especially important for wide-baseline stereo in 3D reconstruction and drone aerial imagery with varying altitude.

How DSP Pooling Works

1. For each keypoint, compute the standard descriptor patch radius r
2. Recompute the descriptor at 3 scaled radii: r/√2, r, r√2
3. Average the 3 histograms element-wise
4. Apply L2 + clip(0.2) + L2 normalization (then RootSIFT)

The output is still a 128-dimensional descriptor — same size, better quality. No change to matching code required.

Enable DSP-SIFT

# One flag enables DSP pooling
sift = PySIFT(dsp=True)

# Everything else stays the same
kp, desc = sift.detectAndCompute(gray, None)
# desc.shape = (N, 128) -- same dimensionality

Impact on Accuracy

GPU Implementation

PySIFT fuses all DSP scales into a single CUDA kernel launch per octave. Each thread block processes one keypoint, looping over the 3 scale multipliers and accumulating into shared memory. This avoids the naive approach of launching 3 separate descriptor kernels (3x overhead).

Who is PySIFT For?

Every vision pipeline that matches images — from surgical navigation to autonomous flight — starts with feature extraction. PySIFT replaces the CPU bottleneck at the foundation of these pipelines with a GPU-resident alternative that requires no retraining, no C++ compilation, and no domain-specific tuning.

# Register two histopathology slide regions
from pysift import PySIFT
sift = PySIFT(dsp=True)
kp1, d1 = sift.detectAndCompute(slide_a)
kp2, d2 = sift.detectAndCompute(slide_b)

# Stitch two aerial frames on Jetson
from pysift import GPUPyStitch
stitcher = GPUPyStitch()
panorama = stitcher.stitch(frame_left, frame_right)

# GPU-resident features for VO pipeline
kp, desc = sift.detectAndCompute(
    frame, None, gpu_output=True
)
# desc is CuPy array in VRAM -- zero-copy
torch_desc = torch.from_dlpack(desc)

# Replace COLMAP's CPU SIFT extraction
sift = PySIFT(dsp=True)
for img in scene_images:
    kp, desc = sift.detectAndCompute(
        img, None, gpu_output=True
    )

Why PySIFT Over Learned Detectors for Downstream Pipelines?

Try PySIFT

pip install staysift

Full Pipeline Example

from pysift import PySIFT
import cv2
import numpy as np

# 1. Detect with DSP-SIFT
sift = PySIFT(dsp=True)
img_a = cv2.imread("scene_left.jpg", cv2.IMREAD_GRAYSCALE)
img_b = cv2.imread("scene_right.jpg", cv2.IMREAD_GRAYSCALE)

kp_a, desc_a = sift.detectAndCompute(img_a, None)
kp_b, desc_b = sift.detectAndCompute(img_b, None)

# 2. Match with ratio test
bf = cv2.BFMatcher(cv2.NORM_L2)
matches = bf.knnMatch(desc_a, desc_b, k=2)
good = [m for m, n in matches if m.distance < 0.85 * n.distance]

# 3. Estimate fundamental matrix
pts_a = np.float32([kp_a[m.queryIdx].pt for m in good])
pts_b = np.float32([kp_b[m.trainIdx].pt for m in good])

F, mask = cv2.findFundamentalMat(
    pts_a, pts_b, cv2.USAC_MAGSAC, 2.0, 0.9999, 10000
)
inliers = mask.ravel().sum()
print(f"Matches: {len(good)}, Inliers: {inliers}")

Advanced: GPU-Resident Pipeline

import cupy as cp

# Detect -- descriptors stay in VRAM
kp_a, desc_a = sift.detectAndCompute(img_a, None, gpu_output=True)
kp_b, desc_b = sift.detectAndCompute(img_b, None, gpu_output=True)

# Match -- GPU matmul, zero PCIe transfer
q, t = gpu_bf_match_mutual(desc_a, desc_b, ratio=0.85)

# Only small coordinate arrays transfer to CPU for RANSAC
pts_a = cp.asnumpy(kp_a[q, :2])
pts_b = cp.asnumpy(kp_b[t, :2])

Resources