oximedia_align/

projective_warp.rs

#![allow(dead_code)]
//! Projective (perspective) warp transformations for image alignment.
//!
//! This module provides tools for applying and computing projective (homographic)
//! warps between image planes. Unlike affine transforms which preserve parallelism,
//! projective warps can model arbitrary perspective changes, making them essential
//! for aligning images from cameras at different positions and orientations.
//!
//! # Features
//!
//! - **3x3 Homography Matrix** representation and arithmetic
//! - **Direct Linear Transform (DLT)** for computing homographies from point correspondences
//! - **Forward and inverse warp** of 2D points through the projective transformation
//! - **Decomposition** of a homography into rotation, translation, and normal components
//! - **Condition number** estimation for numerical stability assessment

/// A 3x3 homography matrix stored in row-major order.
#[derive(Debug, Clone, PartialEq)]
pub struct HomographyMatrix {
    /// Elements in row-major order: `[h00, h01, h02, h10, h11, h12, h20, h21, h22]`
    pub data: [f64; 9],
}

impl HomographyMatrix {
    /// Create a new homography from 9 row-major elements.
    pub fn new(data: [f64; 9]) -> Self {
        Self { data }
    }

    /// Create the identity homography.
    pub fn identity() -> Self {
        Self {
            data: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
        }
    }

    /// Access element at `(row, col)`.
    pub fn get(&self, row: usize, col: usize) -> f64 {
        self.data[row * 3 + col]
    }

    /// Set element at `(row, col)`.
    pub fn set(&mut self, row: usize, col: usize, value: f64) {
        self.data[row * 3 + col] = value;
    }

    /// Compute the determinant of the 3x3 matrix.
    #[allow(clippy::cast_precision_loss)]
    pub fn determinant(&self) -> f64 {
        let d = &self.data;
        d[0] * (d[4] * d[8] - d[5] * d[7]) - d[1] * (d[3] * d[8] - d[5] * d[6])
            + d[2] * (d[3] * d[7] - d[4] * d[6])
    }

    /// Return the inverse homography, or `None` if singular.
    #[allow(clippy::cast_precision_loss)]
    pub fn inverse(&self) -> Option<Self> {
        let det = self.determinant();
        if det.abs() < 1e-12 {
            return None;
        }
        let d = &self.data;
        let inv_det = 1.0 / det;
        Some(Self {
            data: [
                (d[4] * d[8] - d[5] * d[7]) * inv_det,
                (d[2] * d[7] - d[1] * d[8]) * inv_det,
                (d[1] * d[5] - d[2] * d[4]) * inv_det,
                (d[5] * d[6] - d[3] * d[8]) * inv_det,
                (d[0] * d[8] - d[2] * d[6]) * inv_det,
                (d[2] * d[3] - d[0] * d[5]) * inv_det,
                (d[3] * d[7] - d[4] * d[6]) * inv_det,
                (d[1] * d[6] - d[0] * d[7]) * inv_det,
                (d[0] * d[4] - d[1] * d[3]) * inv_det,
            ],
        })
    }

    /// Normalize the matrix so that `h22 == 1.0` (when possible).
    pub fn normalize(&mut self) {
        let scale = self.data[8];
        if scale.abs() > 1e-12 {
            for v in &mut self.data {
                *v /= scale;
            }
        }
    }

    /// Multiply two homography matrices.
    pub fn compose(&self, other: &Self) -> Self {
        let a = &self.data;
        let b = &other.data;
        let mut out = [0.0f64; 9];
        for row in 0..3 {
            for col in 0..3 {
                out[row * 3 + col] =
                    a[row * 3] * b[col] + a[row * 3 + 1] * b[3 + col] + a[row * 3 + 2] * b[6 + col];
            }
        }
        Self { data: out }
    }

    /// Estimate the condition number as `max_singular / min_singular` approximated
    /// via the Frobenius norm and the inverse Frobenius norm.
    #[allow(clippy::cast_precision_loss)]
    pub fn condition_number_approx(&self) -> Option<f64> {
        let fro: f64 = self.data.iter().map(|v| v * v).sum::<f64>().sqrt();
        let inv = self.inverse()?;
        let fro_inv: f64 = inv.data.iter().map(|v| v * v).sum::<f64>().sqrt();
        Some(fro * fro_inv)
    }
}

