Passive Rangefinding Without LiDAR: Depth-From-Defocus and Optical Physics in Pure NumPy

Part 2 of the Hawk Eye Spatial Qualimetry series. Part 1 covered SVD plane fitting and trimmed OLS for 3D tread profiles. Here we answer the question Part 1 left open: how far away was the camera?

Hawk Eye is the internal codename. These components are coming soon to *Tyre Cam AI** (free scans are limited, and a free account unlocks more). This post covers only the rangefinding layer.*

---

1. The problem: relative depth has no units

In Part 1 we extracted a 3D tread relief from a single photo. The output is relative: a monocular U-Net returns values in [0, 1], and the Poisson-integrated relief gives a delta-Z in arbitrary units. Turning that into millimetres needs a scale, and the scale depends on something we never measured: the distance between the lens and the rubber.

The same 8 mm groove spans very different pixel counts at 30 cm and at 90 cm, and its shading gradients change too. If the pipeline treats both shots identically, the fusion layer ends up trusting a 3D estimate that is physically less reliable at one end of the range.

The obvious fix is a depth sensor, and it does not scale:

• LiDAR / ToF is a small minority of devices: roughly 5% in our target fleet (an internal estimate, not a survey).
• Stereo needs two calibrated cameras on a rigid baseline. A mechanic with a phone has neither.
• Hardcoded tire geometry breaks the moment someone photographs a 175 mm hatchback tire instead of a 245 mm SUV tire.

So the Silver Kernel Rangefinder is a passive monocular estimator that fuses two independent cues:

1. Geometry: the angular scale of the detected tread ROI in the frame (a pinhole relation).
2. Optical physics: how defocused and how contrast-attenuated the texture inside the ROI is (Depth-from-Defocus, DfD).

It works from 25 to 100 cm and needs only NumPy plus one PIL resize. Project constraints: no third-party CV packages in the core, no LiDAR, no stereo, no per-tire dimension tables.

![Figure 1. Panel 3 of the Hawk quad-panel HUD: live AI Overlay showing monocular depth estimation (5.7 mm), rangefinder distance readout (Dist.: 59 cm, Z_geom: 58 cm, OptF: 0.990, C_michel: 0.35, Conf: 93%), and 100% model consensus.]

Figure 1. Panel 3 of the Hawk quad-panel HUD: live AI Overlay showing monocular depth estimation (5.7 mm), rangefinder distance readout (Dist.: 59 cm, Z_geom: 58 cm, OptF: 0.990, C_michel: 0.35, Conf: 93%), and 100% model consensus.

---

2. The math

2.1 Geometry cue: pinhole angular scale

Let r = w_roi / w_frame be the fraction of the frame width occupied by the tread ROI. For a pinhole camera, apparent size is inversely proportional to distance:

$$z = b + \frac{K}{r}, \qquad r = \frac{w_{\text{roi}}}{w_{\text{frame}}}$$

K is one scalar per camera system (focal length, sensor width and the framing behaviour of the ROI detector folded together), and b is a small lens-to-reference offset. Section 5 is explicit about what K does and does not assume.

Because z is proportional to 1/r, relative errors pass through 1:1:

$$\frac{\delta z}{z} = -\frac{\delta r}{r}$$

A 3 px jitter on a 600 px ROI is about 0.5% of distance, far below the other error sources.

2.2 Optical cue: thin-lens defocus

For a thin lens with focal length f and f-number N, focused at z_f, a point at distance z is imaged as a blur circle of diameter

$$c(z) = \frac{f^{2}}{N\,(z_f - f)} \cdot \frac{|z - z_f|}{z}$$

Blur grows with |z - z_f|, so sharpness is a magnitude cue. Nearer and farther than the focal plane look alike. Geometry resolves that sign.

To measure blur without a reference image we use the variance of the discrete Laplacian of the luma channel:

$$\nabla^{2} I_{x,y} = I_{x+1,y} + I_{x-1,y} + I_{x,y+1} + I_{x,y-1} - 4\,I_{x,y}$$

