From a808045360fd7904e0303a9859cb4e06ea61e5df Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Thu, 29 Oct 2026 21:27:58 +0100 Subject: [PATCH] bugfixes in vecmat3 and quaternion --- quaternion.cc | 66 +++++++++++++++++++++++++-------------------------- vecmat3.cc | 1 + 2 files changed, 34 insertions(+), 33 deletions(-) diff --git a/quaternion.cc b/quaternion.cc index e6787aa..29e776f 100644 --- a/quaternion.cc +++ b/quaternion.cc @@ -223,90 +223,90 @@ switch(Euler_case(type[0],type[1],type[2])) case Euler_case('x','z','x'): { q[0] = c2*(c1*c3-s1*s3); - q[1] = c2*(s1*c3+c1*s3); - q[2] = -s2*(s1*c3-c1*s3); - q[3] = s2*(c1*c3+s1*s3); + q[1] = -c2*(s1*c3+c1*s3); + q[2] = s2*(s1*c3-c1*s3); + q[3] = -s2*(c1*c3+s1*s3); } break; case Euler_case('x','y','x'): { q[0] = c2*(c1*c3-s1*s3); - q[1] = c2*(s1*c3+c1*s3); - q[2] = s2*(c1*c3+s1*s3); - q[3] = s2*(s1*c3-c1*s3); + q[1] = -c2*(s1*c3+c1*s3); + q[2] = -s2*(c1*c3+s1*s3); + q[3] = -s2*(s1*c3-c1*s3); } break; case Euler_case('y','x','y'): { q[0] = c2*(c1*c3-s1*s3); - q[1] = s2*(c1*c3+s1*s3); - q[2] = c2*(s1*c3+c1*s3); - q[3] = -s2*(s1*c3-c1*s3); + q[1] = -s2*(c1*c3+s1*s3); + q[2] = -c2*(s1*c3+c1*s3); + q[3] = s2*(s1*c3-c1*s3); } break; case Euler_case('y','z','y'): { q[0] = c2*(c1*c3-s1*s3); - q[1] = s2*(s1*c3-c1*s3); - q[2] = c2*(s1*c3+c1*s3); - q[3] = s2*(c1*c3+s1*s3); + q[1] = -s2*(s1*c3-c1*s3); + q[2] = -c2*(s1*c3+c1*s3); + q[3] = -s2*(c1*c3+s1*s3); } break; case Euler_case('z','y','z'): { q[0] = c2*(c1*c3-s1*s3); - q[1] = -s2*(s1*c3-c1*s3); - q[2] = s2*(c1*c3+s1*s3); - q[3] = c2*(s1*c3+c1*s3); + q[1] = s2*(s1*c3-c1*s3); + q[2] = -s2*(c1*c3+s1*s3); + q[3] = -c2*(s1*c3+c1*s3); } break; case Euler_case('z','x','z'): { q[0] = c2*(c1*c3-s1*s3); - q[1] = s2*(c1*c3+s1*s3); - q[2] = s2*(s1*c3-c1*s3); - q[3] = c2*(s1*c3+c1*s3); + q[1] = -s2*(c1*c3+s1*s3); + q[2] = -s2*(s1*c3-c1*s3); + q[3] = -c2*(s1*c3+c1*s3); } break; case Euler_case('x','z','y'): { q[0] = s1*s2*s3 +c1*c2*c3; - q[1] = s1*c2*c3 - c1*s2*s3; - q[2] = -s1*s2*c3+c1*c2*s3; - q[3] = s1*c2*s3 +c1*s2*c3; + q[1] = -(s1*c2*c3 - c1*s2*s3); + q[2] = -(-s1*s2*c3+c1*c2*s3); + q[3] = -(s1*c2*s3 +c1*s2*c3); } break; case Euler_case('x','y','z'): { q[0] = -s1*s2*s3 + c1*c2*c3; - q[1] = s1*c2*c3 +c1*s2*s3; - q[2] = -s1*c2*s3 + c1*s2*c3; - q[3] = s1*s2*c3 + c1*c2*s3; + q[1] = -(s1*c2*c3 +c1*s2*s3); + q[2] = -(-s1*c2*s3 + c1*s2*c3); + q[3] = -(s1*s2*c3 + c1*c2*s3); } break; case Euler_case('y','x','z'): { q[0] = s1*s2*s3 + c1*c2*c3; - q[1] = s1*c2*s3 + c1*s2*c3; - q[2] = s1*c2*c3 - c1*s2*s3; - q[3] = -s1*s2*c3 + c1*c2*s3; + q[1] = -(s1*c2*s3 + c1*s2*c3); + q[2] = -(s1*c2*c3 - c1*s2*s3); + q[3] = -(-s1*s2*c3 + c1*c2*s3); } break; case Euler_case('y','z','x'): { q[0] = -s1*s2*s3 + c1*c2*c3; - q[1] = s1*s2*c3 + c1*c2*s3; - q[2] = s1*c2*c3 + c1*s2*s3; - q[3] = -s1*c2*s3 + c1*s2*c3; + q[1] = -(s1*s2*c3 + c1*c2*s3); + q[2] = -(s1*c2*c3 + c1*s2*s3); + q[3] = -(-s1*c2*s3 + c1*s2*c3); } break; @@ -322,9 +322,9 @@ case Euler_case('z','y','x'): case Euler_case('z','x','y'): { q[0] = -s1*s2*s3 + c1*c2*c3; - q[1] = -s1*c2*s3 + c1*s2*c3; - q[2] = s1*s2*c3 + c1*c2*s3; - q[3] = s1*c2*c3 + c1*s2*s3; + q[1] = -(-s1*c2*s3 + c1*s2*c3); + q[2] = -(s1*s2*c3 + c1*c2*s3); + q[3] = -(s1*c2*c3 + c1*s2*s3); } break; }//switch diff --git a/vecmat3.cc b/vecmat3.cc index 6a50982..f97b154 100644 --- a/vecmat3.cc +++ b/vecmat3.cc @@ -1205,6 +1205,7 @@ void Mat3::svd(Mat3 &u, Vec3 &w, Mat3 &v, bool proper_rotations) const Mat3 ata= this->Ttimes(*this); Vec3 ww; //actually not needed squares of singular values ata.eivec_sym(ww,v,true); +if(proper_rotations && v.determinant()<0) {v[0][2]= -v[0][2]; v[1][2]= -v[1][2]; v[2][2]= -v[2][2];} //det v = +1, the sign goes to w[2] Mat3 uw = (*this)*v; Mat3 r; uw.qrd(u,r);