diff --git a/.github/workflows/nalgebra-ci-build.yml b/.github/workflows/nalgebra-ci-build.yml index 1acd93e15..b3bf13ebc 100644 --- a/.github/workflows/nalgebra-ci-build.yml +++ b/.github/workflows/nalgebra-ci-build.yml @@ -64,7 +64,7 @@ jobs: override: true - uses: actions/checkout@v4 - name: check - run: cargo check --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,rkyv-safe-deser,rayon; + run: cargo check --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,wide-tests,rkyv-safe-deser,rayon; test-nalgebra: runs-on: ubuntu-latest # env: @@ -81,7 +81,7 @@ jobs: override: true - uses: actions/checkout@v4 - name: test - run: cargo test --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,rkyv-safe-deser,rayon; + run: cargo test --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,wide-tests,rkyv-safe-deser,rayon; test-nalgebra-glm: runs-on: ubuntu-latest steps: diff --git a/Cargo.toml b/Cargo.toml index 6f31a54f5..d1dbb1e23 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -98,6 +98,7 @@ arbitrary = ["quickcheck"] proptest-support = ["proptest"] slow-tests = [] rkyv-safe-deser = ["rkyv-serialize", "rkyv/validation"] +wide-tests = ["simba/wide", "simba/rand"] [dependencies] nalgebra-macros = { version = "0.3.0", path = "nalgebra-macros", optional = true } @@ -168,7 +169,7 @@ required-features = ["compare"] name = "nalgebra_bench" harness = false path = "benches/lib.rs" -required-features = ["rand"] +required-features = ["rand", "wide-tests"] #[profile.bench] #opt-level = 0 @@ -176,6 +177,7 @@ required-features = ["rand"] [profile.bench] lto = true +codegen-units = 1 [package.metadata.docs.rs] # Enable all the features when building the docs on docs.rs diff --git a/benches/core/cg.rs b/benches/core/cg.rs new file mode 100644 index 000000000..e6b690697 --- /dev/null +++ b/benches/core/cg.rs @@ -0,0 +1,85 @@ +use na::{Matrix3, Matrix4, Orthographic3, Point2, Point3, Vector2, Vector3}; +use rand::{Rng, SeedableRng}; +use rand_isaac::IsaacRng; +use simba::simd::WideF32x4; + +#[path = "../common/macros.rs"] +mod macros; + +bench_binop_ref!( + mat3_transform_vector2, + Matrix3, + Vector2, + transform_vector +); +bench_binop_ref!( + mat4_transform_vector3, + Matrix4, + Vector3, + transform_vector +); +bench_binop_ref!( + mat3_transform_point2, + Matrix3, + Point2, + transform_point +); +bench_binop_ref!( + mat4_transform_point3, + Matrix4, + Point3, + transform_point +); + +bench_binop_ref!( + mat3_transform_vector2_x4wide, + Matrix3, + Vector2, + simd_transform_vector +); +bench_binop_ref!( + mat4_transform_vector3_x4wide, + Matrix4, + Vector3, + simd_transform_vector +); +bench_binop_ref!( + mat3_transform_point2_x4wide, + Matrix3, + Point2, + simd_transform_point +); +bench_binop_ref!( + mat4_transform_point3_x4wide, + Matrix4, + Point3, + simd_transform_point +); + +fn mat4_transform_vector3_no_division(bench: &mut criterion::Criterion) { + let mut rng = IsaacRng::seed_from_u64(0); + let orthographic = Orthographic3::from_fov( + rng.random_range(0.5..2.0), + rng.random_range(30.0..90.0), + rng.random_range(0.5..1.5), + rng.random_range(2.0..1000.0), + ) + .to_homogeneous(); + let vector = rng.random(); + bench.bench_function("mat4_transform_vector3_no_division", move |bh| { + bh.iter(|| orthographic.transform_vector(&vector)) + }); +} + +criterion_group!( + cg, + mat3_transform_vector2, + mat4_transform_vector3, + mat3_transform_point2, + mat4_transform_point3, + mat3_transform_vector2_x4wide, + mat4_transform_vector3_x4wide, + mat3_transform_point2_x4wide, + mat4_transform_point3_x4wide, + mat4_transform_vector3_no_division, +); diff --git a/benches/core/mod.rs b/benches/core/mod.rs index eda9ddaac..213442298 100644 --- a/benches/core/mod.rs +++ b/benches/core/mod.rs @@ -1,5 +1,7 @@ +pub use self::cg::cg; pub use self::matrix::matrix; pub use self::vector::vector; +mod cg; mod matrix; mod vector; diff --git a/benches/lib.rs b/benches/lib.rs index df866cca3..db2585fd9 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -21,6 +21,7 @@ fn reproducible_dmatrix(nrows: usize, ncols: usize) -> DMatrix { } criterion_main!( + core::cg, core::matrix, core::vector, geometry::quaternion, diff --git a/src/base/cg.rs b/src/base/cg.rs index 8e961e72c..392f20852 100644 --- a/src/base/cg.rs +++ b/src/base/cg.rs @@ -20,6 +20,7 @@ use crate::geometry::{ }; use simba::scalar::{ClosedAddAssign, ClosedMulAssign, RealField}; +use simba::simd::SimdRealField; /// # Translation and scaling in any dimension impl OMatrix @@ -402,6 +403,9 @@ where + Allocator, DimNameDiff>, { /// Transforms the given vector, assuming the matrix `self` uses homogeneous coordinates. + /// + /// Each component of the resulting vector is divided by the last component of the homogeneous + /// coordinates if it is not zero or returned unchanged otherwise. #[inline] pub fn transform_vector( &self, @@ -423,8 +427,40 @@ where } } +/// # Transformation of vectors and points +impl, S: Storage> SquareMatrix +where + DefaultAllocator: Allocator + + Allocator> + + Allocator, DimNameDiff>, +{ + /// Transforms the given vector, assuming the matrix `self` uses homogeneous coordinates. + /// + /// Each component of the resulting vector is divided by the last component of the homogeneous + /// coordinates if it is not zero or returned unchanged otherwise. + #[inline] + pub fn simd_transform_vector( + &self, + v: &OVector>, + ) -> OVector> { + let transform = self.generic_view( + (0, 0), + (DimNameDiff::::name(), DimNameDiff::::name()), + ); + let normalizer = + self.generic_view((D::DIM - 1, 0), (Const::<1>, DimNameDiff::::name())); + let n = normalizer.tr_dot(v); + + let n = n.clone().select(n.simd_ne(T::zero()), T::one()); + transform * (v / n) + } +} + impl, Const<3>>> SquareMatrix, S> { /// Transforms the given point, assuming the matrix `self` uses homogeneous coordinates. + /// + /// Each component of the resulting point is divided by the `z` component if it is not zero + /// or returned unchanged otherwise. #[inline] pub fn transform_point(&self, pt: &Point) -> Point { let transform = self.fixed_view::<2, 2>(0, 0); @@ -440,8 +476,28 @@ impl, Const<3>>> SquareMatrix, } } +impl, Const<3>>> SquareMatrix, S> { + /// Transforms the given point, assuming the matrix `self` uses homogeneous coordinates. + /// + /// Each component of the resulting point is divided by the `z` component if it is not zero + /// or returned unchanged otherwise. + #[inline] + pub fn simd_transform_point(&self, pt: &Point) -> Point { + let transform = self.fixed_view::<2, 2>(0, 0); + let translation = self.fixed_view::<2, 1>(0, 2); + let normalizer = self.fixed_view::<1, 2>(2, 0); + let n = normalizer.tr_dot(&pt.coords) + unsafe { self.get_unchecked((2, 2)).clone() }; + + let n = n.clone().select(n.simd_ne(T::zero()), T::one()); + (transform * pt + translation) / n + } +} + impl, Const<4>>> SquareMatrix, S> { /// Transforms the given point, assuming the matrix `self` uses homogeneous coordinates. + /// + /// Each component of the resulting vector is divided by the `w` component if it is not zero + /// or returned unchanged otherwise. #[inline] pub fn transform_point(&self, pt: &Point) -> Point { let transform = self.fixed_view::<3, 3>(0, 0); @@ -456,3 +512,20 @@ impl, Const<4>>> SquareMatrix, } } } + +impl, Const<4>>> SquareMatrix, S> { + /// Transforms the given point, assuming the matrix `self` uses homogeneous coordinates. + /// + /// Each component of the resulting vector is divided by the `w` component if it is not zero + /// or returned unchanged otherwise. + #[inline] + pub fn simd_transform_point(&self, pt: &Point) -> Point { + let transform = self.fixed_view::<3, 3>(0, 0); + let translation = self.fixed_view::<3, 1>(0, 3); + let normalizer = self.fixed_view::<1, 3>(3, 0); + let n = normalizer.tr_dot(&pt.coords) + unsafe { self.get_unchecked((3, 3)).clone() }; + + let n = n.clone().select(n.simd_ne(T::zero()), T::one()); + (transform * pt + translation) / n + } +} diff --git a/src/geometry/perspective.rs b/src/geometry/perspective.rs index af59f66e9..bef735542 100644 --- a/src/geometry/perspective.rs +++ b/src/geometry/perspective.rs @@ -233,6 +233,8 @@ impl Perspective3 { // TODO: when we get specialization, specialize the Mul impl instead. /// Projects a point. Faster than matrix multiplication. + /// + /// Each component of the resulting point will be NaN or Inf if `p.z` is zero. #[inline] #[must_use] pub fn project_point(&self, p: &Point3) -> Point3 { @@ -261,6 +263,13 @@ impl Perspective3 { // TODO: when we get specialization, specialize the Mul impl instead. /// Projects a vector. Faster than matrix multiplication. + /// + /// `x` and `y` components of the result is a projection as a 2D vector. + /// + /// Components `x` and `y` of the resulting vector will be NaN or Inf if `p.z` is zero. + /// + /// `z` component of the resulting vector always equals to `-self.m33` and does + /// not depend on the components of `p`. #[inline] #[must_use] pub fn project_vector(&self, p: &Vector) -> Vector3 @@ -271,7 +280,7 @@ impl Perspective3 { Vector3::new( self.matrix[(0, 0)].clone() * p[0].clone() * inverse_denom.clone(), self.matrix[(1, 1)].clone() * p[1].clone() * inverse_denom, - self.matrix[(2, 2)].clone(), + -self.matrix[(2, 2)].clone(), ) } diff --git a/tests/core/cg.rs b/tests/core/cg.rs index b061efbba..2f1233934 100644 --- a/tests/core/cg.rs +++ b/tests/core/cg.rs @@ -1,4 +1,9 @@ -use na::{Matrix3, Matrix4, Point2, Point3, Vector2, Vector3}; +use na::{Matrix3, Matrix4, Orthographic3, Perspective3, Point2, Point3, Vector2, Vector3}; + +#[cfg(feature = "wide-tests")] +use na::{Rotation3, SimdValue, Translation3}; +#[cfg(feature = "wide-tests")] +use simba::simd::WideF32x4; /// See Example 3.4 of "Graphics and Visualization: Principles & Algorithms" /// by Theoharis, Papaioannou, Platis, Patrikalakis. @@ -57,3 +62,148 @@ fn test_scaling_wrt_point_3() { assert!(result == expected); } + +#[test] +fn test_perspective_transform_vector() { + let vector = Vector3::new(1.0, 2.0, 3.0); + let perspective = Perspective3::new(2.0, 45.0, 1.0, 1000.0); + + let transformed = perspective.as_matrix().transform_vector(&vector); + + let multiplied = perspective.as_matrix() * vector.to_homogeneous(); + let multiplied = multiplied / multiplied.w; + + assert_relative_eq!(transformed, perspective.project_vector(&vector)); + assert_relative_eq!(transformed.push(1.0), multiplied); +} + +#[test] +fn test_perspective_simd_transform_vector() { + let vector = Vector3::new(1.0, 2.0, 3.0); + let perspective = Perspective3::new(2.0, 45.0, 1.0, 1000.0); + + assert_relative_eq!( + perspective.as_matrix().simd_transform_vector(&vector), + perspective.as_matrix().transform_vector(&vector) + ); +} + +#[test] +fn test_perspective_transform_point3() { + let point = Point3::new(1.0, 2.0, 3.0); + let perspective = Perspective3::new(2.0, 45.0, 1.0, 1000.0); + + let transformed = perspective.as_matrix().transform_point(&point); + + let multiplied = perspective.as_matrix() * point.to_homogeneous(); + let multiplied = multiplied / multiplied.w; + + assert_relative_eq!(transformed, perspective.project_point(&point)); + assert_relative_eq!(transformed.coords.push(1.0), multiplied); +} + +#[test] +fn test_perspective_simd_transform_point3() { + let point = Point3::new(1.0, 2.0, 3.0); + let perspective = Perspective3::new(2.0, 45.0, 1.0, 1000.0); + + assert_relative_eq!( + perspective.as_matrix().simd_transform_point(&point), + perspective.as_matrix().transform_point(&point) + ); +} + +#[test] +fn test_orthographic_transform_vector() { + let vector = Vector3::new(1.0, 2.0, 3.0); + let orthographic = Orthographic3::from_fov(2.0, 45.0, 1.0, 1000.0); + + let transformed = orthographic.as_matrix().transform_vector(&vector); + + let multiplied = orthographic.as_matrix() * vector.to_homogeneous(); + + assert_relative_eq!(transformed, orthographic.project_vector(&vector)); + assert_relative_eq!(transformed.push(0.0), multiplied); +} + +#[test] +fn test_orthographic_simd_transform_vector() { + let vector = Vector3::new(1.0, 2.0, 3.0); + let orthographic = Orthographic3::from_fov(2.0, 45.0, 1.0, 1000.0); + + assert_relative_eq!( + orthographic.as_matrix().simd_transform_vector(&vector), + orthographic.as_matrix().transform_vector(&vector) + ); +} + +#[test] +fn test_orthographic_transform_point3() { + let point = Point3::new(1.0, 2.0, 3.0); + let orthographic = Orthographic3::from_fov(2.0, 45.0, 1.0, 1000.0); + + let transformed = orthographic.as_matrix().transform_point(&point); + + let multiplied = orthographic.as_matrix() * point.to_homogeneous(); + let multiplied = multiplied / multiplied.w; + + assert_relative_eq!(transformed, orthographic.project_point(&point)); + assert_relative_eq!(transformed.coords.push(1.0), multiplied); +} + +#[test] +fn test_orthographic_simd_transform_point3() { + let point = Point3::new(1.0, 2.0, 3.0); + let orthographic = Orthographic3::from_fov(2.0, 45.0, 1.0, 1000.0); + + assert_relative_eq!( + orthographic.as_matrix().simd_transform_point(&point), + orthographic.as_matrix().transform_point(&point) + ); +} + +#[cfg(feature = "wide-tests")] +#[test] +fn test_transform_vector_x4wide() { + let v1 = Vector3::new(0.0, 1.0, 2.0); + let v2 = Vector3::new(3.0, 4.0, 5.0); + let v3 = Vector3::new(6.0, 7.0, 8.0); + let v4 = Vector3::new(9.0, 10.0, 11.0); + + let wide_v = Vector3::new( + WideF32x4::from_arr([v1.x, v2.x, v3.x, v4.x]), + WideF32x4::from_arr([v1.y, v2.y, v3.y, v4.y]), + WideF32x4::from_arr([v1.z, v2.z, v3.z, v4.z]), + ); + + let m1 = Perspective3::new(2.0, 45.0, 1.0, 1000.0).to_homogeneous(); + let m2 = Orthographic3::from_fov(2.0, 45.0, 1.0, 1000.0).to_homogeneous(); + let m3 = Rotation3::from_axis_angle(&Vector3::y_axis(), 2.5).to_homogeneous(); + let m4 = Translation3::new(1.0, 2.0, 3.0).to_homogeneous(); + + let wide_m = Matrix4::new( + WideF32x4::from_arr([m1.m11, m2.m11, m3.m11, m4.m11]), + WideF32x4::from_arr([m1.m12, m2.m12, m3.m12, m4.m12]), + WideF32x4::from_arr([m1.m13, m2.m13, m3.m13, m4.m13]), + WideF32x4::from_arr([m1.m14, m2.m14, m3.m14, m4.m14]), + WideF32x4::from_arr([m1.m21, m2.m21, m3.m21, m4.m21]), + WideF32x4::from_arr([m1.m22, m2.m22, m3.m22, m4.m22]), + WideF32x4::from_arr([m1.m23, m2.m23, m3.m23, m4.m23]), + WideF32x4::from_arr([m1.m24, m2.m24, m3.m24, m4.m24]), + WideF32x4::from_arr([m1.m31, m2.m31, m3.m31, m4.m31]), + WideF32x4::from_arr([m1.m32, m2.m32, m3.m32, m4.m32]), + WideF32x4::from_arr([m1.m33, m2.m33, m3.m33, m4.m33]), + WideF32x4::from_arr([m1.m34, m2.m34, m3.m34, m4.m34]), + WideF32x4::from_arr([m1.m41, m2.m41, m3.m41, m4.m41]), + WideF32x4::from_arr([m1.m42, m2.m42, m3.m42, m4.m42]), + WideF32x4::from_arr([m1.m43, m2.m43, m3.m43, m4.m43]), + WideF32x4::from_arr([m1.m44, m2.m44, m3.m44, m4.m44]), + ); + + let wide_v = wide_m.simd_transform_vector(&wide_v); + + assert_eq!(wide_v.extract(0), m1.transform_vector(&v1)); + assert_eq!(wide_v.extract(1), m2.transform_vector(&v2)); + assert_eq!(wide_v.extract(2), m3.transform_vector(&v3)); + assert_eq!(wide_v.extract(3), m4.transform_vector(&v4)); +} diff --git a/tests/geometry/projection.rs b/tests/geometry/projection.rs index 1e0c9fd51..487199209 100644 --- a/tests/geometry/projection.rs +++ b/tests/geometry/projection.rs @@ -35,7 +35,7 @@ fn perspective_matrix_point_transformation() { #[cfg(feature = "proptest-support")] mod proptest_tests { - use na::{Orthographic3, Perspective3}; + use na::{Orthographic3, Perspective3, Point3}; use crate::proptest::*; use proptest::{prop_assert, proptest}; @@ -60,5 +60,31 @@ mod proptest_tests { prop_assert!(relative_eq!(pt, unprojected, epsilon = 1.0e-7)) } + + #[test] + fn perspective_project_vector(pt in point3(), vec in vector2()) { + let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0); + + let proj_pt = proj.project_point(&pt); + let proj_vec = proj.project_vector(&vec.push(pt.z)); + let proj_pt2 = proj.project_point(&(pt + vec.push(0.0))); + + let proj_pt_plus_proj_vec = Point3::from((proj_pt.xy() + proj_vec.xy()).coords.push(proj_pt.z)); + + prop_assert!(relative_eq!(proj_pt_plus_proj_vec, proj_pt2, epsilon = 1.0e-7)) + } + + #[test] + fn orthographic_project_vector(pt in point3(), vec in vector3()) { + let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0); + + let proj_pt = proj.project_point(&pt); + let proj_vec = proj.project_vector(&vec); + let proj_pt2 = proj.project_point(&(pt + vec)); + + let proj_pt_plus_proj_vec = proj_pt + proj_vec; + + prop_assert!(relative_eq!(proj_pt_plus_proj_vec, proj_pt2, epsilon = 1.0e-7)) + } } }