Why it works. Natural textures have a power-law power spectrum $P(\omega) \propto \omega^{-\alpha}$ with $\alpha \approx 2$ . A Gaussian blur of width $\sigma_b$ multiplies power by $e^{-\sigma_b^2\omega^2}$ , and the Laplacian multiplies it by $\omega^4$ . Integrating over the 2D frequency plane:

$$\text{Var}(\nabla^{2} I) \;\propto\; \int_0^\infty \omega^{5-\alpha}\, e^{-\sigma_b^{2}\omega^{2}}\, d\omega \;\propto\; \sigma_b^{-(6-\alpha)} \;\overset{\alpha=2}{=}\; \sigma_b^{-4}$$

So $\sigma_{\nabla^2} \propto \sigma_b^{-2}$ , and the reciprocal $0.04/\sigma_{\nabla^2}$ grows like the square of the blur radius. The model holds while blur is resolvable ( $\sigma_b \gtrsim 1$ px). Below that, the pixel grid and sensor noise set a floor. This is why every ROI is resampled to a fixed 224x224 grid before measuring anything.

2.3 Michelson contrast, robust version

$$C_{\text{michel}} = \frac{I_{\max} - I_{\min}}{I_{\max} + I_{\min}} \;\;\longrightarrow\;\; \frac{P_{95}(I) - P_{5}(I)}{P_{95}(I) + P_{5}(I) + \varepsilon}$$

Raw min/max are single-pixel statistics, dominated by specular highlights on wet rubber and by hot pixels. Percentiles make the metric robust.

Contrast carries distance information because a fixed physical texture period maps to a pixel period proportional to 1/z. As the camera backs off, micro-texture approaches the sensor Nyquist limit where lens and pixel MTF attenuate it. Contrast therefore decays with distance even when the lens is in focus. That is a second cue with a failure profile different from Laplacian variance.

2.4 Combining the cues: the optic factor

The project uses "optical factor" in two roles, and it helps to keep them apart:

• optic_factor is a bounded multiplier on distance, defined below.
• k_opt = sigma_luma / mu_luma is a quality gate used by the fusion layer (Section 4).

$$s = \tfrac12\left(\frac{0.04}{\sigma_{\nabla^2} + 10^{-3}} - 1\right) + \tfrac12\, C_{\text{michel}}, \qquad \Phi = \text{clip}\big(1 + w_{\text{opt}}\, s,\; 0.60,\; 1.40\big)$$

$$z_{\text{final}} = b + \frac{K}{r}\,\Phi$$

with $w_{\text{opt}} = 0.15$ . The constants here (0.04, the 50/50 split, w_opt, the clamp) are empirically tuned on our own captures. Treat them as calibrated heuristics, not laws of nature. 0.04 is the Laplacian sigma of in-focus reference captures on the 224 grid, so the first term is about 0 for a reference-sharp ROI. The 50/50 split reflects that we had no evidence to prefer one cue, and with a small w_opt the result is insensitive to it.

The design principle matters more than the constants:

• w_opt is small on purpose. Geometry is the backbone. Optics nudges it and, more importantly, tells us how much to trust it.
• The +/-40% clamp guarantees that a motion-blurred or out-of-focus frame cannot send the estimate to infinity.

One more use of spatial frequency: we split the Laplacian into 16x16 tiles and compute a local-variance matrix. The fraction of tiles carrying texture is a coverage score, which rejects ROIs half-filled with asphalt or fender. The per-tile heatmaps we use for diagnostics are not in the hot path, and we can share that code on request.

---

3. The code

Self-contained: pip install numpy pillow, paste, run.

3.1 Primitives

from __future__ import annotations

from typing import NamedTuple

import numpy as np
from PIL import Image

CANON = 224                                    # same grid as the ONNX input tensors
_LUMA = np.array([0.299, 0.587, 0.114], dtype=np.float32)


