added considerations for 180deg rotations in Rotation::powf

This commit is contained in:
Joshua Smith 2022-03-28 16:07:11 -05:00
parent 04d6f4f39c
commit 6eab8f5175
2 changed files with 84 additions and 21 deletions

View File

@ -1065,19 +1065,21 @@ impl<T:RealField, const D: usize> Rotation<T,D>
// println!("q:{}d:{:.3}", q, d); // println!("q:{}d:{:.3}", q, d);
//go down the diagonal and pow every block //go down the diagonal and pow every block
for i in 0..(D-1) { let mut i = 0;
while i < D-1 {
//we've found a 2x2 block! if
//NOTE: the impl of the schur decomposition always sets the inferior diagonal to 0 //For most 2x2 blocks
if !d[(i+1,i)].is_zero() { //NOTE: we use strict equality since `nalgebra`'s schur decomp sets the infradiagonal to zero
!d[(i+1,i)].is_zero() ||
// println!("{}", i); //for +-180 deg rotations
d[(i,i)]<T::zero() && d[(i+1,i+1)]<T::zero()
{
//convert to a complex num and take the arg() //convert to a complex num and find the arg()
let (c, s) = (d[(i,i)].clone(), d[(i+1,i)].clone()); let (c, s) = (d[(i,i)].clone(), d[(i+1,i)].clone());
let angle = s.atan2(c); let angle = s.atan2(c); //for +-180deg rots, this implicitely takes the +180 branch
// println!("{}", angle);
//scale the arg and exponentiate back //scale the arg and exponentiate back
let angle2 = angle * t.clone(); let angle2 = angle * t.clone();
@ -1089,6 +1091,12 @@ impl<T:RealField, const D: usize> Rotation<T,D>
d[(i+1,i )] = s2; d[(i+1,i )] = s2;
d[(i+1,i+1)] = c2; d[(i+1,i+1)] = c2;
//increase by 2 so we don't accidentally misinterpret the
//next line as a 180deg rotation
i += 2;
} else {
i += 1;
} }
} }

View File

@ -39,30 +39,41 @@ mod proptest_tests {
use crate::proptest::*; use crate::proptest::*;
use proptest::{prop_assert, prop_assert_eq, proptest}; use proptest::{prop_assert, prop_assert_eq, proptest};
//creates N rotation planes and angles
macro_rules! gen_rotation_planes {
($($v1:ident, $v2:ident),*) => {
{
//make an orthonormal basis
let mut basis = [$($v1, $v2),*];
Vector::orthonormalize(&mut basis);
let [$($v1, $v2),*] = basis;
//"wedge" the vectors to make an arrary 2-blades representing rotation planes.
[
//Since we start with vector pairs, each bivector is guaranteed to be simple
$($v1.transpose().kronecker(&$v2) - $v2.transpose().kronecker(&$v1)),*
]
}
};
}
macro_rules! gen_powf_rotation_test { macro_rules! gen_powf_rotation_test {
($( ($(
fn $powf_rot_n:ident($($v1:ident in $vec1:ident(), $v2:ident in $vec2:ident()),*); fn $powf_rot_n:ident($($v:ident in $vec:ident()),*);
)*) => { )*) => {
proptest!{$( proptest!{$(
#[test] #[test]
fn $powf_rot_n( fn $powf_rot_n(
$($v1 in $vec1(), $v2 in $vec2(),)* $($v in $vec(),)*
pow in PROPTEST_F64 pow in PROPTEST_F64
) { ) {
use nalgebra::*; use nalgebra::*;
//make an orthonormal basis
let mut basis = [$($v1, $v2),*];
Vector::orthonormalize(&mut basis);
let [$($v1, $v2),*] = basis;
//"wedge" the vectors to make an arrary 2-blades representing rotation planes. //"wedge" the vectors to make an arrary 2-blades representing rotation planes.
let mut bivectors = [ let mut bivectors = gen_rotation_planes!($($v),*);
//Since we start with vector pairs, each bivector is guaranteed to be simple
$($v1.transpose().kronecker(&$v2) - $v2.transpose().kronecker(&$v1)),*
];
//condition the bivectors //condition the bivectors
for b in &mut bivectors { for b in &mut bivectors {
@ -88,7 +99,7 @@ mod proptest_tests {
} }
)*} )*}
} };
} }
gen_powf_rotation_test!( gen_powf_rotation_test!(
@ -101,6 +112,50 @@ mod proptest_tests {
); );
); );
proptest! {
#[test]
fn powf_180deg_rotation_4d(v1 in vector4(), v2 in vector4(), v3 in vector4(), v4 in vector4()) {
use nalgebra::*;
use std::f64::consts::PI;
let [b1,b2] = gen_rotation_planes!(v1,v2,v3,v4);
if let (Some((b1,a1)), Some((b2,a2))) = (
Unit::try_new_and_get(b1,0.0), Unit::try_new_and_get(b2,0.0)
) {
let (b1, b2) = (b1.into_inner(), b2.into_inner());
{
let b = a1*b1 + PI*b2;
let r1 = Rotation::from_matrix_unchecked(b.exp());
let r2 = Rotation::from_matrix_unchecked((b/2.0).exp());
prop_assert!(relative_eq!(r1.powf(0.5), r2, epsilon=1e-7));
}
{
let b = PI*b1 + a2*b2;
let r1 = Rotation::from_matrix_unchecked(b.exp());
let r2 = Rotation::from_matrix_unchecked((b/2.0).exp());
prop_assert!(relative_eq!(r1.powf(0.5), r2, epsilon=1e-7));
}
{
let b = PI*b1 + PI*b2;
let r1 = Rotation::from_matrix_unchecked(b.exp());
let r2 = Rotation::from_matrix_unchecked((b/2.0).exp());
prop_assert!(relative_eq!(r1.powf(0.5), r2, epsilon=1e-7));
}
}
}
}
proptest! { proptest! {
/* /*
* *