Sub-pixel Edge Detection Techniques for High-Precision Metrology

1. Executive Summary & Problem Formulation

In automated industrial metrology, camera systems measure manufactured parts—like engine pistons or microchip wafers—to verify they conform to CAD specifications. If a piston diameter is off by $10\text{ \mu m}$, the engine seizes.

Software engineers tasked with measuring parts usually drop a camera over the assembly line, run the standard OpenCV Canny edge detector (cv2.Canny), find the contour boundaries, and calculate the diameter in pixels.

This approach fails the gauge repeatability requirement immediately. The Canny edge detector operates on the integer pixel grid. It assigns the location of a physical edge to the exact center of a single sensor pixel ($x=142, y=98$). If the camera's optical magnification maps $1\text{ pixel}$ to $20\text{ \mu m}$ of physical space, the absolute best theoretical measurement error you can achieve is $\pm 10\text{ \mu m}$. When the part vibrates on the conveyor belt and shifts by a fraction of a pixel, the Canny edge detector jumps a full integer pixel, causing massive, unpredictable quantization noise in the measurement data.

To achieve micrometer precision without buying a $50,000$ electron microscope, we must break the integer pixel barrier. We must extract sub-pixel edge locations by analyzing the continuous gradient intensity profile across the image pixels, interpolating the exact geometric point where the physical edge crossed the sensor matrix to a precision of $\frac{1}{10}$ or even $\frac{1}{100}$ of a pixel.

2. Mathematical & Architectural Theory

The Flaw of the Canny Edge Detector

Canny edge detection convolves the image with Sobel filters to calculate the gradient magnitude at every pixel. It then performs non-maximum suppression (NMS): if a pixel's gradient magnitude is strictly greater than its neighbors along the gradient direction, it is preserved as an edge; otherwise, it is forced to zero.

Because NMS operates on discrete integer coordinates, it throws away all the critical gradient intensity information in the neighboring pixels. The physical edge didn't snap cleanly to the center of pixel $142$. The photon boundary hit the sensor array and blurred across pixels $141, 142$, and $143$ due to the optical Point Spread Function (PSF) of the lens.

Sub-pixel Parabolic Interpolation

To locate the true sub-pixel peak of the gradient, we capture the gradient magnitude at the integer peak pixel ($g_0$), the pixel behind it ($g_{-1}$), and the pixel ahead of it ($g_{+1}$) along the gradient direction.

We model the gradient intensity as a continuous 1D quadratic parabola: $$g(x) = a x^2 + b x + c$$

Setting the origin $x=0$ at the integer peak pixel ($g_0$), $x=-1$ at ($g_{-1}$), and $x=1$ at ($g_{+1}$), we solve the system of equations for the coefficients. The sub-pixel peak of the parabola (where the derivative equals zero) gives the precise sub-pixel shift $\delta$: $$\delta = \frac{g_{-1} - g_{+1}}{2(g_{-1} - 2g_0 + g_{+1})}$$

The true sub-pixel location of the edge is then $x_{true} = x_{integer} + \delta$. Since this relies on a symmetric parabolic assumption, it provides a stable approximation up to roughly $0.1$ pixels of accuracy.

Zernike Moments for Sub-pixel Edge Detection

When metrology demands $0.01$ pixel accuracy, the 1D parabolic model fails because the optical blur function is 2D and non-linear. Zernike moments offer a mathematically rigorous alternative by projecting the local circular image patch onto a set of orthogonal complex basis polynomials.

An edge within a unit circle can be completely parameterized by three geometric variables: 1. $k$: The background intensity. 2. $h$: The edge step magnitude (foreground minus background). 3. $l$: The perpendicular distance from the center of the pixel to the actual physical edge line.

By calculating just three Zernike moments ($A_{00}$, $A_{11}$, $A_{20}$) using fixed convolution masks over a $7 \times 7$ pixel window, we can algebraically solve for the precise distance $l$ and the angle of the edge. Because Zernike moments integrate the intensity across a 2D area, they are highly robust to sensor noise compared to 1D gradient interpolation.

3. Concrete Implementation: Sub-pixel 1D Interpolation

Implementing full Zernike moments in pure Python is computationally heavy. For high-speed inspection pipelines, 1D parabolic interpolation applied over the Sobel gradient provides the optimal balance of speed and sub-pixel accuracy.

Below is an implementation that locates straight edges (like the width of a metal block) by taking a 1D cross-section array of pixels, computing the gradient, and applying parabolic interpolation to find the exact sub-pixel boundary width.

import cv2
import numpy as np

def compute_1d_gradient(profile: np.ndarray) -> np.ndarray:
    """
    Computes the 1D spatial derivative using a central difference kernel.
    Input: 1D array of pixel intensities.
    Output: 1D array of gradient magnitudes.
    """
    # Central difference kernel: [-0.5, 0, 0.5]
    kernel = np.array([-0.5, 0.0, 0.5])
    # Use valid padding to avoid edge artifacts
    gradient = np.convolve(profile, kernel, mode='same')
    return np.abs(gradient)

