Introduction

calibration-rs is a Rust library for camera calibration — the process of estimating the internal parameters (intrinsics, lens distortion) and external parameters (pose, rig geometry, hand-eye transforms) of camera systems from observed correspondences between known 3D points and their 2D projections.

API reference

Who This Book Is For

This book targets engineers and researchers working in machine vision, robotics, and 3D reconstruction who want to:

• Calibrate single cameras, stereo rigs, or multi-camera systems
• Perform hand-eye calibration for cameras mounted on robot arms
• Calibrate laser triangulation devices (camera + laser plane)
• Understand the mathematical foundations behind calibration algorithms

We assume familiarity with linear algebra (SVD, eigendecomposition, least squares) and basic projective geometry (homogeneous coordinates, projection matrices). Lie group theory (SO(3), SE(3)) is introduced when needed.

What calibration-rs Provides

The library covers the full calibration pipeline:

1. Linear initialization — closed-form solvers (Zhang's method, DLT, PnP, Tsai-Lenz hand-eye, epipolar geometry) that produce approximate parameter estimates
2. Non-linear refinement — Levenberg-Marquardt bundle adjustment that minimizes reprojection error to sub-pixel accuracy
3. Session framework — a high-level API with step functions, configuration, and JSON checkpointing

Workspace Structure

calibration-rs is organized as a 5-crate Rust workspace with a layered architecture:

vision-calibration (facade)
    → vision-calibration-pipeline (sessions, workflows)
        → vision-calibration-optim (non-linear refinement)
        → vision-calibration-linear (initialization)
            → vision-calibration-core (primitives, camera models, RANSAC)

The key dependency rule: the linear and optimization crates are peers — they both depend on vision-calibration-core but not on each other. This separation keeps initialization algorithms independent of the optimization backend.

Relation to OpenCV

Readers familiar with OpenCV's calibration module will find analogous functionality throughout:

calibration-rs differs from OpenCV in several ways: it is written in pure Rust, uses a composable camera model with generic type parameters, provides a backend-agnostic optimization IR, and offers a session framework with JSON checkpointing for production workflows.

Book Organization

The book is structured in seven parts:

• Part I: Camera Model — the composable projection pipeline (pinhole, distortion, sensor tilt, intrinsics)
• Part II: Geometric Primitives — rigid transforms, RANSAC, robust loss functions
• Part III: Linear Initialization — all closed-form solvers with full mathematical derivations
• Part IV: Non-Linear Optimization — Levenberg-Marquardt, manifold constraints, autodiff, the IR architecture
• Part V: Calibration Workflows — end-to-end pipelines for each problem type, with both synthetic and real data examples
• Part VI: Session Framework — the high-level CalibrationSession API
• Part VII: Extending the Library — adding new problems, backends, and pipeline types

Each algorithm chapter includes a formal problem statement, objective function, assumptions, and a full derivation leading to the implementation.

Quickstart

This chapter walks through a minimal camera calibration using synthetic data. By the end you will have estimated camera intrinsics and lens distortion from simulated observations of a planar calibration board.

Add the Dependency

[dependencies]
vision-calibration = "0.2"

Minimal Example

The following program generates synthetic calibration data, runs the two-step calibration pipeline (linear initialization followed by non-linear refinement), and prints the results:

use anyhow::Result;
use vision_calibration::planar_intrinsics::{step_init, step_optimize};
use vision_calibration::prelude::*;
use vision_calibration::synthetic::planar;

fn main() -> Result<()> {
    // 1. Define ground truth camera
    let k_gt = FxFyCxCySkew {
        fx: 800.0, fy: 780.0, cx: 640.0, cy: 360.0, skew: 0.0,
    };
    let dist_gt = BrownConrady5 {
        k1: 0.05, k2: -0.02, k3: 0.0,
        p1: 0.001, p2: -0.001, iters: 8,
    };
    let camera = make_pinhole_camera(k_gt, dist_gt);

    // 2. Generate synthetic observations
    let board = planar::grid_points(8, 6, 0.04); // 8×6 grid, 40 mm squares
    let poses = planar::poses_yaw_y_z(6, -0.2, 0.08, 0.5, 0.05);
    let views = planar::project_views_all(&camera, &board, &poses)?;

    // 3. Create session and set input
    let dataset = PlanarDataset::new(
        views.into_iter().map(View::without_meta).collect()
    )?;
    let mut session = CalibrationSession::<PlanarIntrinsicsProblem>::new();
    session.set_input(dataset)?;

    // 4. Run calibration
    step_init(&mut session, None)?;      // Linear initialization
    step_optimize(&mut session, None)?;   // Non-linear refinement

    // 5. Inspect results
    let export = session.export()?;
    let k = export.params.intrinsics();
    println!("fx={:.1}, fy={:.1}, cx={:.1}, cy={:.1}",
             k.fx, k.fy, k.cx, k.cy);
    println!("Mean reprojection error: {:.4} px",
             export.mean_reproj_error);

    Ok(())
}

What Happens in Each Step

Step 1: step_init — Linear Initialization

1. Computes a homography from each view's 2D-3D correspondences using DLT (direct linear transform)
2. Estimates intrinsics $K$ from the homographies using Zhang's method
3. Estimates lens distortion coefficients from homography residuals
4. Iterates steps 2-3 to refine the joint $K$ + distortion estimate
5. Recovers each view's camera pose from its homography and $K$

After this step, intrinsics are typically within 10-40% of the true values — good enough to initialize non-linear optimization.

Step 2: step_optimize — Bundle Adjustment

Runs Levenberg-Marquardt optimization minimizing the total reprojection error:

$$\underset{K,\mathbf{d},\{T_i\}}{\min }\sum _{i=1}^N\sum _{j=1}^{M_i}{\Vert \pi (K,\mathbf{d},T_i,{\mathbf{P}}_j)-{\mathbf{p}}_{ij}\Vert }^2$$

where $\pi$ is the camera projection function, $T_i$ is the pose of view $i$, $\mathbf{d}$ are distortion coefficients, ${\mathbf{P}}_j$ are known 3D board points, and ${\mathbf{p}}_{ij}$ are observed pixel coordinates.

After optimization, expect <2% error on intrinsics and <1 px mean reprojection error.

Running the Built-in Example

The library ships with this example (and several others):

cargo run -p vision-calibration --example planar_synthetic

Alternative: One-Line Pipeline

If you don't need to inspect intermediate results:

#![allow(unused)]
fn main() {
use vision_calibration::planar_intrinsics::run_calibration;

run_calibration(&mut session)?;
}

Next Steps

• Architecture Overview — understand the crate structure
• Composable Camera Pipeline — the projection model in detail
• Planar Intrinsics Calibration — the full mathematical treatment

Architecture Overview

calibration-rs is organized as a layered workspace of Rust crates plus Python bindings. This chapter explains the dependency structure, data flow, and key design patterns.

Crate Dependency Graph

vision-calibration          (facade: re-exports everything)
    │
    └── vision-calibration-pipeline    (sessions, workflows, step functions)
            │
            ├── vision-calibration-optim   (non-linear optimization)
            │       │
            │       └── vision-calibration-core
            │
            └── vision-calibration-linear  (closed-form initialization)
                    │
                    └── vision-calibration-core

Key rule: vision-calibration-linear and vision-calibration-optim are peers. They both depend on vision-calibration-core but never on each other. This keeps initialization algorithms free of optimization dependencies and vice versa.

Crate Responsibilities

vision-calibration-core

The foundation layer providing:

• Math types — nalgebra-based aliases: Pt2, Pt3, Vec3, Mat3, Iso3, Real (= f64)
• Camera models — composable trait-based pipeline: ProjectionModel, DistortionModel, SensorModel, IntrinsicsModel
• Observation types — CorrespondenceView (2D-3D point pairs), View<Meta>, PlanarDataset, RigDataset
• RANSAC engine — generic Estimator trait with configurable options
• Synthetic data utilities — grid generation, pose sampling, projection
• Reprojection error computation — single-camera and multi-camera rig

vision-calibration-linear

Closed-form initialization solvers. Each produces an approximate estimate suitable for seeding non-linear optimization:

vision-calibration-optim

Non-linear refinement with a backend-agnostic architecture:

1. IR (Intermediate Representation) — ProblemIR with ParamBlock and ResidualBlock types that describe optimization problems independently of any solver
2. Factors — generic residual functions parameterized over RealField for automatic differentiation
3. Backends — currently TinySolverBackend (Levenberg-Marquardt with sparse linear solvers)
4. Problem builders — domain-specific functions that construct IR from calibration data

vision-calibration-pipeline

The session framework providing production-ready workflows:

• CalibrationSession<P: ProblemType> — generic state container with config, input, state, output, exports
• Step functions — free functions operating on &mut CalibrationSession<P> (e.g., step_init, step_optimize)
• Pipeline functions — convenience wrappers chaining all steps
• JSON checkpointing — full serialization for session persistence
• Six problem types: PlanarIntrinsicsProblem, ScheimpflugIntrinsicsProblem, SingleCamHandeyeProblem, RigExtrinsicsProblem, RigHandeyeProblem, LaserlineDeviceProblem

vision-calibration

Unified facade crate that re-exports everything through a clean module hierarchy:

#![allow(unused)]
fn main() {
use vision_calibration::prelude::*;            // Minimal planar hello-world imports
use vision_calibration::planar_intrinsics::*;  // Planar workflow
use vision_calibration::core::*;               // Math types
use vision_calibration::linear::*;             // Linear solvers
use vision_calibration::optim::*;              // Optimization
}

vision-calibration-py

Python bindings crate (maturin/PyO3) exposing high-level session workflows.

Data Flow

Every calibration workflow follows the same pattern:

Observations (2D-3D correspondences)
    │
    ▼
Linear Initialization (vision-calibration-linear)
    │  Closed-form solvers: ~5-40% accuracy
    ▼
Non-Linear Refinement (vision-calibration-optim)
    │  Levenberg-Marquardt: <2% accuracy, <1 px reprojection
    ▼
Calibrated Parameters (K, distortion, poses, ...)

The session framework wraps this flow with configuration, state tracking, and checkpointing:

CalibrationSession::new()
    │
    ▼
session.set_input(dataset)
    │
    ▼
step_init(&mut session)      ← linear initialization
    │
    ▼
step_optimize(&mut session)  ← non-linear refinement
    │
    ▼
session.export()             ← calibrated parameters

Design Patterns

Composable Camera Model

The camera projection pipeline is built from four independent traits composed via generics:

pixel = K(sensor(distortion(projection(direction))))

Each stage can be mixed and matched. For example, a standard camera uses Pinhole + BrownConrady5 + IdentitySensor + FxFyCxCySkew, while a laser profiler might use Pinhole + BrownConrady5 + ScheimpflugParams + FxFyCxCySkew.

Backend-Agnostic Optimization

Problems are defined as an intermediate representation (IR) that is independent of any specific solver. The IR is then compiled to a solver-specific form:

Problem Builder  →  ProblemIR  →  Backend.compile()  →  Backend.solve()
                  (generic)       (solver-specific)

This allows swapping the optimization backend without changing problem definitions.

Step Functions

Calibration workflows are decomposed into discrete steps implemented as free functions. This allows:

• Inspection of intermediate state between steps
• Resumption from any point (via JSON checkpointing)
• Customization of per-step options
• Composition of steps from different problem types

Composable Camera Pipeline

A camera model maps 3D points in the camera frame to 2D pixel coordinates. In calibration-rs, this mapping is decomposed into four composable stages, each implemented as a trait:

$${\mathbf{p}}_{\text{pixel}}=K\circ S\circ D\circ \Pi (\mathbf{d})$$

where:

• $\Pi$ — projection: maps a 3D direction to normalized coordinates on the image plane
• $D$ — distortion: warps normalized coordinates to model lens imperfections
• $S$ — sensor: applies a homography for tilted sensor planes (identity for standard cameras)
• $K$ — intrinsics: scales and translates to pixel coordinates

The Camera Struct

The camera is a generic struct parameterized over the four model traits:

#![allow(unused)]
fn main() {
pub struct Camera<S, P, D, Sm, K>
where
    S: RealField,
    P: ProjectionModel<S>,
    D: DistortionModel<S>,
    Sm: SensorModel<S>,
    K: IntrinsicsModel<S>,
{
    pub proj: P,
    pub dist: D,
    pub sensor: Sm,
    pub k: K,
}
}

The scalar type S is generic over RealField, enabling the same camera to work with f64 for evaluation and with dual numbers for automatic differentiation.

Forward Projection

Given a 3D point ${\mathbf{P}}_c=[X,Y,Z{]}^T$ in the camera frame, projection to pixel coordinates proceeds through the four stages:

#![allow(unused)]
fn main() {
pub fn project_point_c(&self, p_c: &Vector3<S>) -> Option<Point2<S>>
}

1. Projection: Compute normalized coordinates ${\mathbf{n}}_u=\Pi ({\mathbf{P}}_c)$. For pinhole: ${\mathbf{n}}_u=[X/Z,Y/Z{]}^T$. Returns None if the point is behind the camera ($Z\le 0$).

2. Distortion: Apply lens distortion ${\mathbf{n}}_d=D({\mathbf{n}}_u)$. Warps normalized coordinates according to the distortion model (e.g., Brown-Conrady radial + tangential).

3. Sensor transform: Apply sensor model $\mathbf{s}=S({\mathbf{n}}_d)$. For standard cameras this is identity; for tilted sensors it applies a homography.

4. Intrinsics: Map to pixels $\mathbf{p}=K(\mathbf{s})$. Applies focal lengths, principal point, and optional skew.

Inverse (Back-Projection)

The inverse maps a pixel to a ray in the camera frame:

#![allow(unused)]
fn main() {
pub fn backproject_pixel(&self, px: &Point2<S>) -> Ray<S>
}

1. $\mathbf{s}=K^{-1}(\mathbf{p})$ — pixel to sensor coordinates
2. ${\mathbf{n}}_d=S^{-1}(\mathbf{s})$ — sensor to distorted normalized
3. ${\mathbf{n}}_u=D^{-1}({\mathbf{n}}_d)$ — undistort (iterative)
4. $\mathbf{d}={\Pi }^{-1}({\mathbf{n}}_u)$ — normalized to 3D direction

The result is a point on the $Z=1$ plane in the camera frame, defining a ray from the camera center.

Common Type Aliases

For the most common configuration (pinhole + Brown-Conrady + no tilt):

#![allow(unused)]
fn main() {
type PinholeCamera = Camera<f64, Pinhole, BrownConrady5<f64>,
                            IdentitySensor, FxFyCxCySkew<f64>>;
}

The make_pinhole_camera(k, dist) convenience function constructs this type.

Projecting World Points

To project a world point through a posed camera, first transform from world to camera frame using the extrinsic pose $T_{C,W}\in \text{SE}(3)$:

$${\mathbf{P}}_c=T_{C,W}\cdot {\mathbf{P}}_w$$

Then apply the camera projection:

$$\mathbf{p}=\text{camera.project\_point\_c}({\mathbf{P}}_c)$$

Trait Composition

Each stage is an independent trait. This design allows mixing and matching:

Note: ScheimpflugParams stores the tilt angles and has a compile() method that produces a HomographySensor suitable for the Camera struct. See the Sensor Models chapter.

The subsequent sections detail each stage.

Pinhole Projection

The pinhole model is the simplest and most widely used camera projection. It models a camera as an ideal point through which all light rays pass — no lens, no aperture, just a geometric center of projection.

Mathematical Definition

Given a 3D point ${\mathbf{P}}_c=[X,Y,Z{]}^T$ in the camera frame, the pinhole projection maps it to normalized image coordinates:

$$\mathbf{n}=\Pi ({\mathbf{P}}_c)=\left[\begin{array}{c}X/Z\\ Y/Z\end{array}\right]$$

This is perspective division: the point is projected onto the $Z=1$ plane. The projection is defined only for points in front of the camera ($Z>0$).

The inverse (back-projection) maps normalized coordinates to a 3D direction:

$${\Pi }^{-1}(\mathbf{n})=\left[\begin{array}{c}{n}_x\\ n_y\\ 1\end{array}\right]$$

This defines a ray from the camera center through the image point.

Camera Frame Convention

calibration-rs uses a right-handed camera frame:

• $X$: right
• $Y$: down
• $Z$: forward (optical axis)

This matches the convention used by OpenCV and most computer vision libraries. A point with positive $Z$ is in front of the camera; negative $Z$ is behind it.

The ProjectionModel Trait

#![allow(unused)]
fn main() {
pub trait ProjectionModel<S: RealField> {
    fn project_dir(&self, dir_c: &Vector3<S>) -> Option<Point2<S>>;
    fn unproject_dir(&self, n: &Point2<S>) -> Vector3<S>;
}
}

The Pinhole struct implements this trait:

#![allow(unused)]
fn main() {
pub struct Pinhole;

impl<S: RealField> ProjectionModel<S> for Pinhole {
    fn project_dir(&self, dir_c: &Vector3<S>) -> Option<Point2<S>> {
        if dir_c.z > S::zero() {
            Some(Point2::new(
                dir_c.x.clone() / dir_c.z.clone(),
                dir_c.y.clone() / dir_c.z.clone(),
            ))
        } else {
            None
        }
    }

    fn unproject_dir(&self, n: &Point2<S>) -> Vector3<S> {
        Vector3::new(n.x.clone(), n.y.clone(), S::one())
    }
}
}

Note the use of .clone() — this is required for compatibility with dual numbers used in automatic differentiation (see Autodiff and Generic Residual Functions).

Limitations

The pinhole model assumes:

• No lens distortion — addressed by the distortion stage
• No sensor tilt — addressed by the sensor stage
• Central projection — all rays pass through a single point

For cameras with significant distortion (wide-angle, fisheye), the distortion model corrects for lens effects after pinhole projection to normalized coordinates.

Brown-Conrady Distortion

Real lenses introduce geometric distortion: straight lines in the world project as curves in the image. The Brown-Conrady model captures the two dominant distortion effects — radial and tangential — as polynomial functions of the distance from the optical axis.

Distortion Model

Given undistorted normalized coordinates ${\mathbf{n}}_u=[x,y{]}^T$, the distorted coordinates ${\mathbf{n}}_d=[x_d,y_d{]}^T$ are:

$$r^2=x^2+y^2$$

Radial distortion (barrel/pincushion):

$$\alpha =1+k_1{r}^2+k_2{r}^4+k_3{r}^6$$

Tangential distortion (decentering):

$${\delta }_x=2p_{1}xy+p_2(r^2+2x^2)$$ $${\delta }_y=p_1(r^2+2y^2)+2p_{2}xy$$

Combined:

$$x_d=x\cdot \alpha +{\delta }_x$$ $$y_d=y\cdot \alpha +{\delta }_y$$

The model has five parameters: radial coefficients $(k_1,k_2,k_3)$ and tangential coefficients $(p_1,p_2)$.

OpenCV equivalence: This is identical to OpenCV's distortPoints with the 5-parameter model (k1, k2, p1, p2, k3). Note OpenCV's parameter ordering differs from the mathematical ordering.

Physical Interpretation

• $k_1>0$: barrel distortion (lines bow outward from center)
• $k_1<0$: pincushion distortion (lines bow inward)
• $k_2,k_3$: higher-order radial corrections
• $p_1,p_2$: tangential distortion from imperfect lens-sensor alignment (lens elements not perfectly centered)

Typical values for industrial cameras: $|k_1|\sim 0.01\text{-}0.3$, $|k_2|\sim 0.001\text{-}0.1$, $|p_1|,|p_2|\sim 0.0001\text{-}0.01$.

When to Fix $k_3$

The sixth-order radial term $k_3{r}^6$ is only significant for wide-angle lenses where $r$ reaches large values. For typical industrial cameras with moderate field of view:

• $k_3$ is poorly constrained by the data and often absorbs noise
• Estimating $k_3$ can cause overfitting, leading to worse extrapolation outside the calibration region
• The library defaults to fix_k3: true in most problem configurations

Recommendation: Only estimate $k_3$ with high-quality data covering the full image area, or for lenses with field of view > 90°.

Undistortion (Inverse)

Given distorted coordinates ${\mathbf{n}}_d$, recovering undistorted coordinates ${\mathbf{n}}_u$ requires inverting the distortion model. Since the forward model is a polynomial without a closed-form inverse, calibration-rs uses iterative fixed-point refinement:

$${\mathbf{n}}_u^{(0)}={\mathbf{n}}_d$$

$${\mathbf{n}}_u^{(k+1)}={\mathbf{n}}_d-\left(D({\mathbf{n}}_u^{(k)})-{\mathbf{n}}_u^{(k)}\right)$$

This rearranges the distortion equation to isolate ${\mathbf{n}}_u$ and iterates to convergence. The default is 8 iterations, which provides accuracy well below sensor noise for typical distortion magnitudes.

Convergence

The fixed-point iteration converges when the Jacobian of the distortion residual has spectral radius less than 1, which holds for physically realistic distortion values (small $k_i$, $p_i$ relative to $r$). For extreme distortion, more iterations may be needed via the iters parameter.

OpenCV equivalence: cv::undistortPoints performs the same iterative inversion.

The DistortionModel Trait

#![allow(unused)]
fn main() {
pub trait DistortionModel<S: RealField> {
    fn distort(&self, n_undist: &Point2<S>) -> Point2<S>;
    fn undistort(&self, n_dist: &Point2<S>) -> Point2<S>;
}
}

The BrownConrady5<S> struct:

#![allow(unused)]
fn main() {
pub struct BrownConrady5<S: RealField> {
    pub k1: S, pub k2: S, pub k3: S,
    pub p1: S, pub p2: S,
    pub iters: u32,   // undistortion iterations (default: 8)
}
}

NoDistortion is the identity implementation for distortion-free cameras.

Distortion in the Projection Pipeline

Distortion operates in normalized image coordinates (after pinhole projection, before intrinsics):

$$\text{pixel}=K(\underset{\text{distorted normalized}}{\underbrace{D(\Pi ({\mathbf{P}}_c))}})$$

This is the standard convention: distortion is defined on the $Z=1$ plane, not in pixel space. The advantage is that distortion parameters are independent of image resolution and focal length.

Intrinsics Matrix

The intrinsics matrix $K$ maps from the sensor coordinate system (normalized or sensor-transformed coordinates) to pixel coordinates. It encodes the camera's internal geometric properties: focal length, principal point, and optional skew.

Definition

$$K=\left[\begin{array}{ccc}{f}_x& s& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]$$

where:

• $f_x,f_y$ — focal lengths in pixels (horizontal and vertical). These combine the physical focal length with the pixel pitch: $f_x=f/{\Delta }_x$ where $f$ is the focal length in mm and ${\Delta }_x$ is the pixel width in mm.
• $c_x,c_y$ — principal point in pixels. The projection of the optical axis onto the image plane. Typically near the image center but not exactly at it.
• $s$ — skew coefficient. Non-zero only if pixel axes are not perpendicular. Effectively zero for all modern sensors; calibration-rs defaults to zero_skew: true in most configurations.

Forward and Inverse Transform

Sensor to pixel (forward):

$$\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}{f}_x\cdot s_x+s\cdot s_y+c_x\\ f_y\cdot s_y+c_y\end{array}\right]$$

