From 4f9a3cf1e64f9bc7f4067309a6a178e6a3c5715f Mon Sep 17 00:00:00 2001 From: Nathan Date: Wed, 6 Feb 2019 21:15:33 -0600 Subject: [PATCH 01/22] basic algorithm defined --- examples/convolution.rs | 62 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 examples/convolution.rs diff --git a/examples/convolution.rs b/examples/convolution.rs new file mode 100644 index 00000000..f2ab42cc --- /dev/null +++ b/examples/convolution.rs @@ -0,0 +1,62 @@ +extern crate nalgebra as na; +#[allow(unused_imports)] +use na::{Vector,Dim,Real,Vector4,Vector3,Vector2,U1,Matrix,DVector,Dynamic,VecStorage}; +use na::storage::{Storage}; +use std::cmp; + + + +// evum CovvolveMode{ +// Full, +// Valid, +// Same +// } + + + +#[allow(non_snake_case)] +fn Convolve1D, Q: Storage>( + Vector : Vector, + Kernel : Vector + ) -> Matrix> + { + // + // Vector is the vector, Kervel is the kervel + // C is the returv vector + // + if Kernel.len() > Vector.len(){ + return Convolve1D(Kernel, Vector); + } + + let V = Vector.len(); + let K = Kernel.len(); + let L = V + K - 1; + let v = V as i8; + let k = K as i8; + let l = L as i8; + let mut C = DVector::::zeros(L); + + for i in 0..l{ + let u_i = cmp::max(0, i - k); + let u_f = cmp::min(i, v - 1); + if u_i == u_f{ + C[i as usize] += Vector[u_i as usize] * Kernel[(i - u_i) as usize]; + } + else{ + for u in u_i..(u_f+1){ + if i - u < k{ + C[i as usize] += Vector[u as usize] * Kernel[(i - u ) as usize]; + } + } + } + } + C + } + + +fn main() { + let v1 = Vector2::new(3.0,3.0); + let v2 = Vector4::new(1.0,2.0,5.0,9.0); + let x = Convolve1D(v1,v2); + println!("{:?}",x) +} \ No newline at end of file From fd0c497c90d946085aa53388f8d4f16ba0552119 Mon Sep 17 00:00:00 2001 From: Nathan Date: Thu, 7 Feb 2019 19:58:09 -0600 Subject: [PATCH 02/22] Added valid method --- examples/convolution.rs | 123 ++++++++++++++++++++++++++++++---------- 1 file changed, 93 insertions(+), 30 deletions(-) diff --git a/examples/convolution.rs b/examples/convolution.rs index f2ab42cc..4619d8db 100644 --- a/examples/convolution.rs +++ b/examples/convolution.rs @@ -6,57 +6,120 @@ use std::cmp; -// evum CovvolveMode{ -// Full, -// Valid, -// Same -// } +enum ConvolveMode{ + Full, + Valid, + Same +} - - -#[allow(non_snake_case)] -fn Convolve1D, Q: Storage>( - Vector : Vector, - Kernel : Vector +fn convolve_full, Q: Storage>( + vector : Vector, + kernel : Vector ) -> Matrix> { - // - // Vector is the vector, Kervel is the kervel - // C is the returv vector - // - if Kernel.len() > Vector.len(){ - return Convolve1D(Kernel, Vector); - } - let V = Vector.len(); - let K = Kernel.len(); - let L = V + K - 1; - let v = V as i8; - let k = K as i8; - let l = L as i8; - let mut C = DVector::::zeros(L); + let vec = vector.len(); + let ker = kernel.len(); + let len = vec + ker - 1; + let v = vec as i8; + let k = ker as i8; + let l = len as i8; + let mut conv = DVector::::zeros(len); for i in 0..l{ let u_i = cmp::max(0, i - k); let u_f = cmp::min(i, v - 1); + if u_i == u_f{ - C[i as usize] += Vector[u_i as usize] * Kernel[(i - u_i) as usize]; + conv[i as usize] += vector[u_i as usize] * kernel[(i - u_i) as usize]; } else{ for u in u_i..(u_f+1){ if i - u < k{ - C[i as usize] += Vector[u as usize] * Kernel[(i - u ) as usize]; + conv[i as usize] += vector[u as usize] * kernel[(i - u ) as usize]; } } } } - C + conv + } + +fn convolve_valid, Q: Storage>( + vector : Vector, + kernel : Vector + ) -> Matrix> + { + let vec = vector.len(); + let ker = kernel.len(); + let len = vec - ker + 1; + let mut conv = DVector::::zeros(len); + + for i in 0..len { + for j in 0..ker { + conv[i] += vector[i + j] * kernel[ker - j - 1]; + } + } + + conv + } + +fn convolve_same, Q: Storage>( + vector : Vector, + kernel : Vector + ) -> Matrix> + { + + let vec = vector.len(); + let ker = kernel.len(); + let len = vec + ker - 1; + let v = vec as i8; + let k = ker as i8; + let l = len as i8; + let mut conv = DVector::::zeros(len); + + for i in 0..l { + let u_i = cmp::max(0, i - k); + let u_f = cmp::min(i, v - 1); + + if u_i == u_f { + conv[i as usize] += vector[u_i as usize] * kernel[(i - u_i) as usize]; + } + else{ + for u in u_i..(u_f+1){ + if i - u < k{ + conv[i as usize] += vector[u as usize] * kernel[(i - u ) as usize]; + } + } + } + } + conv + } + +fn convolve, Q: Storage>( + vector : Vector, + kernel : Vector, + mode : Option + ) -> Matrix> + { + // + // vector is the vector, Kervel is the kervel + // C is the returv vector + // + if kernel.len() > vector.len(){ + return convolve(kernel, vector, mode); + } + + match mode.unwrap_or(ConvolveMode::Full) { + ConvolveMode::Full => return convolve_full(vector,kernel), + ConvolveMode::Valid => return convolve_valid(vector,kernel), + ConvolveMode::Same => return convolve_same(vector,kernel) + } } fn main() { - let v1 = Vector2::new(3.0,3.0); + let v1 = Vector2::new(3.0,1.0); let v2 = Vector4::new(1.0,2.0,5.0,9.0); - let x = Convolve1D(v1,v2); + let x = convolve(v1,v2,Some(ConvolveMode::Valid)); println!("{:?}",x) } \ No newline at end of file From 3d83860fad91fdd21ad79aad08fd1bfd863d1324 Mon Sep 17 00:00:00 2001 From: Nathan Date: Sat, 9 Feb 2019 20:19:42 -0600 Subject: [PATCH 03/22] Cleaned up some labels and steps --- examples/convolution.rs | 183 +++++++++++++++++++++------------------- 1 file changed, 95 insertions(+), 88 deletions(-) diff --git a/examples/convolution.rs b/examples/convolution.rs index 4619d8db..c3c61eca 100644 --- a/examples/convolution.rs +++ b/examples/convolution.rs @@ -1,125 +1,132 @@ extern crate nalgebra as na; +use na::storage::Storage; #[allow(unused_imports)] -use na::{Vector,Dim,Real,Vector4,Vector3,Vector2,U1,Matrix,DVector,Dynamic,VecStorage}; -use na::storage::{Storage}; +use na::{ + DMatrix, DVector, Dim, Dynamic, Matrix, Matrix2x3, Real, VecStorage, Vector, Vector2, Vector3, + Vector4, U1, +}; use std::cmp; - - -enum ConvolveMode{ +enum ConvolveMode { Full, Valid, - Same + Same, } fn convolve_full, Q: Storage>( - vector : Vector, - kernel : Vector - ) -> Matrix> - { + vector: Vector, + kernel: Vector, +) -> Matrix> { + let vec = vector.len(); + let ker = kernel.len(); + let newlen = vec + ker - 1; - let vec = vector.len(); - let ker = kernel.len(); - let len = vec + ker - 1; - let v = vec as i8; - let k = ker as i8; - let l = len as i8; - let mut conv = DVector::::zeros(len); + let mut conv = DVector::::zeros(newlen); - for i in 0..l{ - let u_i = cmp::max(0, i - k); - let u_f = cmp::min(i, v - 1); + for i in 0..newlen { + let u_i = if i > ker {i - ker} else {0}; + let u_f = cmp::min(i, vec - 1); - if u_i == u_f{ - conv[i as usize] += vector[u_i as usize] * kernel[(i - u_i) as usize]; - } - else{ - for u in u_i..(u_f+1){ - if i - u < k{ - conv[i as usize] += vector[u as usize] * kernel[(i - u ) as usize]; - } + if u_i == u_f { + conv[i] += vector[u_i] * kernel[(i - u_i)]; + } else { + for u in u_i..(u_f + 1) { + if i - u < ker { + conv[i] += vector[u] * kernel[(i - u)]; } } } - conv } + conv +} fn convolve_valid, Q: Storage>( - vector : Vector, - kernel : Vector - ) -> Matrix> - { - let vec = vector.len(); - let ker = kernel.len(); - let len = vec - ker + 1; - let mut conv = DVector::::zeros(len); + vector: Vector, + kernel: Vector, +) -> Matrix> { + let vec = vector.len(); + let ker = kernel.len(); + let newlen = vec - ker + 1; + let mut conv = DVector::::zeros(newlen); - for i in 0..len { - for j in 0..ker { - conv[i] += vector[i + j] * kernel[ker - j - 1]; - } + for i in 0..newlen { + for j in 0..ker { + conv[i] += vector[i + j] * kernel[ker - j - 1]; } - - conv } + conv +} fn convolve_same, Q: Storage>( - vector : Vector, - kernel : Vector - ) -> Matrix> - { + vector: Vector, + kernel: Vector, +) -> Matrix> { + let vec = vector.len(); + let ker = kernel.len(); + let newlen = vec + ker - 1; - let vec = vector.len(); - let ker = kernel.len(); - let len = vec + ker - 1; - let v = vec as i8; - let k = ker as i8; - let l = len as i8; - let mut conv = DVector::::zeros(len); + let mut conv = DVector::::zeros(newlen); - for i in 0..l { - let u_i = cmp::max(0, i - k); - let u_f = cmp::min(i, v - 1); + for i in 0..newlen { + // let u_i = cmp::max(0, i - k); + // let u_f = cmp::min(i, v - 1); - if u_i == u_f { - conv[i as usize] += vector[u_i as usize] * kernel[(i - u_i) as usize]; - } - else{ - for u in u_i..(u_f+1){ - if i - u < k{ - conv[i as usize] += vector[u as usize] * kernel[(i - u ) as usize]; - } - } - } - } - conv + // if u_i == u_f { + // conv[i as usize] += vector[u_i as usize] * kernel[(i - u_i) as usize]; + // } else { + // for u in u_i..(u_f + 1) { + // if i - u < k { + // conv[i as usize] += vector[u as usize] * kernel[(i - u) as usize]; + // } + // } + // } } + conv +} fn convolve, Q: Storage>( - vector : Vector, - kernel : Vector, - mode : Option - ) -> Matrix> - { - // - // vector is the vector, Kervel is the kervel - // C is the returv vector - // - if kernel.len() > vector.len(){ - return convolve(kernel, vector, mode); - } - - match mode.unwrap_or(ConvolveMode::Full) { - ConvolveMode::Full => return convolve_full(vector,kernel), - ConvolveMode::Valid => return convolve_valid(vector,kernel), - ConvolveMode::Same => return convolve_same(vector,kernel) - } + vector: Vector, + kernel: Vector, + mode: Option, +) -> Matrix> { + // + // vector is the vector, Kervel is the kervel + // C is the returv vector + // + if kernel.len() > vector.len() { + return convolve(kernel, vector, mode); } + match mode.unwrap_or(ConvolveMode::Full) { + ConvolveMode::Full => return convolve_full(vector, kernel), + ConvolveMode::Valid => return convolve_valid(vector, kernel), + ConvolveMode::Same => return convolve_same(vector, kernel), + } +} fn main() { let v1 = Vector2::new(3.0,1.0); let v2 = Vector4::new(1.0,2.0,5.0,9.0); let x = convolve(v1,v2,Some(ConvolveMode::Valid)); println!("{:?}",x) -} \ No newline at end of file + + // let m = Matrix2x3::from_anti_diagonal_element(5.0); + // The two additional arguments represent the matrix dimensions. + // let dm = DMatrix::from_anti_diagonal_element(2, 3, 5.0); + let mut m = Matrix2x3::new(1.1, 1.2, 1.3, + 2.1, 2.2, 2.3); + + // assert!(m.m11 == 0.0 && m.m12 == 0.0 && m.m13 == 5.0 && + // m.m21 == 0.0 && m.m22 == 5.0 && m.m23 == 0.0); + // assert!(dm[(0, 0)] == 0.0 && dm[(0, 1)] == 0.0 && dm[(0, 2)] == 5.0 && + // dm[(1, 0)] == 0.0 && dm[(1, 1)] == 5.0 && dm[(1, 2)] == 0.0); + println!("m={:?}",m); + for i in 0..std::cmp::min(m.nrows(),m.ncols()) { + // for j in 0..3 { + println!("m({:?},{:?})={:?}",i,3-i-1,m[(i,3-i-1)]); +unsafe { println!("m({:?},{:?})={:?}",i,3-i-1,*m.get_unchecked_mut((i, 3-i-1))) } + // } + } + + +} From d038ec7b732293557de67424d33a8952f954507d Mon Sep 17 00:00:00 2001 From: Nathan Date: Sat, 9 Feb 2019 21:51:20 -0600 Subject: [PATCH 04/22] All 3 modes of direct done --- examples/convolution.rs | 65 ++++++++++++----------------------------- 1 file changed, 18 insertions(+), 47 deletions(-) diff --git a/examples/convolution.rs b/examples/convolution.rs index c3c61eca..7a3e8780 100644 --- a/examples/convolution.rs +++ b/examples/convolution.rs @@ -1,9 +1,8 @@ extern crate nalgebra as na; use na::storage::Storage; -#[allow(unused_imports)] use na::{ - DMatrix, DVector, Dim, Dynamic, Matrix, Matrix2x3, Real, VecStorage, Vector, Vector2, Vector3, - Vector4, U1, + convert, zero, DMatrix, DVector, Dim, Dynamic, Matrix, Matrix2x3, Real, VecStorage, Vector, + Vector2, Vector3, Vector4, Vector5, U1, }; use std::cmp; @@ -24,7 +23,7 @@ fn convolve_full, Q: Storage>( let mut conv = DVector::::zeros(newlen); for i in 0..newlen { - let u_i = if i > ker {i - ker} else {0}; + let u_i = if i > ker { i - ker } else { 0 }; let u_f = cmp::min(i, vec - 1); if u_i == u_f { @@ -47,6 +46,7 @@ fn convolve_valid, Q: Storage>( let vec = vector.len(); let ker = kernel.len(); let newlen = vec - ker + 1; + let mut conv = DVector::::zeros(newlen); for i in 0..newlen { @@ -63,23 +63,18 @@ fn convolve_same, Q: Storage>( ) -> Matrix> { let vec = vector.len(); let ker = kernel.len(); - let newlen = vec + ker - 1; - let mut conv = DVector::::zeros(newlen); + let mut conv = DVector::::zeros(vec); - for i in 0..newlen { - // let u_i = cmp::max(0, i - k); - // let u_f = cmp::min(i, v - 1); - - // if u_i == u_f { - // conv[i as usize] += vector[u_i as usize] * kernel[(i - u_i) as usize]; - // } else { - // for u in u_i..(u_f + 1) { - // if i - u < k { - // conv[i as usize] += vector[u as usize] * kernel[(i - u) as usize]; - // } - // } - // } + for i in 0..vec { + for j in 0..ker { + let val = if i + j < 1 || i + j >= vec + 1 { + zero::() + } else { + vector[i + j - 1] + }; + conv[i] += val * kernel[ker - j - 1]; + } } conv } @@ -89,10 +84,6 @@ fn convolve, Q: Storage>( kernel: Vector, mode: Option, ) -> Matrix> { - // - // vector is the vector, Kervel is the kervel - // C is the returv vector - // if kernel.len() > vector.len() { return convolve(kernel, vector, mode); } @@ -105,28 +96,8 @@ fn convolve, Q: Storage>( } fn main() { - let v1 = Vector2::new(3.0,1.0); - let v2 = Vector4::new(1.0,2.0,5.0,9.0); - let x = convolve(v1,v2,Some(ConvolveMode::Valid)); - println!("{:?}",x) - - // let m = Matrix2x3::from_anti_diagonal_element(5.0); - // The two additional arguments represent the matrix dimensions. - // let dm = DMatrix::from_anti_diagonal_element(2, 3, 5.0); - let mut m = Matrix2x3::new(1.1, 1.2, 1.3, - 2.1, 2.2, 2.3); - - // assert!(m.m11 == 0.0 && m.m12 == 0.0 && m.m13 == 5.0 && - // m.m21 == 0.0 && m.m22 == 5.0 && m.m23 == 0.0); - // assert!(dm[(0, 0)] == 0.0 && dm[(0, 1)] == 0.0 && dm[(0, 2)] == 5.0 && - // dm[(1, 0)] == 0.0 && dm[(1, 1)] == 5.0 && dm[(1, 2)] == 0.0); - println!("m={:?}",m); - for i in 0..std::cmp::min(m.nrows(),m.ncols()) { - // for j in 0..3 { - println!("m({:?},{:?})={:?}",i,3-i-1,m[(i,3-i-1)]); -unsafe { println!("m({:?},{:?})={:?}",i,3-i-1,*m.get_unchecked_mut((i, 3-i-1))) } - // } - } - - + let v1 = Vector4::new(1.0, 2.0, 1.0, 0.0); + let v2 = Vector4::new(1.0, 2.0, 5.0, 9.0); + let x = convolve(v1, v2, Some(ConvolveMode::Same)); + println!("{:?}", x); } From 645ca8ad52623f926651277a6391cc75ca207f87 Mon Sep 17 00:00:00 2001 From: Nathan Date: Sat, 9 Feb 2019 21:53:22 -0600 Subject: [PATCH 05/22] All 3 modes of direct done --- examples/convolution.rs | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/examples/convolution.rs b/examples/convolution.rs index 7a3e8780..07b1f974 100644 --- a/examples/convolution.rs +++ b/examples/convolution.rs @@ -1,9 +1,6 @@ extern crate nalgebra as na; use na::storage::Storage; -use na::{ - convert, zero, DMatrix, DVector, Dim, Dynamic, Matrix, Matrix2x3, Real, VecStorage, Vector, - Vector2, Vector3, Vector4, Vector5, U1, -}; +use na::{zero, DVector, Dim, Dynamic, Matrix, Real, VecStorage, Vector, U1}; use std::cmp; enum ConvolveMode { @@ -96,8 +93,5 @@ fn convolve, Q: Storage>( } fn main() { - let v1 = Vector4::new(1.0, 2.0, 1.0, 0.0); - let v2 = Vector4::new(1.0, 2.0, 5.0, 9.0); - let x = convolve(v1, v2, Some(ConvolveMode::Same)); - println!("{:?}", x); + } From b3c6492530cf3e93e197d222465da01d5af3f2a7 Mon Sep 17 00:00:00 2001 From: Nathan Date: Sun, 10 Feb 2019 13:40:32 -0600 Subject: [PATCH 06/22] Moved test file to lingal folder, wrote tests based on github ticket request (scipy reference) --- {examples => src/linalg}/convolution.rs | 74 ++++++++++++++----------- src/linalg/mod.rs | 2 + tests/linalg/convolution.rs | 49 ++++++++++++++++ tests/linalg/mod.rs | 1 + 4 files changed, 93 insertions(+), 33 deletions(-) rename {examples => src/linalg}/convolution.rs (58%) create mode 100644 tests/linalg/convolution.rs diff --git a/examples/convolution.rs b/src/linalg/convolution.rs similarity index 58% rename from examples/convolution.rs rename to src/linalg/convolution.rs index 07b1f974..dba6c97d 100644 --- a/examples/convolution.rs +++ b/src/linalg/convolution.rs @@ -1,20 +1,25 @@ -extern crate nalgebra as na; -use na::storage::Storage; -use na::{zero, DVector, Dim, Dynamic, Matrix, Real, VecStorage, Vector, U1}; +use storage::Storage; +use {zero, DVector, Dim, Dynamic, Matrix, Real, VecStorage, Vector, U1}; use std::cmp; -enum ConvolveMode { - Full, - Valid, - Same, -} - -fn convolve_full, Q: Storage>( +/// +/// The output is the full discrete linear convolution of the inputs +/// +pub fn convolve_full, Q: Storage>( vector: Vector, kernel: Vector, ) -> Matrix> { let vec = vector.len(); let ker = kernel.len(); + + if vec == 0 || ker == 0 { + panic!("Convolve's inputs must not be 0-sized. "); + } + + if ker > vec { + return convolve_full(kernel, vector); + } + let newlen = vec + ker - 1; let mut conv = DVector::::zeros(newlen); @@ -36,12 +41,24 @@ fn convolve_full, Q: Storage>( conv } -fn convolve_valid, Q: Storage>( +/// +/// The output consists only of those elements that do not rely on the zero-padding. +/// +pub fn convolve_valid, Q: Storage>( vector: Vector, kernel: Vector, ) -> Matrix> { let vec = vector.len(); let ker = kernel.len(); + + if vec == 0 || ker == 0 { + panic!("Convolve's inputs must not be 0-sized. "); + } + + if ker > vec { + return convolve_valid(kernel, vector); + } + let newlen = vec - ker + 1; let mut conv = DVector::::zeros(newlen); @@ -54,13 +71,24 @@ fn convolve_valid, Q: Storage>( conv } -fn convolve_same, Q: Storage>( +/// +/// The output is the same size as in1, centered with respect to the ‘full’ output. +/// +pub fn convolve_same, Q: Storage>( vector: Vector, kernel: Vector, ) -> Matrix> { let vec = vector.len(); let ker = kernel.len(); + if vec == 0 || ker == 0 { + panic!("Convolve's inputs must not be 0-sized. "); + } + + if ker > vec { + return convolve_same(kernel, vector); + } + let mut conv = DVector::::zeros(vec); for i in 0..vec { @@ -74,24 +102,4 @@ fn convolve_same, Q: Storage>( } } conv -} - -fn convolve, Q: Storage>( - vector: Vector, - kernel: Vector, - mode: Option, -) -> Matrix> { - if kernel.len() > vector.len() { - return convolve(kernel, vector, mode); - } - - match mode.unwrap_or(ConvolveMode::Full) { - ConvolveMode::Full => return convolve_full(vector, kernel), - ConvolveMode::Valid => return convolve_valid(vector, kernel), - ConvolveMode::Same => return convolve_same(vector, kernel), - } -} - -fn main() { - -} +} \ No newline at end of file diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 4418b283..b6a9e8d8 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -17,6 +17,7 @@ mod solve; mod svd; mod symmetric_eigen; mod symmetric_tridiagonal; +mod convolution; //// FIXME: Not complete enough for publishing. //// This handles only cases where each eigenvalue has multiplicity one. @@ -33,3 +34,4 @@ pub use self::schur::*; pub use self::svd::*; pub use self::symmetric_eigen::*; pub use self::symmetric_tridiagonal::*; +pub use self::convolution::*; diff --git a/tests/linalg/convolution.rs b/tests/linalg/convolution.rs new file mode 100644 index 00000000..ef3a02db --- /dev/null +++ b/tests/linalg/convolution.rs @@ -0,0 +1,49 @@ +use na::linalg::{convolve_full,convolve_valid,convolve_same}; +use na::{Vector2,Vector4,DVector}; + +// +// Should mimic calculations in Python's scipy library +// >>>from scipy.signal import convolve +// + +// >>> convolve([1,2,3,4],[1,2],"same") +// array([ 1, 4, 7, 10]) +#[test] +fn convolve_same_check(){ + let vec = Vector4::new(1.0,2.0,3.0,4.0); + let ker = Vector2::new(1.0,2.0); + + let actual = DVector::from_vec(4, vec![1.0,4.0,7.0,10.0]); + + let expected = convolve_same(vec,ker); + + assert!(relative_eq!(actual, expected, epsilon = 1.0e-7)); +} + +// >>> convolve([1,2,3,4],[1,2],"valid") +// array([ 1, 4, 7, 10, 8]) +#[test] +fn convolve_full_check(){ + let vec = Vector4::new(1.0,2.0,3.0,4.0); + let ker = Vector2::new(1.0,2.0); + + let actual = DVector::from_vec(5, vec![1.0,4.0,7.0,10.0,8.0]); + + let expected = convolve_full(vec,ker); + + assert!(relative_eq!(actual, expected, epsilon = 1.0e-7)); +} + +// >>> convolve([1,2,3,4],[1,2],"valid") +// array([ 4, 7, 10]) +#[test] +fn convolve_valid_check(){ + let vec = Vector4::new(1.0,2.0,3.0,4.0); + let ker = Vector2::new(1.0,2.0); + + let actual = DVector::from_vec(3, vec![4.0,7.0,10.0]); + + let expected = convolve_valid(vec,ker); + + assert!(relative_eq!(actual, expected, epsilon = 1.0e-7)); +} \ No newline at end of file diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs index 74a5e03c..4e0bf2eb 100644 --- a/tests/linalg/mod.rs +++ b/tests/linalg/mod.rs @@ -11,3 +11,4 @@ mod real_schur; mod solve; mod svd; mod tridiagonal; +mod convolution; \ No newline at end of file From bca385ea6b132d8300c476e3637452dbb9f12411 Mon Sep 17 00:00:00 2001 From: Nathan Date: Sun, 10 Feb 2019 13:46:37 -0600 Subject: [PATCH 07/22] Quick fix to documentation --- src/linalg/convolution.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs index dba6c97d..e7e20c21 100644 --- a/src/linalg/convolution.rs +++ b/src/linalg/convolution.rs @@ -42,7 +42,7 @@ pub fn convolve_full, Q: Storage } /// -/// The output consists only of those elements that do not rely on the zero-padding. +/// The output convolution consists only of those elements that do not rely on the zero-padding. /// pub fn convolve_valid, Q: Storage>( vector: Vector, @@ -72,7 +72,7 @@ pub fn convolve_valid, Q: Storage, Q: Storage>( vector: Vector, From b08c2ad70d00fb4322f9ed860bb57a51840b2fac Mon Sep 17 00:00:00 2001 From: Nathan Date: Thu, 14 Feb 2019 20:54:26 -0600 Subject: [PATCH 08/22] Feedback updates round 1 --- examples/convolution.rs | 5 +++ src/linalg/convolution.rs | 77 ++++++++++++++++++++++++--------------- 2 files changed, 53 insertions(+), 29 deletions(-) create mode 100644 examples/convolution.rs diff --git a/examples/convolution.rs b/examples/convolution.rs new file mode 100644 index 00000000..5290440c --- /dev/null +++ b/examples/convolution.rs @@ -0,0 +1,5 @@ +fn main(){ + let (x,y) = (1,2); + + println!("{}", x); +} \ No newline at end of file diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs index e7e20c21..ac0ba71e 100644 --- a/src/linalg/convolution.rs +++ b/src/linalg/convolution.rs @@ -1,45 +1,63 @@ use storage::Storage; -use {zero, DVector, Dim, Dynamic, Matrix, Real, VecStorage, Vector, U1}; +use {zero, DVector, Dim, Dynamic, Matrix, Real, VecStorage, Vector, U1, Add}; use std::cmp; -/// -/// The output is the full discrete linear convolution of the inputs -/// -pub fn convolve_full, Q: Storage>( - vector: Vector, - kernel: Vector, -) -> Matrix> { - let vec = vector.len(); - let ker = kernel.len(); +impl> Vector{ - if vec == 0 || ker == 0 { - panic!("Convolve's inputs must not be 0-sized. "); - } + /// Returns the convolution of the vector and a kernel + /// + /// # Arguments + /// + /// * `self` - A DVector with size D > 0 + /// * `kernel` - A DVector with size D > 0 + /// + /// # Note: + /// This function is commutative. If D_kernel > D_vector, + /// they will swap their roles as in + /// (self, kernel) = (kernel,self) + /// + /// # Example + /// + /// ``` + /// + /// ``` + pub fn convolve_full>(&self, kernel: Vector) -> Vector,Add> + { + let vec = self.len(); + let ker = kernel.len(); - if ker > vec { - return convolve_full(kernel, vector); - } + // if vec == 0 || ker == 0 { + // panic!("Convolve's inputs must not be 0-sized. "); + // } - let newlen = vec + ker - 1; + // if ker > vec { + // return kernel::convolve_full(vector); + // } - let mut conv = DVector::::zeros(newlen); + let newlen = vec + ker - 1; + let mut conv = DVector::::zeros(newlen); - for i in 0..newlen { - let u_i = if i > ker { i - ker } else { 0 }; - let u_f = cmp::min(i, vec - 1); + for i in 0..newlen { + let u_i = if i > ker { i - ker } else { 0 }; + let u_f = cmp::min(i, vec - 1); - if u_i == u_f { - conv[i] += vector[u_i] * kernel[(i - u_i)]; - } else { - for u in u_i..(u_f + 1) { - if i - u < ker { - conv[i] += vector[u] * kernel[(i - u)]; + if u_i == u_f { + conv[i] += self[u_i] * kernel[(i - u_i)]; + } else { + for u in u_i..(u_f + 1) { + if i - u < ker { + conv[i] += self[u] * kernel[(i - u)]; + } } } } + // conv } - conv } +/// +/// The output is the full discrete linear convolution of the inputs +/// + /// /// The output convolution consists only of those elements that do not rely on the zero-padding. @@ -102,4 +120,5 @@ pub fn convolve_same, Q: Storage } } conv -} \ No newline at end of file +} + From 9f5201938594056a80a9b20423fa8d4993f7289a Mon Sep 17 00:00:00 2001 From: Nathan Date: Mon, 18 Feb 2019 19:01:18 -0600 Subject: [PATCH 09/22] Fixing type traits based on feedback, `convolve_full` still broken --- src/linalg/convolution.rs | 196 ++++++++++++++++++++++++------------ tests/linalg/convolution.rs | 72 +++++++++---- 2 files changed, 181 insertions(+), 87 deletions(-) diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs index ac0ba71e..c587f8f2 100644 --- a/src/linalg/convolution.rs +++ b/src/linalg/convolution.rs @@ -1,71 +1,110 @@ -use storage::Storage; -use {zero, DVector, Dim, Dynamic, Matrix, Real, VecStorage, Vector, U1, Add}; +use base::allocator::Allocator; +use base::default_allocator::DefaultAllocator; +use base::dimension::{DimAdd, DimDiff, DimMax, DimMaximum, DimName, DimSub, DimSum,Dim}; use std::cmp; +use storage::Storage; +use {zero, Real, Vector, VectorN, U1}; -impl> Vector{ +/// Returns the convolution of the vector and a kernel +/// +/// # Arguments +/// +/// * `vector` - A Vector with size > 0 +/// * `kernel` - A Vector with size > 0 +/// +/// # Note: +/// This function is commutative. If kernel > vector, +/// they will swap their roles as in +/// (self, kernel) = (kernel,self) +/// +/// # Example +/// +/// ``` +/// let vec = Vector3::new(1.0,2.0,3.0); +/// let ker = Vector2::new(0.4,0.6); +/// let convolve = convolve_full(vec,ker); +/// ``` +pub fn convolve_full( + vector: Vector, + kernel: Vector, +) -> VectorN, U1>> +where + N: Real, + D1: DimAdd, + D2: DimAdd>, + DimSum: DimSub, + S1: Storage, + S2: Storage, + DimSum: Dim, + DefaultAllocator: Allocator, U1>>, +{ + let vec = vector.len(); + let ker = kernel.len(); - /// Returns the convolution of the vector and a kernel - /// - /// # Arguments - /// - /// * `self` - A DVector with size D > 0 - /// * `kernel` - A DVector with size D > 0 - /// - /// # Note: - /// This function is commutative. If D_kernel > D_vector, - /// they will swap their roles as in - /// (self, kernel) = (kernel,self) - /// - /// # Example - /// - /// ``` - /// - /// ``` - pub fn convolve_full>(&self, kernel: Vector) -> Vector,Add> - { - let vec = self.len(); - let ker = kernel.len(); + if vec == 0 || ker == 0 { + panic!("Convolve's inputs must not be 0-sized. "); + } - // if vec == 0 || ker == 0 { - // panic!("Convolve's inputs must not be 0-sized. "); - // } + if ker > vec { + return convolve_full(kernel, vector); + } - // if ker > vec { - // return kernel::convolve_full(vector); - // } + let result_len = vector.data.shape().0.add(kernel.data.shape().0).sub(U1); + let mut conv = VectorN::zeros_generic(result_len, U1); - let newlen = vec + ker - 1; - let mut conv = DVector::::zeros(newlen); + for i in 0..(vec + ker - 1) { + let u_i = if i > vec { i - ker } else { 0 }; + let u_f = cmp::min(i, vec - 1); - for i in 0..newlen { - let u_i = if i > ker { i - ker } else { 0 }; - let u_f = cmp::min(i, vec - 1); - - if u_i == u_f { - conv[i] += self[u_i] * kernel[(i - u_i)]; - } else { - for u in u_i..(u_f + 1) { - if i - u < ker { - conv[i] += self[u] * kernel[(i - u)]; - } + if u_i == u_f { + conv[i] += vector[u_i] * kernel[(i - u_i)]; + } else { + for u in u_i..(u_f + 1) { + if i - u < ker { + conv[i] += vector[u] * kernel[(i - u)]; } } } - // conv } + conv } -/// -/// The output is the full discrete linear convolution of the inputs -/// + +/// Returns the convolution of the vector and a kernel +/// The output convolution consists only of those elements that do not rely on the zero-padding. +/// # Arguments /// -/// The output convolution consists only of those elements that do not rely on the zero-padding. -/// -pub fn convolve_valid, Q: Storage>( - vector: Vector, - kernel: Vector, -) -> Matrix> { +/// * `vector` - A Vector with size > 0 +/// * `kernel` - A Vector with size > 0 +/// +/// # Note: +/// This function is commutative. If kernel > vector, +/// they will swap their roles as in +/// (self, kernel) = (kernel,self) +/// +/// # Example +/// +/// ``` +/// let vec = Vector3::new(1.0,2.0,3.0); +/// let ker = Vector2::new(0.4,0.6); +/// let convolve = convolve_valid(vec,ker); +/// ``` +pub fn convolve_valid( + vector: Vector, + kernel: Vector, +) -> VectorN, U1>> +where + N: Real, + D1: DimSub, + D2: DimSub>, + DimDiff: DimAdd, + S1: Storage, + S2: Storage, + DimDiff: DimName, + DefaultAllocator: Allocator, U1>> +{ + let vec = vector.len(); let ker = kernel.len(); @@ -76,12 +115,10 @@ pub fn convolve_valid, Q: Storage vec { return convolve_valid(kernel, vector); } + let result_len = vector.data.shape().0.sub(kernel.data.shape().0).add(U1); + let mut conv = VectorN::zeros_generic(result_len, U1); - let newlen = vec - ker + 1; - - let mut conv = DVector::::zeros(newlen); - - for i in 0..newlen { + for i in 0..(vec - ker + 1) { for j in 0..ker { conv[i] += vector[i + j] * kernel[ker - j - 1]; } @@ -89,13 +126,38 @@ pub fn convolve_valid, Q: Storage, Q: Storage>( - vector: Vector, - kernel: Vector, -) -> Matrix> { +/// # Arguments +/// +/// * `vector` - A Vector with size > 0 +/// * `kernel` - A Vector with size > 0 +/// +/// # Note: +/// This function is commutative. If kernel > vector, +/// they will swap their roles as in +/// (self, kernel) = (kernel,self) +/// +/// # Example +/// +/// ``` +/// let vec = Vector3::new(1.0,2.0,3.0); +/// let ker = Vector2::new(0.4,0.6); +/// let convolve = convolve_same(vec,ker); +/// ``` +pub fn convolve_same( + vector: Vector, + kernel: Vector, +) -> VectorN> +where + N: Real, + D1: DimMax, + D2: DimMax>, + S1: Storage, + S2: Storage, + DimMaximum: Dim, + DefaultAllocator: Allocator>, +{ let vec = vector.len(); let ker = kernel.len(); @@ -107,12 +169,13 @@ pub fn convolve_same, Q: Storage return convolve_same(kernel, vector); } - let mut conv = DVector::::zeros(vec); + let result_len = vector.data.shape().0.max(kernel.data.shape().0); + let mut conv = VectorN::zeros_generic(result_len, U1); for i in 0..vec { for j in 0..ker { let val = if i + j < 1 || i + j >= vec + 1 { - zero::() + zero::() } else { vector[i + j - 1] }; @@ -121,4 +184,3 @@ pub fn convolve_same, Q: Storage } conv } - diff --git a/tests/linalg/convolution.rs b/tests/linalg/convolution.rs index ef3a02db..454c29c6 100644 --- a/tests/linalg/convolution.rs +++ b/tests/linalg/convolution.rs @@ -1,5 +1,7 @@ +#[allow(unused_imports)] // remove after fixing unit test use na::linalg::{convolve_full,convolve_valid,convolve_same}; -use na::{Vector2,Vector4,DVector}; +#[allow(unused_imports)] +use na::{Vector2,Vector3,Vector4,Vector5,DVector}; // // Should mimic calculations in Python's scipy library @@ -10,40 +12,70 @@ use na::{Vector2,Vector4,DVector}; // array([ 1, 4, 7, 10]) #[test] fn convolve_same_check(){ - let vec = Vector4::new(1.0,2.0,3.0,4.0); - let ker = Vector2::new(1.0,2.0); + let vec_s = Vector4::new(1.0,2.0,3.0,4.0); + let ker_s = Vector2::new(1.0,2.0); - let actual = DVector::from_vec(4, vec![1.0,4.0,7.0,10.0]); + let actual_s = Vector4::from_vec(vec![1.0,4.0,7.0,10.0]); - let expected = convolve_same(vec,ker); + let expected_s = convolve_same(vec_s,ker_s); + let expected_s_r = convolve_same(ker_s,vec_s); - assert!(relative_eq!(actual, expected, epsilon = 1.0e-7)); + assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); + assert!(relative_eq!(actual_s, expected_s_r, epsilon = 1.0e-7)); + + let vec_d = DVector::from_vec(4,vec![1.0,2.0,3.0,4.0]); + let ker_d = DVector::from_vec(2,vec![1.0,2.0]); + + let actual_d = DVector::from_vec(4,vec![1.0,4.0,7.0,10.0]); + + let expected_d = convolve_same(vec_d.clone(),ker_d.clone()); + let expected_d_r = convolve_same(ker_d,vec_d); + + assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); + assert!(relative_eq!(actual_d, expected_d_r, epsilon = 1.0e-7)); } -// >>> convolve([1,2,3,4],[1,2],"valid") +// >>> convolve([1,2,3,4],[1,2],"full") // array([ 1, 4, 7, 10, 8]) #[test] fn convolve_full_check(){ - let vec = Vector4::new(1.0,2.0,3.0,4.0); - let ker = Vector2::new(1.0,2.0); + let vec_s = Vector4::new(1.0,2.0,3.0,4.0); + let ker_s = Vector2::new(1.0,2.0); - let actual = DVector::from_vec(5, vec![1.0,4.0,7.0,10.0,8.0]); + let actual_s = Vector5::new(1.0,4.0,7.0,10.0,8.0); - let expected = convolve_full(vec,ker); + let expected_s = convolve_full(vec_s,ker_s); + let expected_s_r = convolve_full(ker_s,vec_s); - assert!(relative_eq!(actual, expected, epsilon = 1.0e-7)); + assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); + assert!(relative_eq!(actual_s, expected_s_r, epsilon = 1.0e-7)); + + let vec_d = DVector::from_vec(4,vec![1.0,2.0,3.0,4.0]); + let ker_d = DVector::from_vec(2,vec![1.0,2.0]); + + let actual_d = DVector::from_vec(5,vec![1.0,4.0,7.0,10.0,8.0]); + + let expected_d = convolve_full(vec_d.clone(),ker_d.clone()); + let expected_d_r = convolve_full(ker_d,vec_d); + + assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); + assert!(relative_eq!(actual_d, expected_d_r, epsilon = 1.0e-7)); } // >>> convolve([1,2,3,4],[1,2],"valid") // array([ 4, 7, 10]) -#[test] -fn convolve_valid_check(){ - let vec = Vector4::new(1.0,2.0,3.0,4.0); - let ker = Vector2::new(1.0,2.0); +// #[test] +// fn convolve_valid_check(){ +// let vec = Vector4::new(1.0,2.0,3.0,4.0); +// let ker = Vector2::new(1.0,2.0); - let actual = DVector::from_vec(3, vec![4.0,7.0,10.0]); +// let actual = Vector3::from_vec(vec![4.0,7.0,10.0]); - let expected = convolve_valid(vec,ker); +// let expected1 = convolve_valid(vec, ker); +// let expected2 = convolve_valid(ker, vec); - assert!(relative_eq!(actual, expected, epsilon = 1.0e-7)); -} \ No newline at end of file + +// assert!(relative_eq!(actual, expected1, epsilon = 1.0e-7)); +// assert!(relative_eq!(actual, expected2, epsilon = 1.0e-7)); + +// } \ No newline at end of file From 2a2debf58be6d30c043cfce9285238cb69605ee1 Mon Sep 17 00:00:00 2001 From: Nathan Date: Wed, 20 Feb 2019 20:32:09 -0600 Subject: [PATCH 10/22] Fixing documentation --- examples/convolution.rs | 18 +++++++++++++--- src/linalg/convolution.rs | 44 +++++++++------------------------------ 2 files changed, 25 insertions(+), 37 deletions(-) diff --git a/examples/convolution.rs b/examples/convolution.rs index 5290440c..d23ed389 100644 --- a/examples/convolution.rs +++ b/examples/convolution.rs @@ -1,5 +1,17 @@ -fn main(){ - let (x,y) = (1,2); +extern crate nalgebra; +use nalgebra::{Vector2,Vector3,Vector4,Vector5,convolve_full,convolve_same,convolve_valid}; - println!("{}", x); +fn main(){ + let vec = Vector4::new(1.0,2.0,3.0,4.0); + let ker = Vector3::new(1.0,2.0,2.1); + + let actual = Vector5::from_vec(vec![1.0,4.0,7.0,10.0,8.0]); + + let expected = convolve_full(vec,ker); + let expected2 = convolve_same(vec,ker); + // let expected3 = convolve_valid(vec,ker); + println!("{}", actual); + println!("{}", expected); + println!("{}", expected2); + // println!("{}", expected3); } \ No newline at end of file diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs index c587f8f2..69ae3f5e 100644 --- a/src/linalg/convolution.rs +++ b/src/linalg/convolution.rs @@ -12,18 +12,10 @@ use {zero, Real, Vector, VectorN, U1}; /// * `vector` - A Vector with size > 0 /// * `kernel` - A Vector with size > 0 /// -/// # Note: -/// This function is commutative. If kernel > vector, -/// they will swap their roles as in -/// (self, kernel) = (kernel,self) +/// This function is commutative. If kernel > vector, +/// they will swap their roles as in +/// (self, kernel) = (kernel,self) /// -/// # Example -/// -/// ``` -/// let vec = Vector3::new(1.0,2.0,3.0); -/// let ker = Vector2::new(0.4,0.6); -/// let convolve = convolve_full(vec,ker); -/// ``` pub fn convolve_full( vector: Vector, kernel: Vector, @@ -77,19 +69,11 @@ where /// /// * `vector` - A Vector with size > 0 /// * `kernel` - A Vector with size > 0 +/// +/// This function is commutative. If kernel > vector, +/// they will swap their roles as in +/// (self, kernel) = (kernel,self) /// -/// # Note: -/// This function is commutative. If kernel > vector, -/// they will swap their roles as in -/// (self, kernel) = (kernel,self) -/// -/// # Example -/// -/// ``` -/// let vec = Vector3::new(1.0,2.0,3.0); -/// let ker = Vector2::new(0.4,0.6); -/// let convolve = convolve_valid(vec,ker); -/// ``` pub fn convolve_valid( vector: Vector, kernel: Vector, @@ -133,18 +117,10 @@ where /// * `vector` - A Vector with size > 0 /// * `kernel` - A Vector with size > 0 /// -/// # Note: -/// This function is commutative. If kernel > vector, -/// they will swap their roles as in -/// (self, kernel) = (kernel,self) +/// This function is commutative. If kernel > vector, +/// they will swap their roles as in +/// (self, kernel) = (kernel,self) /// -/// # Example -/// -/// ``` -/// let vec = Vector3::new(1.0,2.0,3.0); -/// let ker = Vector2::new(0.4,0.6); -/// let convolve = convolve_same(vec,ker); -/// ``` pub fn convolve_same( vector: Vector, kernel: Vector, From 20e9c6f480430caf3ac263b291313d8c85265bb6 Mon Sep 17 00:00:00 2001 From: Jack Wrenn Date: Sat, 23 Feb 2019 09:02:27 -0500 Subject: [PATCH 11/22] Implement `iter::Sum` for `DMatrix` (#552) Fixes #514. --- src/base/ops.rs | 62 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 61 insertions(+), 1 deletion(-) diff --git a/src/base/ops.rs b/src/base/ops.rs index 96d4626f..301232f2 100644 --- a/src/base/ops.rs +++ b/src/base/ops.rs @@ -11,7 +11,7 @@ use base::allocator::{Allocator, SameShapeAllocator, SameShapeC, SameShapeR}; use base::constraint::{ AreMultipliable, DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint, }; -use base::dimension::{Dim, DimMul, DimName, DimProd}; +use base::dimension::{Dim, DimMul, DimName, DimProd, Dynamic}; use base::storage::{ContiguousStorageMut, Storage, StorageMut}; use base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, MatrixSum, Scalar}; @@ -384,6 +384,36 @@ where } } +impl iter::Sum for MatrixMN +where + N: Scalar + ClosedAdd + Zero, + DefaultAllocator: Allocator, +{ + /// # Example + /// ``` + /// # use nalgebra::DVector; + /// assert_eq!(vec![DVector::repeat(3, 1.0f64), + /// DVector::repeat(3, 1.0f64), + /// DVector::repeat(3, 1.0f64)].into_iter().sum::>(), + /// DVector::repeat(3, 1.0f64) + DVector::repeat(3, 1.0f64) + DVector::repeat(3, 1.0f64)); + /// ``` + /// + /// # Panics + /// Panics if the iterator is empty: + /// ```should_panic + /// # use std::iter; + /// # use nalgebra::DMatrix; + /// iter::empty::>().sum::>(); // panics! + /// ``` + fn sum>>(mut iter: I) -> MatrixMN { + if let Some(first) = iter.next() { + iter.fold(first, |acc, x| acc + x) + } else { + panic!("Cannot compute `sum` of empty iterator.") + } + } +} + impl<'a, N, R: DimName, C: DimName> iter::Sum<&'a MatrixMN> for MatrixMN where N: Scalar + ClosedAdd + Zero, @@ -394,6 +424,36 @@ where } } +impl<'a, N, C: Dim> iter::Sum<&'a MatrixMN> for MatrixMN +where + N: Scalar + ClosedAdd + Zero, + DefaultAllocator: Allocator, +{ + /// # Example + /// ``` + /// # use nalgebra::DVector; + /// let v = &DVector::repeat(3, 1.0f64); + /// + /// assert_eq!(vec![v, v, v].into_iter().sum::>(), + /// v + v + v); + /// ``` + /// + /// # Panics + /// Panics if the iterator is empty: + /// ```should_panic + /// # use std::iter; + /// # use nalgebra::DMatrix; + /// iter::empty::<&DMatrix>().sum::>(); // panics! + /// ``` + fn sum>>(mut iter: I) -> MatrixMN { + if let Some(first) = iter.next() { + iter.fold(first.clone(), |acc, x| acc + x) + } else { + panic!("Cannot compute `sum` of empty iterator.") + } + } +} + /* * * Multiplication From 704331be4fbce372dedbad07798086d4dccf9cc2 Mon Sep 17 00:00:00 2001 From: adamnemecek Date: Sat, 23 Feb 2019 07:03:01 -0800 Subject: [PATCH 12/22] added Quaternion accessors for r,i,j,k, refactored conjugate to use imag (#551) --- src/geometry/quaternion.rs | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index b58d885c..49656b2c 100644 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -124,6 +124,12 @@ impl Quaternion { Self::from(self.coords.normalize()) } + /// The imaginary part of this quaternion. + #[inline] + pub fn imag(&self) -> Vector3 { + self.coords.xyz() + } + /// The conjugate of this quaternion. /// /// # Example @@ -135,13 +141,7 @@ impl Quaternion { /// ``` #[inline] pub fn conjugate(&self) -> Self { - let v = Vector4::new( - -self.coords[0], - -self.coords[1], - -self.coords[2], - self.coords[3], - ); - Self::from(v) + Self::from_parts(self.w, -self.imag()) } /// Inverts this quaternion if it is not zero. From 2a2e9d7f8e57af6383f07e77a50a2a89b2407ee9 Mon Sep 17 00:00:00 2001 From: Terence Date: Sun, 24 Feb 2019 11:29:27 -0500 Subject: [PATCH 13/22] Add the `transform` methods as inherent methods on geometric types This adds `transform_point`, `transform_vector`, `inverse_transform_point` and `inverse_transform_vector` as inherent methods on the applicable geometric transformation structures, such that they can be used without the need to import the `Transformation` and `ProjectiveTransformation` traits from `alga`. --- src/geometry/isometry.rs | 89 ++++++++++++++++++++++++++++++- src/geometry/isometry_alga.rs | 9 ++-- src/geometry/quaternion.rs | 80 ++++++++++++++++++++++++++- src/geometry/quaternion_alga.rs | 10 ++-- src/geometry/rotation.rs | 80 ++++++++++++++++++++++++++- src/geometry/rotation_alga.rs | 8 +-- src/geometry/similarity.rs | 79 ++++++++++++++++++++++++++- src/geometry/similarity_alga.rs | 8 +-- src/geometry/transform.rs | 53 +++++++++++++++++- src/geometry/transform_alga.rs | 8 +-- src/geometry/translation.rs | 38 ++++++++++++- src/geometry/translation_alga.rs | 4 +- src/geometry/unit_complex.rs | 70 +++++++++++++++++++++++- src/geometry/unit_complex_alga.rs | 10 ++-- 14 files changed, 506 insertions(+), 40 deletions(-) mode change 100644 => 100755 src/geometry/isometry.rs mode change 100644 => 100755 src/geometry/isometry_alga.rs mode change 100644 => 100755 src/geometry/quaternion.rs mode change 100644 => 100755 src/geometry/quaternion_alga.rs mode change 100644 => 100755 src/geometry/rotation.rs mode change 100644 => 100755 src/geometry/rotation_alga.rs mode change 100644 => 100755 src/geometry/similarity.rs mode change 100644 => 100755 src/geometry/similarity_alga.rs mode change 100644 => 100755 src/geometry/transform.rs mode change 100644 => 100755 src/geometry/transform_alga.rs mode change 100644 => 100755 src/geometry/translation.rs mode change 100644 => 100755 src/geometry/translation_alga.rs mode change 100644 => 100755 src/geometry/unit_complex.rs mode change 100644 => 100755 src/geometry/unit_complex_alga.rs diff --git a/src/geometry/isometry.rs b/src/geometry/isometry.rs old mode 100644 new mode 100755 index 1097be47..04e77b7a --- a/src/geometry/isometry.rs +++ b/src/geometry/isometry.rs @@ -17,7 +17,7 @@ use alga::linear::Rotation; use base::allocator::Allocator; use base::dimension::{DimName, DimNameAdd, DimNameSum, U1}; use base::storage::Owned; -use base::{DefaultAllocator, MatrixN}; +use base::{DefaultAllocator, MatrixN, VectorN}; use geometry::{Point, Translation}; /// A direct isometry, i.e., a rotation followed by a translation, aka. a rigid-body motion, aka. an element of a Special Euclidean (SE) group. @@ -254,6 +254,93 @@ where DefaultAllocator: Allocator pub fn append_rotation_wrt_center_mut(&mut self, r: &R) { self.rotation = self.rotation.append_rotation(r); } + + /// Transform the given point by this isometry. + /// + /// # Example + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Isometry3, Translation3, UnitQuaternion, Vector3, Point3}; + /// let tra = Translation3::new(0.0, 0.0, 3.0); + /// let rot = UnitQuaternion::from_scaled_axis(Vector3::y() * f32::consts::FRAC_PI_2); + /// let iso = Isometry3::from_parts(tra, rot); + /// + /// let transformed_point = iso.transform_point(&Point3::new(1.0, 2.0, 3.0)); + /// assert_relative_eq!(transformed_point, Point3::new(3.0, 2.0, 2.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn transform_point(&self, pt: &Point) -> Point { + self * pt + } + + /// Transform the given vector by this isometry, ignoring the translation + /// component of the isometry. + /// + /// # Example + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Isometry3, Translation3, UnitQuaternion, Vector3}; + /// let tra = Translation3::new(0.0, 0.0, 3.0); + /// let rot = UnitQuaternion::from_scaled_axis(Vector3::y() * f32::consts::FRAC_PI_2); + /// let iso = Isometry3::from_parts(tra, rot); + /// + /// let transformed_point = iso.transform_vector(&Vector3::new(1.0, 2.0, 3.0)); + /// assert_relative_eq!(transformed_point, Vector3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn transform_vector(&self, v: &VectorN) -> VectorN { + self * v + } + + /// Transform the given point by the inverse of this isometry. This may be + /// less expensive than computing the entire isometry inverse and then + /// transforming the point. + /// + /// # Example + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Isometry3, Translation3, UnitQuaternion, Vector3, Point3}; + /// let tra = Translation3::new(0.0, 0.0, 3.0); + /// let rot = UnitQuaternion::from_scaled_axis(Vector3::y() * f32::consts::FRAC_PI_2); + /// let iso = Isometry3::from_parts(tra, rot); + /// + /// let transformed_point = iso.inverse_transform_point(&Point3::new(1.0, 2.0, 3.0)); + /// assert_relative_eq!(transformed_point, Point3::new(0.0, 2.0, 1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn inverse_transform_point(&self, pt: &Point) -> Point { + self.rotation + .inverse_transform_point(&(pt - &self.translation.vector)) + } + + /// Transform the given vector by the inverse of this isometry, ignoring the + /// translation component of the isometry. This may be + /// less expensive than computing the entire isometry inverse and then + /// transforming the point. + /// + /// # Example + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Isometry3, Translation3, UnitQuaternion, Vector3}; + /// let tra = Translation3::new(0.0, 0.0, 3.0); + /// let rot = UnitQuaternion::from_scaled_axis(Vector3::y() * f32::consts::FRAC_PI_2); + /// let iso = Isometry3::from_parts(tra, rot); + /// + /// let transformed_point = iso.inverse_transform_vector(&Vector3::new(1.0, 2.0, 3.0)); + /// assert_relative_eq!(transformed_point, Vector3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn inverse_transform_vector(&self, v: &VectorN) -> VectorN { + self.rotation.inverse_transform_vector(v) + } } // NOTE: we don't require `R: Rotation<...>` here because this is not useful for the implementation diff --git a/src/geometry/isometry_alga.rs b/src/geometry/isometry_alga.rs old mode 100644 new mode 100755 index fee8b63d..a327f3fb --- a/src/geometry/isometry_alga.rs +++ b/src/geometry/isometry_alga.rs @@ -85,12 +85,12 @@ where { #[inline] fn transform_point(&self, pt: &Point) -> Point { - self * pt + self.transform_point(pt) } #[inline] fn transform_vector(&self, v: &VectorN) -> VectorN { - self * v + self.transform_vector(v) } } @@ -101,13 +101,12 @@ where { #[inline] fn inverse_transform_point(&self, pt: &Point) -> Point { - self.rotation - .inverse_transform_point(&(pt - &self.translation.vector)) + self.inverse_transform_point(pt) } #[inline] fn inverse_transform_vector(&self, v: &VectorN) -> VectorN { - self.rotation.inverse_transform_vector(v) + self.inverse_transform_vector(v) } } diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs old mode 100644 new mode 100755 index 96e4f19a..365e55e2 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -19,7 +19,7 @@ use base::dimension::{U1, U3, U4}; use base::storage::{CStride, RStride}; use base::{Matrix3, MatrixN, MatrixSlice, MatrixSliceMut, Unit, Vector3, Vector4}; -use geometry::Rotation; +use geometry::{Point3, Rotation}; /// A quaternion. See the type alias `UnitQuaternion = Unit` for a quaternion /// that may be used as a rotation. @@ -1005,6 +1005,84 @@ impl UnitQuaternion { pub fn to_homogeneous(&self) -> MatrixN { self.to_rotation_matrix().to_homogeneous() } + + /// Rotate a point by this unit quaternion. + /// + /// # Example + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{UnitQuaternion, Vector3, Point3}; + /// let rot = UnitQuaternion::from_axis_angle(&Vector3::y_axis(), f32::consts::FRAC_PI_2); + /// let transformed_point = rot.transform_point(&Point3::new(1.0, 2.0, 3.0)); + /// + /// assert_relative_eq!(transformed_point, Point3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn transform_point(&self, pt: &Point3) -> Point3 { + self * pt + } + + /// Rotate a vector by this unit quaternion. + /// + /// # Example + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{UnitQuaternion, Vector3}; + /// let rot = UnitQuaternion::from_axis_angle(&Vector3::y_axis(), f32::consts::FRAC_PI_2); + /// let transformed_vector = rot.transform_vector(&Vector3::new(1.0, 2.0, 3.0)); + /// + /// assert_relative_eq!(transformed_vector, Vector3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn transform_vector(&self, v: &Vector3) -> Vector3 { + self * v + } + + /// Rotate a point by the inverse of this unit quaternion. This may be + /// cheaper than inverting the unit quaternion and transforming the + /// point. + /// + /// # Example + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{UnitQuaternion, Vector3, Point3}; + /// let rot = UnitQuaternion::from_axis_angle(&Vector3::y_axis(), f32::consts::FRAC_PI_2); + /// let transformed_point = rot.inverse_transform_point(&Point3::new(1.0, 2.0, 3.0)); + /// + /// assert_relative_eq!(transformed_point, Point3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn inverse_transform_point(&self, pt: &Point3) -> Point3 { + // FIXME: would it be useful performancewise not to call inverse explicitly (i-e. implement + // the inverse transformation explicitly here) ? + self.inverse() * pt + } + + /// Rotate a vector by the inverse of this unit quaternion. This may be + /// cheaper than inverting the unit quaternion and transforming the + /// vector. + /// + /// # Example + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{UnitQuaternion, Vector3}; + /// let rot = UnitQuaternion::from_axis_angle(&Vector3::y_axis(), f32::consts::FRAC_PI_2); + /// let transformed_vector = rot.inverse_transform_vector(&Vector3::new(1.0, 2.0, 3.0)); + /// + /// assert_relative_eq!(transformed_vector, Vector3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn inverse_transform_vector(&self, v: &Vector3) -> Vector3 { + self.inverse() * v + } } impl fmt::Display for UnitQuaternion { diff --git a/src/geometry/quaternion_alga.rs b/src/geometry/quaternion_alga.rs old mode 100644 new mode 100755 index 529d27f4..d931b206 --- a/src/geometry/quaternion_alga.rs +++ b/src/geometry/quaternion_alga.rs @@ -197,26 +197,24 @@ impl_structures!( impl Transformation> for UnitQuaternion { #[inline] fn transform_point(&self, pt: &Point3) -> Point3 { - self * pt + self.transform_point(pt) } #[inline] fn transform_vector(&self, v: &Vector3) -> Vector3 { - self * v + self.transform_vector(v) } } impl ProjectiveTransformation> for UnitQuaternion { #[inline] fn inverse_transform_point(&self, pt: &Point3) -> Point3 { - // FIXME: would it be useful performancewise not to call inverse explicitly (i-e. implement - // the inverse transformation explicitly here) ? - self.inverse() * pt + self.inverse_transform_point(pt) } #[inline] fn inverse_transform_vector(&self, v: &Vector3) -> Vector3 { - self.inverse() * v + self.inverse_transform_vector(v) } } diff --git a/src/geometry/rotation.rs b/src/geometry/rotation.rs old mode 100644 new mode 100755 index 6c50230e..961951f1 --- a/src/geometry/rotation.rs +++ b/src/geometry/rotation.rs @@ -18,7 +18,9 @@ use alga::general::Real; use base::allocator::Allocator; use base::dimension::{DimName, DimNameAdd, DimNameSum, U1}; -use base::{DefaultAllocator, MatrixN, Scalar}; +use base::{DefaultAllocator, MatrixN, Scalar, VectorN}; + +use geometry::Point; /// A rotation matrix. #[repr(C)] @@ -351,6 +353,82 @@ where DefaultAllocator: Allocator } } +impl Rotation +where DefaultAllocator: Allocator + Allocator +{ + /// Rotate the given point. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Point3, Rotation2, Rotation3, UnitQuaternion, Vector3}; + /// let rot = Rotation3::new(Vector3::y() * f32::consts::FRAC_PI_2); + /// let transformed_point = rot.transform_point(&Point3::new(1.0, 2.0, 3.0)); + /// + /// assert_relative_eq!(transformed_point, Point3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn transform_point(&self, pt: &Point) -> Point { + self * pt + } + + /// Rotate the given vector. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Rotation2, Rotation3, UnitQuaternion, Vector3}; + /// let rot = Rotation3::new(Vector3::y() * f32::consts::FRAC_PI_2); + /// let transformed_vector = rot.transform_vector(&Vector3::new(1.0, 2.0, 3.0)); + /// + /// assert_relative_eq!(transformed_vector, Vector3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn transform_vector(&self, v: &VectorN) -> VectorN { + self * v + } + + /// Rotate the given point by the inverse of this rotation. This may be + /// cheaper than inverting the rotation and then transforming the given + /// point. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Point3, Rotation2, Rotation3, UnitQuaternion, Vector3}; + /// let rot = Rotation3::new(Vector3::y() * f32::consts::FRAC_PI_2); + /// let transformed_point = rot.inverse_transform_point(&Point3::new(1.0, 2.0, 3.0)); + /// + /// assert_relative_eq!(transformed_point, Point3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn inverse_transform_point(&self, pt: &Point) -> Point { + Point::from(self.inverse_transform_vector(&pt.coords)) + } + + /// Rotate the given vector by the inverse of this rotation. This may be + /// cheaper than inverting the rotation and then transforming the given + /// vector. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Rotation2, Rotation3, UnitQuaternion, Vector3}; + /// let rot = Rotation3::new(Vector3::y() * f32::consts::FRAC_PI_2); + /// let transformed_vector = rot.inverse_transform_vector(&Vector3::new(1.0, 2.0, 3.0)); + /// + /// assert_relative_eq!(transformed_vector, Vector3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn inverse_transform_vector(&self, v: &VectorN) -> VectorN { + self.matrix().tr_mul(v) + } +} + impl Eq for Rotation where DefaultAllocator: Allocator {} impl PartialEq for Rotation diff --git a/src/geometry/rotation_alga.rs b/src/geometry/rotation_alga.rs old mode 100644 new mode 100755 index 18c47b41..cd411c8c --- a/src/geometry/rotation_alga.rs +++ b/src/geometry/rotation_alga.rs @@ -75,12 +75,12 @@ where DefaultAllocator: Allocator + Allocator { #[inline] fn transform_point(&self, pt: &Point) -> Point { - self * pt + self.transform_point(pt) } #[inline] fn transform_vector(&self, v: &VectorN) -> VectorN { - self * v + self.transform_vector(v) } } @@ -89,12 +89,12 @@ where DefaultAllocator: Allocator + Allocator { #[inline] fn inverse_transform_point(&self, pt: &Point) -> Point { - Point::from(self.inverse_transform_vector(&pt.coords)) + self.inverse_transform_point(pt) } #[inline] fn inverse_transform_vector(&self, v: &VectorN) -> VectorN { - self.matrix().tr_mul(v) + self.inverse_transform_vector(v) } } diff --git a/src/geometry/similarity.rs b/src/geometry/similarity.rs old mode 100644 new mode 100755 index f321d575..dcb0f20a --- a/src/geometry/similarity.rs +++ b/src/geometry/similarity.rs @@ -16,7 +16,7 @@ use alga::linear::Rotation; use base::allocator::Allocator; use base::dimension::{DimName, DimNameAdd, DimNameSum, U1}; use base::storage::Owned; -use base::{DefaultAllocator, MatrixN}; +use base::{DefaultAllocator, MatrixN, VectorN}; use geometry::{Isometry, Point, Translation}; /// A similarity, i.e., an uniform scaling, followed by a rotation, followed by a translation. @@ -238,6 +238,83 @@ where pub fn append_rotation_wrt_center_mut(&mut self, r: &R) { self.isometry.append_rotation_wrt_center_mut(r) } + + /// Transform the given point by this similarity. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Point3, Similarity3, Vector3}; + /// let axisangle = Vector3::y() * f32::consts::FRAC_PI_2; + /// let translation = Vector3::new(1.0, 2.0, 3.0); + /// let sim = Similarity3::new(translation, axisangle, 3.0); + /// let transformed_point = sim.transform_point(&Point3::new(4.0, 5.0, 6.0)); + /// assert_relative_eq!(transformed_point, Point3::new(19.0, 17.0, -9.0), epsilon = 1.0e-5); + /// ``` + #[inline] + pub fn transform_point(&self, pt: &Point) -> Point { + self * pt + } + + /// Transform the given vector by this similarity, ignoring the translational + /// component. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Similarity3, Vector3}; + /// let axisangle = Vector3::y() * f32::consts::FRAC_PI_2; + /// let translation = Vector3::new(1.0, 2.0, 3.0); + /// let sim = Similarity3::new(translation, axisangle, 3.0); + /// let transformed_vector = sim.transform_vector(&Vector3::new(4.0, 5.0, 6.0)); + /// assert_relative_eq!(transformed_vector, Vector3::new(18.0, 15.0, -12.0), epsilon = 1.0e-5); + /// ``` + #[inline] + pub fn transform_vector(&self, v: &VectorN) -> VectorN { + self * v + } + + /// Transform the given point by the inverse of this similarity. This may + /// be cheaper than inverting the similarity and then transforming the + /// given point. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Point3, Similarity3, Vector3}; + /// let axisangle = Vector3::y() * f32::consts::FRAC_PI_2; + /// let translation = Vector3::new(1.0, 2.0, 3.0); + /// let sim = Similarity3::new(translation, axisangle, 2.0); + /// let transformed_point = sim.inverse_transform_point(&Point3::new(4.0, 5.0, 6.0)); + /// assert_relative_eq!(transformed_point, Point3::new(-1.5, 1.5, 1.5), epsilon = 1.0e-5); + /// ``` + #[inline] + pub fn inverse_transform_point(&self, pt: &Point) -> Point { + self.isometry.inverse_transform_point(pt) / self.scaling() + } + + /// Transform the given vector by the inverse of this similarity, + /// ignoring the translational component. This may be cheaper than + /// inverting the similarity and then transforming the given vector. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use std::f32; + /// # use nalgebra::{Similarity3, Vector3}; + /// let axisangle = Vector3::y() * f32::consts::FRAC_PI_2; + /// let translation = Vector3::new(1.0, 2.0, 3.0); + /// let sim = Similarity3::new(translation, axisangle, 2.0); + /// let transformed_vector = sim.inverse_transform_vector(&Vector3::new(4.0, 5.0, 6.0)); + /// assert_relative_eq!(transformed_vector, Vector3::new(-3.0, 2.5, 2.0), epsilon = 1.0e-5); + /// ``` + #[inline] + pub fn inverse_transform_vector(&self, v: &VectorN) -> VectorN { + self.isometry.inverse_transform_vector(v) / self.scaling() + } } // NOTE: we don't require `R: Rotation<...>` here because this is not useful for the implementation diff --git a/src/geometry/similarity_alga.rs b/src/geometry/similarity_alga.rs old mode 100644 new mode 100755 index e8a6b154..9b723029 --- a/src/geometry/similarity_alga.rs +++ b/src/geometry/similarity_alga.rs @@ -82,12 +82,12 @@ where { #[inline] fn transform_point(&self, pt: &Point) -> Point { - self * pt + self.transform_point(pt) } #[inline] fn transform_vector(&self, v: &VectorN) -> VectorN { - self * v + self.transform_vector(v) } } @@ -98,12 +98,12 @@ where { #[inline] fn inverse_transform_point(&self, pt: &Point) -> Point { - self.isometry.inverse_transform_point(pt) / self.scaling() + self.inverse_transform_point(pt) } #[inline] fn inverse_transform_vector(&self, v: &VectorN) -> VectorN { - self.isometry.inverse_transform_vector(v) / self.scaling() + self.inverse_transform_vector(v) } } diff --git a/src/geometry/transform.rs b/src/geometry/transform.rs old mode 100644 new mode 100755 index 0d5d9c4a..56a65992 --- a/src/geometry/transform.rs +++ b/src/geometry/transform.rs @@ -6,12 +6,14 @@ use std::marker::PhantomData; #[cfg(feature = "serde-serialize")] use serde::{Deserialize, Deserializer, Serialize, Serializer}; -use alga::general::Real; +use alga::general::{Real, TwoSidedInverse}; use base::allocator::Allocator; use base::dimension::{DimName, DimNameAdd, DimNameSum, U1}; use base::storage::Owned; -use base::{DefaultAllocator, MatrixN}; +use base::{DefaultAllocator, MatrixN, VectorN}; + +use geometry::Point; /// Trait implemented by phantom types identifying the projective transformation type. /// @@ -452,6 +454,53 @@ where DefaultAllocator: Allocator, DimNameSum> } } +impl, C> Transform +where + N: Real, + C: TCategory, + DefaultAllocator: Allocator, DimNameSum> + + Allocator> + + Allocator + + Allocator, +{ + /// Transform the given point by this transformation. + #[inline] + pub fn transform_point(&self, pt: &Point) -> Point { + self * pt + } + + /// Transform the given vector by this transformation, ignoring the + /// translational component of the transformation. + #[inline] + pub fn transform_vector(&self, v: &VectorN) -> VectorN { + self * v + } +} + +impl, C: TCategory> Transform +where C: SubTCategoryOf, + DefaultAllocator: Allocator, DimNameSum> + + Allocator> + + Allocator + + Allocator, +{ + /// Transform the given point by the inverse of this transformation. + /// This may be cheaper than inverting the transformation and transforming + /// the point. + #[inline] + pub fn inverse_transform_point(&self, pt: &Point) -> Point { + self.two_sided_inverse() * pt + } + + /// Transform the given vector by the inverse of this transformation. + /// This may be cheaper than inverting the transformation and transforming + /// the vector. + #[inline] + pub fn inverse_transform_vector(&self, v: &VectorN) -> VectorN { + self.two_sided_inverse() * v + } +} + impl> Transform where DefaultAllocator: Allocator, DimNameSum> { diff --git a/src/geometry/transform_alga.rs b/src/geometry/transform_alga.rs old mode 100644 new mode 100755 index c5ba675b..d58b1d53 --- a/src/geometry/transform_alga.rs +++ b/src/geometry/transform_alga.rs @@ -96,12 +96,12 @@ where { #[inline] fn transform_point(&self, pt: &Point) -> Point { - self * pt + self.transform_point(pt) } #[inline] fn transform_vector(&self, v: &VectorN) -> VectorN { - self * v + self.transform_vector(v) } } @@ -116,12 +116,12 @@ where { #[inline] fn inverse_transform_point(&self, pt: &Point) -> Point { - self.two_sided_inverse() * pt + self.inverse_transform_point(pt) } #[inline] fn inverse_transform_vector(&self, v: &VectorN) -> VectorN { - self.two_sided_inverse() * v + self.inverse_transform_vector(v) } } diff --git a/src/geometry/translation.rs b/src/geometry/translation.rs old mode 100644 new mode 100755 index 7b1a8957..89960961 --- a/src/geometry/translation.rs +++ b/src/geometry/translation.rs @@ -11,13 +11,15 @@ use serde::{Deserialize, Deserializer, Serialize, Serializer}; #[cfg(feature = "abomonation-serialize")] use abomonation::Abomonation; -use alga::general::{ClosedNeg, Real}; +use alga::general::{ClosedAdd, ClosedNeg, ClosedSub, Real}; use base::allocator::Allocator; use base::dimension::{DimName, DimNameAdd, DimNameSum, U1}; use base::storage::Owned; use base::{DefaultAllocator, MatrixN, Scalar, VectorN}; +use geometry::Point; + /// A translation. #[repr(C)] #[derive(Debug)] @@ -190,6 +192,40 @@ where DefaultAllocator: Allocator } } +impl Translation +where DefaultAllocator: Allocator +{ + /// Translate the given point. + /// + /// # Example + /// ``` + /// # use nalgebra::{Translation3, Point3}; + /// let t = Translation3::new(1.0, 2.0, 3.0); + /// let transformed_point = t.transform_point(&Point3::new(4.0, 5.0, 6.0)); + /// assert_eq!(transformed_point, Point3::new(5.0, 7.0, 9.0)); + #[inline] + pub fn transform_point(&self, pt: &Point) -> Point { + pt + &self.vector + } +} + +impl Translation +where DefaultAllocator: Allocator +{ + /// Translate the given point by the inverse of this translation. + /// + /// # Example + /// ``` + /// # use nalgebra::{Translation3, Point3}; + /// let t = Translation3::new(1.0, 2.0, 3.0); + /// let transformed_point = t.inverse_transform_point(&Point3::new(4.0, 5.0, 6.0)); + /// assert_eq!(transformed_point, Point3::new(3.0, 3.0, 3.0)); + #[inline] + pub fn inverse_transform_point(&self, pt: &Point) -> Point { + pt - &self.vector + } +} + impl Eq for Translation where DefaultAllocator: Allocator {} impl PartialEq for Translation diff --git a/src/geometry/translation_alga.rs b/src/geometry/translation_alga.rs old mode 100644 new mode 100755 index fdd24014..c6157a17 --- a/src/geometry/translation_alga.rs +++ b/src/geometry/translation_alga.rs @@ -76,7 +76,7 @@ where DefaultAllocator: Allocator { #[inline] fn transform_point(&self, pt: &Point) -> Point { - pt + &self.vector + self.transform_point(pt) } #[inline] @@ -90,7 +90,7 @@ where DefaultAllocator: Allocator { #[inline] fn inverse_transform_point(&self, pt: &Point) -> Point { - pt - &self.vector + self.inverse_transform_point(pt) } #[inline] diff --git a/src/geometry/unit_complex.rs b/src/geometry/unit_complex.rs old mode 100644 new mode 100755 index f8d94dc8..bd064ef3 --- a/src/geometry/unit_complex.rs +++ b/src/geometry/unit_complex.rs @@ -3,8 +3,8 @@ use num_complex::Complex; use std::fmt; use alga::general::Real; -use base::{Matrix2, Matrix3, Unit, Vector1}; -use geometry::Rotation2; +use base::{Matrix2, Matrix3, Unit, Vector1, Vector2}; +use geometry::{Rotation2, Point2}; /// A complex number with a norm equal to 1. pub type UnitComplex = Unit>; @@ -251,6 +251,72 @@ impl UnitComplex { pub fn to_homogeneous(&self) -> Matrix3 { self.to_rotation_matrix().to_homogeneous() } + + /// Rotate the given point by this unit complex number. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::{UnitComplex, Point2}; + /// # use std::f32; + /// let rot = UnitComplex::new(f32::consts::FRAC_PI_2); + /// let transformed_point = rot.transform_point(&Point2::new(1.0, 2.0)); + /// assert_relative_eq!(transformed_point, Point2::new(-2.0, 1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn transform_point(&self, pt: &Point2) -> Point2 { + self * pt + } + + /// Rotate the given vector by this unit complex number. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::{UnitComplex, Vector2}; + /// # use std::f32; + /// let rot = UnitComplex::new(f32::consts::FRAC_PI_2); + /// let transformed_vector = rot.transform_vector(&Vector2::new(1.0, 2.0)); + /// assert_relative_eq!(transformed_vector, Vector2::new(-2.0, 1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn transform_vector(&self, v: &Vector2) -> Vector2 { + self * v + } + + /// Rotate the given point by the inverse of this unit complex number. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::{UnitComplex, Point2}; + /// # use std::f32; + /// let rot = UnitComplex::new(f32::consts::FRAC_PI_2); + /// let transformed_point = rot.inverse_transform_point(&Point2::new(1.0, 2.0)); + /// assert_relative_eq!(transformed_point, Point2::new(2.0, -1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn inverse_transform_point(&self, pt: &Point2) -> Point2 { + // FIXME: would it be useful performancewise not to call inverse explicitly (i-e. implement + // the inverse transformation explicitly here) ? + self.inverse() * pt + } + + /// Rotate the given vector by the inverse of this unit complex number. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::{UnitComplex, Vector2}; + /// # use std::f32; + /// let rot = UnitComplex::new(f32::consts::FRAC_PI_2); + /// let transformed_vector = rot.inverse_transform_vector(&Vector2::new(1.0, 2.0)); + /// assert_relative_eq!(transformed_vector, Vector2::new(2.0, -1.0), epsilon = 1.0e-6); + /// ``` + #[inline] + pub fn inverse_transform_vector(&self, v: &Vector2) -> Vector2 { + self.inverse() * v + } } impl fmt::Display for UnitComplex { diff --git a/src/geometry/unit_complex_alga.rs b/src/geometry/unit_complex_alga.rs old mode 100644 new mode 100755 index 21b956d9..a52ea4e8 --- a/src/geometry/unit_complex_alga.rs +++ b/src/geometry/unit_complex_alga.rs @@ -63,12 +63,12 @@ where DefaultAllocator: Allocator { #[inline] fn transform_point(&self, pt: &Point2) -> Point2 { - self * pt + self.transform_point(pt) } #[inline] fn transform_vector(&self, v: &Vector2) -> Vector2 { - self * v + self.transform_vector(v) } } @@ -77,14 +77,12 @@ where DefaultAllocator: Allocator { #[inline] fn inverse_transform_point(&self, pt: &Point2) -> Point2 { - // FIXME: would it be useful performancewise not to call inverse explicitly (i-e. implement - // the inverse transformation explicitly here) ? - self.inverse() * pt + self.inverse_transform_point(pt) } #[inline] fn inverse_transform_vector(&self, v: &Vector2) -> Vector2 { - self.inverse() * v + self.inverse_transform_vector(v) } } From 28525bfc203740fb7f8c159c7c712770a73b4123 Mon Sep 17 00:00:00 2001 From: Nathan Date: Sun, 24 Feb 2019 19:53:09 -0600 Subject: [PATCH 14/22] Restructured usage of convolves, added unit testing. --- examples/convolution.rs | 17 ----- src/linalg/convolution.rs | 64 +++++++------------ tests/linalg/convolution.rs | 123 ++++++++++++++++++++++++------------ 3 files changed, 103 insertions(+), 101 deletions(-) delete mode 100644 examples/convolution.rs diff --git a/examples/convolution.rs b/examples/convolution.rs deleted file mode 100644 index d23ed389..00000000 --- a/examples/convolution.rs +++ /dev/null @@ -1,17 +0,0 @@ -extern crate nalgebra; -use nalgebra::{Vector2,Vector3,Vector4,Vector5,convolve_full,convolve_same,convolve_valid}; - -fn main(){ - let vec = Vector4::new(1.0,2.0,3.0,4.0); - let ker = Vector3::new(1.0,2.0,2.1); - - let actual = Vector5::from_vec(vec![1.0,4.0,7.0,10.0,8.0]); - - let expected = convolve_full(vec,ker); - let expected2 = convolve_same(vec,ker); - // let expected3 = convolve_valid(vec,ker); - println!("{}", actual); - println!("{}", expected); - println!("{}", expected2); - // println!("{}", expected3); -} \ No newline at end of file diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs index 69ae3f5e..49c8b154 100644 --- a/src/linalg/convolution.rs +++ b/src/linalg/convolution.rs @@ -1,20 +1,19 @@ use base::allocator::Allocator; use base::default_allocator::DefaultAllocator; -use base::dimension::{DimAdd, DimDiff, DimMax, DimMaximum, DimName, DimSub, DimSum,Dim}; +use base::dimension::{Dim, DimAdd, DimDiff, DimMax, DimMaximum, DimSub, DimSum}; use std::cmp; use storage::Storage; use {zero, Real, Vector, VectorN, U1}; -/// Returns the convolution of the vector and a kernel +/// Returns the convolution of the target vector and a kernel /// /// # Arguments /// /// * `vector` - A Vector with size > 0 /// * `kernel` - A Vector with size > 0 /// -/// This function is commutative. If kernel > vector, -/// they will swap their roles as in -/// (self, kernel) = (kernel,self) +/// # Errors +/// Inputs must statisfy `vector.len() >= kernel.len() > 0`. /// pub fn convolve_full( vector: Vector, @@ -27,18 +26,13 @@ where DimSum: DimSub, S1: Storage, S2: Storage, - DimSum: Dim, DefaultAllocator: Allocator, U1>>, { let vec = vector.len(); let ker = kernel.len(); - if vec == 0 || ker == 0 { - panic!("Convolve's inputs must not be 0-sized. "); - } - - if ker > vec { - return convolve_full(kernel, vector); + if ker == 0 || ker > vec { + panic!("convolve_full expects `vector.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); } let result_len = vector.data.shape().0.add(kernel.data.shape().0).sub(U1); @@ -61,45 +55,38 @@ where conv } - - /// Returns the convolution of the vector and a kernel /// The output convolution consists only of those elements that do not rely on the zero-padding. /// # Arguments /// /// * `vector` - A Vector with size > 0 /// * `kernel` - A Vector with size > 0 -/// -/// This function is commutative. If kernel > vector, -/// they will swap their roles as in -/// (self, kernel) = (kernel,self) +/// +/// +/// # Errors +/// Inputs must statisfy `vector.len() >= kernel.len() > 0`. /// pub fn convolve_valid( vector: Vector, kernel: Vector, -) -> VectorN, U1>> +) -> VectorN, D2>> where N: Real, - D1: DimSub, - D2: DimSub>, - DimDiff: DimAdd, + D1: DimAdd, + D2: Dim, + DimSum: DimSub, S1: Storage, S2: Storage, - DimDiff: DimName, - DefaultAllocator: Allocator, U1>> + DefaultAllocator: Allocator, D2>>, { - let vec = vector.len(); let ker = kernel.len(); - if vec == 0 || ker == 0 { - panic!("Convolve's inputs must not be 0-sized. "); + if ker == 0 || ker > vec { + panic!("convolve_valid expects `vector.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); } - if ker > vec { - return convolve_valid(kernel, vector); - } - let result_len = vector.data.shape().0.sub(kernel.data.shape().0).add(U1); + let result_len = vector.data.shape().0.add(U1).sub(kernel.data.shape().0); let mut conv = VectorN::zeros_generic(result_len, U1); for i in 0..(vec - ker + 1) { @@ -117,10 +104,8 @@ where /// * `vector` - A Vector with size > 0 /// * `kernel` - A Vector with size > 0 /// -/// This function is commutative. If kernel > vector, -/// they will swap their roles as in -/// (self, kernel) = (kernel,self) -/// +/// # Errors +/// Inputs must statisfy `vector.len() >= kernel.len() > 0`. pub fn convolve_same( vector: Vector, kernel: Vector, @@ -131,18 +116,13 @@ where D2: DimMax>, S1: Storage, S2: Storage, - DimMaximum: Dim, DefaultAllocator: Allocator>, { let vec = vector.len(); let ker = kernel.len(); - if vec == 0 || ker == 0 { - panic!("Convolve's inputs must not be 0-sized. "); - } - - if ker > vec { - return convolve_same(kernel, vector); + if ker == 0 || ker > vec { + panic!("convolve_same expects `vector.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); } let result_len = vector.data.shape().0.max(kernel.data.shape().0); diff --git a/tests/linalg/convolution.rs b/tests/linalg/convolution.rs index 454c29c6..ddcfe9f6 100644 --- a/tests/linalg/convolution.rs +++ b/tests/linalg/convolution.rs @@ -1,7 +1,6 @@ -#[allow(unused_imports)] // remove after fixing unit test use na::linalg::{convolve_full,convolve_valid,convolve_same}; -#[allow(unused_imports)] use na::{Vector2,Vector3,Vector4,Vector5,DVector}; +use std::panic; // // Should mimic calculations in Python's scipy library @@ -12,70 +11,110 @@ use na::{Vector2,Vector3,Vector4,Vector5,DVector}; // array([ 1, 4, 7, 10]) #[test] fn convolve_same_check(){ - let vec_s = Vector4::new(1.0,2.0,3.0,4.0); - let ker_s = Vector2::new(1.0,2.0); - + // Static Tests let actual_s = Vector4::from_vec(vec![1.0,4.0,7.0,10.0]); - - let expected_s = convolve_same(vec_s,ker_s); - let expected_s_r = convolve_same(ker_s,vec_s); + let expected_s = convolve_same(Vector4::new(1.0,2.0,3.0,4.0), Vector2::new(1.0,2.0)); assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); - assert!(relative_eq!(actual_s, expected_s_r, epsilon = 1.0e-7)); - - let vec_d = DVector::from_vec(4,vec![1.0,2.0,3.0,4.0]); - let ker_d = DVector::from_vec(2,vec![1.0,2.0]); - let actual_d = DVector::from_vec(4,vec![1.0,4.0,7.0,10.0]); - - let expected_d = convolve_same(vec_d.clone(),ker_d.clone()); - let expected_d_r = convolve_same(ker_d,vec_d); + // Dynamic Tests + let actual_d = DVector::from_vec(vec![1.0,4.0,7.0,10.0]); + let expected_d = convolve_same(DVector::from_vec(vec![1.0,2.0,3.0,4.0]),DVector::from_vec(vec![1.0,2.0])); assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); - assert!(relative_eq!(actual_d, expected_d_r, epsilon = 1.0e-7)); + + // Panic Tests + // These really only apply to dynamic sized vectors + assert!( + panic::catch_unwind(|| { + convolve_same(DVector::from_vec(vec![1.0,2.0]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + convolve_same(DVector::::from_vec(vec![]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + convolve_same(DVector::from_vec(vec![1.0,2.0,3.0,4.0]),DVector::::from_vec(vec![])); + }).is_err() + ); } // >>> convolve([1,2,3,4],[1,2],"full") // array([ 1, 4, 7, 10, 8]) #[test] fn convolve_full_check(){ - let vec_s = Vector4::new(1.0,2.0,3.0,4.0); - let ker_s = Vector2::new(1.0,2.0); - + // Static Tests let actual_s = Vector5::new(1.0,4.0,7.0,10.0,8.0); - - let expected_s = convolve_full(vec_s,ker_s); - let expected_s_r = convolve_full(ker_s,vec_s); + let expected_s = convolve_full(Vector4::new(1.0,2.0,3.0,4.0), Vector2::new(1.0,2.0)); assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); - assert!(relative_eq!(actual_s, expected_s_r, epsilon = 1.0e-7)); - - let vec_d = DVector::from_vec(4,vec![1.0,2.0,3.0,4.0]); - let ker_d = DVector::from_vec(2,vec![1.0,2.0]); - let actual_d = DVector::from_vec(5,vec![1.0,4.0,7.0,10.0,8.0]); - - let expected_d = convolve_full(vec_d.clone(),ker_d.clone()); - let expected_d_r = convolve_full(ker_d,vec_d); + // Dynamic Tests + let actual_d = DVector::from_vec(vec![1.0,4.0,7.0,10.0,8.0]); + let expected_d = convolve_full(DVector::from_vec(vec![1.0,2.0,3.0,4.0]), DVector::from_vec(vec![1.0,2.0])); assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); - assert!(relative_eq!(actual_d, expected_d_r, epsilon = 1.0e-7)); + + // Panic Tests + // These really only apply to dynamic sized vectors + assert!( + panic::catch_unwind(|| { + convolve_full(DVector::from_vec(vec![1.0,2.0]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + convolve_full(DVector::::from_vec(vec![]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + convolve_full(DVector::from_vec(vec![1.0,2.0,3.0,4.0]),DVector::::from_vec(vec![])); + }).is_err() + ); } // >>> convolve([1,2,3,4],[1,2],"valid") // array([ 4, 7, 10]) -// #[test] -// fn convolve_valid_check(){ -// let vec = Vector4::new(1.0,2.0,3.0,4.0); -// let ker = Vector2::new(1.0,2.0); +#[test] +fn convolve_valid_check(){ + // Static Tests + let actual_s = Vector3::from_vec(vec![4.0,7.0,10.0]); + let expected_s = convolve_valid( Vector4::new(1.0,2.0,3.0,4.0), Vector2::new(1.0,2.0)); -// let actual = Vector3::from_vec(vec![4.0,7.0,10.0]); + assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); -// let expected1 = convolve_valid(vec, ker); -// let expected2 = convolve_valid(ker, vec); + // Dynamic Tests + let actual_d = DVector::from_vec(vec![4.0,7.0,10.0]); + let expected_d = convolve_valid(DVector::from_vec(vec![1.0,2.0,3.0,4.0]), DVector::from_vec(vec![1.0,2.0])); + assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); -// assert!(relative_eq!(actual, expected1, epsilon = 1.0e-7)); -// assert!(relative_eq!(actual, expected2, epsilon = 1.0e-7)); + // Panic Tests + // These really only apply to dynamic sized vectors + assert!( + panic::catch_unwind(|| { + convolve_valid(DVector::from_vec(vec![1.0,2.0]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); -// } \ No newline at end of file + assert!( + panic::catch_unwind(|| { + convolve_valid(DVector::::from_vec(vec![]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + convolve_valid(DVector::from_vec(vec![1.0,2.0,3.0,4.0]),DVector::::from_vec(vec![])); + }).is_err() + ); + +} \ No newline at end of file From aeff67ecbdb7cf51900c38901efe86a1da17052c Mon Sep 17 00:00:00 2001 From: Benjamin Hetherington Date: Mon, 25 Feb 2019 11:34:13 -0500 Subject: [PATCH 15/22] Correct typo in quaternion documentation (#559) Correct "looses" to "loses" in quaternion.rs documentation. --- src/geometry/quaternion.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index 49656b2c..6460a446 100644 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -856,7 +856,7 @@ impl UnitQuaternion { /// Compute the exponential of a quaternion. /// - /// Note that this function yields a `Quaternion` because it looses the unit property. + /// Note that this function yields a `Quaternion` because it loses the unit property. #[inline] pub fn exp(&self) -> Quaternion { self.as_ref().exp() @@ -864,7 +864,7 @@ impl UnitQuaternion { /// Compute the natural logarithm of a quaternion. /// - /// Note that this function yields a `Quaternion` because it looses the unit property. + /// Note that this function yields a `Quaternion` because it loses the unit property. /// The vector part of the return value corresponds to the axis-angle representation (divided /// by 2.0) of this unit quaternion. /// From 36feddb8c2277f3a0ffa83c262761755393f704a Mon Sep 17 00:00:00 2001 From: Nathan Date: Sat, 2 Mar 2019 15:00:40 -0600 Subject: [PATCH 16/22] Moving functions into impl for Vector --- src/linalg/convolution.rs | 227 +++++++++++++++++------------------- tests/linalg/convolution.rs | 31 +++-- 2 files changed, 122 insertions(+), 136 deletions(-) diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs index 49c8b154..b121b34a 100644 --- a/src/linalg/convolution.rs +++ b/src/linalg/convolution.rs @@ -5,138 +5,125 @@ use std::cmp; use storage::Storage; use {zero, Real, Vector, VectorN, U1}; -/// Returns the convolution of the target vector and a kernel -/// -/// # Arguments -/// -/// * `vector` - A Vector with size > 0 -/// * `kernel` - A Vector with size > 0 -/// -/// # Errors -/// Inputs must statisfy `vector.len() >= kernel.len() > 0`. -/// -pub fn convolve_full( - vector: Vector, - kernel: Vector, -) -> VectorN, U1>> -where - N: Real, - D1: DimAdd, - D2: DimAdd>, - DimSum: DimSub, - S1: Storage, - S2: Storage, - DefaultAllocator: Allocator, U1>>, -{ - let vec = vector.len(); - let ker = kernel.len(); +impl> Vector { + /// Returns the convolution of the target vector and a kernel + /// + /// # Arguments + /// + /// * `kernel` - A Vector with size > 0 + /// + /// # Errors + /// Inputs must statisfy `vector.len() >= kernel.len() > 0`. + /// + pub fn convolve_full( + &self, + kernel: Vector, + ) -> VectorN, U1>> + where + D1: DimAdd, + D2: DimAdd>, + DimSum: DimSub, + S2: Storage, + DefaultAllocator: Allocator, U1>>, + { + let vec = self.len(); + let ker = kernel.len(); - if ker == 0 || ker > vec { - panic!("convolve_full expects `vector.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); - } + if ker == 0 || ker > vec { + panic!("convolve_full expects `self.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); + } - let result_len = vector.data.shape().0.add(kernel.data.shape().0).sub(U1); - let mut conv = VectorN::zeros_generic(result_len, U1); + let result_len = self.data.shape().0.add(kernel.data.shape().0).sub(U1); + let mut conv = VectorN::zeros_generic(result_len, U1); - for i in 0..(vec + ker - 1) { - let u_i = if i > vec { i - ker } else { 0 }; - let u_f = cmp::min(i, vec - 1); + for i in 0..(vec + ker - 1) { + let u_i = if i > vec { i - ker } else { 0 }; + let u_f = cmp::min(i, vec - 1); - if u_i == u_f { - conv[i] += vector[u_i] * kernel[(i - u_i)]; - } else { - for u in u_i..(u_f + 1) { - if i - u < ker { - conv[i] += vector[u] * kernel[(i - u)]; + if u_i == u_f { + conv[i] += self[u_i] * kernel[(i - u_i)]; + } else { + for u in u_i..(u_f + 1) { + if i - u < ker { + conv[i] += self[u] * kernel[(i - u)]; + } } } } + conv } - conv -} + /// Returns the convolution of the target vector and a kernel + /// The output convolution consists only of those elements that do not rely on the zero-padding. + /// # Arguments + /// + /// * `kernel` - A Vector with size > 0 + /// + /// + /// # Errors + /// Inputs must statisfy `self.len() >= kernel.len() > 0`. + /// + pub fn convolve_valid(&self, kernel: Vector, + ) -> VectorN, D2>> + where + D1: DimAdd, + D2: Dim, + DimSum: DimSub, + S2: Storage, + DefaultAllocator: Allocator, D2>>, + { + let vec = self.len(); + let ker = kernel.len(); -/// Returns the convolution of the vector and a kernel -/// The output convolution consists only of those elements that do not rely on the zero-padding. -/// # Arguments -/// -/// * `vector` - A Vector with size > 0 -/// * `kernel` - A Vector with size > 0 -/// -/// -/// # Errors -/// Inputs must statisfy `vector.len() >= kernel.len() > 0`. -/// -pub fn convolve_valid( - vector: Vector, - kernel: Vector, -) -> VectorN, D2>> -where - N: Real, - D1: DimAdd, - D2: Dim, - DimSum: DimSub, - S1: Storage, - S2: Storage, - DefaultAllocator: Allocator, D2>>, -{ - let vec = vector.len(); - let ker = kernel.len(); - - if ker == 0 || ker > vec { - panic!("convolve_valid expects `vector.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); - } - - let result_len = vector.data.shape().0.add(U1).sub(kernel.data.shape().0); - let mut conv = VectorN::zeros_generic(result_len, U1); - - for i in 0..(vec - ker + 1) { - for j in 0..ker { - conv[i] += vector[i + j] * kernel[ker - j - 1]; + if ker == 0 || ker > vec { + panic!("convolve_valid expects `self.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); } - } - conv -} -/// Returns the convolution of the vector and a kernel -/// The output convolution is the same size as vector, centered with respect to the ‘full’ output. -/// # Arguments -/// -/// * `vector` - A Vector with size > 0 -/// * `kernel` - A Vector with size > 0 -/// -/// # Errors -/// Inputs must statisfy `vector.len() >= kernel.len() > 0`. -pub fn convolve_same( - vector: Vector, - kernel: Vector, -) -> VectorN> -where - N: Real, - D1: DimMax, - D2: DimMax>, - S1: Storage, - S2: Storage, - DefaultAllocator: Allocator>, -{ - let vec = vector.len(); - let ker = kernel.len(); + let result_len = self.data.shape().0.add(U1).sub(kernel.data.shape().0); + let mut conv = VectorN::zeros_generic(result_len, U1); - if ker == 0 || ker > vec { - panic!("convolve_same expects `vector.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); - } - - let result_len = vector.data.shape().0.max(kernel.data.shape().0); - let mut conv = VectorN::zeros_generic(result_len, U1); - - for i in 0..vec { - for j in 0..ker { - let val = if i + j < 1 || i + j >= vec + 1 { - zero::() - } else { - vector[i + j - 1] - }; - conv[i] += val * kernel[ker - j - 1]; + for i in 0..(vec - ker + 1) { + for j in 0..ker { + conv[i] += self[i + j] * kernel[ker - j - 1]; + } } + conv + } + + /// Returns the convolution of the targetvector and a kernel + /// The output convolution is the same size as vector, centered with respect to the ‘full’ output. + /// # Arguments + /// + /// * `kernel` - A Vector with size > 0 + /// + /// # Errors + /// Inputs must statisfy `self.len() >= kernel.len() > 0`. + pub fn convolve_same(&self, kernel: Vector) -> VectorN> + where + D1: DimMax, + D2: DimMax>, + S2: Storage, + DefaultAllocator: Allocator>, + { + let vec = self.len(); + let ker = kernel.len(); + + if ker == 0 || ker > vec { + panic!("convolve_same expects `self.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); + } + + let result_len = self.data.shape().0.max(kernel.data.shape().0); + let mut conv = VectorN::zeros_generic(result_len, U1); + + for i in 0..vec { + for j in 0..ker { + let val = if i + j < 1 || i + j >= vec + 1 { + zero::() + } else { + self[i + j - 1] + }; + conv[i] += val * kernel[ker - j - 1]; + } + } + conv } - conv } diff --git a/tests/linalg/convolution.rs b/tests/linalg/convolution.rs index ddcfe9f6..b0d57f72 100644 --- a/tests/linalg/convolution.rs +++ b/tests/linalg/convolution.rs @@ -1,4 +1,3 @@ -use na::linalg::{convolve_full,convolve_valid,convolve_same}; use na::{Vector2,Vector3,Vector4,Vector5,DVector}; use std::panic; @@ -13,13 +12,13 @@ use std::panic; fn convolve_same_check(){ // Static Tests let actual_s = Vector4::from_vec(vec![1.0,4.0,7.0,10.0]); - let expected_s = convolve_same(Vector4::new(1.0,2.0,3.0,4.0), Vector2::new(1.0,2.0)); + let expected_s = Vector4::new(1.0,2.0,3.0,4.0).convolve_same(Vector2::new(1.0,2.0)); assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); // Dynamic Tests let actual_d = DVector::from_vec(vec![1.0,4.0,7.0,10.0]); - let expected_d = convolve_same(DVector::from_vec(vec![1.0,2.0,3.0,4.0]),DVector::from_vec(vec![1.0,2.0])); + let expected_d = DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_same(DVector::from_vec(vec![1.0,2.0])); assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); @@ -27,19 +26,19 @@ fn convolve_same_check(){ // These really only apply to dynamic sized vectors assert!( panic::catch_unwind(|| { - convolve_same(DVector::from_vec(vec![1.0,2.0]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + DVector::from_vec(vec![1.0,2.0]).convolve_same(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); }).is_err() ); assert!( panic::catch_unwind(|| { - convolve_same(DVector::::from_vec(vec![]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + DVector::::from_vec(vec![]).convolve_same(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); }).is_err() ); assert!( panic::catch_unwind(|| { - convolve_same(DVector::from_vec(vec![1.0,2.0,3.0,4.0]),DVector::::from_vec(vec![])); + DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_same(DVector::::from_vec(vec![])); }).is_err() ); } @@ -50,13 +49,13 @@ fn convolve_same_check(){ fn convolve_full_check(){ // Static Tests let actual_s = Vector5::new(1.0,4.0,7.0,10.0,8.0); - let expected_s = convolve_full(Vector4::new(1.0,2.0,3.0,4.0), Vector2::new(1.0,2.0)); + let expected_s = Vector4::new(1.0,2.0,3.0,4.0).convolve_full(Vector2::new(1.0,2.0)); assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); // Dynamic Tests let actual_d = DVector::from_vec(vec![1.0,4.0,7.0,10.0,8.0]); - let expected_d = convolve_full(DVector::from_vec(vec![1.0,2.0,3.0,4.0]), DVector::from_vec(vec![1.0,2.0])); + let expected_d = DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_full(DVector::from_vec(vec![1.0,2.0])); assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); @@ -64,19 +63,19 @@ fn convolve_full_check(){ // These really only apply to dynamic sized vectors assert!( panic::catch_unwind(|| { - convolve_full(DVector::from_vec(vec![1.0,2.0]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + DVector::from_vec(vec![1.0,2.0]).convolve_full(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); }).is_err() ); assert!( panic::catch_unwind(|| { - convolve_full(DVector::::from_vec(vec![]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + DVector::::from_vec(vec![]).convolve_full(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); }).is_err() ); assert!( panic::catch_unwind(|| { - convolve_full(DVector::from_vec(vec![1.0,2.0,3.0,4.0]),DVector::::from_vec(vec![])); + DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_full(DVector::::from_vec(vec![])); }).is_err() ); } @@ -87,13 +86,13 @@ fn convolve_full_check(){ fn convolve_valid_check(){ // Static Tests let actual_s = Vector3::from_vec(vec![4.0,7.0,10.0]); - let expected_s = convolve_valid( Vector4::new(1.0,2.0,3.0,4.0), Vector2::new(1.0,2.0)); + let expected_s = Vector4::new(1.0,2.0,3.0,4.0).convolve_valid( Vector2::new(1.0,2.0)); assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); // Dynamic Tests let actual_d = DVector::from_vec(vec![4.0,7.0,10.0]); - let expected_d = convolve_valid(DVector::from_vec(vec![1.0,2.0,3.0,4.0]), DVector::from_vec(vec![1.0,2.0])); + let expected_d = DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_valid(DVector::from_vec(vec![1.0,2.0])); assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); @@ -101,19 +100,19 @@ fn convolve_valid_check(){ // These really only apply to dynamic sized vectors assert!( panic::catch_unwind(|| { - convolve_valid(DVector::from_vec(vec![1.0,2.0]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + DVector::from_vec(vec![1.0,2.0]).convolve_valid(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); }).is_err() ); assert!( panic::catch_unwind(|| { - convolve_valid(DVector::::from_vec(vec![]), DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + DVector::::from_vec(vec![]).convolve_valid(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); }).is_err() ); assert!( panic::catch_unwind(|| { - convolve_valid(DVector::from_vec(vec![1.0,2.0,3.0,4.0]),DVector::::from_vec(vec![])); + DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_valid(DVector::::from_vec(vec![])); }).is_err() ); From edb08cd900fd50a2234bc81b1b6b73b6e10cf821 Mon Sep 17 00:00:00 2001 From: Adam Nemecek Date: Tue, 26 Feb 2019 18:12:30 -0800 Subject: [PATCH 17/22] quaternion trigonometry --- Cargo.toml | 3 + src/geometry/quaternion.rs | 252 +++++++++++++++++++++++- src/geometry/quaternion_construction.rs | 21 +- 3 files changed, 266 insertions(+), 10 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 43323c94..0c53df56 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -47,6 +47,9 @@ quickcheck = { version = "0.8", optional = true } pest = { version = "2.0", optional = true } pest_derive = { version = "2.0", optional = true } +[patch.crates-io] +alga = { git = "https://github.com/rustsim/alga", branch = "dev" } + [dev-dependencies] serde_json = "1.0" rand_xorshift = "0.1" diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index 6460a446..a2cc6130 100644 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -506,6 +506,255 @@ impl Quaternion { pub fn normalize_mut(&mut self) -> N { self.coords.normalize_mut() } + + /// Calculates square of a quaternion. + #[inline] + pub fn squared(&self) -> Self { + self * self + } + + /// Divides quaternion into two. + #[inline] + pub fn half(&self) -> Self { + self / ::convert(2.0f64) + } + + /// Calculates square root. + #[inline] + pub fn sqrt(&self) -> Self { + self.powf(::convert(0.5)) + } + + /// Check if the quaternion is pure. + #[inline] + pub fn is_pure(&self) -> bool { + self.w == N::zero() + } + + /// Convert quaternion to pure quaternion. + #[inline] + pub fn pure(&self) -> Self { + Self::from_imag(self.imag()) + } + + /// Calculates the quaternionic cosinus. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let expected = Quaternion::new(58.93364616794395, -34.086183690465596, -51.1292755356984, -68.17236738093119); + /// let result = input.cos(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn cos(&self) -> Self { + let z = self.imag().magnitude(); + let w = -self.w.sin() * z.sinhc(); + Self::from_parts(self.w.cos() * z.cosh(), self.imag() * w) + } + + /// Calculates the quaternionic arccosinus. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let result = input.cos().acos(); + /// assert_relative_eq!(input, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn acos(&self) -> Self { + let u = Self::from_imag(self.imag().normalize()); + let identity = Self::identity(); + + let z = (self + (self.squared() - identity).sqrt()).ln(); + + -(u * z) + } + + /// Calculates the quaternionic sinus. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let expected = Quaternion::new(91.78371578403467, 21.886486853029176, 32.82973027954377, 43.77297370605835); + /// let result = input.sin(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn sin(&self) -> Self { + let z = self.imag().magnitude(); + let w = self.w.cos() * z.sinhc(); + Self::from_parts(self.w.sin() * z.cosh(), self.imag() * w) + } + + /// Calculates the quaternionic arcsinus. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let result = input.sin().asin(); + /// assert_relative_eq!(input, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn asin(&self) -> Self { + let u = Self::from_imag(self.imag().normalize()); + let identity = Self::identity(); + + let z = ((u * self) + (identity - self.squared()).sqrt()).ln(); + + -(u * z) + } + + /// Calculates the quaternionic tangent. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let expected = Quaternion::new(0.00003821631725009489, 0.3713971716439371, 0.5570957574659058, 0.7427943432878743); + /// let result = input.tan(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn tan(&self) -> Self { + let s = self.sin(); + let c = self.cos(); + + let ci = c.try_inverse().unwrap(); + s * ci + } + + /// Calculates the quaternionic arctangent. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let result = input.tan().atan(); + /// assert_relative_eq!(input, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn atan(&self) -> Self { + let u = Self::from_imag(self.imag().normalize()); + let num = u + self; + let den = u - self; + let fr = num * den.try_inverse().unwrap(); + let ln = fr.ln(); + (u.half()) * ln + } + + /// Calculates the hyperbolic quaternionic sinus. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let expected = Quaternion::new(0.7323376060463428, -0.4482074499805421, -0.6723111749708133, -0.8964148999610843); + /// let result = input.sinh(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn sinh(&self) -> Self { + (self.exp() - (-self).exp()).half() + } + + /// Calculates the hyperbolic quaternionic arcsinus. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let expected = Quaternion::new(2.385889902585242, 0.514052600662788, 0.7710789009941821, 1.028105201325576); + /// let result = input.asinh(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn asinh(&self) -> Self { + let identity = Self::identity(); + (self + (identity + self.squared()).sqrt()).ln() + } + + /// Calculates the hyperbolic quaternionic cosinus. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let expected = Quaternion::new(0.9615851176369566, -0.3413521745610167, -0.5120282618415251, -0.6827043491220334); + /// let result = input.cosh(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn cosh(&self) -> Self { + (self.exp() + (-self).exp()).half() + } + + /// Calculates the hyperbolic quaternionic arccosinus. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let expected = Quaternion::new(2.4014472020074007, 0.5162761016176176, 0.7744141524264264, 1.0325522032352352); + /// let result = input.acosh(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn acosh(&self) -> Self { + let identity = Self::identity(); + (self + (self + identity).sqrt() * (self - identity).sqrt()).ln() + } + + /// Calculates the hyperbolic quaternionic tangent. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let expected = Quaternion::new(1.0248695360556623, -0.10229568178876419, -0.1534435226831464, -0.20459136357752844); + /// let result = input.tanh(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn tanh(&self) -> Self { + let s = self.sinh(); + let c = self.cosh(); + + let ci = c.try_inverse().unwrap(); + s * ci + } + + /// Calculates the hyperbolic quaternionic arctangent. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let input = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// let expected = Quaternion::new(0.03230293287000163, 0.5173453683196951, 0.7760180524795426, 1.0346907366393903); + /// let result = input.atanh(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn atanh(&self) -> Self { + let identity = Self::identity(); + ((identity + self).ln() - (identity - self).ln()).half() + } } impl> AbsDiffEq for Quaternion { @@ -879,7 +1128,7 @@ impl UnitQuaternion { #[inline] pub fn ln(&self) -> Quaternion { if let Some(v) = self.axis() { - Quaternion::from_parts(N::zero(), v.into_inner() * self.angle()) + Quaternion::from_imag(v.into_inner() * self.angle()) } else { Quaternion::zero() } @@ -1073,3 +1322,4 @@ impl> UlpsEq for UnitQuaternion { 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 0a30fabe..8eb16b1b 100644 --- a/src/geometry/quaternion_construction.rs +++ b/src/geometry/quaternion_construction.rs @@ -13,9 +13,7 @@ use alga::general::Real; use base::dimension::U3; use base::storage::Storage; -#[cfg(feature = "arbitrary")] -use base::Vector3; -use base::{Unit, Vector, Vector4, Matrix3}; +use base::{Unit, Vector, Vector3, Vector4, Matrix3}; use geometry::{Quaternion, Rotation3, UnitQuaternion}; @@ -43,8 +41,13 @@ impl Quaternion { /// ``` #[inline] pub fn new(w: N, i: N, j: N, k: N) -> Self { - let v = Vector4::::new(i, j, k, w); - Self::from(v) + Self::from(Vector4::new(i, j, k, w)) + } + + /// Constructs a pure quaternion. + #[inline] + pub fn from_imag(vector: Vector3) -> Self { + Self::from_parts(N::zero(), vector) } /// Creates a new quaternion from its scalar and vector parts. Note that the arguments order does @@ -92,7 +95,7 @@ impl Quaternion { /// ``` #[inline] pub fn identity() -> Self { - Self::new(N::one(), N::zero(), N::zero(), N::zero()) + Self::from_parts(N::one(), Vector3::zero()) } } @@ -106,7 +109,7 @@ impl One for Quaternion { impl Zero for Quaternion { #[inline] fn zero() -> Self { - Self::new(N::zero(), N::zero(), N::zero(), N::zero()) + Self::from(Vector4::zero()) } #[inline] @@ -579,7 +582,7 @@ impl UnitQuaternion { pub fn new(axisangle: Vector) -> Self where SB: Storage { let two: N = ::convert(2.0f64); - let q = Quaternion::::from_parts(N::zero(), axisangle / two).exp(); + let q = Quaternion::::from_imag(axisangle / two).exp(); Self::new_unchecked(q) } @@ -608,7 +611,7 @@ impl UnitQuaternion { pub fn new_eps(axisangle: Vector, eps: N) -> Self where SB: Storage { let two: N = ::convert(2.0f64); - let q = Quaternion::::from_parts(N::zero(), axisangle / two).exp_eps(eps); + let q = Quaternion::::from_imag(axisangle / two).exp_eps(eps); Self::new_unchecked(q) } From a2c0a453d35827cb40e62eaa9ff49b78b3f98598 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Fri, 15 Mar 2019 10:50:47 -0400 Subject: [PATCH 18/22] Add operator explanation to docs Co-Authored-By: tpdickso --- src/geometry/isometry.rs | 4 ++++ src/geometry/quaternion.rs | 4 ++++ src/geometry/rotation.rs | 4 ++++ src/geometry/similarity.rs | 4 ++++ src/geometry/transform.rs | 4 ++++ src/geometry/translation.rs | 2 ++ src/geometry/unit_complex.rs | 4 ++++ 7 files changed, 26 insertions(+) diff --git a/src/geometry/isometry.rs b/src/geometry/isometry.rs index 04e77b7a..9519a347 100755 --- a/src/geometry/isometry.rs +++ b/src/geometry/isometry.rs @@ -257,6 +257,8 @@ where DefaultAllocator: Allocator /// Transform the given point by this isometry. /// + /// This is the same as the multiplication `self * pt`. + /// /// # Example /// /// ``` @@ -278,6 +280,8 @@ where DefaultAllocator: Allocator /// Transform the given vector by this isometry, ignoring the translation /// component of the isometry. /// + /// This is the same as the multiplication `self * v`. + /// /// # Example /// /// ``` diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index 365e55e2..c8617f1d 100755 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -1008,6 +1008,8 @@ impl UnitQuaternion { /// Rotate a point by this unit quaternion. /// + /// This is the same as the multiplication `self * pt`. + /// /// # Example /// /// ``` @@ -1026,6 +1028,8 @@ impl UnitQuaternion { /// Rotate a vector by this unit quaternion. /// + /// This is the same as the multiplication `self * v`. + /// /// # Example /// /// ``` diff --git a/src/geometry/rotation.rs b/src/geometry/rotation.rs index 961951f1..10348eb4 100755 --- a/src/geometry/rotation.rs +++ b/src/geometry/rotation.rs @@ -358,6 +358,8 @@ where DefaultAllocator: Allocator + Allocator { /// Rotate the given point. /// + /// This is the same as the multiplication `self * pt`. + /// /// # Example /// ``` /// # #[macro_use] extern crate approx; @@ -375,6 +377,8 @@ where DefaultAllocator: Allocator + Allocator /// Rotate the given vector. /// + /// This is the same as the multiplication `self * v`. + /// /// # Example /// ``` /// # #[macro_use] extern crate approx; diff --git a/src/geometry/similarity.rs b/src/geometry/similarity.rs index dcb0f20a..b890dc67 100755 --- a/src/geometry/similarity.rs +++ b/src/geometry/similarity.rs @@ -241,6 +241,8 @@ where /// Transform the given point by this similarity. /// + /// This is the same as the multiplication `self * pt`. + /// /// # Example /// ``` /// # #[macro_use] extern crate approx; @@ -260,6 +262,8 @@ where /// Transform the given vector by this similarity, ignoring the translational /// component. /// + /// This is the same as the multiplication `self * t`. + /// /// # Example /// ``` /// # #[macro_use] extern crate approx; diff --git a/src/geometry/transform.rs b/src/geometry/transform.rs index 56a65992..e0eb1bb9 100755 --- a/src/geometry/transform.rs +++ b/src/geometry/transform.rs @@ -464,6 +464,8 @@ where + Allocator, { /// Transform the given point by this transformation. + /// + /// This is the same as the multiplication `self * pt`. #[inline] pub fn transform_point(&self, pt: &Point) -> Point { self * pt @@ -471,6 +473,8 @@ where /// Transform the given vector by this transformation, ignoring the /// translational component of the transformation. + /// + /// This is the same as the multiplication `self * v`. #[inline] pub fn transform_vector(&self, v: &VectorN) -> VectorN { self * v diff --git a/src/geometry/translation.rs b/src/geometry/translation.rs index 89960961..a49b2706 100755 --- a/src/geometry/translation.rs +++ b/src/geometry/translation.rs @@ -197,6 +197,8 @@ where DefaultAllocator: Allocator { /// Translate the given point. /// + /// This is the same as the multiplication `self * pt`. + /// /// # Example /// ``` /// # use nalgebra::{Translation3, Point3}; diff --git a/src/geometry/unit_complex.rs b/src/geometry/unit_complex.rs index bd064ef3..6e8cf81d 100755 --- a/src/geometry/unit_complex.rs +++ b/src/geometry/unit_complex.rs @@ -254,6 +254,8 @@ impl UnitComplex { /// Rotate the given point by this unit complex number. /// + /// This is the same as the multiplication `self * pt`. + /// /// # Example /// ``` /// # #[macro_use] extern crate approx; @@ -270,6 +272,8 @@ impl UnitComplex { /// Rotate the given vector by this unit complex number. /// + /// This is the same as the multiplication `self * v`. + /// /// # Example /// ``` /// # #[macro_use] extern crate approx; From 0f09f2a58cd8bee572ff1b4f6ffb14f30dc7221e Mon Sep 17 00:00:00 2001 From: Greizgh Date: Tue, 5 Feb 2019 19:22:17 +0100 Subject: [PATCH 19/22] Fix typo in axpy documentation --- src/base/blas.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/base/blas.rs b/src/base/blas.rs index c16496d4..be71a107 100644 --- a/src/base/blas.rs +++ b/src/base/blas.rs @@ -378,7 +378,7 @@ where { /// Computes `self = a * x + b * self`. /// - /// If be is zero, `self` is never read from. + /// If `b` is zero, `self` is never read from. /// /// # Examples: /// From 1e614db22744b77f982728c11685b85797727173 Mon Sep 17 00:00:00 2001 From: adamnemecek Date: Mon, 18 Mar 2019 01:08:42 -0700 Subject: [PATCH 20/22] Quaternionic division + refactoring (#563) --- src/base/unit.rs | 2 +- src/geometry/quaternion.rs | 49 ++++++++++++++++--------- src/geometry/quaternion_construction.rs | 8 +++- src/geometry/quaternion_ops.rs | 4 +- src/linalg/inverse.rs | 6 +-- 5 files changed, 45 insertions(+), 24 deletions(-) diff --git a/src/base/unit.rs b/src/base/unit.rs index 3d0a9c03..20f15956 100644 --- a/src/base/unit.rs +++ b/src/base/unit.rs @@ -225,7 +225,7 @@ impl Neg for Unit { #[inline] fn neg(self) -> Self::Output { - Unit::new_unchecked(-self.value) + Self::Output::new_unchecked(-self.value) } } diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index 8a56412b..cb7d252c 100755 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -528,7 +528,7 @@ impl Quaternion { /// Check if the quaternion is pure. #[inline] pub fn is_pure(&self) -> bool { - self.w == N::zero() + self.w.is_zero() } /// Convert quaternion to pure quaternion. @@ -537,6 +537,33 @@ impl Quaternion { Self::from_imag(self.imag()) } + /// Left quaternionic division. + /// + /// Calculates B-1 * A where A = self, B = other. + #[inline] + pub fn left_div(&self, other: &Self) -> Option { + other.try_inverse().map(|inv| inv * self) + } + + /// Right quaternionic division. + /// + /// Calculates A * B-1 where A = self, B = other. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let a = Quaternion::new(0.0, 1.0, 2.0, 3.0); + /// let b = Quaternion::new(0.0, 5.0, 2.0, 1.0); + /// let result = a.right_div(&b).unwrap(); + /// let expected = Quaternion::new(0.4, 0.13333333333333336, -0.4666666666666667, 0.26666666666666666); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-7); + /// ``` + #[inline] + pub fn right_div(&self, other: &Self) -> Option { + other.try_inverse().map(|inv| self * inv) + } + /// Calculates the quaternionic cosinus. /// /// # Example @@ -626,11 +653,7 @@ impl Quaternion { /// ``` #[inline] pub fn tan(&self) -> Self { - let s = self.sin(); - let c = self.cos(); - - let ci = c.try_inverse().unwrap(); - s * ci + self.sin().right_div(&self.cos()).unwrap() } /// Calculates the quaternionic arctangent. @@ -648,7 +671,7 @@ impl Quaternion { let u = Self::from_imag(self.imag().normalize()); let num = u + self; let den = u - self; - let fr = num * den.try_inverse().unwrap(); + let fr = num.right_div(&den).unwrap(); let ln = fr.ln(); (u.half()) * ln } @@ -732,11 +755,7 @@ impl Quaternion { /// ``` #[inline] pub fn tanh(&self) -> Self { - let s = self.sinh(); - let c = self.cosh(); - - let ci = c.try_inverse().unwrap(); - s * ci + self.sinh().right_div(&self.cosh()).unwrap() } /// Calculates the hyperbolic quaternionic arctangent. @@ -1096,11 +1115,7 @@ impl UnitQuaternion { /// ``` #[inline] pub fn axis_angle(&self) -> Option<(Unit>, N)> { - if let Some(axis) = self.axis() { - Some((axis, self.angle())) - } else { - None - } + self.axis().map(|axis| (axis, self.angle())) } /// Compute the exponential of a quaternion. diff --git a/src/geometry/quaternion_construction.rs b/src/geometry/quaternion_construction.rs index 8eb16b1b..4b7976bd 100644 --- a/src/geometry/quaternion_construction.rs +++ b/src/geometry/quaternion_construction.rs @@ -71,6 +71,12 @@ impl Quaternion { Self::new(scalar, vector[0], vector[1], vector[2]) } + /// Constructs a real quaternion. + #[inline] + pub fn from_real(r: N) -> Self { + Self::from_parts(r, Vector3::zero()) + } + /// Creates a new quaternion from its polar decomposition. /// /// Note that `axis` is assumed to be a unit vector. @@ -95,7 +101,7 @@ impl Quaternion { /// ``` #[inline] pub fn identity() -> Self { - Self::from_parts(N::one(), Vector3::zero()) + Self::from_real(N::one()) } } diff --git a/src/geometry/quaternion_ops.rs b/src/geometry/quaternion_ops.rs index 2ed72453..d43f104f 100644 --- a/src/geometry/quaternion_ops.rs +++ b/src/geometry/quaternion_ops.rs @@ -552,7 +552,7 @@ impl Neg for Quaternion { #[inline] fn neg(self) -> Self::Output { - Quaternion::from(-self.coords) + Self::Output::from(-self.coords) } } @@ -561,7 +561,7 @@ impl<'a, N: Real> Neg for &'a Quaternion { #[inline] fn neg(self) -> Self::Output { - Quaternion::from(-&self.coords) + Self::Output::from(-&self.coords) } } diff --git a/src/linalg/inverse.rs b/src/linalg/inverse.rs index 5748900f..69f1463d 100644 --- a/src/linalg/inverse.rs +++ b/src/linalg/inverse.rs @@ -36,7 +36,7 @@ impl> SquareMatrix { 0 => true, 1 => { let determinant = self.get_unchecked((0, 0)).clone(); - if determinant == N::zero() { + if determinant.is_zero() { false } else { *self.get_unchecked_mut((0, 0)) = N::one() / determinant; @@ -51,7 +51,7 @@ impl> SquareMatrix { let determinant = m11 * m22 - m21 * m12; - if determinant == N::zero() { + if determinant.is_zero() { false } else { *self.get_unchecked_mut((0, 0)) = m22 / determinant; @@ -83,7 +83,7 @@ impl> SquareMatrix { let determinant = m11 * minor_m12_m23 - m12 * minor_m11_m23 + m13 * minor_m11_m22; - if determinant == N::zero() { + if determinant.is_zero() { false } else { *self.get_unchecked_mut((0, 0)) = minor_m12_m23 / determinant; From e416360fc97d26c54229ca74824c92b3efe2882a Mon Sep 17 00:00:00 2001 From: Adam Nemecek Date: Thu, 21 Mar 2019 11:34:06 -0700 Subject: [PATCH 21/22] geometric operations --- src/geometry/quaternion.rs | 75 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index cb7d252c..ca6f7771 100755 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -307,6 +307,81 @@ impl Quaternion { self.coords.dot(&rhs.coords) } + /// Calculates the inner product (also known as the dot product). + /// See "Foundations of Game Engine Development, Volume 1: Mathematics" by Lengyel + /// Formula 4.89. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let a = Quaternion::new(0.0, 2.0, 3.0, 4.0); + /// let b = Quaternion::new(0.0, 5.0, 2.0, 1.0); + /// let expected = Quaternion::new(-20.0, 0.0, 0.0, 0.0); + /// let result = a.inner(&b); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-5); + #[inline] + pub fn inner(&self, other: &Self) -> Self { + (self * other + other * self).half() + } + + /// Calculates the outer product (also known as the wedge product). + /// See "Foundations of Game Engine Development, Volume 1: Mathematics" by Lengyel + /// Formula 4.89. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let a = Quaternion::new(0.0, 2.0, 3.0, 4.0); + /// let b = Quaternion::new(0.0, 5.0, 2.0, 1.0); + /// let expected = Quaternion::new(0.0, -5.0, 18.0, -11.0); + /// let result = a.outer(&b); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-5); + /// ``` + #[inline] + pub fn outer(&self, other: &Self) -> Self { + (self * other - other * self).half() + } + + /// Calculates the projection of `self` onto `other` (also known as the parallel). + /// See "Foundations of Game Engine Development, Volume 1: Mathematics" by Lengyel + /// Formula 4.94. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let a = Quaternion::new(0.0, 2.0, 3.0, 4.0); + /// let b = Quaternion::new(0.0, 5.0, 2.0, 1.0); + /// let expected = Quaternion::new(0.0, 3.333333333333333, 1.3333333333333333, 0.6666666666666666); + /// let result = a.project(&b).unwrap(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-5); + /// ``` + #[inline] + pub fn project(&self, other: &Self) -> Option { + self.inner(other).right_div(other) + } + + /// Calculates the rejection of `self` from `other` (also known as the perpendicular). + /// See "Foundations of Game Engine Development, Volume 1: Mathematics" by Lengyel + /// Formula 4.94. + /// + /// # Example + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::Quaternion; + /// let a = Quaternion::new(0.0, 2.0, 3.0, 4.0); + /// let b = Quaternion::new(0.0, 5.0, 2.0, 1.0); + /// let expected = Quaternion::new(0.0, -1.3333333333333333, 1.6666666666666665, 3.3333333333333335); + /// let result = a.reject(&b).unwrap(); + /// assert_relative_eq!(expected, result, epsilon = 1.0e-5); + /// ``` + #[inline] + pub fn reject(&self, other: &Self) -> Option { + self.outer(other).right_div(other) + } + /// The polar decomposition of this quaternion. /// /// Returns, from left to right: the quaternion norm, the half rotation angle, the rotation From 1e04053a21aa8cd8136c91a26ae88ee5241217b5 Mon Sep 17 00:00:00 2001 From: Adam Nemecek Date: Thu, 21 Mar 2019 15:51:09 -0700 Subject: [PATCH 22/22] refactoring --- src/geometry/quaternion.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index ca6f7771..d691f48c 100755 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -17,7 +17,7 @@ use alga::general::Real; use base::dimension::{U1, U3, U4}; use base::storage::{CStride, RStride}; -use base::{Matrix3, MatrixN, MatrixSlice, MatrixSliceMut, Unit, Vector3, Vector4}; +use base::{Matrix3, Matrix4, MatrixSlice, MatrixSliceMut, Unit, Vector3, Vector4}; use geometry::{Point3, Rotation}; @@ -1341,7 +1341,7 @@ impl UnitQuaternion { /// assert_relative_eq!(rot.to_homogeneous(), expected, epsilon = 1.0e-6); /// ``` #[inline] - pub fn to_homogeneous(&self) -> MatrixN { + pub fn to_homogeneous(&self) -> Matrix4 { self.to_rotation_matrix().to_homogeneous() }