def canonicalize(roi_rgb: np.ndarray, size: int = CANON) -> np.ndarray:
    """Resize ROI to size x size uint8 RGB so thresholds are resolution-invariant."""
    a = np.asarray(roi_rgb)
    if a.dtype != np.uint8:                    # assume float in [0, 1]
        a = (np.clip(a, 0.0, 1.0) * 255.0 + 0.5).astype(np.uint8)
    if a.shape[:2] != (size, size):
        a = np.asarray(Image.fromarray(a).resize((size, size), Image.Resampling.BILINEAR))
    return a


def to_luma(rgb_u8: np.ndarray) -> np.ndarray:
    """(H, W, 3) uint8 -> (H, W) float32 in [0, 1]."""
    return (rgb_u8.astype(np.float32) @ _LUMA) * (1.0 / 255.0)


def laplacian4(y: np.ndarray) -> np.ndarray:
    """4-neighbour discrete Laplacian, valid region -> (H-2, W-2)."""
    return (y[1:-1, 2:] + y[1:-1, :-2] + y[2:, 1:-1] + y[:-2, 1:-1]
            - 4.0 * y[1:-1, 1:-1])


def p05_p95(y: np.ndarray) -> tuple[float, float]:
    """5th/95th percentiles from one np.partition (cheaper than np.percentile)."""
    f = y.ravel()
    k0, k1 = int(0.05 * (f.size - 1)), int(0.95 * (f.size - 1))
    part = np.partition(f, (k0, k1))
    return float(part[k0]), float(part[k1])

3.2 The rangefinder

class RangeResult(NamedTuple):
    distance_cm: float
    confidence: float          # heuristic score in [0, 1], not a calibrated probability
    info: dict


