From 83c3121c14e0f5d78a87ba843a9c5843adfff24b Mon Sep 17 00:00:00 2001 From: David Date: Mon, 9 Jan 2017 16:43:34 -0600 Subject: [PATCH] NaN value for |AxB| > 1 in some cases --- stim/biomodels/centerline.h | 5 ++++- stim/math/plane.h | 6 +++--- stim/math/quaternion.h | 4 +++- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/stim/biomodels/centerline.h b/stim/biomodels/centerline.h index cfefee0..7a4ca45 100644 --- a/stim/biomodels/centerline.h +++ b/stim/biomodels/centerline.h @@ -25,7 +25,10 @@ protected: if (i == size() - 1) return at(size() - 1) - at(size() - 2); //the last direction vector is oriented towards the last line segment //all other direction vectors are the average direction of the two joined line segments - return ((at(i) - at(i - 1)).norm() + (at(i + 1) - at(i)).norm()).norm(); + vec3 a = at(i) - at(i - 1); + vec3 b = at(i + 1) - at(i); + vec3 ab = a.norm() + b.norm(); + return ab.norm(); } //initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct) void update_L(size_t start = 0) { diff --git a/stim/math/plane.h b/stim/math/plane.h index c85a1cc..6dd54f6 100644 --- a/stim/math/plane.h +++ b/stim/math/plane.h @@ -188,9 +188,9 @@ class plane { quaternion q; q.CreateRotation(N, n); - - N = q.toMatrix3() * N; - U = q.toMatrix3() * U; + matrix_sq M = q.toMatrix3(); + N = M * N; + U = M * U; } diff --git a/stim/math/quaternion.h b/stim/math/quaternion.h index 4170873..80a9115 100644 --- a/stim/math/quaternion.h +++ b/stim/math/quaternion.h @@ -46,7 +46,9 @@ public: from = from.norm(); to = to.norm(); vec3 r = from.cross(to); //compute the rotation vector - T theta = asin(r.len()); //compute the angle of the rotation about r + T l = r.len(); + if (l > 1) l = 1; //we have seen degenerate cases where |r| > 1 (probably due to loss of precision in the cross product) + T theta = asin(l); //compute the angle of the rotation about r //deal with a zero vector (both k and kn point in the same direction) if(theta == (T)0){ return; -- libgit2 0.21.4