where $[s_x,s_y{]}^T$ are sensor coordinates (output of the distortion + sensor stages).

Pixel to sensor (inverse):

$$s_y=\frac{v-c_y}{f_y},s_x=\frac{u-c_x-s\cdot s_y}{f_x}$$

The IntrinsicsModel Trait

#![allow(unused)]
fn main() {
pub trait IntrinsicsModel<S: RealField> {
    fn sensor_to_pixel(&self, sensor: &Point2<S>) -> Point2<S>;
    fn pixel_to_sensor(&self, pixel: &Point2<S>) -> Point2<S>;
}
}

The FxFyCxCySkew<S> struct implements this trait:

#![allow(unused)]
fn main() {
pub struct FxFyCxCySkew<S: RealField> {
    pub fx: S, pub fy: S,
    pub cx: S, pub cy: S,
    pub skew: S,
}
}

Coordinate Utilities

The coordinate_utils module provides convenience functions that combine intrinsics with distortion:

• pixel_to_normalized(pixel, K) — apply $K^{-1}$ to get normalized coordinates on the $Z=1$ plane
• normalized_to_pixel(normalized, K) — apply $K$ to get pixel coordinates
• undistort_pixel(pixel, K, distortion) — pixel → normalized → undistort → return normalized
• distort_to_pixel(normalized, K, distortion) — distort → apply $K$ → return pixel

