oximedia_align/

stereo_rectify.rs

#![allow(dead_code)]
//! Stereo image pair rectification for epipolar geometry alignment.
//!
//! This module implements stereo rectification to transform two camera views so that
//! corresponding epipolar lines become horizontally aligned, simplifying stereo matching.
//!
//! # Features
//!
//! - **Fundamental matrix estimation** from point correspondences
//! - **Essential matrix decomposition** for calibrated cameras
//! - **Hartley rectification** - uncalibrated stereo rectification
//! - **Bouguet rectification** - calibrated stereo rectification splitting rotation
//! - **Epipolar distance computation** for correspondence validation

use crate::{AlignError, AlignResult, Point2D};

/// A 3x3 matrix stored in row-major order.
#[derive(Debug, Clone, PartialEq)]
pub struct Matrix3x3 {
    /// The 9 elements of the matrix in row-major order.
    pub data: [f64; 9],
}

impl Matrix3x3 {
    /// Create a new 3x3 matrix from row-major data.
    #[must_use]
    pub fn new(data: [f64; 9]) -> Self {
        Self { data }
    }

    /// Create the 3x3 identity matrix.
    #[must_use]
    pub fn identity() -> Self {
        Self {
            data: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
        }
    }

    /// Create a zero 3x3 matrix.
    #[must_use]
    pub fn zero() -> Self {
        Self { data: [0.0; 9] }
    }

    /// Get the element at row `r` and column `c`.
    #[must_use]
    pub fn at(&self, r: usize, c: usize) -> f64 {
        self.data[r * 3 + c]
    }

    /// Set the element at row `r` and column `c`.
    pub fn set(&mut self, r: usize, c: usize, val: f64) {
        self.data[r * 3 + c] = val;
    }

    /// Compute the transpose of this matrix.
    #[must_use]
    pub fn transpose(&self) -> Self {
        let d = &self.data;
        Self {
            data: [d[0], d[3], d[6], d[1], d[4], d[7], d[2], d[5], d[8]],
        }
    }

    /// Compute the determinant of this matrix.
    #[must_use]
    #[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])
    }

    /// Multiply this matrix by another 3x3 matrix.
    #[must_use]
    pub fn multiply(&self, other: &Self) -> Self {
        let a = &self.data;
        let b = &other.data;
        let mut result = [0.0; 9];
        for i in 0..3 {
            for j in 0..3 {
                result[i * 3 + j] =
                    a[i * 3] * b[j] + a[i * 3 + 1] * b[3 + j] + a[i * 3 + 2] * b[6 + j];
            }
        }
        Self { data: result }
    }

    /// Multiply this matrix by a 3-vector, returning a 3-vector.
    #[must_use]
    pub fn multiply_vec(&self, v: &[f64; 3]) -> [f64; 3] {
        let d = &self.data;
        [
            d[0] * v[0] + d[1] * v[1] + d[2] * v[2],
            d[3] * v[0] + d[4] * v[1] + d[5] * v[2],
            d[6] * v[0] + d[7] * v[1] + d[8] * v[2],
        ]
    }

    /// Compute the Frobenius norm of the matrix.
    #[must_use]
    pub fn frobenius_norm(&self) -> f64 {
        self.data.iter().map(|x| x * x).sum::<f64>().sqrt()
    }
}

/// A pair of corresponding points in stereo images.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct StereoCorrespondence {
    /// Point in the left image.
    pub left: Point2D,
    /// Point in the right image.
    pub right: Point2D,
}

impl StereoCorrespondence {
    /// Create a new stereo correspondence.
    #[must_use]
    pub fn new(left: Point2D, right: Point2D) -> Self {
        Self { left, right }
    }
}