def estimate_distance_pure_numpy(
    roi_rgb: np.ndarray,
    frame_width_px: int,
    roi_width_px: int | None = None,
    *,
    k_scale_mm: float = 256.25,       # camera-system constant, fit with calibrate_scale()
    offset_mm: float = 15.0,          # lens-to-reference-plane offset
    optic_weight: float = 0.15,       # 0.0 == pure pinhole
    z_range_cm: tuple[float, float] = (25.0, 100.0),
    z_focus_mm: float | None = None,  # optional, see note below
    tile: int = 16,
) -> RangeResult:
    """
    Passive monocular distance from one RGB tread ROI.

    roi_rgb        : (H, W, 3) uint8 (or float in [0, 1]) crop of the tread
    frame_width_px : width of the full frame the ROI was cropped from
    roi_width_px   : ROI width in the ORIGINAL frame (default: roi_rgb.shape[1])
    Returns        : distance in CENTIMETRES (internally computed in mm) and a
                     confidence in [0, 1].
    """
    if frame_width_px <= 0:
        raise ValueError("frame_width_px must be positive")
    roi = np.asarray(roi_rgb)
    if roi.ndim != 3 or roi.shape[2] != 3 or min(roi.shape[:2]) < 32:
        raise ValueError("roi_rgb must be (H, W, 3) with H, W >= 32")
    w_roi = float(roi_width_px if roi_width_px is not None else roi.shape[1])

    # 1. geometry cue (pinhole angular scale)
    r = float(np.clip(w_roi / frame_width_px, 1e-3, 1.0))
    z_geom = k_scale_mm / r

    # 2. optical cues on the canonical grid
    y = to_luma(canonicalize(roi))
    lap = laplacian4(y)
    l_var = float(lap.var())                       # focus measure
    lap_std = l_var ** 0.5
    # edge / lamella density; 0.05 is an empirically tuned threshold on the normalized Laplacian
    s_freq = float((np.abs(lap) > 0.05).mean())
    # stride 2: 4x cheaper, same stats on natural images (synthetic high-freq patterns may alias)
    p05, p95 = p05_p95(y[::2, ::2])
    c_michel = (p95 - p05) / (p95 + p05 + 1e-6)
    k_opt = float(y.std()) / (float(y.mean()) + 1e-6)

    h, w = (lap.shape[0] // tile) * tile, (lap.shape[1] // tile) * tile
    vmap = lap[:h, :w].reshape(h // tile, tile, w // tile, tile).var(axis=(1, 3))
    # a tile counts as textured if it carries at least a fifth of the global Laplacian energy
    coverage = float((vmap > 0.2 * l_var + 1e-9).mean())

    # 3. DfD signal -> bounded multiplicative correction
    s = 0.5 * (0.04 / (lap_std + 1e-3) - 1.0) + 0.5 * c_michel
    sign = -1.0 if (z_focus_mm is not None and z_geom < z_focus_mm) else 1.0
    optic_factor = float(np.clip(1.0 + optic_weight * sign * s, 0.60, 1.40))
    z_cm = (offset_mm + z_geom * optic_factor) / 10.0

    direction = "weak" if abs(s) < 0.05 else ("inverted" if sign < 0 else "normal")

    # 4. confidence: geometric mean, so ANY failing factor drags it down
    # sqrt on the texture term compresses the low end: moderately sharp frames are not over-penalized
    c_tex = float(np.sqrt(np.clip(l_var / 0.005, 0.0, 1.0)))
    c_con = float(np.clip(c_michel / 0.20, 0.0, 1.0))
    c_cov = float(np.clip(coverage / 0.50, 0.0, 1.0))
    outside = max(z_range_cm[0] - z_cm, z_cm - z_range_cm[1], 0.0)
    c_rng = float(np.exp(-outside / 20.0))
    confidence = float((c_tex * c_con * c_cov * c_rng) ** 0.25)

    info = dict(
        z_geom_mm=z_geom, r_percent=100.0 * r, optic_factor=optic_factor,
        dfd_signal=float(s), dfd_direction=direction, l_var=l_var, lap_std=lap_std,
        s_freq=s_freq, c_michel=float(c_michel), k_opt=k_opt, coverage=coverage,
        is_sharp=l_var >= 0.005,
        conf_parts=dict(tex=c_tex, contrast=c_con, coverage=c_cov, range=c_rng),
    )
    return RangeResult(z_cm, confidence, info)

The thresholds in the confidence block (0.005, 0.20, 0.50, the 20 cm range decay) are, like the optic-factor constants, tuned on our captures. The 0.2 factor in coverage and the sqrt on the texture term are empirical too, and the comments in the code say what each one does. The geometric mean is a deliberate choice. With an arithmetic mean, excellent contrast could mask a dead texture channel. With a geometric mean, one near-zero factor vetoes the score.

Optional z_focus_mm. Defocus is symmetric around the focal plane. If the focus distance is known or locked (fixed-focus module, AF lock) and the geometric estimate says the object is nearer than the focal plane, the sign of the optical correction flips and the diagnostic label becomes "inverted". With None, the function reproduces the production formula. One caveat: this resolver trusts z_geom, which is itself an estimate. If the geometric error is large (a 20% miss is not unusual), the sign can flip wrongly near the focal plane. That is why optic_weight is small and z_focus_mm defaults to None. With free-running autofocus, leave it at None.

3.3 Try it

rng = np.random.default_rng(42)
x = np.linspace(0, 8 * np.pi, 300)
ribs = 0.5 + 0.25 * np.sign(np.sin(x))[None, :].repeat(240, axis=0)   # synthetic tread
roi = np.clip(ribs + rng.normal(0, 0.05, ribs.shape), 0, 1)
roi = np.repeat(roi[..., None], 3, axis=2)

# Pretend this crop was 900 px wide in a 1920 px frame:
dist_cm, conf, info = estimate_distance_pure_numpy(roi, 1920, roi_width_px=900)
print(f"{dist_cm:.1f} cm, confidence {conf:.2f}, DfD: {info['dfd_direction']}")

Replace roi with np.asarray(pil_crop) for real photos.

3.4 Calibrating the scale: the same lstsq as Part 1

K and b should be fitted, not typed in. Shoot any tread at 5-10 known ruler distances and solve, in the least-squares sense:

$$\begin{bmatrix} 1/r_1 & 1 \ \vdots & \vdots \ 1/r_n & 1 \end{bmatrix} \begin{bmatrix} K \ b \end{bmatrix} \approx \begin{bmatrix} z_1 \ \vdots \ z_n \end{bmatrix}$$

def calibrate_scale(roi_w_px, frame_w_px, z_true_mm):
    """Fit z = K/r + b. Returns (k_scale_mm, offset_mm, condition_number)."""
    r = np.asarray(roi_w_px, dtype=np.float64) / float(frame_w_px)
    A = np.stack([1.0 / r, np.ones_like(r)], axis=1)
    (k, b), *_ = np.linalg.lstsq(A, np.asarray(z_true_mm, dtype=np.float64), rcond=None)
    return float(k), float(b), float(np.linalg.cond(A))

Use in-focus shots so that optic_factor is about 1 during calibration. As in Part 1, watch the condition number. It stays small when the shots span the working range (30, 45, 60, 80, 100 cm) and blows up if everything is shot at one distance.

---

4. What the distance is actually used for

The distance never converts pixels to millimetres. That job belongs to the TWI anchor, a physical 1.6 mm step covered in Part 3. The range feeds the fusion layer, which decides how much to trust the 3D relief branch versus the CNN wear classifier:

$$w_{\text{base}}(z) = 0.85 - \text{clip}!\left(\frac{z - 800}{700},\,0,\,1\right)\cdot 0.55, \qquad m_{\text{opt}} = \min!\left(1, \tfrac{k_{\text{opt}}}{0.25}\right)\cdot\min!\left(1,\tfrac{l_{\text{var}}}{0.005}\right)$$

with z in mm. The 3D weight is w_base * m_opt * m_twi, clipped to [0.15, 0.85]. The TWI multiplier m_twi and the conflict arbiter are the subject of Part 3. All of these constants are tuned heuristics too.

Two properties matter here. First, m_opt is where the optics acts as a quality gate: a soft, flat frame makes gradient-based Poisson integration unreliable, so the weight shifts to the classifier. Second, the sensitivity to distance is low. In the current build the ramp starts at 800 mm, so the 3D branch sits on a high-trust plateau for most of the 25-100 cm range. Where it ramps, the slope is $0.55/700\ \text{mm} \approx 0.008$ weight per cm, so even a 15 cm error moves $w_{3D}$ by about 0.12. A TWI verdict (x1.25) or a conflict penalty (x0.20) moves it far more.

The fusion layer needs distance to roughly +/-15%, not +/-1%. That is why a passive estimator is good enough.

---

5. Why We Don't Need Hardcoded Tire Dimensions

"Scale from a single image" is where most monocular systems quietly cheat, so let's be precise.

5.1 What the physics gives us for free

We use no tire-size database, no sidewall OCR, no reference object in the frame and no per-vehicle configuration. Every quantity the rangefinder uses is a ratio or a statistic of the image itself:

Cue	What is measured	Needs a tire dimension?
Angular scale r = w_roi / w_frame	Ratio of two pixel lengths	No, only one camera constant K
Defocus (Var of Laplacian)	Spatial-frequency energy of the texture	No, depends on lens PSF, not object size
Michelson contrast	Percentile ratio of luma	No, dimensionless
k_opt	sigma / mu of luma	No, dimensionless

Defocus and contrast are properties of the optical system (aperture, focal plane, MTF). They tie blur to the object-to-focal-plane distance whatever the object is.

5.2 The honest part: one metric anchor is unavoidable

A single image cannot recover absolute scale from geometry alone. Something has to anchor metres. In our design that anchor is the single constant K, and it is:

• a camera-system property, fitted once per device model with calibrate_scale();
• not looked up per tire.

Our shipped default, 256.25 mm, happens to equal 1.25 x 205 mm, because that is how the first calibration was derived. It is stored and treated as one number, and the code never reads a tire width.

What if the photographed tread is wider or narrower than what the calibration implicitly averaged over? The geometric distance scales by the same ratio. For common passenger widths (175-245 mm around a 205 mm baseline) that is roughly -15% / +20% on z_geom. Two things keep this from turning into a wrong answer:

1. Fusion tolerance. As shown in Section 4, the downstream weights are insensitive to +/-15% distance error.
2. Confidence instead of false precision. When the optics say a frame is unusable (soft, flat, half-covered, out of range), confidence drops and fusion leans on the classifier.

If you need tighter absolute accuracy for a vehicle class, recalibrate K for it. It is a five-minute ruler session, not a code change.

5.3 Known failure modes

Failure	Symptom	Mitigation
Motion blur	Looks like defocus, l_var collapses	Low c_tex lowers confidence. Next step: directional Laplacian anisotropy test
Autofocus hunting	Focal plane unknown, DfD sign flips	Small optic_weight, z_focus_mm=None
Wet rubber / specular glare	Inflated l_var, spiky I_max	Percentile Michelson, +/-40% clamp
ROI partly off-tread	Misleading r, texture gaps	coverage term in confidence
Dust / mud film	Low contrast at any distance	Low c_con, TWI arbiter in fusion
Flat uniform lighting	Weak DfD signal	dfd_direction = "weak", geometry dominates

![Figure 2. Panel 2 of the Hawk HUD: Poisson-integrated 3D surface reconstruction relief with detected TWI anchors (Delta Z: 0.638, Approx Depth: 5.42 mm).]

Figure 2. Panel 2 of the Hawk HUD: Poisson-integrated 3D surface reconstruction relief with detected TWI anchors (Delta Z: 0.638, Approx Depth: 5.42 mm).

---

6. Benchmark: under 1 ms on one CPU core

Setup: Intel(R) Core(TM) i7-1165G7 @ 2.80GHz, single thread, NumPy 1.26.4, Python 3.11. Median and p95 over 2000 calls after 200 warm-ups. Input is uniform noise, which makes the texture statistics a worst case.

import time

roi = (np.random.default_rng(0).random((224, 224, 3)) * 255).astype(np.uint8)
for _ in range(200):
    estimate_distance_pure_numpy(roi, 1920)
t = []
for _ in range(2000):
    t0 = time.perf_counter_ns()
    estimate_distance_pure_numpy(roi, 1920)
    t.append(time.perf_counter_ns() - t0)
print(f"median {np.median(t)/1e6:.3f} ms | p95 {np.percentile(t, 95)/1e6:.3f} ms")

Input	Median	p95	Notes
224x224 RGB (already canonical)	0.38 ms	0.48 ms	Hot path: luma, Laplacian, tile variance, partition
128x128 (upsampled to canonical)	0.42 ms	0.52 ms	Resize of a tiny input is cheap
640x480 ROI (incl. PIL resize)	0.85 ms	1.05 ms	Resize dominates, and ONNX preprocessing pays it anyway

Why it is fast:

• No pixel loops. Everything is slicing, reshapes and reductions.
• One np.partition serves both percentiles, on a stride-2 subsample (about 12.5k elements instead of 50k).
• A single (H, W, 3) @ (3,) matmul for luma, and a fixed 224x224 grid, so cost is constant and thresholds are resolution-invariant.

This is orders of magnitude below any CNN stage in the cascade. The rangefinder can run on every frame of a live preview without touching the GUI-thread budget, which is the topic of Part 4.

---

7. Takeaways

1. Relative depth needs a scale, and the scale is the camera-to-rubber distance. It can be estimated passively.
2. Two independent cues beat one. Geometry is the backbone. Laplacian variance and Michelson contrast give a bounded correction and a trust signal.
3. Optics is for confidence as much as for distance. A frame that is too soft or too flat should lower the weight of the 3D branch, not produce a confident wrong number.
4. One calibrated constant, not a tire database. K is fitted with the same np.linalg.lstsq machinery as in Part 1.
5. Fusion needs +/-15%, not +/-1%. A passive monocular estimator delivers that in well under 1 ms.

Next: Terminator Vision: Dual-Stage ONNX Cascade for Tread Wear Indicator (TWI) Grid & Heatmap Detection. We'll look at how the 1.6 mm TWI anchor is detected, how the grid and heatmap detectors reach consensus through NumPy local-maxima search, and how the PyQt6 HUD renders it in real time.

---