These functions are used extensively in the linear initialization algorithms, which need to convert between pixel and normalized coordinate spaces.

Aspect Ratio

For most cameras, $f_x\ne f_y$ because pixels are not perfectly square. The ratio $f_x/f_y$ reflects the pixel aspect ratio. Typical industrial cameras have aspect ratios very close to 1.0 (within 1-3%).

OpenCV equivalence: cv::initCameraMatrix2D provides initial intrinsics estimates. The $K$ matrix format is identical to OpenCV's cameraMatrix.

Sensor Models and Scheimpflug Tilt

In standard cameras, the sensor plane is perpendicular to the optical axis. Some industrial cameras — particularly laser triangulation profilers — deliberately tilt the sensor to achieve a larger depth of field along the laser plane (the Scheimpflug condition). The sensor model stage accounts for this tilt.

The SensorModel Trait

#![allow(unused)]
fn main() {
pub trait SensorModel<S: RealField> {
    fn normalized_to_sensor(&self, n: &Point2<S>) -> Point2<S>;
    fn sensor_to_normalized(&self, s: &Point2<S>) -> Point2<S>;
}
}

The sensor model sits between distortion and intrinsics in the pipeline:

$$\mathbf{p}=K(\underset{\text{sensor coords}}{\underbrace{S(D(\Pi ({\mathbf{P}}_c)))}})$$

IdentitySensor

For standard cameras with untilted sensors:

$$S(\mathbf{n})=\mathbf{n},S^{-1}(\mathbf{s})=\mathbf{s}$$

This is the default and most common case.

Scheimpflug Tilt Model

The Scheimpflug condition states that when the lens plane, image plane, and object plane intersect along a common line, the entire object plane is in sharp focus. Laser profilers exploit this by tilting the sensor to keep the laser line in focus.

Tilt Projection Matrix

The tilt is parameterized by two rotation angles:

• ${\tau }_x$ — tilt around the horizontal (X) axis
• ${\tau }_y$ — tilt around the vertical (Y) axis

The tilt homography is constructed as:

$$R_x=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\tau }_x& -\sin {\tau }_x\\ 0& \sin {\tau }_x& \cos {\tau }_x\end{array}\right],R_y=\left[\begin{array}{ccc}\cos {\tau }_y& 0& \sin {\tau }_y\\ 0& 1& 0\\ -\sin {\tau }_y& 0& \cos {\tau }_y\end{array}\right]$$

$$R=R_y\cdot R_x$$

The projection onto the tilted sensor plane normalizes by the third row of $R$ applied to the point in homogeneous coordinates:

$$H_{\text{tilt}}={\text{proj}}_z\cdot R$$

where ${\text{proj}}_z$ normalizes by the $Z$-component (perspective division through the tilted plane). In practice, the first two rows of $R$ divided by the third row give a $2\times 2$ rational transform in the image plane.

OpenCV equivalence: This matches OpenCV's computeTiltProjectionMatrix used in the 14-parameter rational model (CALIB_TILTED_MODEL).

ScheimpflugParams