/// Configuration for stereo rectification.
#[derive(Debug, Clone)]
pub struct StereoRectifyConfig {
    /// Image width in pixels.
    pub image_width: u32,
    /// Image height in pixels.
    pub image_height: u32,
    /// Maximum number of RANSAC iterations for fundamental matrix estimation.
    pub max_iterations: u32,
    /// Inlier threshold in pixels for epipolar distance.
    pub inlier_threshold: f64,
    /// Minimum number of inliers required.
    pub min_inliers: usize,
}

impl Default for StereoRectifyConfig {
    fn default() -> Self {
        Self {
            image_width: 1920,
            image_height: 1080,
            max_iterations: 2000,
            inlier_threshold: 2.0,
            min_inliers: 8,
        }
    }
}

/// Result of stereo rectification containing the two rectifying homographies.
#[derive(Debug, Clone)]
pub struct RectificationResult {
    /// Rectifying homography for the left image.
    pub h_left: Matrix3x3,
    /// Rectifying homography for the right image.
    pub h_right: Matrix3x3,
    /// The estimated fundamental matrix.
    pub fundamental: Matrix3x3,
    /// Number of inlier correspondences used.
    pub num_inliers: usize,
    /// Mean epipolar error after rectification in pixels.
    pub mean_error: f64,
}

/// Stereo rectification engine.
#[derive(Debug, Clone)]
pub struct StereoRectifier {
    /// Configuration for the rectification process.
    config: StereoRectifyConfig,
}

impl StereoRectifier {
    /// Create a new stereo rectifier with the given configuration.
    #[must_use]
    pub fn new(config: StereoRectifyConfig) -> Self {
        Self { config }
    }

    /// Create a stereo rectifier with default configuration.
    #[must_use]
    pub fn with_defaults() -> Self {
        Self {
            config: StereoRectifyConfig::default(),
        }
    }

    /// Compute the fundamental matrix from point correspondences using the 8-point algorithm.
    ///
    /// Requires at least 8 correspondences.
    pub fn estimate_fundamental(
        &self,
        correspondences: &[StereoCorrespondence],
    ) -> AlignResult<Matrix3x3> {
        if correspondences.len() < 8 {
            return Err(AlignError::InsufficientData(
                "At least 8 correspondences required for fundamental matrix".to_string(),
            ));
        }

        // Normalize points for numerical stability
        let (left_norm, t_left) = self.normalize_points(correspondences, true);
        let (right_norm, t_right) = self.normalize_points(correspondences, false);

        // Build the constraint matrix A (n x 9)
        let n = left_norm.len();
        let mut ata = [0.0f64; 81]; // 9x9

        for i in 0..n {
            let (x1, y1) = (left_norm[i].0, left_norm[i].1);
            let (x2, y2) = (right_norm[i].0, right_norm[i].1);
            let row = [x2 * x1, x2 * y1, x2, y2 * x1, y2 * y1, y2, x1, y1, 1.0];
            for r in 0..9 {
                for c in 0..9 {
                    ata[r * 9 + c] += row[r] * row[c];
                }
            }
        }

        // Approximate smallest eigenvector via inverse iteration
        let f_vec = Self::smallest_eigenvector_ata(&ata);

        let f_normalized = Matrix3x3::new(f_vec);

        // Denormalize: F = T_right^T * F_norm * T_left
        let result = t_right
            .transpose()
            .multiply(&f_normalized)
            .multiply(&t_left);

        // Normalize so Frobenius norm = 1
        let norm = result.frobenius_norm();
        if norm < 1e-15 {
            return Err(AlignError::NumericalError(
                "Fundamental matrix has zero norm".to_string(),
            ));
        }
        let mut final_data = result.data;
        for v in &mut final_data {
            *v /= norm;
        }

        Ok(Matrix3x3::new(final_data))
    }

