diff --git a/Cargo.toml b/Cargo.toml index 37a43717..398f0db9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -19,7 +19,6 @@ path = "src/lib.rs" arbitrary = [ "quickcheck" ] [dependencies] -rustc-serialize = "0.3" typenum = "1.4" generic-array = "0.2" rand = "0.3" @@ -27,8 +26,13 @@ num-traits = "0.1" num-complex = "0.1" approx = "0.1" alga = "0.4" +serde = "0.9" +serde_derive = "0.9" # clippy = "*" [dependencies.quickcheck] optional = true version = "0.3" + +[dev-dependencies] +serde_json = "0.9" diff --git a/Makefile b/Makefile index f9586c02..00457c11 100644 --- a/Makefile +++ b/Makefile @@ -1,11 +1,11 @@ all: - cargo build --features "arbitrary" + CARGO_INCREMENTAL=1 cargo build --features "arbitrary" doc: - cargo doc + CARGO_INCREMENTAL=1 cargo doc bench: cargo bench test: - cargo test --features "arbitrary" + CARGO_INCREMENTAL=1 cargo test --features "arbitrary" diff --git a/README.md b/README.md index c860659c..7d6a99fa 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,6 @@ +

+ crates.io +

crates.io @@ -8,75 +11,6 @@

- Documentation | Forum + Users guide | Documentation | Forum