#![allow(unused)]
fn main() {
pub struct ScheimpflugParams {
    pub tilt_x: f64,   // tau_x in radians
    pub tilt_y: f64,   // tau_y in radians
}
}

The compile() method generates a HomographySensor containing the $3\times 3$ homography and its precomputed inverse. ScheimpflugParams::default() returns zero tilt (identity sensor).

HomographySensor

The general homography sensor applies an arbitrary $3\times 3$ projective transform:

$$\mathbf{s}=\text{dehomogenize}(H\cdot [\mathbf{n},1{]}^T)$$

$$\mathbf{n}=\text{dehomogenize}(H^{-1}\cdot [\mathbf{s},1{]}^T)$$

This is used as the runtime representation compiled from ScheimpflugParams.

When to Use Scheimpflug

• Laser triangulation profilers: The laser plane is at an angle to the camera axis. Tilting the sensor brings the entire laser line into focus.
• Machine vision with tilted object planes: When the depth of field is insufficient to keep the entire object in focus.

For standard cameras (no tilt), use IdentitySensor (the default).

Optimization of Tilt Parameters

In the non-linear optimization stage, Scheimpflug parameters $({\tau }_x,{\tau }_y)$ are treated as additional optimization variables. The tilt_projection_matrix_generic<T>() function in the factor module computes the tilt homography using dual numbers for automatic differentiation, enabling joint optimization of tilt with intrinsics, distortion, and poses.

Serialization and Runtime-Dynamic Types

The generic Camera<S, P, D, Sm, K> type provides compile-time composition of camera stages. However, for JSON serialization (session checkpointing, data exchange) and runtime-dynamic camera construction, calibration-rs provides an enum-based parameter system.

CameraParams

The CameraParams struct holds camera parameters as serializable enums:

#![allow(unused)]
fn main() {
pub struct CameraParams {
    pub projection: ProjectionParams,
    pub distortion: DistortionParams,
    pub sensor: SensorParams,
    pub intrinsics: IntrinsicsParams,
}
}

Each component is an enum of supported variants:

#![allow(unused)]
fn main() {
pub enum ProjectionParams { Pinhole }

pub enum DistortionParams {
    None,
    BrownConrady5 { k1, k2, k3, p1, p2 },
}

pub enum SensorParams {
    Identity,
    Scheimpflug { tilt_x, tilt_y },
}

pub enum IntrinsicsParams {
    FxFyCxCySkew { fx, fy, cx, cy, skew },
}
}

Building a Camera from Parameters

The build() method constructs a concrete Camera from serialized parameters:

#![allow(unused)]
fn main() {
let params = CameraParams { /* ... */ };
let camera = params.build();
}

This is used internally by the session framework to reconstruct cameras from JSON-checkpointed state.

Convenience Accessors

CameraParams provides typed accessors:

#![allow(unused)]
fn main() {
params.intrinsics()   // → &FxFyCxCySkew<f64>
params.distortion()   // → &BrownConrady5<f64> or zero distortion
}

Use in Session Framework

CameraParams appears throughout the pipeline:

• Export records store calibrated parameters as CameraParams
• Session JSON serializes intermediate and final camera parameters
• Config types reference CameraParams for initial guesses

This provides a stable, version-safe serialization format decoupled from the generic type system.

Rigid Transformations and SE(3)

Camera calibration is fundamentally about estimating rigid transformations — rotations and translations that relate different coordinate frames: camera, world, robot gripper, calibration target, rig, etc. This chapter covers the rotation representations and transformation algebra used throughout calibration-rs.

Rotation Representations

Rotation Matrix $R\in \text{SO}(3)$

A $3\times 3$ orthogonal matrix with determinant $+1$:

$$R^{T}R=I,\det (R)=1$$

The set of all such matrices forms the special orthogonal group $\text{SO}(3)$. It has 3 degrees of freedom (9 entries minus 6 constraints).

Unit Quaternion

A 4-element vector $\mathbf{q}=[q_x,q_y,q_z,q_w{]}^T$ with unit norm $\Vert \mathbf{q}\Vert =1$. The rotation of a vector $\mathbf{v}$ is:

$${\mathbf{v}}^{\prime }=\mathbf{q}\otimes \mathbf{v}\otimes {\mathbf{q}}^{\ast }$$

where $\otimes$ is quaternion multiplication and ${\mathbf{q}}^{\ast }$ is the conjugate.

Quaternions are preferred for optimization because they avoid gimbal lock and have simpler interpolation properties than Euler angles. calibration-rs (via nalgebra) stores quaternions as [i, j, k, w] internally, which maps to $[q_x,q_y,q_z,q_w]$.

Axis-Angle

A rotation of angle $\theta$ around unit axis $\hat{\mathbf{a}}$ is represented as:

$$\mathbf{\omega }=\theta \hat{\mathbf{a}}\in {\mathbb{R}}^3$$

with $\theta =\Vert \mathbf{\omega }\Vert$ and $\hat{\mathbf{a}}=\mathbf{\omega }/\theta$.

The Rodrigues formula converts axis-angle to rotation matrix:

$$R=I+\frac{\sin \theta }{\theta }[\mathbf{\omega }{]}_{\times }+\frac{1-\cos \theta }{{\theta }^2}[\mathbf{\omega }{]}_{\times }^2$$

where $[\mathbf{\omega }{]}_{\times }$ is the skew-symmetric matrix of $\mathbf{\omega }$:

$$[\mathbf{\omega }{]}_{\times }=\left[\begin{array}{ccc}0& -{\omega }_z& {\omega }_y\\ {\omega }_z& 0& -{\omega }_x\\ -{\omega }_y& {\omega }_x& 0\end{array}\right]$$

Rigid Transformations: SE(3)

A rigid body transformation combines a rotation $R$ and translation $\mathbf{t}$:

$$T=\left[\begin{array}{cc}R& \mathbf{t}\\ 0^T& 1\end{array}\right]\in \text{SE}(3)$$

The special Euclidean group $\text{SE}(3)$ has 6 degrees of freedom (3 rotation + 3 translation).

Applying a Transform

A point $\mathbf{P}$ transformed by $T$:

$${\mathbf{P}}^{\prime }=R\mathbf{P}+\mathbf{t}$$

or in homogeneous coordinates:

$$\left[\begin{array}{c}{\mathbf{P}}^{\prime }\\ 1\end{array}\right]=T\left[\begin{array}{c}\mathbf{P}\\ 1\end{array}\right]$$

Composition

Transforms compose by matrix multiplication:

$$T_{A,C}=T_{A,B}\cdot T_{B,C}$$

Inverse

$$T^{-1}=\left[\begin{array}{cc}{R}^T& -R^T\mathbf{t}\\ 0^T& 1\end{array}\right]$$

Naming Convention

calibration-rs uses the notation $T_{A,B}$ to mean "the transform from frame $B$ to frame $A$":

$${\mathbf{P}}_A=T_{A,B}\cdot {\mathbf{P}}_B$$

Common transforms:

nalgebra Type: Iso3

calibration-rs uses nalgebra's Isometry3<f64> (aliased as Iso3) for SE(3) transforms:

#![allow(unused)]
fn main() {
pub type Iso3 = nalgebra::Isometry3<f64>;
}

This stores a UnitQuaternion (rotation) and a Translation3 (translation) separately, avoiding the redundancy of a full $4\times 4$ matrix.

SE(3) Storage for Optimization

In the optimization IR, SE(3) parameters are stored as a 7-element vector:

$$[q_x,q_y,q_z,q_w,t_x,t_y,t_z]$$

The quaternion $[q_x,q_y,q_z,q_w]$ is unit-constrained. The optimizer uses the SE(3) manifold (see Manifold Optimization) to maintain this constraint during parameter updates.

The Exponential Map

The Lie algebra $\mathfrak{se}(3)$ provides a 6-dimensional tangent space at the identity. A tangent vector $\mathbf{\xi }=[\mathbf{\omega },\mathbf{v}{]}^T\in {\mathbb{R}}^6$ maps to an SE(3) element via the exponential map:

$$T=\exp ({\mathbf{\xi }}^{\wedge })=\left[\begin{array}{cc}\exp ([\mathbf{\omega }{]}_{\times })& V\mathbf{v}\\ 0^T& 1\end{array}\right]$$

where $V$ is the left Jacobian of SO(3):

$$V=I+\frac{1-\cos \theta }{{\theta }^2}[\mathbf{\omega }{]}_{\times }+\frac{\theta -\sin \theta }{{\theta }^3}[\mathbf{\omega }{]}_{\times }^2$$

The exponential map is central to manifold optimization: parameter updates are computed in the tangent space and then retracted to the manifold (see Manifold Optimization).

RANSAC

RANSAC (Random Sample Consensus) is a robust estimation method that fits a model to data containing outliers. Unlike least-squares methods that use all data points and can be corrupted by even a single outlier, RANSAC repeatedly samples minimal subsets, fits a model to each, and selects the model with the most support (inliers).

Algorithm

Given $N$ data points and a model requiring $m$ points to fit:

1. Repeat for up to $k_{\max }$ iterations:

   1. Sample $m$ random data points
   2. Fit a model to the $m$-point sample (skip if degenerate)
   3. Score: count inliers — points with residual $<ϵ$
   4. If this model has more inliers than the current best, update the best
   5. Optionally refit the model on all inliers for a tighter fit
   6. Update the iteration bound dynamically based on the current inlier ratio