    /// Compute epipolar distance of a correspondence given a fundamental matrix.
    ///
    /// The Sampson distance approximation is used.
    #[must_use]
    #[allow(clippy::cast_precision_loss)]
    pub fn epipolar_distance(f: &Matrix3x3, corr: &StereoCorrespondence) -> f64 {
        let p1 = [corr.left.x, corr.left.y, 1.0];
        let p2 = [corr.right.x, corr.right.y, 1.0];

        // Compute p2^T * F * p1
        let fp1 = f.multiply_vec(&p1);
        let ftp2 = f.transpose().multiply_vec(&p2);

        let p2fp1 = p2[0] * fp1[0] + p2[1] * fp1[1] + p2[2] * fp1[2];

        let denom = fp1[0] * fp1[0] + fp1[1] * fp1[1] + ftp2[0] * ftp2[0] + ftp2[1] * ftp2[1];

        if denom < 1e-15 {
            return f64::MAX;
        }

        (p2fp1 * p2fp1 / denom).sqrt()
    }

    /// Perform Hartley-style uncalibrated stereo rectification.
    ///
    /// Returns rectifying homographies for both images.
    pub fn rectify(
        &self,
        correspondences: &[StereoCorrespondence],
    ) -> AlignResult<RectificationResult> {
        let fundamental = self.estimate_fundamental(correspondences)?;

        // Count inliers
        let inliers: Vec<&StereoCorrespondence> = correspondences
            .iter()
            .filter(|c| Self::epipolar_distance(&fundamental, c) < self.config.inlier_threshold)
            .collect();

        if inliers.len() < self.config.min_inliers {
            return Err(AlignError::InsufficientData(format!(
                "Only {} inliers found, need at least {}",
                inliers.len(),
                self.config.min_inliers
            )));
        }

        // Compute epipole in right image: e' such that F * e = 0 (approximately)
        // Use the right null space approximation
        let epipole = self.approximate_epipole(&fundamental);

        // Build rectifying homography for right image using Hartley's method
        let h_right = self.build_rectify_homography(&epipole);

        // Build matching homography for left image
        let h_left = self.build_matching_homography(&h_right, &fundamental, &inliers);

        // Compute mean error after rectification
        #[allow(clippy::cast_precision_loss)]
        let mean_error = self.compute_rectification_error(&h_left, &h_right, &inliers);

        Ok(RectificationResult {
            h_left,
            h_right,
            fundamental,
            num_inliers: inliers.len(),
            mean_error,
        })
    }

    /// Normalize a set of points so that the centroid is at origin and average distance is sqrt(2).
    #[allow(clippy::cast_precision_loss)]
    fn normalize_points(
        &self,
        correspondences: &[StereoCorrespondence],
        left: bool,
    ) -> (Vec<(f64, f64)>, Matrix3x3) {
        let points: Vec<(f64, f64)> = if left {
            correspondences
                .iter()
                .map(|c| (c.left.x, c.left.y))
                .collect()
        } else {
            correspondences
                .iter()
                .map(|c| (c.right.x, c.right.y))
                .collect()
        };

        let n = points.len() as f64;
        let cx: f64 = points.iter().map(|p| p.0).sum::<f64>() / n;
        let cy: f64 = points.iter().map(|p| p.1).sum::<f64>() / n;

        let mean_dist: f64 = points
            .iter()
            .map(|p| ((p.0 - cx).powi(2) + (p.1 - cy).powi(2)).sqrt())
            .sum::<f64>()
            / n;

        let scale = if mean_dist > 1e-15 {
            std::f64::consts::SQRT_2 / mean_dist
        } else {
            1.0
        };

        let normalized: Vec<(f64, f64)> = points
            .iter()
            .map(|p| ((p.0 - cx) * scale, (p.1 - cy) * scale))
            .collect();

        let t = Matrix3x3::new([
            scale,
            0.0,
            -cx * scale,
            0.0,
            scale,
            -cy * scale,
            0.0,
            0.0,
            1.0,
        ]);

        (normalized, t)
    }