- -nalgebra -======== - -**nalgebra** is a low-dimensional linear algebra library written for Rust targeting: - -* General-purpose linear algebra (still lacks a lot of features…) -* Real time computer graphics. -* Real time computer physics. - -## Using **nalgebra** -You will need the last stable build of the [rust compiler](http://www.rust-lang.org) -and the official package manager: [cargo](https://github.com/rust-lang/cargo). - -Simply add the following to your `Cargo.toml` file: - -```.ignore -[dependencies] -nalgebra = "0.10.*" -``` - - -All the functionality of **nalgebra** is grouped in one place: the root module `nalgebra::`. This -module re-exports everything and includes free functions for all traits methods performing -out-of-place operations. - -Thus, you can import the whole prelude using: - -```.ignore -use nalgebra::*; -``` - -However, the recommended way to use **nalgebra** is to import types and traits -explicitly, and call free-functions using the `na::` prefix: - -```.rust -extern crate nalgebra as na; -use na::{Vector3, Rotation3, Rotation}; - -fn main() { - let a = Vector3::new(1.0f64, 1.0, 1.0); - let mut b = Rotation3::new(na::zero()); - - b.append_rotation_mut(&a); - - assert!(na::approx_eq(&na::rotation(&b), &a)); -} -``` - - -## Features -**nalgebra** is meant to be a general-purpose, low-dimensional, linear algebra library, with -an optimized set of tools for computer graphics and physics. Those features include: - -* Vectors with predefined static sizes: `Vector1`, `Vector2`, `Vector3`, `Vector4`, `Vector5`, `Vector6`. -* Vector with a user-defined static size: `VectorN` (available only with the `generic_sizes` feature). -* Points with static sizes: `Point1`, `Point2`, `Point3`, `Point4`, `Point5`, `Point6`. -* Square matrices with static sizes: `Matrix1`, `Matrix2`, `Matrix3`, `Matrix4`, `Matrix5`, `Matrix6 `. -* Rotation matrices: `Rotation2`, `Rotation3` -* Quaternions: `Quaternion`, `Unit`. -* Unit-sized values (unit vectors, unit quaternions, etc.): `Unit`, e.g., `Unit>`. -* Isometries (translation ⨯ rotation): `Isometry2`, `Isometry3` -* Similarity transformations (translation ⨯ rotation ⨯ uniform scale): `Similarity2`, `Similarity3`. -* 3D projections for computer graphics: `Persp3`, `PerspMatrix3`, `Ortho3`, `OrthoMatrix3`. -* Dynamically sized heap-allocated vector: `DVector`. -* Dynamically sized stack-allocated vectors with a maximum size: `DVector1` to `DVector6`. -* Dynamically sized heap-allocated (square or rectangular) matrix: `DMatrix`. -* Linear algebra and data analysis operators: `Covariance`, `Mean`, `qr`, `cholesky`. -* Almost one trait per functionality: useful for generic programming. diff --git a/examples/dimensional_genericity.rs b/examples/dimensional_genericity.rs new file mode 100644 index 00000000..df22b945 --- /dev/null +++ b/examples/dimensional_genericity.rs @@ -0,0 +1,66 @@ +extern crate alga; +extern crate nalgebra as na; + +use alga::general::Real; +use alga::linear::FiniteDimInnerSpace; +use na::{Unit, ColumnVector, OwnedColumnVector, Vector2, Vector3}; +use na::storage::Storage; +use na::dimension::{DimName, U1}; + +/// Reflects a vector wrt. the hyperplane with normal `plane_normal`. +fn reflect_wrt_hyperplane_with_algebraic_genericity(plane_normal: &Unit, vector: &V) -> V + where V: FiniteDimInnerSpace + Copy { + let n = plane_normal.as_ref(); // Get the underlying vector of type `V`. + *vector - *n * (n.dot(vector) * na::convert(2.0)) +} + + +/// Reflects a vector wrt. the hyperplane with normal `plane_normal`. +fn reflect_wrt_hyperplane_with_structural_genericity(plane_normal: &Unit>, + vector: &ColumnVector) + -> OwnedColumnVector + where N: Real, + D: DimName, + S: Storage { + let n = plane_normal.as_ref(); // Get the underlying V. + vector - n * (n.dot(vector) * na::convert(2.0)) +} + +/// Reflects a 2D vector wrt. the 2D line with normal `plane_normal`. +fn reflect_wrt_hyperplane2(plane_normal: &Unit>, + vector: &Vector2) + -> Vector2 + where N: Real { + let n = plane_normal.as_ref(); // Get the underlying Vector2 + vector - n * (n.dot(vector) * na::convert(2.0)) +} + +/// Reflects a 3D vector wrt. the 3D plane with normal `plane_normal`. +/// /!\ This is an exact replicate of `reflect_wrt_hyperplane2, but for 3D. +fn reflect_wrt_hyperplane3(plane_normal: &Unit>, + vector: &Vector3) + -> Vector3 + where N: Real { + let n = plane_normal.as_ref(); // Get the underlying Vector3 + vector - n * (n.dot(vector) * na::convert(2.0)) +} + + +fn main() { + let plane2 = Vector2::y_axis(); // 2D plane normal. + let plane3 = Vector3::y_axis(); // 3D plane normal. + + let v2 = Vector2::new(1.0, 2.0); // 2D vector to be reflected. + let v3 = Vector3::new(1.0, 2.0, 3.0); // 3D vector to be reflected. + + // We can call the same function for 2D and 3D. + assert_eq!(reflect_wrt_hyperplane_with_algebraic_genericity(&plane2, &v2).y, -2.0); + assert_eq!(reflect_wrt_hyperplane_with_algebraic_genericity(&plane3, &v3).y, -2.0); + + assert_eq!(reflect_wrt_hyperplane_with_structural_genericity(&plane2, &v2).y, -2.0); + assert_eq!(reflect_wrt_hyperplane_with_structural_genericity(&plane3, &v3).y, -2.0); + + // Call each specific implementation depending on the dimension. + assert_eq!(reflect_wrt_hyperplane2(&plane2, &v2).y, -2.0); + assert_eq!(reflect_wrt_hyperplane3(&plane3, &v3).y, -2.0); +} diff --git a/examples/homogeneous_coordinates.rs b/examples/homogeneous_coordinates.rs new file mode 100644 index 00000000..a18b5a72 --- /dev/null +++ b/examples/homogeneous_coordinates.rs @@ -0,0 +1,45 @@ +#[macro_use] +extern crate approx; +extern crate nalgebra as na; + +use std::f32; +use na::{Vector2, Point2, Isometry2}; + + +fn use_dedicated_types() { + let iso = Isometry2::new(Vector2::new(1.0, 1.0), f32::consts::PI); + let pt = Point2::new(1.0, 0.0); + let vec = Vector2::x(); + + let transformed_pt = iso * pt; + let transformed_vec = iso * vec; + + assert_relative_eq!(transformed_pt, Point2::new(0.0, 1.0)); + assert_relative_eq!(transformed_vec, Vector2::new(-1.0, 0.0)); +} + +fn use_homogeneous_coordinates() { + let iso = Isometry2::new(Vector2::new(1.0, 1.0), f32::consts::PI); + let pt = Point2::new(1.0, 0.0); + let vec = Vector2::x(); + + // Compute using homogeneous coordinates. + let hom_iso = iso.to_homogeneous(); + let hom_pt = pt.to_homogeneous(); + let hom_vec = vec.to_homogeneous(); + + let hom_transformed_pt = hom_iso * hom_pt; + let hom_transformed_vec = hom_iso * hom_vec; + + // Convert back to the cartesian coordinates. + let transformed_pt = Point2::from_homogeneous(hom_transformed_pt).unwrap(); + let transformed_vec = Vector2::from_homogeneous(hom_transformed_vec).unwrap(); + + assert_relative_eq!(transformed_pt, Point2::new(0.0, 1.0)); + assert_relative_eq!(transformed_vec, Vector2::new(-1.0, 0.0)); +} + +fn main() { + use_dedicated_types(); + use_homogeneous_coordinates(); +} diff --git a/examples/identity.rs b/examples/identity.rs new file mode 100644 index 00000000..1805bf1e --- /dev/null +++ b/examples/identity.rs @@ -0,0 +1,42 @@ +extern crate alga; +extern crate nalgebra as na; + + +use alga::linear::Transformation; +use na::{Id, Vector3, Point3, Isometry3}; + +/* + * Applies `n` times the transformation `t` to the vector `v` and sum each + * intermediate value. + */ +fn complicated_algorithm(v: &Vector3, t: &T, n: usize) -> Vector3 + where T: Transformation> { + + let mut result = *v; + + // Do lots of operations involving t. + for _ in 0 .. n { + result = v + t.transform_vector(&result); + } + + result +} + + +/* + * The two following calls are equivalent in term of result. + */ +fn main() { + let v = Vector3::new(1.0, 2.0, 3.0); + + // The specialization generated by the compiler will do vector additions only. + let result1 = complicated_algorithm(&v, &Id::new(), 100000); + + // The specialization generated by the compiler will also include matrix multiplications. + let iso = Isometry3::identity(); + let result2 = complicated_algorithm(&v, &iso, 100000); + + // They both return the same result. + assert!(result1 == Vector3::new(100001.0, 200002.0, 300003.0)); + assert!(result2 == Vector3::new(100001.0, 200002.0, 300003.0)); +} diff --git a/examples/matrix_construction.rs b/examples/matrix_construction.rs new file mode 100644 index 00000000..4eb4ad06 --- /dev/null +++ b/examples/matrix_construction.rs @@ -0,0 +1,62 @@ +extern crate nalgebra as na; + +use na::{Vector2, RowVector3, Matrix2x3, DMatrix}; + + +fn main() { + // All the following matrices are equal but constructed in different ways. + let m = Matrix2x3::new(1.1, 1.2, 1.3, + 2.1, 2.2, 2.3); + + let m1 = Matrix2x3::from_rows(&[ + RowVector3::new(1.1, 1.2, 1.3), + RowVector3::new(2.1, 2.2, 2.3) + ]); + + let m2 = Matrix2x3::from_columns(&[ + Vector2::new(1.1, 2.1), + Vector2::new(1.2, 2.2), + Vector2::new(1.3, 2.3) + ]); + + let m3 = Matrix2x3::from_row_slice(&[ + 1.1, 1.2, 1.3, + 2.1, 2.2, 2.3 + ]); + + let m4 = Matrix2x3::from_column_slice(&[ + 1.1, 2.1, + 1.2, 2.2, + 1.3, 2.3 + ]); + + let m5 = Matrix2x3::from_fn(|r, c| (r + 1) as f32 + (c + 1) as f32 / 10.0); + + let m6 = Matrix2x3::from_iterator([ 1.1f32, 2.1, 1.2, 2.2, 1.3, 2.3 ].iter().cloned()); + + assert_eq!(m, m1); assert_eq!(m, m2); assert_eq!(m, m3); + assert_eq!(m, m4); assert_eq!(m, m5); assert_eq!(m, m6); + + // All the following matrices are equal but constructed in different ways. + // This time, we used a dynamically-sized matrix to show the extra arguments + // for the matrix shape. + let dm = DMatrix::from_row_slice(4, 3, &[ + 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0, + 0.0, 0.0, 0.0 + ]); + + let dm1 = DMatrix::from_diagonal_element(4, 3, 1.0); + let dm2 = DMatrix::identity(4, 3); + let dm3 = DMatrix::from_fn(4, 3, |r, c| if r == c { 1.0 } else { 0.0 }); + let dm4 = DMatrix::from_iterator(4, 3, [ + // Components listed column-by-column. + 1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0 + ].iter().cloned()); + + assert_eq!(dm, dm1); assert_eq!(dm, dm2); + assert_eq!(dm, dm3); assert_eq!(dm, dm4); +} diff --git a/examples/mvp.rs b/examples/mvp.rs new file mode 100644 index 00000000..7a38c671 --- /dev/null +++ b/examples/mvp.rs @@ -0,0 +1,28 @@ +#![allow(unused_variables)] + +extern crate nalgebra as na; + +use na::{Vector3, Point3, Isometry3, Perspective3}; + +fn main() { + // Our object is translated along the x axis. + let model = Isometry3::new(Vector3::x(), na::zero()); + + // Our camera looks toward the point (1.0, 0.0, 0.0). + // It is located at (0.0, 0.0, 1.0). + let eye = Point3::new(0.0, 0.0, 1.0); + let target = Point3::new(1.0, 0.0, 0.0); + let view = Isometry3::look_at_rh(&eye, &target, &Vector3::y()); + + // A perspective projection. + let projection = Perspective3::new(16.0 / 9.0, 3.14 / 2.0, 1.0, 1000.0); + + // The combination of the model with the view is still an isometry. + let model_view = model * view; + + // Convert everything to a `Matrix4` so that they can be combined. + let mat_model_view = model_view.to_homogeneous(); + + // Combine everything. + let model_view_projection = projection.as_matrix() * mat_model_view; +} diff --git a/examples/raw_pointer.rs b/examples/raw_pointer.rs new file mode 100644 index 00000000..ec960da4 --- /dev/null +++ b/examples/raw_pointer.rs @@ -0,0 +1,35 @@ +extern crate nalgebra as na; + +use na::{Vector3, Point3, Matrix3}; + +fn main() { + let v = Vector3::new(1.0f32, 0.0, 1.0); + let p = Point3::new(1.0f32, 0.0, 1.0); + let m = na::one::>(); + + // Convert to arrays. + let v_array = v.as_slice(); + let p_array = p.coords.as_slice(); + let m_array = m.as_slice(); + + // Get data pointers. + let v_pointer = v_array.as_ptr(); + let p_pointer = p_array.as_ptr(); + let m_pointer = m_array.as_ptr(); + + /* Then pass the raw pointers to some graphics API. */ + + unsafe { + assert_eq!(*v_pointer, 1.0); + assert_eq!(*v_pointer.offset(1), 0.0); + assert_eq!(*v_pointer.offset(2), 1.0); + + assert_eq!(*p_pointer, 1.0); + assert_eq!(*p_pointer.offset(1), 0.0); + assert_eq!(*p_pointer.offset(2), 1.0); + + assert_eq!(*m_pointer, 1.0); + assert_eq!(*m_pointer.offset(4), 1.0); + assert_eq!(*m_pointer.offset(8), 1.0); + } +} diff --git a/examples/scalar_genericity.rs b/examples/scalar_genericity.rs new file mode 100644 index 00000000..156d93ab --- /dev/null +++ b/examples/scalar_genericity.rs @@ -0,0 +1,34 @@ +extern crate alga; +extern crate nalgebra as na; + +use alga::general::{RingCommutative, Real}; +use na::{Vector3, Scalar}; + +fn print_vector(m: &Vector3) { + println!("{:?}", m) +} + +fn print_squared_norm(v: &Vector3) { + // NOTE: alternatively, nalgebra already defines `v.squared_norm()`. + let sqnorm = v.dot(v); + + println!("{:?}", sqnorm); +} + +fn print_norm(v: &Vector3) { + // NOTE: alternatively, nalgebra already defines `v.norm()`. + let norm = v.dot(v).sqrt(); + + // The Real bound implies that N is Display so we can + // use "{}" instead of "{:?}" for the format string. + println!("{}", norm) +} + +fn main() { + let v1 = Vector3::new(1, 2, 3); + let v2 = Vector3::new(1.0, 2.0, 3.0); + + print_vector(&v1); + print_squared_norm(&v1); + print_norm(&v2); +} diff --git a/examples/screen_to_view_coords.rs b/examples/screen_to_view_coords.rs new file mode 100644 index 00000000..4d6bc351 --- /dev/null +++ b/examples/screen_to_view_coords.rs @@ -0,0 +1,24 @@ +#![allow(unused_variables)] + +extern crate nalgebra as na; + +use na::{Point2, Point3, Perspective3, Unit}; + + +fn main() { + let projection = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0); + let screen_point = Point2::new(10.0f32, 20.0); + + // Compute two points in clip-space. + // "ndc" = normalized device coordinates. + let near_ndc_point = Point3::new(screen_point.x / 800.0, screen_point.y / 600.0, -1.0); + let far_ndc_point = Point3::new(screen_point.x / 800.0, screen_point.y / 600.0, 1.0); + + // Unproject them to view-space. + let near_view_point = projection.unproject_point(&near_ndc_point); + let far_view_point = projection.unproject_point(&far_ndc_point); + + // Compute the view-space line parameters. + let line_location = near_view_point; + let line_direction = Unit::new_normalize(far_view_point - near_view_point); +} diff --git a/examples/transform_conversion.rs b/examples/transform_conversion.rs new file mode 100644 index 00000000..37b266c4 --- /dev/null +++ b/examples/transform_conversion.rs @@ -0,0 +1,23 @@ +extern crate nalgebra as na; + +use na::{Vector2, Isometry2, Similarity2}; + +fn main() { + // Isometry -> Similarity conversion always succeeds. + let iso = Isometry2::new(Vector2::new(1.0f32, 2.0), na::zero()); + let _: Similarity2 = na::convert(iso); + + // Similarity -> Isometry conversion fails if the scaling factor is not 1.0. + let sim_without_scaling = Similarity2::new(Vector2::new(1.0f32, 2.0), 3.14, 1.0); + let sim_with_scaling = Similarity2::new(Vector2::new(1.0f32, 2.0), 3.14, 2.0); + + let iso_success: Option> = na::try_convert(sim_without_scaling); + let iso_fail: Option> = na::try_convert(sim_with_scaling); + + assert!(iso_success.is_some()); + assert!(iso_fail.is_none()); + + // Similarity -> Isometry conversion can be forced at your own risks. + let iso_forced: Isometry2 = unsafe { na::convert_unchecked(sim_with_scaling) }; + assert_eq!(iso_success.unwrap(), iso_forced); +} diff --git a/examples/transform_matrix4.rs b/examples/transform_matrix4.rs new file mode 100644 index 00000000..7d473b6b --- /dev/null +++ b/examples/transform_matrix4.rs @@ -0,0 +1,39 @@ +#[macro_use] +extern crate approx; +extern crate alga; +extern crate nalgebra as na; + +use alga::linear::Transformation; +use na::{Vector3, Point3, Matrix4}; + + +fn main() { + // Create a uniform scaling matrix with scaling factor 2. + let mut m = Matrix4::new_scaling(2.0); + + assert_eq!(m.transform_vector(&Vector3::x()), Vector3::x() * 2.0); + assert_eq!(m.transform_vector(&Vector3::y()), Vector3::y() * 2.0); + assert_eq!(m.transform_vector(&Vector3::z()), Vector3::z() * 2.0); + + // Append a nonuniform scaling in-place. + m.append_nonuniform_scaling_mut(&Vector3::new(1.0, 2.0, 3.0)); + + assert_eq!(m.transform_vector(&Vector3::x()), Vector3::x() * 2.0); + assert_eq!(m.transform_vector(&Vector3::y()), Vector3::y() * 4.0); + assert_eq!(m.transform_vector(&Vector3::z()), Vector3::z() * 6.0); + + // Append a translation out-of-place. + let m2 = m.append_translation(&Vector3::new(42.0, 0.0, 0.0)); + + assert_eq!(m2.transform_point(&Point3::new(1.0, 1.0, 1.0)), Point3::new(42.0 + 2.0, 4.0, 6.0)); + + // Create rotation. + let rot = Matrix4::from_scaled_axis(&Vector3::x() * 3.14); + let rot_then_m = m * rot; // Right-multiplication is equivalent to prepending `rot` to `m`. + let m_then_rot = rot * m; // Left-multiplication is equivalent to appending `rot` to `m`. + + let pt = Point3::new(1.0, 2.0, 3.0); + + assert_relative_eq!(m.transform_point(&rot.transform_point(&pt)), rot_then_m.transform_point(&pt)); + assert_relative_eq!(rot.transform_point(&m.transform_point(&pt)), m_then_rot.transform_point(&pt)); +} diff --git a/examples/transform_vector_point.rs b/examples/transform_vector_point.rs new file mode 100644 index 00000000..f94b88e3 --- /dev/null +++ b/examples/transform_vector_point.rs @@ -0,0 +1,20 @@ +#[macro_use] +extern crate approx; +extern crate nalgebra as na; + +use std::f32; +use na::{Vector2, Point2, Isometry2}; + +fn main() { + let t = Isometry2::new(Vector2::new(1.0, 1.0), f32::consts::PI); + let p = Point2::new(1.0, 0.0); // Will be affected by te rotation and the translation. + let v = Vector2::x(); // Will *not* be affected by the translation. + + assert_relative_eq!(t * p, Point2::new(-1.0 + 1.0, 1.0)); + // ^^^^ │ ^^^^^^^^ + // rotated │ translated + + assert_relative_eq!(t * v, Vector2::new(-1.0, 0.0)); + // ^^^^^ + // rotated only +} diff --git a/examples/transform_vector_point3.rs b/examples/transform_vector_point3.rs new file mode 100644 index 00000000..f1a010bf --- /dev/null +++ b/examples/transform_vector_point3.rs @@ -0,0 +1,27 @@ +extern crate alga; +extern crate nalgebra as na; + +use alga::linear::Transformation; +use na::{Vector3, Vector4, Point3, Matrix4}; + + +fn main() { + let mut m = Matrix4::new_rotation_wrt_point(Vector3::x() * 1.57, Point3::new(1.0, 2.0, 1.0)); + m.append_scaling_mut(2.0); + + let point1 = Point3::new(2.0, 3.0, 4.0); + let homogeneous_point2 = Vector4::new(2.0, 3.0, 4.0, 1.0); + + // First option: use the dedicated `.transform_point(...)` method. + let transformed_point1 = m.transform_point(&point1); + // Second option: use the homogeneous coordinates of the point. + let transformed_homogeneous_point2 = m * homogeneous_point2; + + // Recover the 3D point from its 4D homogeneous coordinates. + let transformed_point2 = Point3::from_homogeneous(transformed_homogeneous_point2); + + // Check that transforming the 3D point with the `.transform_point` method is + // indeed equivalent to multiplying its 4D homogeneous coordinates by the 4x4 + // matrix. + assert_eq!(transformed_point1, transformed_point2.unwrap()); +} diff --git a/examples/transformation_pointer.rs b/examples/transformation_pointer.rs new file mode 100644 index 00000000..68d35e09 --- /dev/null +++ b/examples/transformation_pointer.rs @@ -0,0 +1,25 @@ +extern crate nalgebra as na; + +use na::{Vector3, Isometry3}; + +fn main() { + let iso = Isometry3::new(Vector3::new(1.0f32, 0.0, 1.0), na::zero()); + + // Compute the homogeneous coordinates first. + let iso_matrix = iso.to_homogeneous(); + let iso_array = iso_matrix.as_slice(); + let iso_pointer = iso_array.as_ptr(); + + /* Then pass the raw pointer to some graphics API. */ + + unsafe { + assert_eq!(*iso_pointer, 1.0); + assert_eq!(*iso_pointer.offset(5), 1.0); + assert_eq!(*iso_pointer.offset(10), 1.0); + assert_eq!(*iso_pointer.offset(15), 1.0); + + assert_eq!(*iso_pointer.offset(12), 1.0); + assert_eq!(*iso_pointer.offset(13), 0.0); + assert_eq!(*iso_pointer.offset(14), 1.0); + } +} diff --git a/examples/unit_wrapper.rs b/examples/unit_wrapper.rs new file mode 100644 index 00000000..dc898c50 --- /dev/null +++ b/examples/unit_wrapper.rs @@ -0,0 +1,30 @@ +extern crate nalgebra as na; + +use na::{Unit, Vector3}; + +fn length_on_direction_with_unit(v: &Vector3, dir: &Unit>) -> f32 { + // No need to normalize `dir`: we know that it is non-zero and normalized. + v.dot(dir.as_ref()) +} + + + +fn length_on_direction_without_unit(v: &Vector3, dir: &Vector3) -> f32 { + // Obligatory normalization of the direction vector (and test, for robustness). + if let Some(unit_dir) = dir.try_normalize(1.0e-6) { + v.dot(&unit_dir) + } + else { + // Normalization failed because the norm was too small. + panic!("Invalid input direction.") + } +} + +fn main() { + let v = Vector3::new(1.0, 2.0, 3.0); + + let l1 = length_on_direction_with_unit(&v, &Vector3::y_axis()); + let l2 = length_on_direction_without_unit(&v, &Vector3::y()); + + assert_eq!(l1, l2) +} diff --git a/src/core/alias.rs b/src/core/alias.rs index e4395dcd..df296257 100644 --- a/src/core/alias.rs +++ b/src/core/alias.rs @@ -19,52 +19,88 @@ pub type MatrixNM = Matrix>; /// A staticaly sized column-major square matrix with `D` rows and columns. pub type MatrixN = MatrixNM; +/// A stack-allocated, column-major, 1x1 square matrix. pub type Matrix1 = MatrixN; +/// A stack-allocated, column-major, 2x2 square matrix. pub type Matrix2 = MatrixN; +/// A stack-allocated, column-major, 3x3 square matrix. pub type Matrix3 = MatrixN; +/// A stack-allocated, column-major, 4x4 square matrix. pub type Matrix4 = MatrixN; +/// A stack-allocated, column-major, 5x5 square matrix. pub type Matrix5 = MatrixN; +/// A stack-allocated, column-major, 6x6 square matrix. pub type Matrix6 = MatrixN; +/// A stack-allocated, column-major, 1x2 square matrix. pub type Matrix1x2 = MatrixNM; +/// A stack-allocated, column-major, 1x3 square matrix. pub type Matrix1x3 = MatrixNM; +/// A stack-allocated, column-major, 1x4 square matrix. pub type Matrix1x4 = MatrixNM; +/// A stack-allocated, column-major, 1x5 square matrix. pub type Matrix1x5 = MatrixNM; +/// A stack-allocated, column-major, 1x6 square matrix. pub type Matrix1x6 = MatrixNM; +/// A stack-allocated, column-major, 2x3 square matrix. pub type Matrix2x3 = MatrixNM; +/// A stack-allocated, column-major, 2x4 square matrix. pub type Matrix2x4 = MatrixNM; +/// A stack-allocated, column-major, 2x5 square matrix. pub type Matrix2x5 = MatrixNM; +/// A stack-allocated, column-major, 2x6 square matrix. pub type Matrix2x6 = MatrixNM; +/// A stack-allocated, column-major, 3x4 square matrix. pub type Matrix3x4 = MatrixNM; +/// A stack-allocated, column-major, 3x5 square matrix. pub type Matrix3x5 = MatrixNM; +/// A stack-allocated, column-major, 3x6 square matrix. pub type Matrix3x6 = MatrixNM; +/// A stack-allocated, column-major, 4x5 square matrix. pub type Matrix4x5 = MatrixNM; +/// A stack-allocated, column-major, 4x6 square matrix. pub type Matrix4x6 = MatrixNM; +/// A stack-allocated, column-major, 5x6 square matrix. pub type Matrix5x6 = MatrixNM; +/// A stack-allocated, column-major, 2x1 square matrix. pub type Matrix2x1 = MatrixNM; +/// A stack-allocated, column-major, 3x1 square matrix. pub type Matrix3x1 = MatrixNM; +/// A stack-allocated, column-major, 4x1 square matrix. pub type Matrix4x1 = MatrixNM; +/// A stack-allocated, column-major, 5x1 square matrix. pub type Matrix5x1 = MatrixNM; +/// A stack-allocated, column-major, 6x1 square matrix. pub type Matrix6x1 = MatrixNM; +/// A stack-allocated, column-major, 3x2 square matrix. pub type Matrix3x2 = MatrixNM; +/// A stack-allocated, column-major, 4x2 square matrix. pub type Matrix4x2 = MatrixNM; +/// A stack-allocated, column-major, 5x2 square matrix. pub type Matrix5x2 = MatrixNM; +/// A stack-allocated, column-major, 6x2 square matrix. pub type Matrix6x2 = MatrixNM; +/// A stack-allocated, column-major, 4x3 square matrix. pub type Matrix4x3 = MatrixNM; +/// A stack-allocated, column-major, 5x3 square matrix. pub type Matrix5x3 = MatrixNM; +/// A stack-allocated, column-major, 6x3 square matrix. pub type Matrix6x3 = MatrixNM; +/// A stack-allocated, column-major, 5x4 square matrix. pub type Matrix5x4 = MatrixNM; +/// A stack-allocated, column-major, 6x4 square matrix. pub type Matrix6x4 = MatrixNM; +/// A stack-allocated, column-major, 6x5 square matrix. pub type Matrix6x5 = MatrixNM; @@ -81,11 +117,17 @@ pub type DVector = Matrix>; /// A statically sized D-dimensional column vector. pub type VectorN = MatrixNM; +/// A stack-allocated, 1-dimensional column vector. pub type Vector1 = VectorN; +/// A stack-allocated, 2-dimensional column vector. pub type Vector2 = VectorN; +/// A stack-allocated, 3-dimensional column vector. pub type Vector3 = VectorN; +/// A stack-allocated, 4-dimensional column vector. pub type Vector4 = VectorN; +/// A stack-allocated, 5-dimensional column vector. pub type Vector5 = VectorN; +/// A stack-allocated, 6-dimensional column vector. pub type Vector6 = VectorN; @@ -102,9 +144,15 @@ pub type RowDVector = Matrix>; /// A statically sized D-dimensional row vector. pub type RowVectorN = MatrixNM; +/// A stack-allocated, 1-dimensional row vector. pub type RowVector1 = RowVectorN; +/// A stack-allocated, 2-dimensional row vector. pub type RowVector2 = RowVectorN; +/// A stack-allocated, 3-dimensional row vector. pub type RowVector3 = RowVectorN; +/// A stack-allocated, 4-dimensional row vector. pub type RowVector4 = RowVectorN; +/// A stack-allocated, 5-dimensional row vector. pub type RowVector5 = RowVectorN; +/// A stack-allocated, 6-dimensional row vector. pub type RowVector6 = RowVectorN; diff --git a/src/core/allocator.rs b/src/core/allocator.rs index 23bfd356..d0c0a416 100644 --- a/src/core/allocator.rs +++ b/src/core/allocator.rs @@ -1,3 +1,5 @@ +//! Abstract definition of a matrix data storage allocator. + use std::any::Any; use core::Scalar; @@ -15,6 +17,7 @@ use core::storage::{Storage, OwnedStorage}; /// Every allocator must be both static and dynamic. Though not all implementations may share the /// same `Buffer` type. pub trait Allocator: Any + Sized { + /// The type of buffer this allocator can instanciate. type Buffer: OwnedStorage; /// Allocates a buffer with the given number of rows and columns without initializing its content. diff --git a/src/core/cg.rs b/src/core/cg.rs index cfd40ac3..e409d006 100644 --- a/src/core/cg.rs +++ b/src/core/cg.rs @@ -25,7 +25,7 @@ impl SquareMatrix #[inline] pub fn new_scaling(scaling: N) -> Self { let mut res = Self::from_diagonal_element(scaling); - res[(D::dim(), D::dim())] = N::one(); + res[(D::dim() - 1, D::dim() - 1)] = N::one(); res } @@ -50,7 +50,7 @@ impl SquareMatrix SB: Storage, U1>, S::Alloc: Allocator, U1> { let mut res = Self::one(); - res.fixed_slice_mut::, U1>(0, D::dim()).copy_from(translation); + res.fixed_slice_mut::, U1>(0, D::dim() - 1).copy_from(translation); res } @@ -311,7 +311,7 @@ impl SquareMatrix S::Alloc: Allocator, U1> { for i in 0 .. D::dim() { for j in 0 .. D::dim() - 1 { - self[(j, i)] += shift[i] * self[(D::dim(), j)]; + self[(j, i)] += shift[j] * self[(D::dim() - 1, i)]; } } } @@ -324,12 +324,12 @@ impl SquareMatrix S::Alloc: Allocator, U1> + Allocator, DimNameDiff> + Allocator> { - let scale = self.fixed_slice::>(D::dim(), 0).tr_dot(&shift); + let scale = self.fixed_slice::>(D::dim() - 1, 0).tr_dot(&shift); let post_translation = self.fixed_slice::, DimNameDiff>(0, 0) * shift; - self[(D::dim(), D::dim())] += scale; + self[(D::dim() - 1, D::dim() - 1)] += scale; - let mut translation = self.fixed_slice_mut::, U1>(0, D::dim()); + let mut translation = self.fixed_slice_mut::, U1>(0, D::dim() - 1); translation += post_translation; } } @@ -349,7 +349,7 @@ impl Transformation, SB>> for Squa fn transform_vector(&self, v: &ColumnVector, SB>) -> ColumnVector, SB> { let transform = self.fixed_slice::, DimNameDiff>(0, 0); - let normalizer = self.fixed_slice::>(D::dim(), 0); + let normalizer = self.fixed_slice::>(D::dim() - 1, 0); let n = normalizer.tr_dot(&v); if !n.is_zero() { @@ -363,9 +363,9 @@ impl Transformation, SB>> for Squa fn transform_point(&self, pt: &PointBase, SB>) -> PointBase, SB> { let transform = self.fixed_slice::, DimNameDiff>(0, 0); - let translation = self.fixed_slice::, U1>(0, D::dim()); - let normalizer = self.fixed_slice::>(D::dim(), 0); - let n = normalizer.tr_dot(&pt.coords) + unsafe { *self.get_unchecked(D::dim(), D::dim()) }; + let translation = self.fixed_slice::, U1>(0, D::dim() - 1); + let normalizer = self.fixed_slice::>(D::dim() - 1, 0); + let n = normalizer.tr_dot(&pt.coords) + unsafe { *self.get_unchecked(D::dim() - 1, D::dim() - 1) }; if !n.is_zero() { return transform * (pt / n) + translation; diff --git a/src/core/constraint.rs b/src/core/constraint.rs index 182bf2ae..0f017e5b 100644 --- a/src/core/constraint.rs +++ b/src/core/constraint.rs @@ -1,6 +1,8 @@ +//! Compatibility constraints between matrix shapes, e.g., for addition or multiplication. + use core::dimension::{Dim, DimName, Dynamic}; -/// A type for enforcing constraints. +/// A type used in `where` clauses for enforcing constraints. pub struct ShapeConstraint; /// Constraints `C1` and `R2` to be equivalent. @@ -14,9 +16,12 @@ where ShapeConstraint: DimEq { } macro_rules! equality_trait_decl( - ($($Trait: ident),* $(,)*) => {$( + ($($doc: expr, $Trait: ident),* $(,)*) => {$( // XXX: we can't do something like `DimEq for D2` because we would require a blancket impl… + #[doc = $doc] pub trait $Trait { + /// This is either equal to `D1` or `D2`, always choosing the one (if any) which is a type-level + /// constant. type Representative: Dim; } @@ -34,11 +39,26 @@ macro_rules! equality_trait_decl( )*} ); -equality_trait_decl!(DimEq, SameNumberOfRows, SameNumberOfColumns); +equality_trait_decl!( + "Constraints `D1` and `D2` to be equivalent.", + DimEq, -/// Constraints D1 and D2 to be equivalent, where the both designates dimensions of algebraic + "Constraints `D1` and `D2` to be equivalent. \ + They are both assumed to be the number of \ + rows of a matrix.", + SameNumberOfRows, + + "Constraints `D1` and `D2` to be equivalent. \ + They are both assumed to be the number of \ + columns of a matrix.", + SameNumberOfColumns +); + +/// Constraints D1 and D2 to be equivalent, where they both designate dimensions of algebraic /// entities (e.g. square matrices). pub trait SameDimension: SameNumberOfRows + SameNumberOfColumns { + /// This is either equal to `D1` or `D2`, always choosing the one (if any) which is a type-level + /// constant. type Representative: Dim; } diff --git a/src/core/construction.rs b/src/core/construction.rs index db586a6c..606f0a2f 100644 --- a/src/core/construction.rs +++ b/src/core/construction.rs @@ -116,6 +116,10 @@ impl> Matrix res } + /// Builds a new matrix from its rows. + /// + /// Panics if not enough rows are provided (for statically-sized matrices), or if all rows do + /// not have the same dimensions. #[inline] pub fn from_rows(rows: &[Matrix]) -> Matrix where SB: Storage { @@ -127,13 +131,18 @@ impl> Matrix if C::try_to_usize().is_none() { assert!(rows.iter().all(|r| r.len() == ncols), - "The rows provided must all have the same dimension."); + "The provided rows must all have the same dimension."); } // FIXME: optimize that. Self::from_fn_generic(R::from_usize(nrows), C::from_usize(ncols), |i, j| rows[i][(0, j)]) } + + /// Builds a new matrix from its columns. + /// + /// Panics if not enough columns are provided (for statically-sized matrices), or if all + /// columns do not have the same dimensions. #[inline] pub fn from_columns(columns: &[ColumnVector]) -> Matrix where SB: Storage { @@ -239,12 +248,17 @@ macro_rules! impl_constructors( Self::from_fn_generic($($gargs, )* f) } + /// Creates an identity matrix. If the matrix is not square, the largest square + /// submatrix (starting at the first row and column) is set to the identity while all + /// other entries are set to zero. #[inline] pub fn identity($($args: usize,)*) -> Matrix where N: Zero + One { Self::identity_generic($($gargs),* ) } + /// Creates a matrix filled with its diagonal filled with `elt` and all other + /// components set to zero. #[inline] pub fn from_diagonal_element($($args: usize,)* elt: N) -> Matrix where N: Zero + One { diff --git a/src/core/coordinates.rs b/src/core/coordinates.rs index 6fff9193..3a8f371f 100644 --- a/src/core/coordinates.rs +++ b/src/core/coordinates.rs @@ -1,3 +1,9 @@ +#![allow(missing_docs)] + +//! Structures to which matrices and vector can be auto-dereferenced (through `Deref`) to access +//! components using their names. For example, if `v` is a 3D vector, one can write `v.z` instead +//! of `v[2]`. + use std::mem; use std::ops::{Deref, DerefMut}; @@ -14,8 +20,10 @@ use core::allocator::OwnedAllocator; macro_rules! coords_impl( ($T: ident; $($comps: ident),*) => { + /// Data structure used to provide access to matrix and vector coordinates with the dot + /// notation, e.g., `v.x` is the same as `v[0]` for a vector. #[repr(C)] - #[derive(Eq, PartialEq, RustcEncodable, RustcDecodable, Clone, Hash, Debug, Copy)] + #[derive(Eq, PartialEq, Clone, Hash, Debug, Copy, Serialize, Deserialize)] pub struct $T { $(pub $comps: N),* } diff --git a/src/core/decompositions.rs b/src/core/decompositions.rs index 7c19bd21..41183fcb 100644 --- a/src/core/decompositions.rs +++ b/src/core/decompositions.rs @@ -272,14 +272,7 @@ impl> SquareMatrix { /// /// Matrix symmetricness is not checked. Returns `None` if `self` is not definite positive. #[inline] - pub fn cholesky(&self) -> Option> { - let out = self.transpose(); - self.do_cholesky(out).ok() - } - - /// Cholesky decomposition G of a square symmetric positive definite matrix A, such that A = G * G^T - #[inline] - pub fn cholesky_unchecked(&self) -> Result, &'static str> { + pub fn cholesky(&self) -> Result, &'static str> { let out = self.transpose(); if !out.relative_eq(self, N::default_epsilon(), N::default_max_relative()) { @@ -289,6 +282,13 @@ impl> SquareMatrix { self.do_cholesky(out) } + /// Cholesky decomposition G of a square symmetric positive definite matrix A, such that A = G * G^T + #[inline] + pub fn cholesky_unchecked(&self) -> Result, &'static str> { + let out = self.transpose(); + self.do_cholesky(out) + } + #[inline(always)] fn do_cholesky(&self, mut out: OwnedSquareMatrix) -> Result, &'static str> { diff --git a/src/core/default_allocator.rs b/src/core/default_allocator.rs index 41cdf50d..acfbe97c 100644 --- a/src/core/default_allocator.rs +++ b/src/core/default_allocator.rs @@ -1,3 +1,8 @@ +//! The default matrix data storage allocator. +//! +//! This will use stack-allocated buffers for matrices with dimensions known at compile-time, and +//! heap-allocated buffers for matrices with at least one dimension unknown at compile-time. + use std::mem; use std::ops::Mul; diff --git a/src/core/dimension.rs b/src/core/dimension.rs index 1f723ad4..20819b68 100644 --- a/src/core/dimension.rs +++ b/src/core/dimension.rs @@ -1,10 +1,14 @@ +#![allow(missing_docs)] + +//! Traits and tags for identifying the dimension of all algebraic entities. + use std::fmt::Debug; use std::any::Any; use std::ops::{Add, Sub, Mul, Div}; use typenum::{self, Unsigned, UInt, B1, Bit, UTerm, Sum, Prod, Diff, Quot}; /// Dim of dynamically-sized algebraic entities. -#[derive(Clone, Copy, Eq, PartialEq, Debug)] +#[derive(Clone, Copy, Eq, PartialEq, Debug, Serialize, Deserialize)] pub struct Dynamic { value: usize } @@ -19,17 +23,27 @@ impl Dynamic { } } -/// Trait implemented by dynamic dimensions. +/// Trait implemented by `Dynamic`. pub trait IsDynamic { } -/// Trait implemented by dimensions that are not equal to U1. +/// Trait implemented by `Dynamic` and type-level integers different from `U1`. pub trait IsNotStaticOne { } impl IsDynamic for Dynamic { } impl IsNotStaticOne for Dynamic { } +/// Trait implemented by any type that can be used as a dimension. This includes type-level +/// integers and `Dynamic` (for dimensions not known at compile-time). pub trait Dim: Any + Debug + Copy + PartialEq + Send { + /// Gets the compile-time value of `Self`. Returns `None` if it is not known, i.e., if `Self = + /// Dynamic`. fn try_to_usize() -> Option; + + /// Gets the run-time value of `self`. For type-level integers, this is the same as + /// `Self::try_to_usize().unwrap()`. fn value(&self) -> usize; + + /// Builds an instance of `Self` from a run-time value. Panics if `Self` is a type-level + /// integer and `dim != Self::try_to_usize().unwrap()`. fn from_usize(dim: usize) -> Self; } @@ -127,6 +141,7 @@ dim_ops!( ); +/// Trait implemented exclusively by type-level integers. pub trait DimName: Dim { type Value: NamedDim; @@ -146,7 +161,7 @@ pub trait NamedDim: Sized + Any + Unsigned { type Name: DimName; } -#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)] +#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)] pub struct U1; impl Dim for U1 { @@ -182,7 +197,7 @@ impl NamedDim for typenum::U1{ macro_rules! named_dimension( ($($D: ident),* $(,)*) => {$( - #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)] + #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)] pub struct $D; impl Dim for $D { diff --git a/src/core/iter.rs b/src/core/iter.rs index 50af9da5..3c99a4ae 100644 --- a/src/core/iter.rs +++ b/src/core/iter.rs @@ -1,3 +1,5 @@ +//! Matrix iterators. + use std::marker::PhantomData; use std::mem; diff --git a/src/core/matrix.rs b/src/core/matrix.rs index cde94f38..c92b76a8 100644 --- a/src/core/matrix.rs +++ b/src/core/matrix.rs @@ -7,9 +7,9 @@ use std::any::TypeId; use std::mem; use approx::ApproxEq; -use alga::general::{Field, Real}; +use alga::general::{Ring, Real}; -use core::Scalar; +use core::{Scalar, Unit}; use core::dimension::{Dim, DimAdd, DimSum, U1, U2}; use core::constraint::{ShapeConstraint, SameNumberOfRows, SameNumberOfColumns}; use core::iter::{MatrixIter, MatrixIterMut}; @@ -55,10 +55,37 @@ pub type MatrixTrMul = Matrix = Matrix>::Alloc as Allocator>::Buffer>; +/// The most generic column-major matrix (and vector) type. +/// +/// It combines four type parameters: +/// - `N`: for the matrix components scalar type. +/// - `R`: for the matrix number of rows. +/// - `C`: for the matrix number of columns. +/// - `S`: for the matrix data storage, i.e., the buffer that actually contains the matrix +/// components. +/// +/// The matrix dimensions parameters `R` and `C` can either be: +/// - type-level unsigned integer contants (e.g. `U1`, `U124`) from the `nalgebra::` root module. +/// All numbers from 0 to 127 are defined that way. +/// - type-level unsigned integer constants (e.g. `U1024`, `U10000`) from the `typenum::` crate. +/// Using those, you will not get error messages as nice as for numbers smaller than 128 defined on +/// the `nalgebra::` module. +/// - the special value `Dynamic` from the `nalgebra::` root module. This indicates that the +/// specified dimension is not known at compile-time. Note that this will generally imply that the +/// matrix data storage `S` performs a dynamic allocation and contains extra metadata for the +/// matrix shape. +/// +/// Note that mixing `Dynamic` with type-level unsigned integers is allowed. Actually, a +/// dynamically-sized column vector should be represented as a `Matrix` (given +/// some concrete types for `N` and a compatible data storage type `S`). #[repr(C)] -#[derive(RustcEncodable, RustcDecodable, Hash, Debug, Clone, Copy)] +#[derive(Serialize, Deserialize, Hash, Debug, Clone, Copy)] pub struct Matrix { + /// The data storage that contains all the matrix components and informations about its number + /// of rows and column (if needed). pub data: S, + + #[serde(skip_serializing, skip_deserializing)] _phantoms: PhantomData<(N, R, C)> } @@ -149,7 +176,7 @@ impl> Matrix { nrows * ncols } - /// The shape (number of rows, number of columns) of this matrix. + /// The shape of this matrix returned as the tuple (number of rows, number of columns). #[inline] pub fn shape(&self) -> (usize, usize) { let (nrows, ncols) = self.data.shape(); @@ -405,6 +432,20 @@ impl ColumnVector res } + + /// Constructs a vector from coordinates in projective space, i.e., removes a `0` at the end of + /// `self`. Returns `None` if this last component is not zero. + #[inline] + pub fn from_homogeneous(v: ColumnVector, SB>) -> Option> + where SB: Storage, U1, Alloc = S::Alloc> { + if v[v.len() - 1].is_zero() { + let nrows = D::from_usize(v.len() - 1); + Some(v.generic_slice((0, 0), (nrows, U1)).into_owned()) + } + else { + None + } + } } @@ -677,7 +718,7 @@ impl PartialOrd for Matrix if let Some(mut first_ord) = first_ord { let mut it = self.iter().zip(other.iter()); - it.next(); // Drop the first elements (we already tested it). + let _ = it.next(); // Drop the first elements (we already tested it). for (left, right) in it { if let Some(ord) = left.partial_cmp(right) { @@ -806,7 +847,7 @@ impl fmt::Display for Matrix impl Matrix - where N: Scalar + Field, + where N: Scalar + Ring, S: Storage { /// The dot product between two matrices (seen as vectors). #[inline] @@ -1011,3 +1052,35 @@ impl Matrix } } } + +impl ApproxEq for Unit> + where N: Scalar + ApproxEq, + S: Storage, + N::Epsilon: Copy { + type Epsilon = N::Epsilon; + + #[inline] + fn default_epsilon() -> Self::Epsilon { + N::default_epsilon() + } + + #[inline] + fn default_max_relative() -> Self::Epsilon { + N::default_max_relative() + } + + #[inline] + fn default_max_ulps() -> u32 { + N::default_max_ulps() + } + + #[inline] + fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool { + self.as_ref().relative_eq(other.as_ref(), epsilon, max_relative) + } + + #[inline] + fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { + self.as_ref().ulps_eq(other.as_ref(), epsilon, max_ulps) + } +} diff --git a/src/core/matrix_alga.rs b/src/core/matrix_alga.rs index 0c4ae645..ed0bad14 100644 --- a/src/core/matrix_alga.rs +++ b/src/core/matrix_alga.rs @@ -263,7 +263,7 @@ impl FiniteDimInnerSpace for Matrix a = Self::from_column_slice(&[N::zero(), -v[2], v[1]]); }; - a.normalize_mut(); + let _ = a.normalize_mut(); if f(&a.cross(v)) { f(&a); diff --git a/src/core/matrix_array.rs b/src/core/matrix_array.rs index 6d2e7154..5ef36198 100644 --- a/src/core/matrix_array.rs +++ b/src/core/matrix_array.rs @@ -1,6 +1,11 @@ +use std::mem; +use std::marker::PhantomData; use std::ops::{Deref, DerefMut, Mul}; -use std::fmt::{Debug, Formatter, Result}; +use std::fmt::{self, Debug, Formatter}; use std::hash::{Hash, Hasher}; +use serde::{Serialize, Serializer, Deserialize, Deserializer}; +use serde::ser::SerializeSeq; +use serde::de::{SeqVisitor, Visitor}; use typenum::Prod; use generic_array::{ArrayLength, GenericArray}; @@ -28,6 +33,7 @@ where R: DimName, data: GenericArray> } + impl Hash for MatrixArray where N: Hash, R: DimName, @@ -70,7 +76,7 @@ where N: Debug, R::Value: Mul, Prod: ArrayLength { #[inline] - fn fmt(&self, fmt: &mut Formatter) -> Result { + fn fmt(&self, fmt: &mut Formatter) -> fmt::Result { self.data.fmt(fmt) } } @@ -185,3 +191,97 @@ unsafe impl OwnedStorage for MatrixArray &mut self[..] } } + + +/* + * + * Allocation-less serde impls. + * + */ +// XXX: open an issue for GenericArray so that it implements serde traits? +impl Serialize for MatrixArray +where N: Scalar + Serialize, + R: DimName, + C: DimName, + R::Value: Mul, + Prod: ArrayLength { + + + fn serialize(&self, serializer: S) -> Result + where S: Serializer { + let mut serializer = serializer.serialize_seq_fixed_size(R::dim() * C::dim())?; + + for e in self.iter() { + serializer.serialize_element(e)?; + } + + serializer.end() + } +} + + +impl Deserialize for MatrixArray +where N: Scalar + Deserialize, + R: DimName, + C: DimName, + R::Value: Mul, + Prod: ArrayLength { + + + fn deserialize(deserializer: D) -> Result + where D: Deserializer { + + let len = R::dim() * C::dim(); + deserializer.deserialize_seq_fixed_size(len, MatrixArrayVisitor::new()) + } +} + + +/// A visitor that produces a matrix array. +struct MatrixArrayVisitor { + marker: PhantomData<(N, R, C)> +} + +impl MatrixArrayVisitor +where N: Scalar, + R: DimName, + C: DimName, + R::Value: Mul, + Prod: ArrayLength { + + /// Construct a new sequence visitor. + pub fn new() -> Self { + MatrixArrayVisitor { + marker: PhantomData, + } + } +} + +impl Visitor for MatrixArrayVisitor +where N: Scalar + Deserialize, + R: DimName, + C: DimName, + R::Value: Mul, + Prod: ArrayLength { + + type Value = MatrixArray; + + fn expecting(&self, formatter: &mut Formatter) -> fmt::Result { + formatter.write_str("a matrix array") + } + + #[inline] + fn visit_seq(self, mut visitor: V) -> Result, V::Error> + where V: SeqVisitor { + + let mut out: Self::Value = unsafe { mem::uninitialized() }; + let mut curr = 0; + + while let Some(value) = try!(visitor.visit()) { + out[curr] = value; + curr += 1; + } + + Ok(out) + } +} diff --git a/src/core/matrix_slice.rs b/src/core/matrix_slice.rs index 0008c6bc..8a09fb5f 100644 --- a/src/core/matrix_slice.rs +++ b/src/core/matrix_slice.rs @@ -7,7 +7,8 @@ use core::storage::{Storage, StorageMut, Owned}; use core::allocator::Allocator; macro_rules! slice_storage_impl( - ($Storage: ident as $SRef: ty; $T: ident.$get_addr: ident ($Ptr: ty as $Ref: ty)) => { + ($doc: expr; $Storage: ident as $SRef: ty; $T: ident.$get_addr: ident ($Ptr: ty as $Ref: ty)) => { + #[doc = $doc] pub struct $T<'a, N: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim, Alloc> { ptr: $Ptr, shape: (R, C), @@ -54,8 +55,14 @@ macro_rules! slice_storage_impl( } ); -slice_storage_impl!(Storage as &'a S; SliceStorage.get_address_unchecked(*const N as &'a N)); -slice_storage_impl!(StorageMut as &'a mut S; SliceStorageMut.get_address_unchecked_mut(*mut N as &'a mut N)); +slice_storage_impl!("A matrix data storage for a matrix slice. Only contains an internal reference \ + to another matrix data storage."; + Storage as &'a S; SliceStorage.get_address_unchecked(*const N as &'a N)); + +slice_storage_impl!("A mutable matrix data storage for mutable matrix slice. Only contains an \ + internal mutable reference to another matrix data storage."; + StorageMut as &'a mut S; SliceStorageMut.get_address_unchecked_mut(*mut N as &'a mut N) +); impl<'a, N: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim, Alloc> Copy @@ -345,6 +352,11 @@ macro_rules! matrix_slice_impl( } } + + /// Slices this matrix starting at its component `(start.0, start.1)` and with + /// `(shape.0, shape.1)` components. Each row (resp. column) of the sliced matrix is + /// separated by `steps.0` (resp. `steps.1`) ignored rows (resp. columns) of the + /// original matrix. #[inline] pub fn $slice_with_steps($me: $Me, start: (usize, usize), shape: (usize, usize), steps: (usize, usize)) -> $MatrixSlice { @@ -371,6 +383,10 @@ macro_rules! matrix_slice_impl( } } + /// Slices this matrix starting at its component `(start.0, start.1)` and with + /// `(R::dim(), CSlice::dim())` components. Each row (resp. column) of the sliced + /// matrix is separated by `steps.0` (resp. `steps.1`) ignored rows (resp. columns) of + /// the original matrix. #[inline] pub fn $fixed_slice_with_steps($me: $Me, start: (usize, usize), steps: (usize, usize)) -> $MatrixSlice diff --git a/src/core/matrix_vec.rs b/src/core/matrix_vec.rs index 2393a15d..b492ec47 100644 --- a/src/core/matrix_vec.rs +++ b/src/core/matrix_vec.rs @@ -5,7 +5,6 @@ use core::dimension::{Dim, DimName, Dynamic, U1}; use core::storage::{Storage, StorageMut, Owned, OwnedStorage}; use core::default_allocator::DefaultAllocator; - /* * * Storage. @@ -13,7 +12,7 @@ use core::default_allocator::DefaultAllocator; */ /// A Vec-based matrix data storage. It may be dynamically-sized. #[repr(C)] -#[derive(Eq, Debug, Clone, PartialEq)] +#[derive(Eq, Debug, Clone, PartialEq, Serialize, Deserialize)] pub struct MatrixVec { data: Vec, nrows: R, @@ -21,8 +20,10 @@ pub struct MatrixVec { } impl MatrixVec { + /// Creates a new dynamic matrix data storage from the given vector and shape. #[inline] pub fn new(nrows: R, ncols: C, data: Vec) -> MatrixVec { + assert!(nrows.value() * ncols.value() == data.len(), "Data storage buffer dimension mismatch."); MatrixVec { data: data, nrows: nrows, diff --git a/src/core/mod.rs b/src/core/mod.rs index e56db888..f37ec9b8 100644 --- a/src/core/mod.rs +++ b/src/core/mod.rs @@ -1,9 +1,11 @@ +//! Data structures for vector and matrix computations. + pub mod dimension; pub mod constraint; pub mod allocator; pub mod storage; pub mod coordinates; -pub mod ops; +mod ops; pub mod iter; pub mod default_allocator; diff --git a/src/core/ops.rs b/src/core/ops.rs index ebc444cf..f310f68b 100644 --- a/src/core/ops.rs +++ b/src/core/ops.rs @@ -33,7 +33,9 @@ impl Index<(usize, usize)> for Matrix #[inline] fn index(&self, ij: (usize, usize)) -> &N { - assert!(ij < self.shape(), "Matrix index out of bounds."); + let shape = self.shape(); + assert!(ij.0 < shape.0 && ij.1 < shape.1, "Matrix index out of bounds."); + unsafe { self.get_unchecked(ij.0, ij.1) } } } @@ -53,7 +55,9 @@ impl IndexMut<(usize, usize)> for Matrix #[inline] fn index_mut(&mut self, ij: (usize, usize)) -> &mut N { - assert!(ij < self.shape(), "Matrix index out of bounds."); + let shape = self.shape(); + assert!(ij.0 < shape.0 && ij.1 < shape.1, "Matrix index out of bounds."); + unsafe { self.get_unchecked_mut(ij.0, ij.1) } } } @@ -121,7 +125,13 @@ macro_rules! componentwise_binop_impl( assert!(self.shape() == right.shape(), "Matrix addition/subtraction dimensions mismatch."); let mut res = self.into_owned_sum::(); - for (left, right) in res.iter_mut().zip(right.iter()) { + // XXX: optimize our iterator! + // + // Using our own iterator prevents loop unrolling, wich breaks some optimization + // (like SIMD). On the other hand, using the slice iterator is 4x faster. + + // for (left, right) in res.iter_mut().zip(right.iter()) { + for (left, right) in res.as_mut_slice().iter_mut().zip(right.iter()) { *left = left.$method(*right) } @@ -143,7 +153,13 @@ macro_rules! componentwise_binop_impl( assert!(self.shape() == right.shape(), "Matrix addition/subtraction dimensions mismatch."); let mut res = right.into_owned_sum::(); - for (left, right) in self.iter().zip(res.iter_mut()) { + // XXX: optimize our iterator! + // + // Using our own iterator prevents loop unrolling, wich breaks some optimization + // (like SIMD). On the other hand, using the slice iterator is 4x faster. + + // for (left, right) in self.iter().zip(res.iter_mut()) { + for (left, right) in self.iter().zip(res.as_mut_slice().iter_mut()) { *right = left.$method(*right) } @@ -237,7 +253,13 @@ macro_rules! componentwise_scalarop_impl( fn $method(self, rhs: N) -> Self::Output { let mut res = self.into_owned(); - for left in res.iter_mut() { + // XXX: optimize our iterator! + // + // Using our own iterator prevents loop unrolling, wich breaks some optimization + // (like SIMD). On the other hand, using the slice iterator is 4x faster. + + // for left in res.iter_mut() { + for left in res.as_mut_slice().iter_mut() { *left = left.$method(rhs) } @@ -282,7 +304,13 @@ macro_rules! left_scalar_mul_impl( fn mul(self, right: Matrix<$T, R, C, S>) -> Self::Output { let mut res = right.into_owned(); - for right in res.iter_mut() { + // XXX: optimize our iterator! + // + // Using our own iterator prevents loop unrolling, wich breaks some optimization + // (like SIMD). On the other hand, using the slice iterator is 4x faster. + + // for right in res.iter_mut() { + for right in res.as_mut_slice().iter_mut() { *right = self * *right } diff --git a/src/core/storage.rs b/src/core/storage.rs index a3f98a43..d70e416c 100644 --- a/src/core/storage.rs +++ b/src/core/storage.rs @@ -1,3 +1,5 @@ +//! Abstract definition of a matrix data storage. + use std::mem; use std::any::Any; @@ -110,6 +112,12 @@ pub unsafe trait Storage: Sized { } } + +/// Trait implemented by matrix data storage that can provide a mutable access to its elements. +/// +/// Note that a mutable access does not mean that the matrix owns its data. For example, a mutable +/// matrix slice can provide mutable access to its elements even if it does not own its data (it +/// contains only an internal reference to them). pub unsafe trait StorageMut: Storage { /// The matrix mutable data pointer. fn ptr_mut(&mut self) -> *mut N; diff --git a/src/core/unit.rs b/src/core/unit.rs index 9c35bed0..c8e8c578 100644 --- a/src/core/unit.rs +++ b/src/core/unit.rs @@ -1,4 +1,5 @@ -use std::ops::Neg; +use std::mem; +use std::ops::{Neg, Deref}; use approx::ApproxEq; use alga::general::SubsetOf; @@ -9,7 +10,7 @@ use alga::linear::NormedSpace; /// /// Use `.as_ref()` or `.unwrap()` to obtain the undelying value by-reference or by-move. #[repr(C)] -#[derive(Eq, PartialEq, Clone, Hash, Debug, Copy)] +#[derive(Eq, PartialEq, Clone, Hash, Debug, Copy, Serialize, Deserialize)] pub struct Unit { value: T } @@ -114,34 +115,34 @@ where T::Field: ApproxEq { } } -impl ApproxEq for Unit { - type Epsilon = T::Epsilon; - - #[inline] - fn default_epsilon() -> Self::Epsilon { - T::default_epsilon() - } - - #[inline] - fn default_max_relative() -> Self::Epsilon { - T::default_max_relative() - } - - #[inline] - fn default_max_ulps() -> u32 { - T::default_max_ulps() - } - - #[inline] - fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool { - self.value.relative_eq(&other.value, epsilon, max_relative) - } - - #[inline] - fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { - self.value.ulps_eq(&other.value, epsilon, max_ulps) - } -} +// impl ApproxEq for Unit { +// type Epsilon = T::Epsilon; +// +// #[inline] +// fn default_epsilon() -> Self::Epsilon { +// T::default_epsilon() +// } +// +// #[inline] +// fn default_max_relative() -> Self::Epsilon { +// T::default_max_relative() +// } +// +// #[inline] +// fn default_max_ulps() -> u32 { +// T::default_max_ulps() +// } +// +// #[inline] +// fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool { +// self.value.relative_eq(&other.value, epsilon, max_relative) +// } +// +// #[inline] +// fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { +// self.value.ulps_eq(&other.value, epsilon, max_ulps) +// } +// } // FIXME:re-enable this impl when spacialization is possible. @@ -163,3 +164,12 @@ impl Neg for Unit { Unit::new_unchecked(-self.value) } } + +impl Deref for Unit { + type Target = T; + + #[inline] + fn deref(&self) -> &T { + unsafe { mem::transmute(self) } + } +} diff --git a/src/geometry/isometry.rs b/src/geometry/isometry.rs index 7b5fc6ee..5be41990 100644 --- a/src/geometry/isometry.rs +++ b/src/geometry/isometry.rs @@ -11,14 +11,24 @@ use core::storage::{Storage, OwnedStorage}; use core::allocator::{Allocator, OwnedAllocator}; use geometry::{TranslationBase, PointBase}; + +/// An isometry that uses a data storage deduced from the allocator `A`. +pub type OwnedIsometryBase = + IsometryBase>::Buffer, R>; + /// A direct isometry, i.e., a rotation followed by a translation. #[repr(C)] -#[derive(Hash, Debug, Clone, Copy)] +#[derive(Hash, Debug, Clone, Copy, Serialize, Deserialize)] pub struct IsometryBase { + /// The pure rotational part of this isometry. pub rotation: R, + /// The pure translational part of this isometry. pub translation: TranslationBase, - // One private field just to prevent explicit construction. - _noconstruct: PhantomData + + + // One dummy private field just to prevent explicit construction. + #[serde(skip_serializing, skip_deserializing)] + _noconstruct: PhantomData } impl IsometryBase diff --git a/src/geometry/isometry_ops.rs b/src/geometry/isometry_ops.rs index 02f924ed..0ca1bca8 100644 --- a/src/geometry/isometry_ops.rs +++ b/src/geometry/isometry_ops.rs @@ -214,8 +214,8 @@ isometry_binop_assign_impl_all!( ); -// IsometryBase × RotationBase -// IsometryBase ÷ RotationBase +// IsometryBase × R +// IsometryBase ÷ R isometry_binop_impl_all!( Mul, mul; self: IsometryBase, rhs: R, Output = IsometryBase; @@ -357,10 +357,10 @@ isometry_from_composition_impl_all!( Mul, mul; (D, D), (D, U1) for D: DimName; self: RotationBase, right: TranslationBase, Output = IsometryBase>; - [val val] => IsometryBase::from_parts(TranslationBase::from_vector(&self * right.vector), self); - [ref val] => IsometryBase::from_parts(TranslationBase::from_vector(self * right.vector), self.clone()); + [val val] => IsometryBase::from_parts(TranslationBase::from_vector(&self * right.vector), self); + [ref val] => IsometryBase::from_parts(TranslationBase::from_vector(self * right.vector), self.clone()); [val ref] => IsometryBase::from_parts(TranslationBase::from_vector(&self * &right.vector), self); - [ref ref] => IsometryBase::from_parts(TranslationBase::from_vector(self * &right.vector), self.clone()); + [ref ref] => IsometryBase::from_parts(TranslationBase::from_vector(self * &right.vector), self.clone()); ); @@ -370,10 +370,10 @@ isometry_from_composition_impl_all!( (U4, U1), (U3, U1); self: UnitQuaternionBase, right: TranslationBase, Output = IsometryBase>; - [val val] => IsometryBase::from_parts(TranslationBase::from_vector(&self * right.vector), self); - [ref val] => IsometryBase::from_parts(TranslationBase::from_vector(self * right.vector), self.clone()); + [val val] => IsometryBase::from_parts(TranslationBase::from_vector(&self * right.vector), self); + [ref val] => IsometryBase::from_parts(TranslationBase::from_vector( self * right.vector), self.clone()); [val ref] => IsometryBase::from_parts(TranslationBase::from_vector(&self * &right.vector), self); - [ref ref] => IsometryBase::from_parts(TranslationBase::from_vector(self * &right.vector), self.clone()); + [ref ref] => IsometryBase::from_parts(TranslationBase::from_vector( self * &right.vector), self.clone()); ); // RotationBase × IsometryBase diff --git a/src/geometry/mod.rs b/src/geometry/mod.rs index 21b51c61..0ee4c9e0 100644 --- a/src/geometry/mod.rs +++ b/src/geometry/mod.rs @@ -1,3 +1,5 @@ +//! Data structures for points and usual transformations (rotations, isometries, etc.) + mod op_macros; mod point; @@ -28,6 +30,7 @@ mod unit_complex; mod unit_complex_construction; mod unit_complex_ops; mod unit_complex_alga; +mod unit_complex_conversion; mod translation; mod translation_construction; diff --git a/src/geometry/orthographic.rs b/src/geometry/orthographic.rs index 18ebd9ed..900fd74a 100644 --- a/src/geometry/orthographic.rs +++ b/src/geometry/orthographic.rs @@ -13,7 +13,7 @@ use core::helper; use geometry::{PointBase, OwnedPoint}; /// A 3D orthographic projection stored as an homogeneous 4x4 matrix. -#[derive(Debug, Clone, Copy)] // FIXME: Hash +#[derive(Debug, Clone, Copy, Serialize, Deserialize)] // FIXME: Hash pub struct OrthographicBase> { matrix: SquareMatrix } @@ -25,9 +25,7 @@ impl Eq for OrthographicBase where N: Scalar + Eq, S: Storage { } -impl PartialEq for OrthographicBase - where N: Scalar, - S: Storage { +impl> PartialEq for OrthographicBase { #[inline] fn eq(&self, right: &Self) -> bool { self.matrix == right.matrix @@ -80,10 +78,7 @@ impl OrthographicBase } } -impl OrthographicBase - where N: Real, - S: Storage { - +impl> OrthographicBase { /// A reference to the underlying homogeneous transformation matrix. #[inline] pub fn as_matrix(&self) -> &SquareMatrix { @@ -96,6 +91,26 @@ impl OrthographicBase self.matrix } + /// Retrieves the inverse of the underlying homogeneous matrix. + #[inline] + pub fn inverse(&self) -> OwnedSquareMatrix { + let mut res = self.to_homogeneous(); + + let inv_m11 = N::one() / self.matrix[(0, 0)]; + let inv_m22 = N::one() / self.matrix[(1, 1)]; + let inv_m33 = N::one() / self.matrix[(2, 2)]; + + res[(0, 0)] = inv_m11; + res[(1, 1)] = inv_m22; + res[(2, 2)] = inv_m33; + + res[(0, 3)] = -self.matrix[(0, 3)] * inv_m11; + res[(1, 3)] = -self.matrix[(1, 3)] * inv_m22; + res[(2, 3)] = -self.matrix[(2, 3)] * inv_m33; + + res + } + /// Computes the corresponding homogeneous matrix. #[inline] pub fn to_homogeneous(&self) -> OwnedSquareMatrix { @@ -151,6 +166,18 @@ impl OrthographicBase ) } + /// Un-projects a point. Faster than multiplication by the underlying matrix inverse. + #[inline] + pub fn unproject_point(&self, p: &PointBase) -> OwnedPoint + where SB: Storage { + + OwnedPoint::::new( + (p[0] - self.matrix[(0, 3)]) / self.matrix[(0, 0)], + (p[1] - self.matrix[(1, 3)]) / self.matrix[(1, 1)], + (p[2] - self.matrix[(2, 3)]) / self.matrix[(2, 2)] + ) + } + // FIXME: when we get specialization, specialize the Mul impl instead. /// Projects a vector. Faster than matrix multiplication. #[inline] @@ -165,9 +192,7 @@ impl OrthographicBase } } -impl OrthographicBase - where N: Real, - S: StorageMut { +impl> OrthographicBase { /// Sets the smallest x-coordinate of the view cuboid. #[inline] pub fn set_left(&mut self, left: N) { diff --git a/src/geometry/perspective.rs b/src/geometry/perspective.rs index 55bf59b5..c9110ea9 100644 --- a/src/geometry/perspective.rs +++ b/src/geometry/perspective.rs @@ -13,7 +13,7 @@ use core::helper; use geometry::{PointBase, OwnedPoint}; /// A 3D perspective projection stored as an homogeneous 4x4 matrix. -#[derive(Debug, Clone, Copy)] // FIXME: Hash +#[derive(Debug, Clone, Copy, Serialize, Deserialize)] // FIXME: Hash pub struct PerspectiveBase> { matrix: SquareMatrix } @@ -85,6 +85,25 @@ impl PerspectiveBase self.matrix } + /// Retrieves the inverse of the underlying homogeneous matrix. + #[inline] + pub fn inverse(&self) -> OwnedSquareMatrix { + let mut res = self.to_homogeneous(); + + res[(0, 0)] = N::one() / self.matrix[(0, 0)]; + res[(1, 1)] = N::one() / self.matrix[(1, 1)]; + res[(2, 2)] = N::zero(); + + let m23 = self.matrix[(2, 3)]; + let m32 = self.matrix[(3, 2)]; + + res[(2, 3)] = N::one() / m32; + res[(3, 2)] = N::one() / m23; + res[(3, 3)] = -self.matrix[(2, 2)] / (m23 * m32); + + res + } + /// Computes the corresponding homogeneous matrix. #[inline] pub fn to_homogeneous(&self) -> OwnedSquareMatrix { @@ -137,6 +156,20 @@ impl PerspectiveBase ) } + /// Un-projects a point. Faster than multiplication by the matrix inverse. + #[inline] + pub fn unproject_point(&self, p: &PointBase) -> OwnedPoint + where SB: Storage { + + let inverse_denom = self.matrix[(2, 3)] / (p[2] + self.matrix[(2, 2)]); + + OwnedPoint::::new( + p[0] * inverse_denom / self.matrix[(0, 0)], + p[1] * inverse_denom / self.matrix[(1, 1)], + -inverse_denom + ) + } + // FIXME: when we get specialization, specialize the Mul impl instead. /// Projects a vector. Faster than matrix multiplication. #[inline] diff --git a/src/geometry/point.rs b/src/geometry/point.rs index e29158ee..70cba65f 100644 --- a/src/geometry/point.rs +++ b/src/geometry/point.rs @@ -9,13 +9,11 @@ use core::dimension::{DimName, DimNameSum, DimNameAdd, U1}; use core::storage::{Storage, StorageMut, MulStorage}; use core::allocator::{Allocator, SameShapeR}; -pub type PointSumStorage = - <>::Alloc as Allocator, U1>>::Buffer; - // XXX Bad name: we can't even add points… /// The type of the result of the sum of a point with a vector. pub type PointSum = - PointBase, PointSumStorage>; + PointBase, + <>::Alloc as Allocator, U1>>::Buffer>; /// The type of the result of the multiplication of a point by a matrix. pub type PointMul = PointBase>; @@ -23,10 +21,11 @@ pub type PointMul = PointBase = PointBase>::Buffer>; -/// A point in n-dimensional euclidean space. +/// A point in a n-dimensional euclidean space. #[repr(C)] -#[derive(Hash, Debug)] +#[derive(Hash, Debug, Serialize, Deserialize)] pub struct PointBase> { + /// The coordinates of this point, i.e., the shift from the origin. pub coords: ColumnVector } diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index 9332dfa2..e970ae15 100644 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -21,8 +21,9 @@ pub type OwnedUnitQuaternionBase = UnitQuaternionBase` for a quaternion /// that may be used as a rotation. #[repr(C)] -#[derive(Hash, Debug, Copy, Clone)] +#[derive(Hash, Debug, Copy, Clone, Serialize, Deserialize)] pub struct QuaternionBase> { + /// This quaternion as a 4D vector of coordinates in the `[ x, y, z, w ]` storage order. pub coords: ColumnVector } @@ -86,6 +87,12 @@ impl QuaternionBase self.coords.norm_squared() } + /// Normalizes this quaternion. + #[inline] + pub fn normalize(&self) -> OwnedQuaternionBase { + QuaternionBase::from_vector(self.coords.normalize()) + } + /// Compute the conjugate of this quaternion. #[inline] pub fn conjugate(&self) -> OwnedQuaternionBase { @@ -108,6 +115,13 @@ impl QuaternionBase None } } + + /// Linear interpolation between two quaternion. + #[inline] + pub fn lerp(&self, other: &QuaternionBase, t: N) -> OwnedQuaternionBase + where S2: Storage { + self * (N::one() - t) + other * t + } } @@ -208,6 +222,12 @@ impl QuaternionBase true } } + + /// Normalizes this quaternion. + #[inline] + pub fn normalize_mut(&mut self) -> N { + self.coords.normalize_mut() + } } impl ApproxEq for QuaternionBase @@ -262,6 +282,18 @@ pub type UnitQuaternionBase = Unit>; impl UnitQuaternionBase where N: Real, S: Storage { + /// Moves this unit quaternion into one that owns its data. + #[inline] + pub fn into_owned(self) -> OwnedUnitQuaternionBase { + UnitQuaternionBase::new_unchecked(self.unwrap().into_owned()) + } + + /// Clones this unit quaternion into one that owns its data. + #[inline] + pub fn clone_owned(&self) -> OwnedUnitQuaternionBase { + UnitQuaternionBase::new_unchecked(self.as_ref().clone_owned()) + } + /// The rotation angle in [0; pi] of this unit quaternion. #[inline] pub fn angle(&self) -> N { @@ -312,6 +344,75 @@ impl UnitQuaternionBase where S2: Storage { other / self } + + /// Linear interpolation between two unit quaternions. + /// + /// The result is not normalized. + #[inline] + pub fn lerp(&self, other: &UnitQuaternionBase, t: N) -> OwnedQuaternionBase + where S2: Storage { + self.as_ref().lerp(other.as_ref(), t) + } + + /// Normalized linear interpolation between two unit quaternions. + #[inline] + pub fn nlerp(&self, other: &UnitQuaternionBase, t: N) -> OwnedUnitQuaternionBase + where S2: Storage { + let mut res = self.lerp(other, t); + let _ = res.normalize_mut(); + + UnitQuaternionBase::new_unchecked(res) + } + + /// Spherical linear interpolation between two unit quaternions. + /// + /// Panics if the angle between both quaternion is 180 degrees (in which case the interpolation + /// is not well-defined). + #[inline] + pub fn slerp(&self, other: &UnitQuaternionBase, t: N) -> OwnedUnitQuaternionBase + where S2: Storage { + self.try_slerp(other, t, N::zero()).expect( + "Unable to perform a spherical quaternion interpolation when they \ + are 180 degree apart (the result is not unique).") + } + + /// Computes the spherical linear interpolation between two unit quaternions or returns `None` + /// if both quaternions are approximately 180 degrees apart (in which case the interpolation is + /// not well-defined). + /// + /// # Arguments + /// * `self`: the first quaternion to interpolate from. + /// * `other`: the second quaternion to interpolate toward. + /// * `t`: the interpolation parameter. Should be between 0 and 1. + /// * `epsilon`: the value bellow which the sinus of the angle separating both quaternion + /// must be to return `None`. + #[inline] + pub fn try_slerp(&self, other: &UnitQuaternionBase, t: N, epsilon: N) + -> Option> + where S2: Storage { + + let c_hang = self.coords.dot(&other.coords); + + // self == other + if c_hang.abs() >= N::one() { + return Some(self.clone_owned()) + } + + let hang = c_hang.acos(); + let s_hang = (N::one() - c_hang * c_hang).sqrt(); + + // FIXME: what if s_hang is 0.0 ? The result is not well-defined. + if relative_eq!(s_hang, N::zero(), epsilon = epsilon) { + None + } + else { + let ta = ((N::one() - t) * hang).sin() / s_hang; + let tb = (t * hang).sin() / s_hang; + let res = self.as_ref() * ta + other.as_ref() * tb; + + Some(UnitQuaternionBase::new_unchecked(res)) + } + } } impl UnitQuaternionBase @@ -453,3 +554,34 @@ impl fmt::Display for UnitQuaternionBase } } } + +impl ApproxEq for UnitQuaternionBase + where N: Real + ApproxEq, + S: Storage { + type Epsilon = N; + + #[inline] + fn default_epsilon() -> Self::Epsilon { + N::default_epsilon() + } + + #[inline] + fn default_max_relative() -> Self::Epsilon { + N::default_max_relative() + } + + #[inline] + fn default_max_ulps() -> u32 { + N::default_max_ulps() + } + + #[inline] + fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool { + self.as_ref().relative_eq(other.as_ref(), epsilon, max_relative) + } + + #[inline] + fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { + self.as_ref().ulps_eq(other.as_ref(), epsilon, max_ulps) + } +} diff --git a/src/geometry/quaternion_construction.rs b/src/geometry/quaternion_construction.rs index c83372be..2a627d93 100644 --- a/src/geometry/quaternion_construction.rs +++ b/src/geometry/quaternion_construction.rs @@ -33,7 +33,7 @@ impl QuaternionBase /// Creates a new quaternion from its individual components. Note that the arguments order does /// **not** follow the storage order. /// - /// The storage order is [ x, y, z, w ]. + /// The storage order is `[ x, y, z, w ]`. #[inline] pub fn new(w: N, x: N, y: N, z: N) -> Self { let v = ColumnVector::::new(x, y, z, w); @@ -169,19 +169,44 @@ impl UnitQuaternionBase pub fn from_rotation_matrix(rotmat: &RotationBase) -> Self where SB: Storage, SB::Alloc: Allocator { - let angle = rotmat.angle(); - if let Some(axis) = rotmat.axis() { - Self::from_axis_angle(&axis, angle) + // Robust matrix to quaternion transformation. + // See http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion + let tr = rotmat[(0, 0)] + rotmat[(1, 1)] + rotmat[(2, 2)]; + let res; + + let _0_25: N = ::convert(0.25); + + if tr > N::zero() { + let denom = (tr + N::one()).sqrt() * ::convert(2.0); + res = QuaternionBase::new(_0_25 * denom, + (rotmat[(2, 1)] - rotmat[(1, 2)]) / denom, + (rotmat[(0, 2)] - rotmat[(2, 0)]) / denom, + (rotmat[(1, 0)] - rotmat[(0, 1)]) / denom); } - else if angle > ::convert(1.0f64) { - // The angle is 3.14. - -Self::identity() + else if rotmat[(0, 0)] > rotmat[(1, 1)] && rotmat[(0, 0)] > rotmat[(2, 2)] { + let denom = (N::one() + rotmat[(0, 0)] - rotmat[(1, 1)] - rotmat[(2, 2)]).sqrt() * ::convert(2.0); + res = QuaternionBase::new((rotmat[(2, 1)] - rotmat[(1, 2)]) / denom, + _0_25 * denom, + (rotmat[(0, 1)] + rotmat[(1, 0)]) / denom, + (rotmat[(0, 2)] + rotmat[(2, 0)]) / denom); + } + else if rotmat[(1, 1)] > rotmat[(2, 2)] { + let denom = (N::one() + rotmat[(1, 1)] - rotmat[(0, 0)] - rotmat[(2, 2)]).sqrt() * ::convert(2.0); + res = QuaternionBase::new((rotmat[(0, 2)] - rotmat[(2, 0)]) / denom, + (rotmat[(0, 1)] + rotmat[(1, 0)]) / denom, + _0_25 * denom, + (rotmat[(1, 2)] + rotmat[(2, 1)]) / denom); } else { - // The angle is 0. - Self::identity() + let denom = (N::one() + rotmat[(2, 2)] - rotmat[(0, 0)] - rotmat[(1, 1)]).sqrt() * ::convert(2.0); + res = QuaternionBase::new((rotmat[(1, 0)] - rotmat[(0, 1)]) / denom, + (rotmat[(0, 2)] + rotmat[(2, 0)]) / denom, + (rotmat[(1, 2)] + rotmat[(2, 1)]) / denom, + _0_25 * denom); } + + Self::new_unchecked(res) } /// The unit quaternion needed to make `a` and `b` be collinear and point toward the same diff --git a/src/geometry/quaternion_coordinates.rs b/src/geometry/quaternion_coordinates.rs index aeb5beca..fdf4169c 100644 --- a/src/geometry/quaternion_coordinates.rs +++ b/src/geometry/quaternion_coordinates.rs @@ -1,5 +1,5 @@ -use std::ops::{Deref, DerefMut}; use std::mem; +use std::ops::{Deref, DerefMut}; use alga::general::Real; diff --git a/src/geometry/rotation.rs b/src/geometry/rotation.rs index 651ade5d..d447d585 100644 --- a/src/geometry/rotation.rs +++ b/src/geometry/rotation.rs @@ -15,8 +15,8 @@ pub type OwnedRotation = RotationBase>:: /// A rotation matrix. #[repr(C)] -#[derive(Hash, Debug, Clone, Copy)] -pub struct RotationBase> { +#[derive(Hash, Debug, Clone, Copy, Serialize, Deserialize)] +pub struct RotationBase { matrix: SquareMatrix } diff --git a/src/geometry/rotation_alga.rs b/src/geometry/rotation_alga.rs index e53c53f3..de226095 100644 --- a/src/geometry/rotation_alga.rs +++ b/src/geometry/rotation_alga.rs @@ -201,20 +201,23 @@ impl Rotation> for RotationBase, SB::Alloc: OwnedAllocator { #[inline] - fn powf(&self, n: N) -> Option { - // XXX: add the general case. + fn powf(&self, _: N) -> Option { + // XXX: Add the general case. + // XXX: Use specialization for 2D and 3D. unimplemented!() } #[inline] - fn rotation_between(a: &ColumnVector, b: &ColumnVector) -> Option { - // XXX: add the general case. + fn rotation_between(_: &ColumnVector, _: &ColumnVector) -> Option { + // XXX: Add the general case. + // XXX: Use specialization for 2D and 3D. unimplemented!() } #[inline] - fn scaled_rotation_between(a: &ColumnVector, b: &ColumnVector, n: N) -> Option { - // XXX: add the general case. + fn scaled_rotation_between(_: &ColumnVector, _: &ColumnVector, _: N) -> Option { + // XXX: Add the general case. + // XXX: Use specialization for 2D and 3D. unimplemented!() } } diff --git a/src/geometry/rotation_ops.rs b/src/geometry/rotation_ops.rs index c3643149..ca9b907b 100644 --- a/src/geometry/rotation_ops.rs +++ b/src/geometry/rotation_ops.rs @@ -118,7 +118,7 @@ md_impl_all!( ); -// RotationBase *= RotationBase +// RotationBase ×= RotationBase // FIXME: try not to call `inverse()` explicitly. md_assign_impl_all!( diff --git a/src/geometry/rotation_specialization.rs b/src/geometry/rotation_specialization.rs index 1a1190c4..9c0f13bb 100644 --- a/src/geometry/rotation_specialization.rs +++ b/src/geometry/rotation_specialization.rs @@ -44,7 +44,7 @@ where N: Real, pub fn rotation_between(a: &ColumnVector, b: &ColumnVector) -> Self where SB: Storage, SC: Storage { - UnitComplex::rotation_between(a, b).to_rotation_matrix() + ::convert(UnitComplex::rotation_between(a, b).to_rotation_matrix()) } /// The smallest rotation needed to make `a` and `b` collinear and point toward the same @@ -53,7 +53,7 @@ where N: Real, pub fn scaled_rotation_between(a: &ColumnVector, b: &ColumnVector, s: N) -> Self where SB: Storage, SC: Storage { - UnitComplex::scaled_rotation_between(a, b, s).to_rotation_matrix() + ::convert(UnitComplex::scaled_rotation_between(a, b, s).to_rotation_matrix()) } } @@ -81,6 +81,8 @@ where N: Real, other * self.inverse() } + /// Raise the quaternion to a given floating power, i.e., returns the rotation with the angle + /// of `self` multiplied by `n`. #[inline] pub fn powf(&self, n: N) -> OwnedRotation { OwnedRotation::<_, _, S::Alloc>::new(self.angle() * n) @@ -338,6 +340,8 @@ where N: Real, other * self.inverse() } + /// Raise the quaternion to a given floating power, i.e., returns the rotation with the same + /// axis as `self` and an angle equal to `self.angle()` multiplied by `n`. #[inline] pub fn powf(&self, n: N) -> OwnedRotation { if let Some(axis) = self.axis() { diff --git a/src/geometry/similarity.rs b/src/geometry/similarity.rs index 7b4d76f8..2a072bd6 100644 --- a/src/geometry/similarity.rs +++ b/src/geometry/similarity.rs @@ -10,10 +10,16 @@ use core::storage::{Storage, OwnedStorage}; use core::allocator::{Allocator, OwnedAllocator}; use geometry::{PointBase, TranslationBase, IsometryBase}; + +/// A similarity that uses a data storage deduced from the allocator `A`. +pub type OwnedSimilarityBase = + SimilarityBase>::Buffer, R>; + /// A similarity, i.e., an uniform scaling, followed by a rotation, followed by a translation. #[repr(C)] -#[derive(Hash, Debug, Clone, Copy)] +#[derive(Hash, Debug, Clone, Copy, Serialize, Deserialize)] pub struct SimilarityBase { + /// The part of this similarity that does not include the scaling factor. pub isometry: IsometryBase, scaling: N } @@ -62,12 +68,6 @@ impl SimilarityBase self.isometry.translation.vector *= self.scaling; } - /// The scaling factor of this similarity transformation. - #[inline] - pub fn scaling(&self) -> N { - self.scaling - } - /// The scaling factor of this similarity transformation. #[inline] pub fn set_scaling(&mut self, scaling: N) { @@ -162,6 +162,12 @@ impl SimilarityBase res } + + /// The scaling factor of this similarity transformation. + #[inline] + pub fn scaling(&self) -> N { + self.scaling + } } diff --git a/src/geometry/similarity_alias.rs b/src/geometry/similarity_alias.rs index 870d3a88..3bf96e75 100644 --- a/src/geometry/similarity_alias.rs +++ b/src/geometry/similarity_alias.rs @@ -12,7 +12,7 @@ pub type Similarity2 = SimilarityBase, UnitComp /// A 3-dimensional similarity. pub type Similarity3 = SimilarityBase, UnitQuaternion>; -/// A 3-dimensional similarity using a rotation matrix for its rotation part. +/// A 2-dimensional similarity using a rotation matrix for its rotation part. pub type SimilarityMatrix2 = Similarity; /// A 3-dimensional similarity using a rotation matrix for its rotation part. diff --git a/src/geometry/similarity_construction.rs b/src/geometry/similarity_construction.rs index 443de37a..d8f039f1 100644 --- a/src/geometry/similarity_construction.rs +++ b/src/geometry/similarity_construction.rs @@ -12,7 +12,8 @@ use core::dimension::{DimName, U1, U2, U3, U4}; use core::allocator::{OwnedAllocator, Allocator}; use core::storage::OwnedStorage; -use geometry::{PointBase, TranslationBase, RotationBase, SimilarityBase, UnitQuaternionBase, IsometryBase}; +use geometry::{PointBase, TranslationBase, RotationBase, SimilarityBase, + UnitComplex, UnitQuaternionBase, IsometryBase}; impl SimilarityBase @@ -106,6 +107,17 @@ impl SimilarityBase> } } +impl SimilarityBase> + where N: Real, + S: OwnedStorage, + S::Alloc: OwnedAllocator { + /// Creates a new similarity from a translation and a rotation angle. + #[inline] + pub fn new(translation: ColumnVector, angle: N, scaling: N) -> Self { + Self::from_parts(TranslationBase::from_vector(translation), UnitComplex::new(angle), scaling) + } +} + // 3D rotation. macro_rules! similarity_construction_impl( ($Rot: ty, $RotId: ident, $RRDim: ty, $RCDim: ty) => { diff --git a/src/geometry/transform.rs b/src/geometry/transform.rs index 831f6cab..726e0342 100644 --- a/src/geometry/transform.rs +++ b/src/geometry/transform.rs @@ -14,6 +14,8 @@ use core::allocator::Allocator; /// /// NOTE: this trait is not intended to be implementable outside of the `nalgebra` crate. pub trait TCategory: Any + Debug + Copy + PartialEq + Send { + /// Indicates whether a `Transform` with the category `Self` has a bottom-row different from + /// `0 0 .. 1`. #[inline] fn has_normalizer() -> bool { true @@ -28,16 +30,19 @@ pub trait TCategory: Any + Debug + Copy + PartialEq + Send { N::Epsilon: Copy; } -/// Traits that gives the transformation category that is compatible with the result of the +/// Traits that gives the `Transform` category that is compatible with the result of the /// multiplication of transformations with categories `Self` and `Other`. pub trait TCategoryMul: TCategory { + /// The transform category that results from the multiplication of a `Transform` to a + /// `Transform`. This is usually equal to `Self` or `Other`, whichever is the most + /// general category. type Representative: TCategory; } -/// Indicates that `Self` is a more general transformation category than `Other`. +/// Indicates that `Self` is a more general `Transform` category than `Other`. pub trait SuperTCategoryOf: TCategory { } -/// Indicates that `Self` is a more specific transformation category than `Other`. +/// Indicates that `Self` is a more specific `Transform` category than `Other`. /// /// Automatically implemented based on `SuperTCategoryOf`. pub trait SubTCategoryOf: TCategory { } @@ -46,17 +51,17 @@ where T1: TCategory, T2: SuperTCategoryOf { } -/// Tag representing the most general (not necessarily inversible) transformation type. -#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)] -pub struct TGeneral; +/// Tag representing the most general (not necessarily inversible) `Transform` type. +#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)] +pub enum TGeneral { } -/// Tag representing the most general inversible transformation type. -#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)] -pub struct TProjective; +/// Tag representing the most general inversible `Transform` type. +#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)] +pub enum TProjective { } -/// Tag representing an affine transformation. Its bottom-row is equal to `(0, 0 ... 0, 1)`. -#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)] -pub struct TAffine; +/// Tag representing an affine `Transform`. Its bottom-row is equal to `(0, 0 ... 0, 1)`. +#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq, Serialize, Deserialize)] +pub enum TAffine { } impl TCategory for TGeneral { #[inline] @@ -150,9 +155,11 @@ pub type OwnedTransform /// It is stored as a matrix with dimensions `(D + 1, D + 1)`, e.g., it stores a 4x4 matrix for a /// 3D transformation. #[repr(C)] -#[derive(Debug, Clone, Copy)] // FIXME: Hash +#[derive(Debug, Clone, Copy, Serialize, Deserialize)] // FIXME: Hash pub struct TransformBase, S, C: TCategory> { matrix: SquareMatrix, S>, + + #[serde(skip_serializing, skip_deserializing)] _phantom: PhantomData } diff --git a/src/geometry/transform_alias.rs b/src/geometry/transform_alias.rs index 1919916c..c387fb4f 100644 --- a/src/geometry/transform_alias.rs +++ b/src/geometry/transform_alias.rs @@ -3,14 +3,27 @@ use core::dimension::{U1, U2, U3, DimNameSum}; use geometry::{TransformBase, TGeneral, TProjective, TAffine}; -pub type Transform = TransformBase, DimNameSum>, TGeneral>; +/// A `D`-dimensional general transformation that may not be inversible. Stored as an homogeneous +/// `(D + 1) × (D + 1)` matrix. +pub type Transform = TransformBase, DimNameSum>, TGeneral>; + +/// An inversible `D`-dimensional general transformation. Stored as an homogeneous +/// `(D + 1) × (D + 1)` matrix. pub type Projective = TransformBase, DimNameSum>, TProjective>; -pub type Affine = TransformBase, DimNameSum>, TAffine>; -pub type Transform2 = Transform; +/// A `D`-dimensional affine transformation. Stored as an homogeneous `(D + 1) × (D + 1)` matrix. +pub type Affine = TransformBase, DimNameSum>, TAffine>; + +/// A 2D general transformation that may not be inversible. Stored as an homogeneous 3x3 matrix. +pub type Transform2 = Transform; +/// An inversible 2D general transformation. Stored as an homogeneous 3x3 matrix. pub type Projective2 = Projective; -pub type Affine2 = Affine; +/// A 2D affine transformation. Stored as an homogeneous 3x3 matrix. +pub type Affine2 = Affine; -pub type Transform3 = Transform; +/// A 3D general transformation that may not be inversible. Stored as an homogeneous 4x4 matrix. +pub type Transform3 = Transform; +/// An inversible 3D general transformation. Stored as an homogeneous 4x4 matrix. pub type Projective3 = Projective; -pub type Affine3 = Affine; +/// A 3D affine transformation. Stored as an homogeneous 4x4 matrix. +pub type Affine3 = Affine; diff --git a/src/geometry/transform_ops.rs b/src/geometry/transform_ops.rs index 21a8cf81..bacc883e 100644 --- a/src/geometry/transform_ops.rs +++ b/src/geometry/transform_ops.rs @@ -13,6 +13,68 @@ use geometry::{PointBase, OwnedPoint, TransformBase, OwnedTransform, TCategory, SubTCategoryOf, SuperTCategoryOf, TGeneral, TProjective, TAffine, RotationBase, UnitQuaternionBase, IsometryBase, SimilarityBase, TranslationBase}; +/* + * + * In the following, we provide: + * ========================= + * + * Index<(usize, usize)> + * IndexMut<(usize, usize)> (where TCategory == TGeneral) + * + * (Operators) + * + * TransformBase × IsometryBase + * TransformBase × RotationBase + * TransformBase × SimilarityBase + * TransformBase × TransformBase + * TransformBase × UnitQuaternion + * FIXME: TransformBase × UnitComplex + * TransformBase × TranslationBase + * TransformBase × ColumnVector + * TransformBase × PointBase + * + * IsometryBase × TransformBase + * RotationBase × TransformBase + * SimilarityBase × TransformBase + * TranslationBase × TransformBase + * UnitQuaternionBase × TransformBase + * FIXME: UnitComplex × TransformBase + * + * FIXME: TransformBase ÷ IsometryBase + * TransformBase ÷ RotationBase + * FIXME: TransformBase ÷ SimilarityBase + * TransformBase ÷ TransformBase + * TransformBase ÷ UnitQuaternion + * TransformBase ÷ TranslationBase + * + * FIXME: IsometryBase ÷ TransformBase + * RotationBase ÷ TransformBase + * FIXME: SimilarityBase ÷ TransformBase + * TranslationBase ÷ TransformBase + * UnitQuaternionBase ÷ TransformBase + * FIXME: UnitComplex ÷ TransformBase + * + * + * (Assignment Operators) + * + * + * TransformBase ×= TransformBase + * TransformBase ×= SimilarityBase + * TransformBase ×= IsometryBase + * TransformBase ×= RotationBase + * TransformBase ×= UnitQuaternionBase + * FIXME: TransformBase ×= UnitComplex + * TransformBase ×= TranslationBase + * + * TransformBase ÷= TransformBase + * FIXME: TransformBase ÷= SimilarityBase + * FIXME: TransformBase ÷= IsometryBase + * TransformBase ÷= RotationBase + * TransformBase ÷= UnitQuaternionBase + * FIXME: TransformBase ÷= UnitComplex + * + */ + /* * * Indexing. @@ -41,61 +103,6 @@ impl IndexMut<(usize, usize)> for TransformBase } } -/* - * - * In the following, we provide: - * ========================= - * - * - * (Operators) - * - * [OK] TransformBase × IsometryBase - * [OK] TransformBase × RotationBase - * [OK] TransformBase × SimilarityBase - * [OK] TransformBase × TransformBase - * [OK] TransformBase × UnitQuaternion - * [OK] TransformBase × TranslationBase - * [OK] TransformBase × ColumnVector - * [OK] TransformBase × PointBase - * - * [OK] IsometryBase × TransformBase - * [OK] RotationBase × TransformBase - * [OK] SimilarityBase × TransformBase - * [OK] TranslationBase × TransformBase - * [OK] UnitQuaternionBase × TransformBase - * - * TransformBase ÷ IsometryBase - * [OK] TransformBase ÷ RotationBase - * TransformBase ÷ SimilarityBase - * [OK] TransformBase ÷ TransformBase - * [OK] TransformBase ÷ UnitQuaternion - * [OK] TransformBase ÷ TranslationBase - * - * IsometryBase ÷ TransformBase - * [OK] RotationBase ÷ TransformBase - * SimilarityBase ÷ TransformBase - * [OK] TranslationBase ÷ TransformBase - * [OK] UnitQuaternionBase ÷ TransformBase - * - * - * (Assignment Operators) - * - * - * [OK] TransformBase ×= TransformBase - * [OK] TransformBase ×= SimilarityBase - * [OK] TransformBase ×= IsometryBase - * [OK] TransformBase ×= RotationBase - * [OK] TransformBase ×= UnitQuaternionBase - * [OK] TransformBase ×= TranslationBase - * - * [OK] TransformBase ÷= TransformBase - * TransformBase ÷= SimilarityBase - * TransformBase ÷= IsometryBase - * [OK] TransformBase ÷= RotationBase - * [OK] TransformBase ÷= UnitQuaternionBase - * - */ - // TransformBase × ColumnVector md_impl_all!( diff --git a/src/geometry/translation.rs b/src/geometry/translation.rs index 182ca2b0..2fb28671 100644 --- a/src/geometry/translation.rs +++ b/src/geometry/translation.rs @@ -14,8 +14,10 @@ pub type OwnedTranslation = TranslationBase*/> { + /// The translation coordinates, i.e., how much is added to a point's coordinates when it is + /// translated. pub vector: ColumnVector } diff --git a/src/geometry/unit_complex.rs b/src/geometry/unit_complex.rs index 29d07ab6..0d4eb876 100644 --- a/src/geometry/unit_complex.rs +++ b/src/geometry/unit_complex.rs @@ -1,12 +1,11 @@ use std::fmt; +use approx::ApproxEq; use num_complex::Complex; use alga::general::Real; -use core::{Unit, SquareMatrix, Vector1}; -use core::dimension::{U2, U3}; -use core::allocator::{OwnedAllocator, Allocator}; -use core::storage::OwnedStorage; -use geometry::{RotationBase, OwnedRotation}; +use core::{Unit, SquareMatrix, Vector1, Matrix3}; +use core::dimension::U2; +use geometry::Rotation2; /// A complex number with a norm equal to 1. pub type UnitComplex = Unit>; @@ -15,7 +14,7 @@ impl UnitComplex { /// The rotation angle in `]-pi; pi]` of this unit complex number. #[inline] pub fn angle(&self) -> N { - self.complex().im.atan2(self.complex().re) + self.im.atan2(self.re) } /// The rotation angle returned as a 1-dimensional vector. @@ -35,7 +34,7 @@ impl UnitComplex { /// Compute the conjugate of this unit complex number. #[inline] pub fn conjugate(&self) -> Self { - UnitComplex::new_unchecked(self.as_ref().conj()) + UnitComplex::new_unchecked(self.conj()) } /// Inverts this complex number if it is not zero. @@ -83,13 +82,11 @@ impl UnitComplex { /// Builds the rotation matrix corresponding to this unit complex number. #[inline] - pub fn to_rotation_matrix(&self) -> RotationBase - where S: OwnedStorage, - S::Alloc: OwnedAllocator { - let r = self.complex().re; - let i = self.complex().im; + pub fn to_rotation_matrix(&self) -> Rotation2 { + let r = self.re; + let i = self.im; - RotationBase::from_matrix_unchecked( + Rotation2::from_matrix_unchecked( SquareMatrix::<_, U2, _>::new( r, -i, i, r @@ -99,12 +96,8 @@ impl UnitComplex { /// Converts this unit complex number into its equivalent homogeneous transformation matrix. #[inline] - pub fn to_homogeneous(&self) -> SquareMatrix - where S: OwnedStorage, - S::Alloc: OwnedAllocator + - Allocator { - let r: OwnedRotation = self.to_rotation_matrix(); - r.to_homogeneous() + pub fn to_homogeneous(&self) -> Matrix3 { + self.to_rotation_matrix().to_homogeneous() } } @@ -113,3 +106,34 @@ impl fmt::Display for UnitComplex { write!(f, "UnitComplex angle: {}", self.angle()) } } + +impl ApproxEq for UnitComplex { + type Epsilon = N; + + #[inline] + fn default_epsilon() -> Self::Epsilon { + N::default_epsilon() + } + + #[inline] + fn default_max_relative() -> Self::Epsilon { + N::default_max_relative() + } + + #[inline] + fn default_max_ulps() -> u32 { + N::default_max_ulps() + } + + #[inline] + fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool { + self.re.relative_eq(&other.re, epsilon, max_relative) && + self.im.relative_eq(&other.im, epsilon, max_relative) + } + + #[inline] + fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { + self.re.ulps_eq(&other.re, epsilon, max_ulps) && + self.im.ulps_eq(&other.im, epsilon, max_ulps) + } +} diff --git a/src/geometry/unit_complex_conversion.rs b/src/geometry/unit_complex_conversion.rs new file mode 100644 index 00000000..b5d3bdf2 --- /dev/null +++ b/src/geometry/unit_complex_conversion.rs @@ -0,0 +1,173 @@ +use num::Zero; +use num_complex::Complex; + +use alga::general::{SubsetOf, SupersetOf, Real}; +use alga::linear::Rotation as AlgaRotation; + +use core::SquareMatrix; +use core::dimension::{U1, U2, U3}; +use core::storage::OwnedStorage; +use core::allocator::{Allocator, OwnedAllocator}; +use geometry::{PointBase, UnitComplex, RotationBase, OwnedRotation, IsometryBase, + SimilarityBase, TransformBase, SuperTCategoryOf, TAffine, TranslationBase}; + +/* + * This file provides the following conversions: + * ============================================= + * + * UnitComplex -> UnitComplex + * UnitComplex -> RotationBase + * UnitComplex -> IsometryBase + * UnitComplex -> SimilarityBase + * UnitComplex -> TransformBase + * UnitComplex -> Matrix (homogeneous) + * + * NOTE: + * UnitComplex -> Complex is already provided by: Unit -> T + */ + +impl SubsetOf> for UnitComplex + where N1: Real, + N2: Real + SupersetOf { + #[inline] + fn to_superset(&self) -> UnitComplex { + UnitComplex::new_unchecked(self.as_ref().to_superset()) + } + + #[inline] + fn is_in_subset(uq: &UnitComplex) -> bool { + ::is_convertible::<_, Complex>(uq.as_ref()) + } + + #[inline] + unsafe fn from_superset_unchecked(uq: &UnitComplex) -> Self { + Self::new_unchecked(::convert_ref_unchecked(uq.as_ref())) + } +} + +impl SubsetOf> for UnitComplex + where N1: Real, + N2: Real + SupersetOf, + S: OwnedStorage, + S::Alloc: OwnedAllocator + + Allocator + + Allocator + + Allocator { + #[inline] + fn to_superset(&self) -> RotationBase { + let q: UnitComplex = self.to_superset(); + q.to_rotation_matrix().to_superset() + } + + #[inline] + fn is_in_subset(rot: &RotationBase) -> bool { + ::is_convertible::<_, OwnedRotation>(rot) + } + + #[inline] + unsafe fn from_superset_unchecked(rot: &RotationBase) -> Self { + let q = UnitComplex::::from_rotation_matrix(rot); + ::convert_unchecked(q) + } +} + + +impl SubsetOf> for UnitComplex + where N1: Real, + N2: Real + SupersetOf, + S: OwnedStorage, + R: AlgaRotation> + SupersetOf>, + S::Alloc: OwnedAllocator { + #[inline] + fn to_superset(&self) -> IsometryBase { + IsometryBase::from_parts(TranslationBase::identity(), ::convert_ref(self)) + } + + #[inline] + fn is_in_subset(iso: &IsometryBase) -> bool { + iso.translation.vector.is_zero() + } + + #[inline] + unsafe fn from_superset_unchecked(iso: &IsometryBase) -> Self { + ::convert_ref_unchecked(&iso.rotation) + } +} + + +impl SubsetOf> for UnitComplex + where N1: Real, + N2: Real + SupersetOf, + S: OwnedStorage, + R: AlgaRotation> + SupersetOf>, + S::Alloc: OwnedAllocator { + #[inline] + fn to_superset(&self) -> SimilarityBase { + SimilarityBase::from_isometry(::convert_ref(self), N2::one()) + } + + #[inline] + fn is_in_subset(sim: &SimilarityBase) -> bool { + sim.isometry.translation.vector.is_zero() && + sim.scaling() == N2::one() + } + + #[inline] + unsafe fn from_superset_unchecked(sim: &SimilarityBase) -> Self { + ::convert_ref_unchecked(&sim.isometry) + } +} + + +impl SubsetOf> for UnitComplex + where N1: Real, + N2: Real + SupersetOf, + S: OwnedStorage, + C: SuperTCategoryOf, + S::Alloc: OwnedAllocator + + Allocator + + Allocator + + Allocator + + Allocator { + #[inline] + fn to_superset(&self) -> TransformBase { + TransformBase::from_matrix_unchecked(self.to_homogeneous().to_superset()) + } + + #[inline] + fn is_in_subset(t: &TransformBase) -> bool { + >::is_in_subset(t.matrix()) + } + + #[inline] + unsafe fn from_superset_unchecked(t: &TransformBase) -> Self { + Self::from_superset_unchecked(t.matrix()) + } +} + + +impl SubsetOf> for UnitComplex + where N1: Real, + N2: Real + SupersetOf, + S: OwnedStorage, + S::Alloc: OwnedAllocator + + Allocator + + Allocator + + Allocator + + Allocator { + #[inline] + fn to_superset(&self) -> SquareMatrix { + self.to_homogeneous().to_superset() + } + + #[inline] + fn is_in_subset(m: &SquareMatrix) -> bool { + ::is_convertible::<_, OwnedRotation>(m) + } + + #[inline] + unsafe fn from_superset_unchecked(m: &SquareMatrix) -> Self { + let rot: OwnedRotation = ::convert_ref_unchecked(m); + Self::from_rotation_matrix(&rot) + } +} diff --git a/src/geometry/unit_complex_ops.rs b/src/geometry/unit_complex_ops.rs index db09fe0f..f7331903 100644 --- a/src/geometry/unit_complex_ops.rs +++ b/src/geometry/unit_complex_ops.rs @@ -3,8 +3,13 @@ use std::ops::{Mul, MulAssign, Div, DivAssign}; use alga::general::Real; use core::{Unit, ColumnVector, OwnedColumnVector}; use core::dimension::{U1, U2}; -use core::storage::Storage; -use geometry::{UnitComplex, RotationBase, PointBase, OwnedPoint}; +use core::storage::{Storage, OwnedStorage}; +use core::allocator::OwnedAllocator; +use geometry::{UnitComplex, RotationBase, + PointBase, OwnedPoint, + IsometryBase, OwnedIsometryBase, + SimilarityBase, OwnedSimilarityBase, + TranslationBase}; /* * This file provides: @@ -23,8 +28,11 @@ use geometry::{UnitComplex, RotationBase, PointBase, OwnedPoint}; * UnitComplex × ColumnVector * UnitComplex × Unit * - * NOTE: -UnitComplex is already provided by `Unit`. + * UnitComplex × IsometryBase + * UnitComplex × SimilarityBase + * UnitComplex × TranslationBase -> IsometryBase * + * NOTE: -UnitComplex is already provided by `Unit`. * * (Assignment Operators) * @@ -34,8 +42,8 @@ use geometry::{UnitComplex, RotationBase, PointBase, OwnedPoint}; * UnitComplex ÷= UnitComplex * UnitComplex ÷= RotationBase * - * FIXME: RotationBase ×= UnitComplex - * FIXME: RotationBase ÷= UnitComplex + * RotationBase ×= UnitComplex + * RotationBase ÷= UnitComplex * */ @@ -232,7 +240,7 @@ complex_op_impl_all!( [ref ref] => { let i = self.as_ref().im; let r = self.as_ref().re; - OwnedColumnVector::<_, U2, S::Alloc>::new(r * rhs[0] - i * rhs[0], i * rhs[1] + r * rhs[1]) + OwnedColumnVector::<_, U2, S::Alloc>::new(r * rhs[0] - i * rhs[1], i * rhs[0] + r * rhs[1]) }; ); @@ -247,6 +255,45 @@ complex_op_impl_all!( [ref ref] => Unit::new_unchecked(self * rhs.as_ref()); ); +// UnitComplex × IsometryBase +complex_op_impl_all!( + Mul, mul; + (U2, U1); + self: UnitComplex, rhs: IsometryBase>, + Output = OwnedIsometryBase>; + [val val] => &self * &rhs; + [ref val] => self * &rhs; + [val ref] => &self * rhs; + [ref ref] => { + let shift = self * &rhs.translation.vector; + IsometryBase::from_parts(TranslationBase::from_vector(shift), self * &rhs.rotation) + }; +); + +// UnitComplex × SimilarityBase +complex_op_impl_all!( + Mul, mul; + (U2, U1); + self: UnitComplex, rhs: SimilarityBase>, + Output = OwnedSimilarityBase>; + [val val] => &self * &rhs; + [ref val] => self * &rhs; + [val ref] => &self * rhs; + [ref ref] => SimilarityBase::from_isometry(self * &rhs.isometry, rhs.scaling()); +); + +// UnitComplex × TranslationBase +complex_op_impl_all!( + Mul, mul; + (U2, U1); + self: UnitComplex, rhs: TranslationBase, + Output = OwnedIsometryBase>; + [val val] => IsometryBase::from_parts(TranslationBase::from_vector(&self * rhs.vector), self); + [ref val] => IsometryBase::from_parts(TranslationBase::from_vector( self * rhs.vector), self.clone()); + [val ref] => IsometryBase::from_parts(TranslationBase::from_vector(&self * &rhs.vector), self); + [ref ref] => IsometryBase::from_parts(TranslationBase::from_vector( self * &rhs.vector), self.clone()); +); + // UnitComplex ×= UnitComplex impl MulAssign> for UnitComplex { #[inline] @@ -307,3 +354,42 @@ impl<'b, N: Real, S: Storage> DivAssign<&'b RotationBase> f *self = &*self / rhs } } + + +// RotationBase ×= UnitComplex +impl MulAssign> for RotationBase + where S: OwnedStorage, + S::Alloc: OwnedAllocator { + #[inline] + fn mul_assign(&mut self, rhs: UnitComplex) { + self.mul_assign(rhs.to_rotation_matrix()) + } +} + +impl<'b, N: Real, S: Storage> MulAssign<&'b UnitComplex> for RotationBase + where S: OwnedStorage, + S::Alloc: OwnedAllocator { + #[inline] + fn mul_assign(&mut self, rhs: &'b UnitComplex) { + self.mul_assign(rhs.to_rotation_matrix()) + } +} + +// RotationBase ÷= UnitComplex +impl DivAssign> for RotationBase + where S: OwnedStorage, + S::Alloc: OwnedAllocator { + #[inline] + fn div_assign(&mut self, rhs: UnitComplex) { + self.div_assign(rhs.to_rotation_matrix()) + } +} + +impl<'b, N: Real, S: Storage> DivAssign<&'b UnitComplex> for RotationBase + where S: OwnedStorage, + S::Alloc: OwnedAllocator { + #[inline] + fn div_assign(&mut self, rhs: &'b UnitComplex) { + self.div_assign(rhs.to_rotation_matrix()) + } +} diff --git a/src/lib.rs b/src/lib.rs index 9d269e25..f713768d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,14 +1,97 @@ +/*! +# nalgebra + +**nalgebra** is a linear algebra library written for Rust targeting: + +* General-purpose linear algebra (still lacks a lot of features…) +* Real time computer graphics. +* Real time computer physics. + +## Using **nalgebra** +You will need the last stable build of the [rust compiler](http://www.rust-lang.org) +and the official package manager: [cargo](https://github.com/rust-lang/cargo). + +Simply add the following to your `Cargo.toml` file: + +```.ignore +[dependencies] +nalgebra = "0.11" +``` + + +Most useful functionalities of **nalgebra** are grouped in the root module `nalgebra::`. + +However, the recommended way to use **nalgebra** is to import types and traits +explicitly, and call free-functions using the `na::` prefix: + +```.rust +#[macro_use] +extern crate approx; // For the macro relative_eq! +extern crate nalgebra as na; +use na::{Vector3, Rotation3}; + +fn main() { + let axis = Vector3::x_axis(); + let angle = 1.57; + let b = Rotation3::from_axis_angle(&axis, angle); + + relative_eq!(b.axis().unwrap(), axis); + relative_eq!(b.angle(), angle); +} +``` + + +## Features +**nalgebra** is meant to be a general-purpose, low-dimensional, linear algebra library, with +an optimized set of tools for computer graphics and physics. Those features include: + +* A single parametrizable type `Matrix` for vectors, (square or rectangular) matrices, and slices + with dimensions known either at compile-time (using type-level integers) or at runtime. +* Matrices and vectors with compile-time sizes are statically allocated while dynamic ones are + allocated on the heap. +* Convenient aliases for low-dimensional matrices and vectors: `Vector1` to `Vector6` and + `Matrix1x1` to `Matrix6x6` (including rectangular matrices like `Matrix2x5`. +* Points sizes known at compile time, and convenience aliases: `Point1` to `Point6`. +* Translation (seen as a transformation that composes by multiplication): `Translation2`, + `Translation3`. +* Rotation matrices: `Rotation2`, `Rotation3`. +* Quaternions: `Quaternion`, `UnitQuaternion` (for 3D rotation). +* Unit complex numbers can be used for 2D rotation: `UnitComplex`. +* Algebraic entities with a norm equal to one: `Unit`, e.g., `Unit>`. +* Isometries (translation ⨯ rotation): `Isometry2`, `Isometry3` +* Similarity transformations (translation ⨯ rotation ⨯ uniform scale): `Similarity2`, `Similarity3`. +* Affine transformations stored as an homogeneous matrix: `Affine2`, `Affine3`. +* Projective (i.e. invertible) transformations stored as an homogeneous matrix: `Projective2`, + `Projective3`. +* General transformations that does not have to be invertible, stored as an homogeneous matrix: + `Transform2`, `Transform3`. +* 3D projections for computer graphics: `Perspective3`, `Orthographic3`. +* Linear algebra and data analysis operators: QR decomposition, eigen-decomposition. +* Implements all meaningful traits from the [alga](https://crates.io/crates/alga) crate for + generic programming. +*/ + + + + + // #![feature(plugin)] // // #![plugin(clippy)] -// #![allow(too_many_arguments)] -// #![allow(derive_hash_xor_eq)] -// #![allow(len_without_is_empty)] -// #![allow(transmute_ptr_to_ref)] + +#![deny(non_camel_case_types)] +#![deny(unused_parens)] +#![deny(non_upper_case_globals)] +#![deny(unused_qualifications)] +#![deny(unused_results)] +#![warn(missing_docs)] +#![doc(html_root_url = "http://nalgebra.org/doc")] #[cfg(feature = "arbitrary")] extern crate quickcheck; -extern crate rustc_serialize; +extern crate serde; +#[macro_use] +extern crate serde_derive; extern crate num_traits as num; extern crate num_complex; extern crate rand; @@ -32,11 +115,14 @@ pub use traits::*; use std::cmp::{self, PartialOrd, Ordering}; use num::Signed; -use alga::general::{Id, Identity, SupersetOf, MeetSemilattice, JoinSemilattice, Lattice, Inverse, +use alga::general::{Identity, SupersetOf, MeetSemilattice, JoinSemilattice, Lattice, Inverse, Multiplicative, Additive, AdditiveGroup}; use alga::linear::SquareMatrix as AlgaSquareMatrix; use alga::linear::{InnerSpace, NormedSpace, FiniteDimVectorSpace, EuclideanSpace}; +pub use alga::general::Id; + + /* * * Multiplicative identity. @@ -279,7 +365,7 @@ pub fn try_inverse(m: &M) -> Option { m.try_inverse() } -/// Computes the the multiplicative inverse (transformation, unit quaternion, etc.) +/// Computes the multiplicative inverse of an (always invertible) algebraic entity. #[inline] pub fn inverse>(m: &M) -> M { m.inverse() @@ -295,7 +381,7 @@ pub fn dot(a: &V, b: &V) -> V::Field { a.dot(b) } -/// Computes the angle between two vectors. +/// Computes the smallest angle between two vectors. #[inline] pub fn angle(a: &V, b: &V) -> V::Real { a.angle(b) @@ -305,25 +391,25 @@ pub fn angle(a: &V, b: &V) -> V::Real { * Normed space */ -/// Computes the L2 norm of a vector. +/// Computes the L2 (euclidean) norm of a vector. #[inline] pub fn norm(v: &V) -> V::Field { v.norm() } -/// Computes the squared L2 norm of a vector. +/// Computes the squared L2 (euclidean) norm of the vector `v`. #[inline] pub fn norm_squared(v: &V) -> V::Field { v.norm_squared() } -/// Gets the normalized version of a vector. +/// Computes the normalized version of the vector `v`. #[inline] pub fn normalize(v: &V) -> V { v.normalize() } -/// Gets the normalized version of a vector or `None` if its norm is smaller than `min_norm`. +/// Computes the normalized version of the vector `v` or returns `None` if its norm is smaller than `min_norm`. #[inline] pub fn try_normalize(v: &V, min_norm: V::Field) -> Option { v.try_normalize(min_norm) diff --git a/tests/isometry.rs b/tests/isometry.rs index 7c04fac8..2a5990b3 100644 --- a/tests/isometry.rs +++ b/tests/isometry.rs @@ -10,7 +10,10 @@ extern crate alga; extern crate nalgebra as na; use alga::linear::{Transformation, ProjectiveTransformation}; -use na::{Vector3, Point3, Rotation3, Isometry3, Translation3, UnitQuaternion}; +use na::{ + Vector3, Point3, Rotation3, Isometry3, Translation3, UnitQuaternion, + Vector2, Point2, Rotation2, Isometry2, Translation2, UnitComplex +}; quickcheck!( fn append_rotation_wrt_point_to_id(r: UnitQuaternion, p: Point3) -> bool { @@ -65,8 +68,39 @@ quickcheck!( relative_eq!(i.inverse() * p, i.inverse_transform_point(&p), epsilon = 1.0e-7) } - fn composition(i: Isometry3, uq: UnitQuaternion, r: Rotation3, - t: Translation3, v: Vector3, p: Point3) -> bool { + fn composition2(i: Isometry2, uc: UnitComplex, r: Rotation2, + t: Translation2, v: Vector2, p: Point2) -> bool { + // (rotation × translation) * point = rotation × (translation * point) + relative_eq!((uc * t) * v, uc * v, epsilon = 1.0e-7) && + relative_eq!((r * t) * v, r * v, epsilon = 1.0e-7) && + relative_eq!((uc * t) * p, uc * (t * p), epsilon = 1.0e-7) && + relative_eq!((r * t) * p, r * (t * p), epsilon = 1.0e-7) && + + // (translation × rotation) * point = translation × (rotation * point) + (t * uc) * v == uc * v && + (t * r) * v == r * v && + (t * uc) * p == t * (uc * p) && + (t * r) * p == t * (r * p) && + + // (rotation × isometry) * point = rotation × (isometry * point) + relative_eq!((uc * i) * v, uc * (i * v), epsilon = 1.0e-7) && + relative_eq!((uc * i) * p, uc * (i * p), epsilon = 1.0e-7) && + + // (isometry × rotation) * point = isometry × (rotation * point) + relative_eq!((i * uc) * v, i * (uc * v), epsilon = 1.0e-7) && + relative_eq!((i * uc) * p, i * (uc * p), epsilon = 1.0e-7) && + + // (translation × isometry) * point = translation × (isometry * point) + relative_eq!((t * i) * v, (i * v), epsilon = 1.0e-7) && + relative_eq!((t * i) * p, t * (i * p), epsilon = 1.0e-7) && + + // (isometry × translation) * point = isometry × (translation * point) + relative_eq!((i * t) * v, i * v, epsilon = 1.0e-7) && + relative_eq!((i * t) * p, i * (t * p), epsilon = 1.0e-7) + } + + fn composition3(i: Isometry3, uq: UnitQuaternion, r: Rotation3, + t: Translation3, v: Vector3, p: Point3) -> bool { // (rotation × translation) * point = rotation × (translation * point) relative_eq!((uq * t) * v, uq * v, epsilon = 1.0e-7) && relative_eq!((r * t) * v, r * v, epsilon = 1.0e-7) && diff --git a/tests/matrix_slice.rs b/tests/matrix_slice.rs index 10bdc9f4..6587fecb 100644 --- a/tests/matrix_slice.rs +++ b/tests/matrix_slice.rs @@ -1,8 +1,3 @@ -#[cfg(feature = "arbitrary")] -#[macro_use] -extern crate quickcheck; -#[macro_use] -extern crate approx; extern crate num_traits as num; extern crate nalgebra as na; diff --git a/tests/projection.rs b/tests/projection.rs new file mode 100644 index 00000000..6e25e5a3 --- /dev/null +++ b/tests/projection.rs @@ -0,0 +1,51 @@ +#[cfg(feature = "arbitrary")] +#[macro_use] +extern crate quickcheck; +#[macro_use] +extern crate approx; + +extern crate nalgebra as na; + +use na::{Point3, Perspective3, Orthographic3}; + +#[test] +fn perspective_inverse() { + let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0); + let inv = proj.inverse(); + + let id = inv * proj.unwrap(); + + assert!(id.is_identity(1.0e-7)); +} + +#[test] +fn orthographic_inverse() { + let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0); + let inv = proj.inverse(); + + let id = inv * proj.unwrap(); + + assert!(id.is_identity(1.0e-7)); +} + + +#[cfg(feature = "arbitrary")] +quickcheck!{ + fn perspective_project_unproject(pt: Point3) -> bool { + let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0); + + let projected = proj.project_point(&pt); + let unprojected = proj.unproject_point(&projected); + + relative_eq!(pt, unprojected, epsilon = 1.0e-7) + } + + fn orthographic_project_unproject(pt: Point3) -> bool { + let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0); + + let projected = proj.project_point(&pt); + let unprojected = proj.unproject_point(&projected); + + relative_eq!(pt, unprojected, epsilon = 1.0e-7) + } +} diff --git a/tests/rotation.rs b/tests/rotation.rs index 2a2e0af3..d683c3a9 100644 --- a/tests/rotation.rs +++ b/tests/rotation.rs @@ -7,6 +7,7 @@ extern crate num_traits as num; extern crate alga; extern crate nalgebra as na; +use std::f64; use alga::general::Real; use na::{Vector2, Vector3, Rotation2, Rotation3, Unit}; @@ -67,7 +68,6 @@ quickcheck!( a.angle(&b) == b.angle(&a) } - /* * * RotationBase matrix between vectors. diff --git a/tests/serde.rs b/tests/serde.rs new file mode 100644 index 00000000..f55512c5 --- /dev/null +++ b/tests/serde.rs @@ -0,0 +1,47 @@ +extern crate rand; +extern crate nalgebra as na; + +extern crate serde_json; + +use na::{ + DMatrix, + Matrix3x4, + Point3, + Translation3, + Rotation3, + Isometry3, + IsometryMatrix3, + Similarity3, + SimilarityMatrix3, + Quaternion +}; + +macro_rules! test_serde( + ($($test: ident, $ty: ident);* $(;)*) => {$( + #[test] + fn $test() { + let v: $ty = rand::random(); + let serialized = serde_json::to_string(&v).unwrap(); + assert_eq!(v, serde_json::from_str(&serialized).unwrap()); + } + )*} +); + +#[test] +fn serde_dmatrix() { + let v: DMatrix = DMatrix::new_random(3, 4); + let serialized = serde_json::to_string(&v).unwrap(); + assert_eq!(v, serde_json::from_str(&serialized).unwrap()); +} + +test_serde!( + serde_matrix3x4, Matrix3x4; + serde_point3, Point3; + serde_translation3, Translation3; + serde_rotation3, Rotation3; + serde_isometry3, Isometry3; + serde_isometry_matrix3, IsometryMatrix3; + serde_similarity3, Similarity3; + serde_similarity_matrix3, SimilarityMatrix3; + serde_quaternion, Quaternion; +); diff --git a/tests/unit_complex.rs b/tests/unit_complex.rs new file mode 100644 index 00000000..3d578276 --- /dev/null +++ b/tests/unit_complex.rs @@ -0,0 +1,171 @@ +#![allow(non_snake_case)] + +#[cfg(feature = "arbitrary")] +#[macro_use] +extern crate quickcheck; +#[macro_use] +extern crate approx; +extern crate num_traits as num; +extern crate alga; +extern crate nalgebra as na; + +use na::{Unit, UnitComplex, Vector2, Point2, Rotation2}; + + +quickcheck!( + + /* + * + * From/to rotation matrix. + * + */ + fn unit_complex_rotation_conversion(c: UnitComplex) -> bool { + let r = c.to_rotation_matrix(); + let cc = UnitComplex::from_rotation_matrix(&r); + let rr = cc.to_rotation_matrix(); + + relative_eq!(c, cc, epsilon = 1.0e-7) && + relative_eq!(r, rr, epsilon = 1.0e-7) + } + + /* + * + * PointBase/Vector transformation. + * + */ + fn unit_complex_transformation(c: UnitComplex, v: Vector2, p: Point2) -> bool { + let r = c.to_rotation_matrix(); + let rv = r * v; + let rp = r * p; + + relative_eq!( c * v, rv, epsilon = 1.0e-7) && + relative_eq!( c * &v, rv, epsilon = 1.0e-7) && + relative_eq!(&c * v, rv, epsilon = 1.0e-7) && + relative_eq!(&c * &v, rv, epsilon = 1.0e-7) && + + relative_eq!( c * p, rp, epsilon = 1.0e-7) && + relative_eq!( c * &p, rp, epsilon = 1.0e-7) && + relative_eq!(&c * p, rp, epsilon = 1.0e-7) && + relative_eq!(&c * &p, rp, epsilon = 1.0e-7) + } + + /* + * + * Inversion. + * + */ + fn unit_complex_inv(c: UnitComplex) -> bool { + let iq = c.inverse(); + relative_eq!(&iq * &c, UnitComplex::identity(), epsilon = 1.0e-7) && + relative_eq!( iq * &c, UnitComplex::identity(), epsilon = 1.0e-7) && + relative_eq!(&iq * c, UnitComplex::identity(), epsilon = 1.0e-7) && + relative_eq!( iq * c, UnitComplex::identity(), epsilon = 1.0e-7) && + + relative_eq!(&c * &iq, UnitComplex::identity(), epsilon = 1.0e-7) && + relative_eq!( c * &iq, UnitComplex::identity(), epsilon = 1.0e-7) && + relative_eq!(&c * iq, UnitComplex::identity(), epsilon = 1.0e-7) && + relative_eq!( c * iq, UnitComplex::identity(), epsilon = 1.0e-7) + } + + /* + * + * Quaterion * Vector == RotationBase * Vector + * + */ + fn unit_complex_mul_vector(c: UnitComplex, v: Vector2, p: Point2) -> bool { + let r = c.to_rotation_matrix(); + + relative_eq!(c * v, r * v, epsilon = 1.0e-7) && + relative_eq!(c * p, r * p, epsilon = 1.0e-7) + } + + // Test that all operators (incl. all combinations of references) work. + // See the top comment on `geometry/quaternion_ops.rs` for details on which operations are + // supported. + fn all_op_exist(uc: UnitComplex, v: Vector2, p: Point2, r: Rotation2) -> bool { + let uv = Unit::new_normalize(v); + + let ucMuc = uc * uc; + let ucMr = uc * r; + let rMuc = r * uc; + let ucDuc = uc / uc; + let ucDr = uc / r; + let rDuc = r / uc; + + let ucMp = uc * p; + let ucMv = uc * v; + let ucMuv = uc * uv; + + let mut ucMuc1 = uc; + let mut ucMuc2 = uc; + + let mut ucMr1 = uc; + let mut ucMr2 = uc; + + let mut ucDuc1 = uc; + let mut ucDuc2 = uc; + + let mut ucDr1 = uc; + let mut ucDr2 = uc; + + ucMuc1 *= uc; + ucMuc2 *= &uc; + + ucMr1 *= r; + ucMr2 *= &r; + + ucDuc1 /= uc; + ucDuc2 /= &uc; + + ucDr1 /= r; + ucDr2 /= &r; + + ucMuc1 == ucMuc && + ucMuc1 == ucMuc2 && + + ucMr1 == ucMr && + ucMr1 == ucMr2 && + + ucDuc1 == ucDuc && + ucDuc1 == ucDuc2 && + + ucDr1 == ucDr && + ucDr1 == ucDr2 && + + ucMuc == &uc * &uc && + ucMuc == uc * &uc && + ucMuc == &uc * uc && + + ucMr == &uc * &r && + ucMr == uc * &r && + ucMr == &uc * r && + + rMuc == &r * &uc && + rMuc == r * &uc && + rMuc == &r * uc && + + ucDuc == &uc / &uc && + ucDuc == uc / &uc && + ucDuc == &uc / uc && + + ucDr == &uc / &r && + ucDr == uc / &r && + ucDr == &uc / r && + + rDuc == &r / &uc && + rDuc == r / &uc && + rDuc == &r / uc && + + ucMp == &uc * &p && + ucMp == uc * &p && + ucMp == &uc * p && + + ucMv == &uc * &v && + ucMv == uc * &v && + ucMv == &uc * v && + + ucMuv == &uc * &uv && + ucMuv == uc * &uv && + ucMuv == &uc * uv + } +);