affine-rs 0.1.0

//! Two-dimensional affine transformations.
//!
//! An [`Affine`] stores the six effective coefficients of the augmented
//! matrix below. The final row is fixed and cannot be put into an invalid
//! state.
//!
//! ```text
//! | a  b  c |
//! | d  e  f |
//! | 0  0  1 |
//! ```
//!
//! # Example
//!
//! ```
//! use affine_rs::Affine;
//!
//! let transform = Affine::translation(10.0, 20.0)
//!     .compose(Affine::scale(2.0, 2.0));
//! assert_eq!(transform.transform_point([3.0, 4.0]), [16.0, 28.0]);
//! ```

mod error;

use core::fmt;
use core::ops::Mul;

pub use error::AffineError;

/// Default tolerance for approximate comparisons.
pub const DEFAULT_EPSILON: f64 = 1.0e-5;
const DEFAULT_EPSILON_SQUARED: f64 = DEFAULT_EPSILON * DEFAULT_EPSILON;

/// A two-dimensional affine transform.
///
/// The coefficients map an input point `[x, y]` to:
///
/// ```text
/// x' = a*x + b*y + c
/// y' = d*x + e*y + f
/// ```
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Affine {
    pub a: f64,
    pub b: f64,
    pub c: f64,
    pub d: f64,
    pub e: f64,
    pub f: f64,
}

impl Affine {
    /// The identity transform.
    pub const IDENTITY: Self = Self::new(1.0, 0.0, 0.0, 0.0, 1.0, 0.0);

    /// Creates a transform from its six effective coefficients.
    #[must_use]
    #[allow(clippy::many_single_char_names)]
    pub const fn new(a: f64, b: f64, c: f64, d: f64, e: f64, f: f64) -> Self {
        Self { a, b, c, d, e, f }
    }

    /// Creates a translation transform.
    #[must_use]
    pub const fn translation(x_offset: f64, y_offset: f64) -> Self {
        Self::new(1.0, 0.0, x_offset, 0.0, 1.0, y_offset)
    }

    /// Creates independent x- and y-axis scaling.
    #[must_use]
    pub const fn scale(x_scale: f64, y_scale: f64) -> Self {
        Self::new(x_scale, 0.0, 0.0, 0.0, y_scale, 0.0)
    }

    /// Creates uniform scaling on both axes.
    #[must_use]
    pub const fn uniform_scale(scale: f64) -> Self {
        Self::scale(scale, scale)
    }

    /// Creates x- and y-axis shear from angles in degrees.
    #[must_use]
    pub fn shear(x_angle_degrees: f64, y_angle_degrees: f64) -> Self {
        let x_shear = x_angle_degrees.to_radians().tan();
        let y_shear = y_angle_degrees.to_radians().tan();
        Self::new(1.0, x_shear, 0.0, y_shear, 1.0, 0.0)
    }

    /// Creates a counter-clockwise rotation around the origin.
    #[must_use]
    pub fn rotation(angle_degrees: f64) -> Self {
        let (cosine, sine) = cos_sin_degrees(angle_degrees);
        Self::new(cosine, -sine, 0.0, sine, cosine, 0.0)
    }

    /// Creates a counter-clockwise rotation around `pivot`.
    #[must_use]
    pub fn rotation_around(angle_degrees: f64, pivot: [f64; 2]) -> Self {
        let (cosine, sine) = cos_sin_degrees(angle_degrees);
        let [pivot_x, pivot_y] = pivot;
        Self::new(
            cosine,
            -sine,
            pivot_x - pivot_x * cosine + pivot_y * sine,
            sine,
            cosine,
            pivot_y - pivot_x * sine - pivot_y * cosine,
        )
    }

    /// Creates the non-identity permutation matrix for two dimensions.
    #[must_use]
    pub const fn permutation() -> Self {
        Self::new(0.0, 1.0, 0.0, 1.0, 0.0, 0.0)
    }

    /// Creates a transform from GDAL `GeoTransform` coefficient order.
    #[must_use]
    #[allow(clippy::many_single_char_names)]
    pub const fn from_gdal(coefficients: [f64; 6]) -> Self {
        let [c, a, b, f, d, e] = coefficients;
        Self::new(a, b, c, d, e, f)
    }

    /// Returns coefficients in GDAL `GeoTransform` order.
    #[must_use]
    pub const fn to_gdal(self) -> [f64; 6] {
        [self.c, self.a, self.b, self.f, self.d, self.e]
    }