2. Return the best model and its inlier set

Dynamic Iteration Bound

After finding a model with inlier ratio $w=n_{\text{inlier}}/N$, the number of iterations needed to find an all-inlier sample with probability $p$ is:

$$k=\frac{\log (1-p)}{\log (1-w^m)}$$

where $m$ is the minimum sample size. As more inliers are found, $w$ increases and $k$ decreases, allowing early termination. The iteration count is clamped to $[k_{\text{current}},k_{\max }]$ — it can only decrease.

Example: For $m=4$ (homography), 50% inliers, and 99% confidence: $k=\frac{\log (0.01)}{\log (1-0.5^4)}=\frac{-4.6}{-0.065}\approx 72$ iterations.

The Estimator Trait

calibration-rs provides a generic RANSAC engine that works with any model implementing the Estimator trait:

#![allow(unused)]
fn main() {
pub trait Estimator {
    type Datum;
    type Model;

    const MIN_SAMPLES: usize;

    fn fit(data: &[Self::Datum], sample_indices: &[usize])
        -> Option<Self::Model>;

    fn residual(model: &Self::Model, datum: &Self::Datum) -> f64;

    fn is_degenerate(data: &[Self::Datum], sample_indices: &[usize])
        -> bool { false }

    fn refit(data: &[Self::Datum], inliers: &[usize])
        -> Option<Self::Model> {
        Self::fit(data, inliers)
    }
}
}

Configuration

#![allow(unused)]
fn main() {
pub struct RansacOptions {
    pub max_iters: usize,      // Maximum iterations
    pub thresh: f64,           // Inlier distance threshold
    pub min_inliers: usize,    // Minimum consensus set size
    pub confidence: f64,       // Desired probability (0-1) for dynamic bound
    pub seed: u64,             // RNG seed for reproducibility
    pub refit_on_inliers: bool, // Refit model on full inlier set
}
}

Choosing the threshold $ϵ$: The threshold should reflect the expected noise level. For pixel-space residuals, $ϵ=2\text{-}5$ pixels is typical. For normalized-coordinate residuals, scale accordingly by the focal length.

Result

#![allow(unused)]
fn main() {
pub struct RansacResult<M> {
    pub success: bool,
    pub model: Option<M>,
    pub inliers: Vec<usize>,
    pub inlier_rms: f64,
    pub iters: usize,
}
}

Instantiated Models

RANSAC is used with several models in calibration-rs:

Deterministic Seeding

All RANSAC calls in calibration-rs use a deterministic seed for the random number generator. This ensures:

• Reproducible results across runs
• Deterministic tests that do not flake
• Debuggable behavior — the same input always produces the same output

The seed can be configured via RansacOptions::seed.

Best-Model Selection

When multiple models achieve the same inlier count, RANSAC selects the model with the lowest inlier RMS (root mean square of inlier residuals). This breaks ties in favor of more accurate models.

Robust Loss Functions

Standard least squares minimizes the sum of squared residuals $\sum r_i^2$. This objective is highly sensitive to outliers: a single point with a large residual can dominate the entire cost function and corrupt the solution. Robust loss functions (M-estimators) reduce the influence of large residuals, making optimization tolerant to outliers in the data.

Problem Setup

In non-linear least squares, we minimize:

$$\underset{\theta }{\min }\sum _{i=1}^N\rho (r_i(\theta ))$$

where $\rho$ is the loss function applied to each residual $r_i$. The standard (non-robust) case uses $\rho (r)=\frac{1}{2}{r}^2$.

Available Loss Functions

calibration-rs provides three robust loss functions, each parameterized by a scale $c>0$ that controls the transition from quadratic (inlier) to robust (outlier) behavior.

Huber Loss