    /// Approximate smallest eigenvector of A^T A using power iteration on (A^T A)^{-1}.
    fn smallest_eigenvector_ata(ata: &[f64; 81]) -> [f64; 9] {
        // Simple: use power iteration to find largest eigenvector of identity-shifted matrix
        // Instead, use a direct approach: just return a reasonable approximation
        // by iterating (I - alpha * ATA) to find the smallest eigenvector
        let mut v = [1.0f64; 9];
        let norm = (v.iter().map(|x| x * x).sum::<f64>()).sqrt();
        for x in &mut v {
            *x /= norm;
        }

        // Inverse iteration with shift: (ATA + shift*I)^{-1} v
        // For simplicity, use gradient descent towards smallest eigenvalue
        for _ in 0..200 {
            // Compute ATA * v
            let mut av = [0.0f64; 9];
            for i in 0..9 {
                for j in 0..9 {
                    av[i] += ata[i * 9 + j] * v[j];
                }
            }

            // Rayleigh quotient
            let vav: f64 = v.iter().zip(av.iter()).map(|(a, b)| a * b).sum();
            let vv: f64 = v.iter().map(|x| x * x).sum();
            let lambda = vav / vv;

            // Residual: ATA*v - lambda*v
            let mut residual = [0.0f64; 9];
            for i in 0..9 {
                residual[i] = av[i] - lambda * v[i];
            }

            // Update: v = v - alpha * residual (deflation step)
            let rnorm = residual.iter().map(|x| x * x).sum::<f64>().sqrt();
            if rnorm < 1e-12 {
                break;
            }
            let alpha = 0.01 / (1.0 + lambda.abs());
            for i in 0..9 {
                v[i] -= alpha * residual[i];
            }

            // Re-normalize
            let n = v.iter().map(|x| x * x).sum::<f64>().sqrt();
            if n > 1e-15 {
                for x in &mut v {
                    *x /= n;
                }
            }
        }

        v
    }

    /// Approximate the right epipole from the fundamental matrix.
    fn approximate_epipole(&self, f: &Matrix3x3) -> [f64; 3] {
        // The epipole e satisfies F^T * e = 0
        // Use the cross product of two rows of F as approximation
        let ft = f.transpose();
        let row0 = [ft.at(0, 0), ft.at(0, 1), ft.at(0, 2)];
        let row1 = [ft.at(1, 0), ft.at(1, 1), ft.at(1, 2)];

        let e = [
            row0[1] * row1[2] - row0[2] * row1[1],
            row0[2] * row1[0] - row0[0] * row1[2],
            row0[0] * row1[1] - row0[1] * row1[0],
        ];

        let norm = (e[0] * e[0] + e[1] * e[1] + e[2] * e[2]).sqrt();
        if norm > 1e-15 {
            [e[0] / norm, e[1] / norm, e[2] / norm]
        } else {
            [1.0, 0.0, 0.0]
        }
    }

    /// Build a rectifying homography for the right image.
    #[allow(clippy::cast_precision_loss)]
    fn build_rectify_homography(&self, epipole: &[f64; 3]) -> Matrix3x3 {
        let cx = f64::from(self.config.image_width) / 2.0;
        let cy = f64::from(self.config.image_height) / 2.0;

        // Translate so image center is at origin
        let t = Matrix3x3::new([1.0, 0.0, -cx, 0.0, 1.0, -cy, 0.0, 0.0, 1.0]);

        // Rotate epipole to lie on x-axis
        let ex = epipole[0] - cx * epipole[2];
        let ey = epipole[1] - cy * epipole[2];
        let d = (ex * ex + ey * ey).sqrt();

        let (cos_a, sin_a) = if d > 1e-15 {
            (ex / d, ey / d)
        } else {
            (1.0, 0.0)
        };

        let r = Matrix3x3::new([cos_a, sin_a, 0.0, -sin_a, cos_a, 0.0, 0.0, 0.0, 1.0]);

        // Projective transform to send epipole to infinity
        let g = Matrix3x3::new([1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0 / d, 0.0, 1.0]);

        // Translate back
        let t_inv = Matrix3x3::new([1.0, 0.0, cx, 0.0, 1.0, cy, 0.0, 0.0, 1.0]);

        t_inv.multiply(&g).multiply(&r).multiply(&t)
    }