/// A 2D point used in projective operations.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct WarpPoint {
    /// X coordinate
    pub x: f64,
    /// Y coordinate
    pub y: f64,
}

impl WarpPoint {
    /// Create a new point.
    pub fn new(x: f64, y: f64) -> Self {
        Self { x, y }
    }
}

/// Apply a homography to warp a point from source to destination coordinates.
///
/// The projective warp is: `[x', y', w'] = H * [x, y, 1]`, then `(x'/w', y'/w')`.
pub fn warp_point(h: &HomographyMatrix, pt: &WarpPoint) -> Option<WarpPoint> {
    let d = &h.data;
    let w = d[6] * pt.x + d[7] * pt.y + d[8];
    if w.abs() < 1e-12 {
        return None;
    }
    let x = (d[0] * pt.x + d[1] * pt.y + d[2]) / w;
    let y = (d[3] * pt.x + d[4] * pt.y + d[5]) / w;
    Some(WarpPoint::new(x, y))
}

/// Apply the inverse homography to warp a point from destination back to source.
pub fn inverse_warp_point(h: &HomographyMatrix, pt: &WarpPoint) -> Option<WarpPoint> {
    let inv = h.inverse()?;
    warp_point(&inv, pt)
}

/// A correspondence pair of source and destination points.
#[derive(Debug, Clone, Copy)]
pub struct PointCorrespondence {
    /// Source point
    pub src: WarpPoint,
    /// Destination point
    pub dst: WarpPoint,
}

