pub use num_traits::{Bounded, Float, One, Zero};
pub(crate) use std::ops::{Add, Div, Mul, Neg, Sub, AddAssign, SubAssign, MulAssign, DivAssign};
use crate::rgb::RGBf;
use crate::xyz::XYZ;
use std::ops::{Index, IndexMut};
pub fn clamp<T>(x: T, a: T, b: T) -> T
where
T: PartialOrd,
{
if x < a {
a
} else if x > b {
b
} else {
x
}
}
pub fn lerp(a: f32, b: f32, t: f32) -> f32 {
(1.0 - t) * a + t * b
}
#[inline(always)]
pub fn sqrt<T>(x: T) -> T where T: Real {
x.sqrt()
}
#[inline(always)]
pub fn sqr<T>(x: T) -> T where T: Real {
x * x
}
#[inline(always)]
pub fn abs<T>(x: T) -> T where T: Real {
x.abs()
}
#[inline(always)]
pub fn sin<T>(x: T) -> T where T: Real {
x.sin()
}
#[inline(always)]
pub fn asin<T>(x: T) -> T where T: Real {
x.asin()
}
#[inline(always)]
pub fn cos<T>(x: T) -> T where T: Real {
x.cos()
}
#[inline(always)]
pub fn acos<T>(x: T) -> T where T: Real {
x.acos()
}
#[inline(always)]
pub fn tan<T>(x: T) -> T where T: Real {
x.tan()
}
#[inline(always)]
pub fn atan2<T>(x: T, y: T) -> T where T: Real {
x.atan2(y)
}
#[inline(always)]
pub fn exp<T>(x: T) -> T where T: Real {
x.exp()
}
#[inline(always)]
pub fn pow<T>(x: T, y: T) -> T where T: Real {
x.powf(y)
}
#[inline(always)]
pub fn hypot<T>(x: T, y: T) -> T where T: Real {
x.hypot(y)
}
#[inline(always)]
pub fn powi<T>(x: T, i: i32) -> T where T: Real {
x.powi(i)
}
#[repr(C)]
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct Matrix33<T> where T: Real {
pub x: [T; 9],
}
pub type M3f64 = Matrix33<f64>;
pub type M3f32 = Matrix33<f32>;
impl<T> Matrix33<T> where T: Real {
pub fn make_identity() -> Matrix33<T> {
Matrix33::<T> {
x: [T::one(), T::zero(), T::zero(),
T::zero(), T::one(), T::zero(),
T::zero(), T::zero(), T::one()],
}
}
pub fn new(values: [T; 9]) -> Matrix33<T> {
Matrix33 { x: values }
}
pub fn transposed(&self) -> Matrix33<T> {
Matrix33 {
x: [
self[0][0], self[1][0], self[2][0], self[0][1], self[1][1],
self[2][1], self[0][2], self[1][2], self[2][2],
],
}
}
pub fn determinant(&self) -> T {
self[0][0] * (self[1][1] * self[2][2] - self[1][2] * self[2][1])
+ self[0][1] * (self[1][2] * self[2][0] - self[1][0] * self[2][2])
+ self[0][2] * (self[1][0] * self[2][1] - self[1][1] * self[2][0])
}
pub fn gj_inverse(&self) -> Option<Matrix33<T>> {
let mut mtx_t = self.clone();
let mut mtx_s = Matrix33::make_identity();
for i in 0..2 {
let mut pivot = i;
let mut pivot_size: T = self[i][i].abs();
#[allow(clippy::needless_range_loop)]
for j in (i + 1)..3 {
let tmp = self[j][i].abs();
if tmp > pivot_size {
pivot = j;
pivot_size = tmp;
}
}
if pivot_size == T::zero() {
return None;
}
if pivot != i {
#[allow(clippy::manual_swap)]
for j in 0..3 {
let tmp1 = mtx_t[i][j];
mtx_t[i][j] = mtx_t[pivot][j];
mtx_t[pivot][j] = tmp1;
let tmp2 = mtx_s[i][j];
mtx_s[i][j] = mtx_s[pivot][j];
mtx_s[pivot][j] = tmp2;
}
}
for j in (i + 1)..3 {
let f = mtx_t[j][i] / mtx_t[i][i];
for k in 0..3 {
let t = mtx_t[i][k];
let s = mtx_s[i][k];
mtx_t[j][k] -= f * t;
mtx_s[j][k] -= f * s;
}
}
}
for i in (0..3).rev() {
let f = mtx_t[i][i];
if f == T::zero() {
return None;
}
for j in 0..3 {
mtx_t[i][j] /= f;
mtx_s[i][j] /= f;
}
for j in 0..i {
let f = mtx_t[j][i];
for k in 0..3 {
let t = mtx_t[i][k];
let s = mtx_s[i][k];
mtx_t[j][k] -= f * t;
mtx_s[j][k] -= f * s;
}
}
}
Some(mtx_s)
}
pub fn inverse(self) -> Option<Matrix33<T>> {
if self[0][2] > T::epsilon()
|| self[1][2] > T::epsilon()
|| (self[2][2] - T::one()).abs() > T::epsilon()
{
let mut mtx_s = Matrix33::new([
self[1][1] * self[2][2] - self[2][1] * self[1][2],
self[2][1] * self[0][2] - self[0][1] * self[2][2],
self[0][1] * self[1][2] - self[1][1] * self[0][2],
self[2][0] * self[1][2] - self[1][0] * self[2][2],
self[0][0] * self[2][2] - self[2][0] * self[0][2],
self[1][0] * self[0][2] - self[0][0] * self[1][2],
self[1][0] * self[2][1] - self[2][0] * self[1][1],
self[2][0] * self[0][1] - self[0][0] * self[2][1],
self[0][0] * self[1][1] - self[1][0] * self[0][1],
]);
let r = self[0][0] * mtx_s[0][0]
+ self[0][1] * mtx_s[1][0]
+ self[0][2] * mtx_s[2][0];
if r.abs() >= T::one() {
for s in mtx_s.x.iter_mut() {
*s /= r;
}
} else {
let mr = r.abs() / T::min_positive_value();
for s in mtx_s.x.iter_mut() {
if mr > s.abs() {
*s /= r;
} else {
return None;
}
}
}
Some(mtx_s)
} else {
let mut mtx_s = Matrix33::new([
self[1][1],
-self[0][1],
T::zero(),
-self[1][0],
self[0][0],
T::zero(),
T::zero(),
T::zero(),
T::one(),
]);
let r = self[0][0] * self[1][1] - self[1][0] * self[0][1];
if r.abs() >= T::one() {
for s in mtx_s.x.iter_mut() {
*s /= r;
}
} else {
let mr = r.abs() / T::min_positive_value();
for s in mtx_s.x.iter_mut() {
if mr > s.abs() {
*s /= r;
} else {
return None;
}
}
}
mtx_s[2][0] = -self[2][0] * mtx_s[0][0] - self[2][1] * mtx_s[1][0];
mtx_s[2][1] = -self[2][0] * mtx_s[0][1] - self[2][1] * mtx_s[1][1];
Some(mtx_s)
}
}
}
impl From<M3f64> for M3f32 {
fn from(m: M3f64) -> M3f32 {
M3f32 {
x: [
m.x[0] as f32,
m.x[1] as f32,
m.x[2] as f32,
m.x[3] as f32,
m.x[4] as f32,
m.x[5] as f32,
m.x[6] as f32,
m.x[7] as f32,
m.x[8] as f32,
]
}
}
}
impl<T> Index<usize> for Matrix33<T> where T: Real {
type Output = [T];
fn index(&self, index: usize) -> &[T] {
let offset = index * 3;
&self.x[offset..(offset + 3)]
}
}
impl<T> IndexMut<usize> for Matrix33<T> where T: Real {
fn index_mut(&mut self, index: usize) -> &mut [T] {
let offset = index * 3;
&mut self.x[offset..(offset + 3)]
}
}
impl<T> Mul for Matrix33<T> where T: Real {
type Output = Self;
fn mul(self, rhs: Self) -> Self {
let mut m = Matrix33::<T>::new([T::zero(); 9]);
#[allow(clippy::needless_range_loop)]
for i in 0..3 {
for j in 0..3 {
for k in 0..3 {
m[i][j] += self[i][k] * rhs[k][j];
}
}
}
m
}
}
impl<T> Add<T> for Matrix33<T> where T: Real {
type Output = Self;
fn add(self, rhs: T) -> Self {
Matrix33::new([
self.x[0] + rhs,
self.x[1] + rhs,
self.x[2] + rhs,
self.x[3] + rhs,
self.x[4] + rhs,
self.x[5] + rhs,
self.x[6] + rhs,
self.x[7] + rhs,
self.x[8] + rhs,
])
}
}
impl<T> Sub<T> for Matrix33<T> where T: Real {
type Output = Self;
fn sub(self, rhs: T) -> Self {
Matrix33::new([
self.x[0] - rhs,
self.x[1] - rhs,
self.x[2] - rhs,
self.x[3] - rhs,
self.x[4] - rhs,
self.x[5] - rhs,
self.x[6] - rhs,
self.x[7] - rhs,
self.x[8] - rhs,
])
}
}
impl<T> Mul<T> for Matrix33<T> where T: Real {
type Output = Self;
fn mul(self, rhs: T) -> Self {
Matrix33::new([
self.x[0] * rhs,
self.x[1] * rhs,
self.x[2] * rhs,
self.x[3] * rhs,
self.x[4] * rhs,
self.x[5] * rhs,
self.x[6] * rhs,
self.x[7] * rhs,
self.x[8] * rhs,
])
}
}
impl<T> Div<T> for Matrix33<T> where T: Real {
type Output = Self;
fn div(self, rhs: T) -> Self {
Matrix33::new([
self.x[0] / rhs,
self.x[1] / rhs,
self.x[2] / rhs,
self.x[3] / rhs,
self.x[4] / rhs,
self.x[5] / rhs,
self.x[6] / rhs,
self.x[7] / rhs,
self.x[8] / rhs,
])
}
}
impl<T> Neg for Matrix33<T> where T: Real {
type Output = Self;
fn neg(self) -> Self {
Matrix33::new([
-self.x[0], -self.x[1], -self.x[2], -self.x[3], -self.x[4],
-self.x[5], -self.x[6], -self.x[7], -self.x[8],
])
}
}
impl<T> Mul<XYZ<T>> for Matrix33<T> where T: Real {
type Output = XYZ<T>;
fn mul(self, xyz: XYZ<T>) -> XYZ<T> {
XYZ::new(
self.x[0] * xyz.x + self.x[1] * xyz.y + self.x[2] * xyz.z,
self.x[3] * xyz.x + self.x[4] * xyz.y + self.x[5] * xyz.z,
self.x[6] * xyz.x + self.x[7] * xyz.y + self.x[8] * xyz.z,
)
}
}
impl<T> Mul<RGBf<T>> for Matrix33<T> where T: Real {
type Output = RGBf<T>;
fn mul(self, rgb: RGBf<T>) -> RGBf<T> {
RGBf::new(
self.x[0] * rgb.r + self.x[1] * rgb.g + self.x[2] * rgb.b,
self.x[3] * rgb.r + self.x[4] * rgb.g + self.x[5] * rgb.b,
self.x[6] * rgb.r + self.x[7] * rgb.g + self.x[8] * rgb.b,
)
}
}
pub trait Scalar:
Copy
+ Zero
+ One
+ Bounded
+ Neg<Output = Self>
+ PartialOrd
+ PartialEq
+ Add<Output = Self>
+ Sub<Output = Self>
+ Mul<Output = Self>
+ Div<Output = Self>
+ AddAssign
+ SubAssign
+ MulAssign
+ DivAssign
+ std::fmt::Display
+ std::fmt::Debug
{
}
pub trait Real: Scalar + Float {}
impl<T> Real for T where T: Scalar + Float {}
impl Scalar for f32 {}
impl Scalar for f64 {}