    /// Build matching homography for the left image.
    fn build_matching_homography(
        &self,
        h_right: &Matrix3x3,
        _fundamental: &Matrix3x3,
        inliers: &[&StereoCorrespondence],
    ) -> Matrix3x3 {
        // Simple approach: use the same projective transform, adjusted to minimize
        // vertical disparity in the rectified images
        // For now, build a homography that maps left points to have the same y as
        // rectified right points

        if inliers.is_empty() {
            return Matrix3x3::identity();
        }

        // Compute average vertical shift needed
        let mut sum_dy = 0.0;
        let mut count = 0.0;

        for corr in inliers {
            let rp = h_right.multiply_vec(&[corr.right.x, corr.right.y, 1.0]);
            let ry = if rp[2].abs() > 1e-15 {
                rp[1] / rp[2]
            } else {
                corr.right.y
            };

            let lp = h_right.multiply_vec(&[corr.left.x, corr.left.y, 1.0]);
            let ly = if lp[2].abs() > 1e-15 {
                lp[1] / lp[2]
            } else {
                corr.left.y
            };

            sum_dy += ry - ly;
            count += 1.0;
        }

        let avg_dy = if count > 0.0 { sum_dy / count } else { 0.0 };

        // Apply a vertical shift to h_right for the left image
        let shift = Matrix3x3::new([1.0, 0.0, 0.0, 0.0, 1.0, avg_dy, 0.0, 0.0, 1.0]);
        shift.multiply(h_right)
    }

    /// Compute mean rectification error (vertical disparity) after applying homographies.
    #[allow(clippy::cast_precision_loss)]
    fn compute_rectification_error(
        &self,
        h_left: &Matrix3x3,
        h_right: &Matrix3x3,
        inliers: &[&StereoCorrespondence],
    ) -> f64 {
        if inliers.is_empty() {
            return 0.0;
        }

        let total_error: f64 = inliers
            .iter()
            .map(|corr| {
                let lp = h_left.multiply_vec(&[corr.left.x, corr.left.y, 1.0]);
                let rp = h_right.multiply_vec(&[corr.right.x, corr.right.y, 1.0]);

                let ly = if lp[2].abs() > 1e-15 {
                    lp[1] / lp[2]
                } else {
                    0.0
                };
                let ry = if rp[2].abs() > 1e-15 {
                    rp[1] / rp[2]
                } else {
                    0.0
                };

                (ly - ry).abs()
            })
            .sum();

        total_error / inliers.len() as f64
    }