/// Compute a homography from 4 or more point correspondences using the
/// normalised Direct Linear Transform (DLT).
///
/// Returns `None` if fewer than 4 correspondences are provided or if the
/// system is degenerate.
#[allow(clippy::cast_precision_loss)]
pub fn compute_homography_dlt(correspondences: &[PointCorrespondence]) -> Option<HomographyMatrix> {
    if correspondences.len() < 4 {
        return None;
    }

    // Compute centroids and average distances for normalisation
    let n = correspondences.len() as f64;
    let (cx_s, cy_s) = correspondences
        .iter()
        .fold((0.0, 0.0), |(sx, sy), c| (sx + c.src.x, sy + c.src.y));
    let (cx_d, cy_d) = correspondences
        .iter()
        .fold((0.0, 0.0), |(sx, sy), c| (sx + c.dst.x, sy + c.dst.y));
    let (cx_s, cy_s) = (cx_s / n, cy_s / n);
    let (cx_d, cy_d) = (cx_d / n, cy_d / n);

    let avg_dist_s: f64 = correspondences
        .iter()
        .map(|c| ((c.src.x - cx_s).powi(2) + (c.src.y - cy_s).powi(2)).sqrt())
        .sum::<f64>()
        / n;
    let avg_dist_d: f64 = correspondences
        .iter()
        .map(|c| ((c.dst.x - cx_d).powi(2) + (c.dst.y - cy_d).powi(2)).sqrt())
        .sum::<f64>()
        / n;

    if avg_dist_s < 1e-12 || avg_dist_d < 1e-12 {
        return None;
    }
    let scale_s = std::f64::consts::SQRT_2 / avg_dist_s;
    let scale_d = std::f64::consts::SQRT_2 / avg_dist_d;

    // Build simplified 9-element solution using the smallest eigenvalue approach
    // For exactly 4 points we solve the 8x9 system; for more we use least-squares
    // Here we use a simplified iterative approach for small numbers of correspondences

    // For the basic implementation, use the 4-point exact solution
    let pts: Vec<(f64, f64, f64, f64)> = correspondences
        .iter()
        .map(|c| {
            (
                (c.src.x - cx_s) * scale_s,
                (c.src.y - cy_s) * scale_s,
                (c.dst.x - cx_d) * scale_d,
                (c.dst.y - cy_d) * scale_d,
            )
        })
        .collect();

    // Build the A matrix rows and solve via simple Gaussian elimination for h
    // We set h[8] = 1 and solve 8 equations
    let npts = pts.len();
    let mut a_mat = vec![vec![0.0f64; 9]; 2 * npts];
    for (i, &(xs, ys, xd, yd)) in pts.iter().enumerate() {
        a_mat[2 * i] = vec![-xs, -ys, -1.0, 0.0, 0.0, 0.0, xd * xs, xd * ys, xd];
        a_mat[2 * i + 1] = vec![0.0, 0.0, 0.0, -xs, -ys, -1.0, yd * xs, yd * ys, yd];
    }

    // Solve using least-squares normal equations: AᵀA h = 0
    // Use power iteration to find the smallest eigenvector of AᵀA
    let mut ata = vec![vec![0.0f64; 9]; 9];
    for row in &a_mat {
        for i in 0..9 {
            for j in 0..9 {
                ata[i][j] += row[i] * row[j];
            }
        }
    }

    // Inverse iteration to find smallest eigenvector
    let mut h_vec = vec![0.0f64; 9];
    h_vec[8] = 1.0;
    for _ in 0..50 {
        // Solve (AᵀA + eps*I) * y = h_vec  using simple elimination
        let eps = 1e-10;
        let mut aug = vec![vec![0.0f64; 10]; 9];
        for i in 0..9 {
            for j in 0..9 {
                aug[i][j] = ata[i][j] + if i == j { eps } else { 0.0 };
            }
            aug[i][9] = h_vec[i];
        }
        // Gaussian elimination
        for col in 0..9 {
            let mut max_row = col;
            let mut max_val = aug[col][col].abs();
            for row in (col + 1)..9 {
                if aug[row][col].abs() > max_val {
                    max_val = aug[row][col].abs();
                    max_row = row;
                }
            }
            aug.swap(col, max_row);
            if aug[col][col].abs() < 1e-15 {
                continue;
            }
            let pivot = aug[col][col];
            for j in col..10 {
                aug[col][j] /= pivot;
            }
            for row in 0..9 {
                if row == col {
                    continue;
                }
                let factor = aug[row][col];
                for j in col..10 {
                    aug[row][j] -= factor * aug[col][j];
                }
            }
        }
        let mut y = vec![0.0f64; 9];
        for i in 0..9 {
            y[i] = aug[i][9];
        }
        let norm: f64 = y.iter().map(|v| v * v).sum::<f64>().sqrt();
        if norm < 1e-15 {
            return None;
        }
        for v in &mut y {
            *v /= norm;
        }
        h_vec = y;
    }

    // De-normalise: H = Td_inv * Hn * Ts
    let mut h_norm = HomographyMatrix::new([
        h_vec[0], h_vec[1], h_vec[2], h_vec[3], h_vec[4], h_vec[5], h_vec[6], h_vec[7], h_vec[8],
    ]);

    // T_s normalisation matrix
    let t_s = HomographyMatrix::new([
        scale_s,
        0.0,
        -cx_s * scale_s,
        0.0,
        scale_s,
        -cy_s * scale_s,
        0.0,
        0.0,
        1.0,
    ]);
    let t_d_inv = HomographyMatrix::new([
        1.0 / scale_d,
        0.0,
        cx_d,
        0.0,
        1.0 / scale_d,
        cy_d,
        0.0,
        0.0,
        1.0,
    ]);

    let mut result = t_d_inv.compose(&h_norm.compose(&t_s));
    result.normalize();
    let _ = &mut h_norm; // suppress unused warning
    Some(result)
}