    /// Returns coefficients in Shapely affine-transform order.
    #[must_use]
    pub const fn to_shapely(self) -> [f64; 6] {
        [self.a, self.b, self.d, self.e, self.c, self.f]
    }

    /// Returns the area scaling factor of the linear part.
    #[must_use]
    pub fn determinant(self) -> f64 {
        self.a.mul_add(self.e, -(self.b * self.d))
    }

    /// Returns `true` when all stored coefficients are finite.
    #[must_use]
    pub fn is_finite(self) -> bool {
        [self.a, self.b, self.c, self.d, self.e, self.f]
            .into_iter()
            .all(f64::is_finite)
    }

    /// Returns `true` when the determinant is exactly zero.
    #[must_use]
    pub fn is_degenerate(self) -> bool {
        self.determinant() == 0.0
    }

    /// Returns `true` when the transform preserves orientation.
    #[must_use]
    pub fn is_proper(self) -> bool {
        self.determinant() > 0.0
    }

    /// Returns the matrix as its three two-dimensional column vectors.
    #[must_use]
    pub const fn column_vectors(self) -> [[f64; 2]; 3] {
        [[self.a, self.d], [self.b, self.e], [self.c, self.f]]
    }

    /// Returns `true` when the transform is approximately the identity.
    #[must_use]
    pub fn is_identity(self) -> bool {
        self.approx_eq(Self::IDENTITY, DEFAULT_EPSILON)
    }

    /// Returns `true` when axis-aligned shapes remain axis-aligned.
    #[must_use]
    pub fn is_rectilinear(self) -> bool {
        (self.a.abs() < DEFAULT_EPSILON && self.e.abs() < DEFAULT_EPSILON)
            || (self.d.abs() < DEFAULT_EPSILON && self.b.abs() < DEFAULT_EPSILON)
    }

    /// Returns `true` when the linear transform preserves angles.
    #[must_use]
    pub fn is_conformal(self) -> bool {
        self.a.mul_add(self.b, self.d * self.e).abs() < DEFAULT_EPSILON
    }

    /// Returns `true` when the linear transform is orthonormal.
    #[must_use]
    pub fn is_orthonormal(self) -> bool {
        self.is_conformal()
            && (1.0 - self.a.mul_add(self.a, self.d * self.d)).abs() < DEFAULT_EPSILON
            && (1.0 - self.b.mul_add(self.b, self.e * self.e)).abs() < DEFAULT_EPSILON
    }

    /// Returns the two singular values, ordered from greatest to least.
    #[must_use]
    pub fn scaling(self) -> [f64; 2] {
        let trace = self.a * self.a + self.b * self.b + self.d * self.d + self.e * self.e;
        let determinant_squared = self.determinant().powi(2);
        let mut delta = trace.mul_add(trace / 4.0, -determinant_squared);
        if delta < DEFAULT_EPSILON_SQUARED {
            delta = 0.0;
        }
        let root = delta.sqrt();
        [(trace / 2.0 + root).sqrt(), (trace / 2.0 - root).sqrt()]
    }

    /// Returns the eccentricity induced by the linear transform.
    #[must_use]
    pub fn eccentricity(self) -> f64 {
        let [major, minor] = self.scaling();
        major.mul_add(major, -(minor * minor)).sqrt() / major
    }

    /// Returns the counter-clockwise rotation angle in degrees.
    ///
    /// # Errors
    ///
    /// Returns [`AffineError::UndefinedRotation`] for improper or degenerate
    /// transforms, and [`AffineError::NonFiniteTransform`] for non-finite input.
    pub fn rotation_angle(self) -> Result<f64, AffineError> {
        if !self.is_finite() {
            return Err(AffineError::NonFiniteTransform);
        }
        if !self.is_proper() || self.is_degenerate() {
            return Err(AffineError::UndefinedRotation);
        }
        let [major, _] = self.scaling();
        Ok((self.d / major).atan2(self.a / major).to_degrees())
    }

    /// Applies the transform to an `[x, y]` point.
    #[must_use]
    pub fn transform_point(self, [x, y]: [f64; 2]) -> [f64; 2] {
        [
            x.mul_add(self.a, y.mul_add(self.b, self.c)),
            x.mul_add(self.d, y.mul_add(self.e, self.f)),
        ]
    }

    /// Applies the transform to every point in a mutable slice.
    pub fn transform_points_in_place(self, points: &mut [[f64; 2]]) {
        for point in points {
            *point = self.transform_point(*point);
        }
    }