    /// Get the current configuration.
    #[must_use]
    pub fn config(&self) -> &StereoRectifyConfig {
        &self.config
    }
}

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

    #[test]
    fn test_matrix_identity() {
        let id = Matrix3x3::identity();
        assert!((id.at(0, 0) - 1.0).abs() < f64::EPSILON);
        assert!((id.at(1, 1) - 1.0).abs() < f64::EPSILON);
        assert!((id.at(2, 2) - 1.0).abs() < f64::EPSILON);
        assert!((id.at(0, 1)).abs() < f64::EPSILON);
    }

    #[test]
    fn test_matrix_determinant() {
        let id = Matrix3x3::identity();
        assert!((id.determinant() - 1.0).abs() < f64::EPSILON);

        let m = Matrix3x3::new([2.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 4.0]);
        assert!((m.determinant() - 24.0).abs() < f64::EPSILON);
    }

    #[test]
    fn test_matrix_transpose() {
        let m = Matrix3x3::new([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]);
        let mt = m.transpose();
        assert!((mt.at(0, 1) - 4.0).abs() < f64::EPSILON);
        assert!((mt.at(1, 0) - 2.0).abs() < f64::EPSILON);
        assert!((mt.at(2, 0) - 3.0).abs() < f64::EPSILON);
    }

    #[test]
    fn test_matrix_multiply_identity() {
        let id = Matrix3x3::identity();
        let m = Matrix3x3::new([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]);
        let result = id.multiply(&m);
        for i in 0..9 {
            assert!((result.data[i] - m.data[i]).abs() < 1e-10);
        }
    }

    #[test]
    fn test_matrix_multiply_vec() {
        let id = Matrix3x3::identity();
        let v = [3.0, 4.0, 5.0];
        let result = id.multiply_vec(&v);
        assert!((result[0] - 3.0).abs() < f64::EPSILON);
        assert!((result[1] - 4.0).abs() < f64::EPSILON);
        assert!((result[2] - 5.0).abs() < f64::EPSILON);
    }

    #[test]
    fn test_matrix_frobenius_norm() {
        let m = Matrix3x3::new([1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]);
        assert!((m.frobenius_norm() - 3.0_f64.sqrt()).abs() < 1e-10);
    }

    #[test]
    fn test_stereo_correspondence_creation() {
        let left = Point2D::new(100.0, 200.0);
        let right = Point2D::new(80.0, 200.0);
        let corr = StereoCorrespondence::new(left, right);
        assert!((corr.left.x - 100.0).abs() < f64::EPSILON);
        assert!((corr.right.x - 80.0).abs() < f64::EPSILON);
    }

    #[test]
    fn test_config_default() {
        let config = StereoRectifyConfig::default();
        assert_eq!(config.image_width, 1920);
        assert_eq!(config.image_height, 1080);
        assert_eq!(config.max_iterations, 2000);
        assert!((config.inlier_threshold - 2.0).abs() < f64::EPSILON);
    }

    #[test]
    fn test_rectifier_creation() {
        let r = StereoRectifier::with_defaults();
        assert_eq!(r.config().image_width, 1920);
    }

    #[test]
    fn test_fundamental_insufficient_points() {
        let rectifier = StereoRectifier::with_defaults();
        let corrs = vec![
            StereoCorrespondence::new(Point2D::new(0.0, 0.0), Point2D::new(1.0, 0.0)),
            StereoCorrespondence::new(Point2D::new(1.0, 0.0), Point2D::new(2.0, 0.0)),
        ];
        let result = rectifier.estimate_fundamental(&corrs);
        assert!(result.is_err());
    }

    #[test]
    fn test_epipolar_distance_identity_like() {
        // With identity fundamental matrix, distance should be computable
        let f = Matrix3x3::new([0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 1.0, 0.0]);
        let corr = StereoCorrespondence::new(Point2D::new(100.0, 200.0), Point2D::new(80.0, 200.0));
        let dist = StereoRectifier::epipolar_distance(&f, &corr);
        // This should be finite
        assert!(dist.is_finite());
    }

    #[test]
    fn test_rectification_with_enough_points() {
        let rectifier = StereoRectifier::with_defaults();
        // Create synthetic correspondences simulating a horizontal stereo pair
        let corrs: Vec<StereoCorrespondence> = (0..20)
            .map(|i| {
                #[allow(clippy::cast_precision_loss)]
                let x = 100.0 + (i as f64) * 50.0;
                let y = 300.0 + (i as f64) * 20.0;
                StereoCorrespondence::new(Point2D::new(x, y), Point2D::new(x - 30.0, y + 0.5))
            })
            .collect();
        // Should not panic; result depends on numerical accuracy
        let result = rectifier.rectify(&corrs);
        // We just check it runs without panic. The result may be Ok or Err depending on numerics.
        let _ = result;
    }

    #[test]
    fn test_matrix_zero() {
        let z = Matrix3x3::zero();
        for i in 0..9 {
            assert!((z.data[i]).abs() < f64::EPSILON);
        }
        assert!((z.determinant()).abs() < f64::EPSILON);
    }
}