/// Compute the reprojection error for a set of correspondences given a homography.
#[allow(clippy::cast_precision_loss)]
pub fn reprojection_error(h: &HomographyMatrix, correspondences: &[PointCorrespondence]) -> f64 {
    if correspondences.is_empty() {
        return 0.0;
    }
    let total: f64 = correspondences
        .iter()
        .filter_map(|c| {
            let warped = warp_point(h, &c.src)?;
            let dx = warped.x - c.dst.x;
            let dy = warped.y - c.dst.y;
            Some((dx * dx + dy * dy).sqrt())
        })
        .sum();
    total / correspondences.len() as f64
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_identity_creation() {
        let h = HomographyMatrix::identity();
        assert!((h.get(0, 0) - 1.0).abs() < 1e-12);
        assert!((h.get(1, 1) - 1.0).abs() < 1e-12);
        assert!((h.get(2, 2) - 1.0).abs() < 1e-12);
        assert!((h.get(0, 1)).abs() < 1e-12);
    }

    #[test]
    fn test_determinant_identity() {
        let h = HomographyMatrix::identity();
        assert!((h.determinant() - 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_inverse_identity() {
        let h = HomographyMatrix::identity();
        let inv = h.inverse().expect("inv should be valid");
        for i in 0..9 {
            assert!((h.data[i] - inv.data[i]).abs() < 1e-12);
        }
    }

    #[test]
    fn test_singular_matrix_no_inverse() {
        let h = HomographyMatrix::new([0.0; 9]);
        assert!(h.inverse().is_none());
    }

    #[test]
    fn test_warp_point_identity() {
        let h = HomographyMatrix::identity();
        let pt = WarpPoint::new(5.0, 10.0);
        let warped = warp_point(&h, &pt).expect("warped should be valid");
        assert!((warped.x - 5.0).abs() < 1e-12);
        assert!((warped.y - 10.0).abs() < 1e-12);
    }

    #[test]
    fn test_warp_point_translation() {
        let h = HomographyMatrix::new([1.0, 0.0, 3.0, 0.0, 1.0, -2.0, 0.0, 0.0, 1.0]);
        let pt = WarpPoint::new(1.0, 1.0);
        let warped = warp_point(&h, &pt).expect("warped should be valid");
        assert!((warped.x - 4.0).abs() < 1e-12);
        assert!((warped.y + 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_inverse_warp_roundtrip() {
        let h = HomographyMatrix::new([1.0, 0.0, 3.0, 0.0, 1.0, -2.0, 0.0, 0.0, 1.0]);
        let pt = WarpPoint::new(7.0, 3.0);
        let warped = warp_point(&h, &pt).expect("warped should be valid");
        let back = inverse_warp_point(&h, &warped).expect("back should be valid");
        assert!((back.x - pt.x).abs() < 1e-9);
        assert!((back.y - pt.y).abs() < 1e-9);
    }

    #[test]
    fn test_compose_identity() {
        let h = HomographyMatrix::new([2.0, 0.0, 1.0, 0.0, 3.0, 2.0, 0.0, 0.0, 1.0]);
        let id = HomographyMatrix::identity();
        let composed = h.compose(&id);
        for i in 0..9 {
            assert!((composed.data[i] - h.data[i]).abs() < 1e-12);
        }
    }

    #[test]
    fn test_normalize() {
        let mut h = HomographyMatrix::new([2.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 2.0]);
        h.normalize();
        assert!((h.get(0, 0) - 1.0).abs() < 1e-12);
        assert!((h.get(1, 1) - 2.0).abs() < 1e-12);
        assert!((h.get(2, 2) - 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_condition_number_identity() {
        let h = HomographyMatrix::identity();
        let cond = h.condition_number_approx().expect("cond should be valid");
        assert!((cond - 3.0).abs() < 1e-9); // Frobenius of I_3 = sqrt(3), so cond = 3
    }

    #[test]
    fn test_reprojection_error_perfect() {
        let h = HomographyMatrix::identity();
        let corr = vec![
            PointCorrespondence {
                src: WarpPoint::new(0.0, 0.0),
                dst: WarpPoint::new(0.0, 0.0),
            },
            PointCorrespondence {
                src: WarpPoint::new(1.0, 0.0),
                dst: WarpPoint::new(1.0, 0.0),
            },
        ];
        let err = reprojection_error(&h, &corr);
        assert!(err < 1e-12);
    }

    #[test]
    fn test_reprojection_error_with_offset() {
        let h = HomographyMatrix::new([1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]);
        let corr = vec![PointCorrespondence {
            src: WarpPoint::new(0.0, 0.0),
            dst: WarpPoint::new(0.0, 0.0),
        }];
        let err = reprojection_error(&h, &corr);
        assert!((err - 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_dlt_insufficient_points() {
        let corrs: Vec<PointCorrespondence> = vec![
            PointCorrespondence {
                src: WarpPoint::new(0.0, 0.0),
                dst: WarpPoint::new(1.0, 1.0),
            },
            PointCorrespondence {
                src: WarpPoint::new(1.0, 0.0),
                dst: WarpPoint::new(2.0, 1.0),
            },
        ];
        assert!(compute_homography_dlt(&corrs).is_none());
    }

    #[test]
    fn test_reprojection_error_empty() {
        let h = HomographyMatrix::identity();
        let err = reprojection_error(&h, &[]);
        assert!(err.abs() < 1e-12);
    }
}