$$\rho (r)=\{\begin{array}{ll}\frac{1}{2}{r}^2& \text{if }|r|\le c\\ c\left(|r|-\frac{c}{2}\right)& \text{if }|r|>c\end{array}$$

Properties:

• Quadratic for small residuals, linear for large residuals
• Continuous first derivative
• Influence function: bounded at $\pm c$ — outliers contribute constant gradient, not growing

When to use: The default robust loss. Good general-purpose choice when you expect a moderate number of outliers.

Cauchy Loss

$$\rho (r)=\frac{c^2}{2}\ln \left(1+\frac{r^2}{c^2}\right)$$

Properties:

• Grows logarithmically for large residuals (slower than linear)
• Smooth everywhere
• Influence function: $\psi (r)=r/(1+r^2/c^2)$ — decreases to zero for large $|r|$, effectively down-weighting far outliers

When to use: When outliers are far from the bulk of the data and should have near-zero influence.

Arctan Loss

$$\rho (r)=c^2\arctan \left(\frac{r^2}{c^2}\right)$$

Properties:

• Bounded: $\rho (r)\to \frac{\pi }{2}{c}^2$ as $|r|\to \infty$
• Influence function approaches zero for large residuals (redescending)

When to use: When very strong outlier rejection is needed. More aggressive than Cauchy but can make convergence harder.

Comparison

Choosing the Scale Parameter $c$

The scale $c$ sets the boundary between "inlier" and "outlier" behavior:

• Too small: Treats good data as outliers, reducing effective sample size
• Too large: Outliers still dominate (approaches standard least squares)
• Rule of thumb: Set $c$ to the expected residual magnitude for good data points. For reprojection residuals, $c=1\text{-}3$ pixels is typical.

Usage in calibration-rs

Robust losses are specified per-residual block in the optimization IR:

#![allow(unused)]
fn main() {
pub enum RobustLoss {
    None,
    Huber { scale: f64 },
    Cauchy { scale: f64 },
    Arctan { scale: f64 },
}
}

Each problem type exposes the loss function as a configuration option:

#![allow(unused)]
fn main() {
session.update_config(|c| {
    c.robust_loss = RobustLoss::Huber { scale: 2.0 };
})?;
}

The backend applies the loss function during residual evaluation, modifying both the cost and the Jacobian.

Iteratively Reweighted Least Squares (IRLS)

Under the hood, robust loss functions are typically implemented via IRLS: each residual is weighted by $w_i={\rho }^{\prime }(r_i)/r_i$, and the weighted least-squares problem is solved iteratively. The Levenberg-Marquardt backend handles this automatically.

Interaction with RANSAC

RANSAC and robust losses address outliers at different stages:

• RANSAC (linear initialization): Binary inlier/outlier classification. Used during model fitting to reject gross outliers before any optimization.
• Robust losses (non-linear refinement): Soft down-weighting. Used during optimization to reduce the influence of moderate outliers that passed RANSAC.

The two approaches are complementary: RANSAC handles gross outliers during initialization, while robust losses handle smaller outliers during refinement.

Why Linear Initialization Matters

Non-linear optimization is the workhorse of camera calibration: it achieves sub-pixel reprojection error by jointly optimizing all parameters. But non-linear optimizers — Levenberg-Marquardt, Gauss-Newton, and their variants — are local methods. They converge to the nearest local minimum, which is only the global minimum if the starting point is sufficiently close.

The role of linear initialization is to provide that starting point.

The Init-Then-Refine Paradigm

Every calibration workflow in calibration-rs follows a two-phase pattern:

1. Linear initialization: Closed-form solvers estimate approximate parameters. These are fast, deterministic, and require no initial guess — but they achieve only ~5-40% accuracy on camera parameters.

2. Non-linear refinement: Levenberg-Marquardt bundle adjustment minimizes reprojection error starting from the linear estimate. This converges to <2% parameter accuracy and <1 px reprojection error.

The linear estimate does not need to be highly accurate. It needs to be in the basin of convergence of the non-linear optimizer — close enough that gradient-based iteration converges to the correct solution rather than a local minimum.

Why Not Just Optimize Directly?

Without initialization, you would need to guess intrinsics ($f_x,f_y,c_x,c_y$), distortion ($k_1,k_2,p_1,p_2$), and all camera poses ($R,\mathbf{t}$ per view). For a typical 20-view calibration:

• 4 intrinsics + 5 distortion + 20 × 6 pose = 129 parameters
• The cost landscape has many local minima
• A random starting point will almost certainly diverge or converge to a wrong solution

Linear initialization eliminates this problem by computing a reasonable estimate from the data structure alone.

Linear Solvers in calibration-rs

Common Mathematical Tools

Most linear solvers share a common toolkit:

• SVD (Singular Value Decomposition): Solves homogeneous systems $A\mathbf{x}=0$ via the right singular vector corresponding to the smallest singular value. This is the core of DLT methods.
• Hartley normalization: Conditions the DLT system for numerical stability by centering and scaling the input points.
• Polynomial root finding: Minimal solvers (P3P, 7-point F, 5-point E) reduce to polynomial equations.

The following chapters derive each algorithm in detail.

Hartley Normalization

DLT algorithms build large linear systems from point coordinates and solve them via SVD. When the coordinates span vastly different scales (e.g., pixel values in the hundreds vs. homogeneous scale factor of 1), the resulting system is ill-conditioned: small perturbations in the data cause large changes in the solution. Hartley normalization is a simple preprocessing step that dramatically improves numerical stability.

The Problem

Consider the homography DLT: we build a $2N\times 9$ matrix $A$ from pixel coordinates $(u,v)$ and solve $A\mathbf{h}=0$. If $(u,v)\sim (640,480)$, the entries of $A$ range from $\sim 1$ to $\sim 640^2\approx 4\times 10^5$. The condition number of $A$ becomes large, and the SVD solution is sensitive to floating-point errors.

The Solution

Normalize the input points so they are centered at the origin with a controlled scale, solve the system, then denormalize the result.

2D Normalization

Given $N$ points $\{{\mathbf{p}}_i=(x_i,y_i)\}$:

1. Centroid: $\bar{x}=\frac{1}{N}\sum x_i$, $\bar{y}=\frac{1}{N}\sum y_i$

2. Mean distance from centroid: $\bar{d}=\frac{1}{N}\sum \sqrt{(x_i-\bar{x}{)}^2+(y_i-\bar{y}{)}^2}$

3. Scale factor: $s=\sqrt{2}/\bar{d}$

4. Normalization matrix:

$$T=\left[\begin{array}{ccc}s& 0& -s\bar{x}\\ 0& s& -s\bar{y}\\ 0& 0& 1\end{array}\right]$$

5. Normalized points: ${\tilde{\mathbf{p}}}_i=T\left[\begin{array}{c}{x}_i\\ y_i\\ 1\end{array}\right]$

After normalization, the points have centroid at the origin and mean distance $\sqrt{2}$ from the origin.

3D Normalization

For 3D points $\{{\mathbf{P}}_i=(X_i,Y_i,Z_i)\}$, the analogous normalization uses:

• Mean distance target: $\sqrt{3}$ (instead of $\sqrt{2}$)
• A $4\times 4$ normalization matrix $T_{\text{3D}}$

Denormalization

After solving the normalized system, the result must be denormalized. For a homography estimated from normalized points:

$$H=T_{\text{image}}^{-1}\cdot H_{\text{norm}}\cdot T_{\text{world}}$$

For a fundamental matrix:

$$F=T_2^T\cdot F_{\text{norm}}\cdot T_1$$

The denormalization formula depends on the specific problem. Each chapter states the appropriate denormalization.

When It Is Used

Hartley normalization is applied in every DLT-based solver in calibration-rs:

Reference

Hartley, R.I. (1997). "In Defense of the Eight-Point Algorithm." IEEE Transactions on Pattern Analysis and Machine Intelligence, 19(6), 580-593.

The normalization technique is described in Algorithm 4.2 of Hartley & Zisserman (2004).

Homography Estimation (DLT)

A homography is a $3\times 3$ projective transformation between two planes. In calibration, the most common case is the mapping from a planar calibration board to the image plane. Homography estimation is the first step in Zhang's calibration method.

Problem Statement

Given: $N\ge 4$ point correspondences $\{({\mathbf{x}}_i,{\mathbf{x}}_i^{\prime })\}$ where ${\mathbf{x}}_i$ are points on the world plane and ${\mathbf{x}}_i^{\prime }$ are the corresponding image points.

Find: $H\in {\mathbb{R}}^{3\times 3}$ such that ${\mathbf{x}}_i^{\prime }\sim H{\mathbf{x}}_i$ (equality up to scale) for all $i$.

Assumptions:

• The world points are coplanar (the calibration board is flat)
• Correspondences are correct (or RANSAC is used to handle outliers)
• At least 4 non-collinear point correspondences

Derivation

From Projective Relation to Linear System

The relation ${\mathbf{x}}^{\prime }\sim H\mathbf{x}$ in homogeneous coordinates means:

$$\left[\begin{array}{c}{u}^{\prime }\\ v^{\prime }\\ 1\end{array}\right]\sim \left[\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$

Expanding and eliminating the scale factor using cross-product (the constraint ${\mathbf{x}}^{\prime }\times H\mathbf{x}=0$), we get two independent equations per correspondence:

$$-x\cdot h_{11}-y\cdot h_{12}-h_{13}+u^{\prime }(x\cdot h_{31}+y\cdot h_{32}+h_{33})=0$$

$$-x\cdot h_{21}-y\cdot h_{22}-h_{23}+v^{\prime }(x\cdot h_{31}+y\cdot h_{32}+h_{33})=0$$

The DLT System

Stacking all $N$ correspondences into a $2N\times 9$ matrix $A$:

$$A=\left[\begin{array}{ccccccccc}-x_1& -y_1& -1& 0& 0& 0& u_1^{\prime }{x}_1& u_1^{\prime }{y}_1& u_1^{\prime }\\ 0& 0& 0& -x_1& -y_1& -1& v_1^{\prime }{x}_1& v_1^{\prime }{y}_1& v_1^{\prime }\\ ⋮& & & & & & & & ⋮\end{array}\right]$$

where $\mathbf{h}=[h_{11},h_{12},h_{13},h_{21},h_{22},h_{23},h_{31},h_{32},h_{33}{]}^T$.

Solving via SVD

We seek $\mathbf{h}$ minimizing $\Vert A\mathbf{h}{\Vert }^2$ subject to $\Vert \mathbf{h}\Vert =1$ (to avoid the trivial solution $\mathbf{h}=0$). This is the standard homogeneous least squares problem, solved by the SVD of $A$:

$$A=U\Sigma V^T$$

The solution is the last column of $V$ (corresponding to the smallest singular value). Reshaping this 9-vector into a $3\times 3$ matrix gives $H$.

Normalization and Denormalization

The complete algorithm with Hartley normalization:

1. Normalize world and image points: $({\mathbf{x}}_i,T_w)$ and $({\mathbf{x}}_i^{\prime },T_i)$
2. Build the $2N\times 9$ matrix $A$ from normalized points
3. Solve via SVD to get $H_{\text{norm}}$
4. Denormalize: $H=T_i^{-1}\cdot H_{\text{norm}}\cdot T_w$
5. Scale: $H\leftarrow H/H_{33}$

Degrees of Freedom

$H$ is a $3\times 3$ matrix with 9 entries, but it is defined up to scale, giving 8 degrees of freedom. Each correspondence provides 2 equations, so the minimum is $N=4$ correspondences (giving $2\times 4=8$ equations for 8 unknowns).

With $N>4$, the system is overdetermined and the SVD solution minimizes the algebraic error in a least-squares sense.

RANSAC Wrapper

For real data with potential mismatched correspondences, dlt_homography_ransac wraps the DLT solver in RANSAC:

#![allow(unused)]
fn main() {
let (H, inliers) = dlt_homography_ransac(&world_pts, &image_pts, &opts)?;
}

• MIN_SAMPLES = 4
• Residual: Euclidean reprojection error in pixels: $\Vert \text{dehomogenize}(H\mathbf{x})-{\mathbf{x}}^{\prime }\Vert$
• Degeneracy check: Tests if the first 3 world points are collinear (a degenerate configuration for homography estimation)

OpenCV equivalence: cv::findHomography with RANSAC method.

API

#![allow(unused)]
fn main() {
// Direct DLT (all points assumed inliers)
let H = dlt_homography(&world_2d, &image_2d)?;

// With RANSAC
let opts = RansacOptions {
    max_iters: 1000,
    thresh: 3.0,   // pixels
    confidence: 0.99,
    ..Default::default()
};
let (H, inliers) = dlt_homography_ransac(&world_2d, &image_2d, &opts)?;
}

Geometric Interpretation

For a planar calibration board at $Z=0$ in world coordinates, the $3\times 4$ camera projection matrix $P=K[R|\mathbf{t}]$ reduces to a $3\times 3$ homography:

$$H\sim K[{\mathbf{r}}_1{\mathbf{r}}_2\mathbf{t}]$$

where ${\mathbf{r}}_1,{\mathbf{r}}_2$ are the first two columns of $R$. This relationship is the foundation of Zhang's calibration method (next chapter).

Zhang's Intrinsics from Homographies

Zhang's method is the most widely used technique for estimating camera intrinsics from a planar calibration board. Given homographies from multiple views of the board, it recovers the intrinsics matrix $K$ using constraints from the image of the absolute conic (IAC).

Problem Statement

Given: $M\ge 3$ homographies $\{H_k{\}}_{k=1}^M$ from a planar calibration board to the image, computed from different viewpoints.

Find: Camera intrinsics $K=\left[\begin{array}{ccc}{f}_x& s& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]$.

Assumptions:

• The calibration board is planar
• At least 3 homographies from distinct viewpoints
• The camera intrinsics are constant across all views
• Distortion is negligible (or has been accounted for — see Iterative Intrinsics)

Derivation

Homography and Intrinsics

From the previous chapter, the homography from a $Z=0$ board plane is:

$$H=\lambda K[{\mathbf{r}}_1{\mathbf{r}}_2\mathbf{t}]$$

where $\lambda$ is an arbitrary scale and ${\mathbf{r}}_1,{\mathbf{r}}_2$ are columns of the rotation matrix. Let $H=[{\mathbf{h}}_1{\mathbf{h}}_2{\mathbf{h}}_3]$. Then:

$${\mathbf{h}}_1=\lambda K{\mathbf{r}}_1,{\mathbf{h}}_2=\lambda K{\mathbf{r}}_2$$

Orthogonality Constraints

Since ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$ are columns of a rotation matrix, they satisfy:

$${\mathbf{r}}_1^T{\mathbf{r}}_2=0\text{(orthogonality)}$$ $${\mathbf{r}}_1^T{\mathbf{r}}_1={\mathbf{r}}_2^T{\mathbf{r}}_2\text{(equal norm)}$$

Substituting ${\mathbf{r}}_i=\frac{1}{\lambda }{K}^{-1}{\mathbf{h}}_i$:

$${\mathbf{h}}_1^T{K}^{-T}{K}^{-1}{\mathbf{h}}_2=0$$ $${\mathbf{h}}_1^T{K}^{-T}{K}^{-1}{\mathbf{h}}_1={\mathbf{h}}_2^T{K}^{-T}{K}^{-1}{\mathbf{h}}_2$$

The Image of the Absolute Conic (IAC)

Define the symmetric matrix:

$$B=K^{-T}{K}^{-1}=\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{12}& B_{22}& B_{23}\\ B_{13}& B_{23}& B_{33}\end{array}\right]$$

Since $B$ is symmetric, it has 6 independent entries. The constraints become:

$${\mathbf{h}}_1^{T}B{\mathbf{h}}_2=0$$ $${\mathbf{h}}_1^{T}B{\mathbf{h}}_1-{\mathbf{h}}_2^{T}B{\mathbf{h}}_2=0$$

Each homography gives 2 linear equations in the 6 unknowns of $B$.

The ${\mathbf{v}}_{ij}$ Vectors

To write the constraints as a linear system, define the 6-vector:

$${\mathbf{v}}_{ij}=\left[\begin{array}{c}{h}_{1i}{h}_{1j}\\ h_{1i}{h}_{2j}+h_{2i}{h}_{1j}\\ h_{2i}{h}_{2j}\\ h_{3i}{h}_{1j}+h_{1i}{h}_{3j}\\ h_{3i}{h}_{2j}+h_{2i}{h}_{3j}\\ h_{3i}{h}_{3j}\end{array}\right]$$

where $h_{ki}$ is the $(k,i)$-th element of $H$ (1-indexed, column $i$, row $k$). Then:

$${\mathbf{h}}_i^{T}B{\mathbf{h}}_j={\mathbf{v}}_{ij}^T\mathbf{b}$$

where $\mathbf{b}=[B_{11},B_{12},B_{22},B_{13},B_{23},B_{33}{]}^T$.

The Linear System

For each homography $H_k$, the two constraints become:

$$\left[\begin{array}{c}{\mathbf{v}}_{12}^T\\ ({\mathbf{v}}_{11}-{\mathbf{v}}_{22}{)}^T\end{array}\right]\mathbf{b}=0$$

Stacking $M$ homographies gives a $2M\times 6$ system:

$$V\mathbf{b}=0$$

With $M\ge 3$ (giving $\ge 6$ equations), this is solved via SVD: $\mathbf{b}$ is the right singular vector of $V$ corresponding to the smallest singular value.

Extracting Intrinsics from $B$

Given $\mathbf{b}=[B_{11},B_{12},B_{22},B_{13},B_{23},B_{33}]$, the intrinsics are recovered by a Cholesky-like factorization of $B=K^{-T}{K}^{-1}$.

The closed-form extraction:

$$v_0=\frac{B_{12}{B}_{13}-B_{11}{B}_{23}}{B_{11}{B}_{22}-B_{12}^2}$$

$$\lambda =B_{33}-\frac{B_{13}^2+v_0(B_{12}{B}_{13}-B_{11}{B}_{23})}{B_{11}}$$

$$f_x=\sqrt{\frac{\lambda }{B_{11}}}$$

$$f_y=\sqrt{\frac{\lambda B_{11}}{B_{11}{B}_{22}-B_{12}^2}}$$

$$s=-\frac{B_{12}{f}_x^2{f}_y}{\lambda }$$

$$c_x=\frac{s\cdot v_0}{f_y}-\frac{B_{13}{f}_x^2}{\lambda }$$

$$c_y=v_0$$

Validity Checks

The extraction can fail when:

• $B_{11}{B}_{22}-B_{12}^2\approx 0$ (degenerate: insufficient view diversity)
• $|B_{11}|\approx 0$
• $\lambda$ and $B_{11}$ have different signs (negative focal length squared)

These cases indicate that the input homographies do not sufficiently constrain the intrinsics, typically because the views are too similar (e.g., all near-frontal).

Minimum Number of Views

• $M=3$ homographies: 6 equations for 6 unknowns — the minimum for a unique solution (with skew)
• $M=2$: Only 4 equations; requires the additional assumption that $s=0$ (zero skew) to reduce unknowns to 5
• $M>3$: Overdetermined system; SVD gives the least-squares solution

Practical Considerations

View diversity: The views should include significant rotation around both axes. Pure translations or rotations around a single axis lead to degenerate configurations where some intrinsic parameters are undetermined.

Distortion: Zhang's method assumes distortion-free observations. When applied to distorted pixels, the estimated $K$ is biased. The Iterative Intrinsics chapter addresses this with an alternating refinement scheme.

OpenCV equivalence: Zhang's method is the internal initialization step of cv::calibrateCamera. OpenCV does not expose it as a separate API.

API

#![allow(unused)]
fn main() {
let K = estimate_intrinsics_from_homographies(&homographies)?;
}

Returns FxFyCxCySkew with the estimated intrinsics. Typically followed by distortion estimation and non-linear refinement.

Reference

Zhang, Z. (2000). "A Flexible New Technique for Camera Calibration." IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(11), 1330-1334.

Distortion Estimation from Homography Residuals

After estimating intrinsics $K$ via Zhang's method (which assumes distortion-free observations), the residuals between observed and predicted pixel positions encode the lens distortion. This chapter derives a linear method to estimate Brown-Conrady distortion coefficients from these residuals.

Problem Statement

Given:

• Camera intrinsics $K$ (from Zhang's method)
• $M$ homographies $\{H_k\}$ from the calibration board to the image
• Point correspondences: board points $\{{\mathbf{P}}_j\}$ and observed pixels $\{{\mathbf{p}}_{kj}\}$

Find: Distortion coefficients $(k_1,k_2,k_3,p_1,p_2)$.

Assumptions:

• $K$ is a reasonable estimate (possibly biased by distortion)
• Distortion is moderate (the linear approximation holds)
• The homographies were computed from distorted pixel observations

Derivation

Ideal vs. Observed Coordinates

For each correspondence in view $k$, the ideal (undistorted) pixel position predicted by the homography is:

$$\tilde{\mathbf{p}}=\text{dehomogenize}(H_k\cdot [{\mathbf{P}}_j,1{]}^T)$$

Convert both observed and ideal pixels to normalized coordinates:

$${\mathbf{n}}_{\text{obs}}=K^{-1}\left[\begin{array}{c}{\mathbf{p}}_{\text{obs}}\\ 1\end{array}\right],{\mathbf{n}}_{\text{ideal}}=K^{-1}\left[\begin{array}{c}\tilde{\mathbf{p}}\\ 1\end{array}\right]$$

The residual in normalized coordinates encodes the distortion effect:

$$\Delta \mathbf{n}={\mathbf{n}}_{\text{obs}}-{\mathbf{n}}_{\text{ideal}}\approx D({\mathbf{n}}_{\text{ideal}})-{\mathbf{n}}_{\text{ideal}}$$

Linearized Distortion Model

The Brown-Conrady distortion applied to point $\mathbf{n}=[x,y{]}^T$ with $r^2=x^2+y^2$ produces a displacement:

$$\Delta x=x(k_1{r}^2+k_2{r}^4+k_3{r}^6)+2p_{1}xy+p_2(r^2+2x^2)$$ $$\Delta y=y(k_1{r}^2+k_2{r}^4+k_3{r}^6)+p_1(r^2+2y^2)+2p_{2}xy$$

This is linear in the distortion coefficients $(k_1,k_2,k_3,p_1,p_2)$.

Building the Linear System

For each correspondence, write the $x$ and $y$ components:

$$\left[\begin{array}{ccccc}x{r}^2& x{r}^4& x{r}^6& 2xy& r^2+2x^2\\ y{r}^2& y{r}^4& y{r}^6& r^2+2y^2& 2xy\end{array}\right]\left[\begin{array}{c}{k}_1\\ k_2\\ k_3\\ p_1\\ p_2\end{array}\right]=\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]$$

where $x,y$ are the ideal normalized coordinates and $\Delta x,\Delta y$ are the observed residuals.

Stacking all correspondences from all views gives an overdetermined system:

$$A\mathbf{c}=\mathbf{b}$$

where $A$ is $2NM\times 5$ (or smaller if some coefficients are fixed), and $\mathbf{c}=[k_1,k_2,k_3,p_1,p_2{]}^T$.

Solving

The least-squares solution is:

$$\mathbf{c}=(A^{T}A{)}^{-1}{A}^T\mathbf{b}$$

In practice, this is solved via SVD of $A$ for numerical stability.

Options

#![allow(unused)]
fn main() {
pub struct DistortionFitOptions {
    pub fix_tangential: bool,  // Fix p1 = p2 = 0
    pub fix_k3: bool,          // Fix k3 = 0 (default: true)
    pub iters: u32,            // Undistortion iterations
}
}

• fix_k3 (default true): Removes $k_3$ from the system, reducing to a 4-parameter or 2-parameter model. Recommended unless the lens has strong higher-order radial distortion.
• fix_tangential (default false): Removes $p_1,p_2$, estimating only radial distortion. Useful when tangential distortion is known to be negligible.

When parameters are fixed, the corresponding columns are removed from $A$.

Accuracy

This linear estimate is typically within 10-50% of the true distortion values. The accuracy depends on:

• Quality of $K$: If intrinsics are biased (from distorted observations), the estimated distortion absorbs some of the bias
• Number of points: More correspondences improve the overdetermined system
• Point distribution: Points across the full image area constrain distortion better than points clustered near the center (where distortion is small)

This accuracy is sufficient for initializing non-linear refinement.

API

#![allow(unused)]
fn main() {
let dist = estimate_distortion_from_homographies(
    &k_matrix, &views, opts
)?;
}

Returns BrownConrady5 with the estimated coefficients.

OpenCV equivalence: OpenCV estimates distortion jointly with intrinsics inside cv::calibrateCamera. The separate distortion estimation step is specific to calibration-rs's initialization approach.

Iterative Intrinsics + Distortion

Zhang's method assumes distortion-free observations. When applied to images from a camera with significant distortion, the estimated intrinsics are biased because the distorted pixels violate the linear homography model. The iterative intrinsics solver addresses this by alternating between intrinsics estimation and distortion correction.

Problem Statement

Given: $M$ views of a planar calibration board, with observed (distorted) pixel coordinates.

Find: Camera intrinsics $K$ and distortion coefficients $(k_1,k_2,k_3,p_1,p_2)$ jointly.

Assumptions:

• The observations are distorted (raw detector output)
• Distortion is moderate enough that Zhang's method on distorted pixels gives a usable initial $K$
• 1-3 iterations of alternating refinement suffice

Why Alternation Works

The joint estimation of $K$ and distortion is a non-convex problem. However, it decomposes naturally:

• Given $K$: Distortion estimation is a linear problem (see Distortion Fit)
• Given distortion: Undistorting the pixels and re-estimating $K$ is a linear problem (Zhang's method)

Each step solves a convex subproblem. The alternation converges because:

1. The initial Zhang estimate (ignoring distortion) is typically within the basin where the alternation contracts
2. Each step reduces the residual between the model and observations
3. The coupling between $K$ and distortion is relatively weak for moderate distortion

Algorithm

$$\boxed{\text{Iterative Intrinsics Estimation}}$$

Input: Views $\{({\mathbf{P}}_j,{\mathbf{p}}_{kj})\}$ with distorted observations ${\mathbf{p}}_{kj}$

Output: Intrinsics $K$, distortion $\mathbf{d}$

1. Compute homographies $\{H_k\}$ from distorted pixels via DLT

2. Estimate initial $K^{(0)}$ from $\{H_k\}$ via Zhang's method

3. For $i=1,\dots ,n_{\text{iter}}$:

   a. Estimate distortion ${\mathbf{d}}^{(i)}$ from homography residuals using $K^{(i-1)}$

   b. Undistort all observed pixels: ${\hat{\mathbf{p}}}_{kj}=\text{undistort}({\mathbf{p}}_{kj},K^{(i-1)},{\mathbf{d}}^{(i)})$

   c. Recompute homographies $\{H_k^{(i)}\}$ from undistorted pixels

   d. Re-estimate $K^{(i)}$ from $\{H_k^{(i)}\}$ via Zhang's method

   e. (Optional) Enforce zero skew: $K_{12}^{(i)}=0$

4. Return $K^{(n_{\text{iter}})},{\mathbf{d}}^{(n_{\text{iter}})}$

Convergence

Typically 1-2 iterations suffice:

• Iteration 0 (Zhang on distorted pixels): 10-40% intrinsics error
• Iteration 1: Distortion estimate corrects the dominant radial effect; intrinsics error drops to 5-20%
• Iteration 2: Further refinement; diminishing returns

More iterations are safe but rarely improve the estimate significantly. The default is 2 iterations.

Configuration

#![allow(unused)]
fn main() {
pub struct IterativeIntrinsicsOptions {
    pub iterations: usize,                    // 1-3 typical (default: 2)
    pub distortion_opts: DistortionFitOptions, // Controls fix_k3, fix_tangential, iters
    pub zero_skew: bool,                      // Force skew = 0 (default: true)
}
}

The Undistortion Step

Step 3b converts distorted pixels back to "ideal" pixels by inverting the distortion model:

1. Convert distorted pixel to normalized coordinates: ${\mathbf{n}}_d=K^{-1}[\mathbf{p},1{]}^T$
2. Undistort using the current distortion estimate: ${\mathbf{n}}_u=D^{-1}({\mathbf{n}}_d)$ (iterative fixed-point, see Brown-Conrady Distortion)
3. Convert back to pixel: $\hat{\mathbf{p}}=K\cdot [{\mathbf{n}}_u,1{]}^T$

Note that this uses $K$ both for normalization and for converting back — the undistorted pixels are in the same coordinate system as the original distorted pixels, just without the distortion effect.

Accuracy Expectations

The iterative linear estimate is not meant to be highly accurate. Its purpose is to provide a starting point good enough for the non-linear optimizer to converge.

API

#![allow(unused)]
fn main() {
use vision_calibration::linear::iterative_intrinsics::*;
use vision_calibration::linear::DistortionFitOptions;

let opts = IterativeIntrinsicsOptions {
    iterations: 2,
    distortion_opts: DistortionFitOptions {
        fix_k3: true,
        fix_tangential: false,
        iters: 8,
    },
    zero_skew: true,
};

let camera = estimate_intrinsics_iterative(&dataset, opts)?;
// camera is PinholeCamera = Camera<f64, Pinhole, BrownConrady5, IdentitySensor, FxFyCxCySkew>
// camera.k: FxFyCxCySkew (intrinsics)
// camera.dist: BrownConrady5 (distortion)
}

Pose from Homography

Given camera intrinsics $K$ and a homography $H$ from a planar calibration board to the image, we can decompose $H$ to recover the camera pose (rotation $R$ and translation $\mathbf{t}$) relative to the board.

Problem Statement

Given: Intrinsics $K$ and homography $H$ such that $\mathbf{p}\sim H\cdot [{\mathbf{P}}_{xy},1{]}^T$ for board points at $Z=0$.

Find: Rigid transform $T_{C,B}=[R|\mathbf{t}]$ (board to camera).

Assumptions:

• The board lies at $Z=0$ in world coordinates
• $K$ is known (or has been estimated)
• The homography $H$ was computed from correct correspondences

Derivation

Extracting Rotation and Translation

Recall from the Homography chapter:

$$H\sim K[{\mathbf{r}}_1{\mathbf{r}}_2\mathbf{t}]$$

where ${\mathbf{r}}_1,{\mathbf{r}}_2$ are the first two columns of the rotation matrix and $\mathbf{t}$ is the translation.

Removing the intrinsics:

$$K^{-1}H=\lambda [{\mathbf{r}}_1{\mathbf{r}}_2\mathbf{t}]$$

Let $[{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3]=K^{-1}H$. The scale factor $\lambda$ is recovered from the constraint that ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$ have unit norm:

$$\lambda =\frac{2}{\Vert {\mathbf{a}}_1\Vert +\Vert {\mathbf{a}}_2\Vert }$$

Then:

$${\mathbf{r}}_1=\lambda {\mathbf{a}}_1,{\mathbf{r}}_2=\lambda {\mathbf{a}}_2,\mathbf{t}=\lambda {\mathbf{a}}_3$$

The third rotation column is:

$${\mathbf{r}}_3={\mathbf{r}}_1\times {\mathbf{r}}_2$$

Projecting onto SO(3)

Due to noise, the matrix $R_{\text{approx}}=[{\mathbf{r}}_1{\mathbf{r}}_2{\mathbf{r}}_3]$ is not exactly orthonormal. We project it onto $\text{SO}(3)$ using SVD:

$$R_{\text{approx}}=U\Sigma V^T$$

$$R=U{V}^T$$

If $\det (R)=-1$, flip the sign of the third column of $U$ and recompute.

Ensuring Forward-Facing

If $t_z<0$, the board is behind the camera. In this case, flip the sign of both $R$ and $\mathbf{t}$:

$$R\leftarrow -R,\mathbf{t}\leftarrow -\mathbf{t}$$

This resolves the sign ambiguity inherent in the scale factor $\lambda$.

Accuracy

The pose from homography is an approximate estimate because:

1. The homography itself is subject to noise (DLT algebraic error minimization, not geometric)
2. The SVD projection onto SO(3) corrects non-orthogonality but introduces additional error
3. Distortion (if not corrected) biases the homography

Typical rotation error: 1-5°. Typical translation direction error: 5-15%. These estimates are refined in non-linear optimization.

OpenCV equivalence: cv::decomposeHomographyMat provides a similar decomposition, returning up to 4 pose candidates. calibration-rs returns a single pose by resolving ambiguities via the forward-facing constraint.

API

#![allow(unused)]
fn main() {
let pose = estimate_planar_pose_from_h(&K, &H)?;
// pose: Iso3 (T_C_B: board to camera transform)
}

Usage in Calibration Pipeline

In the planar intrinsics calibration pipeline, pose estimation is applied to every view after $K$ has been estimated:

1. Compute homography $H_k$ per view
2. Estimate $K$ from all homographies (Zhang's method)
3. Decompose each $H_k$ to get pose $T_{C,B}^{(k)}$
4. Use poses as initial values for non-linear optimization