    /// Composes two transforms.
    ///
    /// `self.compose(rhs)` applies `rhs` first and `self` second.
    #[must_use]
    pub fn compose(self, rhs: Self) -> Self {
        Self::new(
            self.a.mul_add(rhs.a, self.b * rhs.d),
            self.a.mul_add(rhs.b, self.b * rhs.e),
            self.a.mul_add(rhs.c, self.b.mul_add(rhs.f, self.c)),
            self.d.mul_add(rhs.a, self.e * rhs.d),
            self.d.mul_add(rhs.b, self.e * rhs.e),
            self.d.mul_add(rhs.c, self.e.mul_add(rhs.f, self.f)),
        )
    }

    /// Returns the inverse transform.
    ///
    /// # Errors
    ///
    /// Returns [`AffineError::NonFiniteTransform`] for non-finite
    /// coefficients and [`AffineError::NonInvertibleTransform`] when the
    /// determinant is zero.
    pub fn inverse(self) -> Result<Self, AffineError> {
        if !self.is_finite() {
            return Err(AffineError::NonFiniteTransform);
        }

        let determinant = self.determinant();
        if determinant == 0.0 {
            return Err(AffineError::NonInvertibleTransform);
        }

        let inverse_determinant = determinant.recip();
        let a = self.e * inverse_determinant;
        let b = -self.b * inverse_determinant;
        let d = -self.d * inverse_determinant;
        let e = self.a * inverse_determinant;

        Ok(Self::new(
            a,
            b,
            -self.c.mul_add(a, self.f * b),
            d,
            e,
            -self.c.mul_add(d, self.f * e),
        ))
    }

    /// Compares all six coefficients using an absolute tolerance.
    #[must_use]
    pub fn approx_eq(self, other: Self, epsilon: f64) -> bool {
        let lhs = [self.a, self.b, self.c, self.d, self.e, self.f];
        let rhs = [other.a, other.b, other.c, other.d, other.e, other.f];
        lhs.into_iter()
            .zip(rhs)
            .all(|(left, right)| (left - right).abs() < epsilon)
    }

    /// Parses six coefficients in World File order.
    ///
    /// World Files store the center of the upper-left pixel. The returned
    /// transform uses its upper-left corner, so a half-pixel translation is
    /// applied automatically.
    ///
    /// # Errors
    ///
    /// Returns [`AffineError::InvalidWorldFile`] unless the input contains
    /// exactly six finite floating-point values.
    #[allow(clippy::many_single_char_names)]
    pub fn from_world_file(text: &str) -> Result<Self, AffineError> {
        let values = text
            .split_whitespace()
            .map(str::parse::<f64>)
            .collect::<Result<Vec<_>, _>>()
            .map_err(|_| AffineError::InvalidWorldFile)?;
        let [a, d, b, e, c, f] =
            <[f64; 6]>::try_from(values).map_err(|_| AffineError::InvalidWorldFile)?;
        let center = Self::new(a, b, c, d, e, f);
        if !center.is_finite() {
            return Err(AffineError::InvalidWorldFile);
        }
        Ok(center.compose(Self::translation(-0.5, -0.5)))
    }

    /// Serializes this transform in six-line World File order.
    #[must_use]
    pub fn to_world_file(self) -> String {
        let center = self.compose(Self::translation(0.5, 0.5));
        [center.a, center.d, center.b, center.e, center.c, center.f]
            .into_iter()
            .map(|value| format!("{value:?}"))
            .collect::<Vec<_>>()
            .join("\n")
            + "\n"
    }
}

impl Default for Affine {
    fn default() -> Self {
        Self::IDENTITY
    }
}

impl Mul for Affine {
    type Output = Self;

    fn mul(self, rhs: Self) -> Self::Output {
        self.compose(rhs)
    }
}

impl Mul<[f64; 2]> for Affine {
    type Output = [f64; 2];

    fn mul(self, rhs: [f64; 2]) -> Self::Output {
        self.transform_point(rhs)
    }
}

impl fmt::Display for Affine {
    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            formatter,
            "|{:.2},{:.2},{:.2}|\n|{:.2},{:.2},{:.2}|\n|0.00,0.00,1.00|",
            self.a, self.b, self.c, self.d, self.e, self.f
        )
    }
}

fn cos_sin_degrees(angle: f64) -> (f64, f64) {
    let normalized = angle.rem_euclid(360.0);
    match normalized {
        90.0 => (0.0, 1.0),
        180.0 => (-1.0, 0.0),
        270.0 => (0.0, -1.0),
        _ => {
            let (sine, cosine) = normalized.to_radians().sin_cos();
            (cosine, sine)
        }
    }
}