def find_subpixel_peaks(gradient: np.ndarray, threshold: float) -> list[float]:
    """
    Locates sub-pixel peaks in a 1D gradient array using parabolic interpolation.
    """
    peaks = []
    
    # Iterate through the array, skipping the very first and last pixels
    for x in range(1, len(gradient) - 1):
        g_prev = gradient[x - 1]
        g_curr = gradient[x]
        g_next = gradient[x + 1]
        
        # Check for integer local maximum (Non-maximum suppression)
        if g_curr >= g_prev and g_curr >= g_next and g_curr > threshold:
            
            # Compute parabolic sub-pixel shift delta
            # Prevent division by zero if the gradient is perfectly flat
            denominator = 2.0 * (g_prev - 2.0 * g_curr + g_next)
            if denominator == 0:
                delta = 0.0
            else:
                delta = (g_prev - g_next) / denominator
                
            # Clamp the delta to ensure it stays within the [-0.5, 0.5] pixel range
            delta = max(-0.5, min(0.5, delta))
            
            # Store the absolute sub-pixel coordinate
            subpixel_coord = float(x) + delta
            peaks.append(subpixel_coord)
            
    return peaks

def measure_part_width(image_path: str, scan_y: int) -> float:
    """
    Measures the physical width of a part across a specific horizontal scanline.
    """
    # Load image as grayscale float32 for precise math
    img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
    if img is None:
        raise ValueError("Image not found")
        
    img_f32 = img.astype(np.float32)
    
    # Extract the 1D horizontal intensity profile at the specified Y coordinate
    profile = img_f32[scan_y, :]
    
    # Smooth the profile slightly with a 1D Gaussian to suppress sensor noise
    profile_smooth = cv2.GaussianBlur(profile, (5, 1), 0).flatten()
    
    # Calculate gradients
    gradient = compute_1d_gradient(profile_smooth)
    
    # Extract peaks (Set threshold based on minimum expected contrast)
    peaks = find_subpixel_peaks(gradient, threshold=15.0)
    
    if len(peaks) < 2:
        print("Failed to find two clear part boundaries.")
        return 0.0
        
    # Assuming the part is a dark block on a light background, 
    # we take the distance between the first and last detected edges
    left_edge_px = peaks[0]
    right_edge_px = peaks[-1]
    
    pixel_width = right_edge_px - left_edge_px
    print(f"Left Edge: {left_edge_px:.4f} px")
    print(f"Right Edge: {right_edge_px:.4f} px")
    print(f"Sub-pixel Width: {pixel_width:.4f} px")
    
    return pixel_width

if __name__ == '__main__':
    # Pseudo-execution context
    # Suppose our camera is calibrated to 0.025 mm per pixel
    # width_px = measure_part_width("machined_part.png", scan_y=480)
    # physical_width_mm = width_px * 0.025
    # print(f"Physical Width: {physical_width_mm:.4f} mm")
    pass

4. Edge Cases, Optimization & Memory Considerations

Telecentric Lens Necessity

Sub-pixel algorithms assume the pixel boundaries directly correlate to the physical boundaries of the object. Standard entocentric camera lenses suffer from perspective distortion. If a machined cylinder shifts slightly closer to a standard lens, it will occupy a wider footprint on the sensor array. The algorithm will correctly report a wider pixel measurement, causing the system to falsely fail the part for being physically too large.

High-precision metrology requires a Telecentric Lens. These lenses have an entrance pupil at infinity, ensuring constant magnification regardless of the object's distance from the camera. If the part shifts closer or further away, its size on the sensor array remains mathematically constant.

Sensor Illumination and Chromatic Aberration

If you illuminate a steel part with generic white LEDs, the sub-pixel accuracy collapses. White light consists of multiple wavelengths. Because camera lenses refract different wavelengths at different angles (chromatic aberration), the red channel edge will physically map to a different sub-pixel location than the blue channel edge. If you capture in grayscale, these edges smear together, destroying the sharp gradient peak.

Metrology systems use monochromatic illumination (usually narrow-band red or blue lasers/LED arrays paired with corresponding physical bandpass filters screwed onto the lens). This entirely eliminates chromatic aberration and produces a pristine, mathematically modeled gradient curve.

Sub-pixel Contour Extraction

If you need to measure the roundness of a circular bore, you cannot just scan a 1D line. You must extract a continuous 2D sub-pixel curve. You can achieve this by first finding the standard integer contour using cv2.findContours. For every point on that integer contour, compute the gradient normal vector, sample the intensity pixels along that specific normal vector, and execute the 1D parabolic interpolation to shift that single contour point onto the sub-pixel boundary.

5. Benchmarks & Practical Engineering Takeaways

We benchmarked integer Canny edge extraction vs Parabolic Sub-pixel interpolation on a static $50.000\text{ mm}$ gauge block, capturing 1,000 frames under factory vibration conditions.

Engineering Guidelines

• Disable Image Compression: Never save measurement templates or calibration images as JPEGs. The DCT block quantization algorithms inside JPEG ruin sub-pixel gradient profiles, injecting invisible block artifacts that shift the parabolic curve. Always use uncompressed TIFF or PNG formats.
• Perform Camera Calibration: Lens distortion will warp straight edges. You must undistort the image coordinates before you compute physical measurements, or the sub-pixel precision is wasted on a distorted geometry space.
• Do not upscale images: Some developers attempt "sub-pixel" accuracy by resizing the image by 400% using cubic interpolation and running integer Canny. This is mathematically invalid. Image interpolation algorithms invent fake pixel data based on smoothing splines, which physically alters the location of the actual gradient peak.

6. References & Cross-Links

• Ghosal, S., & Mehrotra, R. (1993). Orthogonal moment operators for subpixel edge detection. Pattern Recognition.
• Steger, C. (1998). An unbiased detector of curvilinear structures. IEEE Transactions on Pattern Analysis and Machine Intelligence.
• Susam, A. (2026). Stereoscopic Vision Depth Estimation: Algorithms and Calibration. Read Article.
• Susam, A. (2026). Hardware-Level Synchronization of Multi-Camera Systems using GPIO Triggers. Read Article.