NaN value for |AxB| > 1 in some cases

David Mayerich
1 parent 90f83371
Showing 3 changed files with 10 additions and 5 deletions Show diff stats
stim/biomodels/centerline.h
stim/math/plane.h
stim/math/quaternion.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<T> a = at(i) - at(i - 1);
+		vec3<T> b = at(i + 1) - at(i);
+		vec3<T> 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) {
@@ -188,9 +188,9 @@ class plane
 		{
 			quaternion<T> q;
 			q.CreateRotation(N, n);
-			
-			N = q.toMatrix3() * N;
-			U = q.toMatrix3() * U;
+			matrix_sq<T, 3> M = q.toMatrix3();
+			N = M * N;
+			U = M * U;
  
 		}
  
@@ -46,7 +46,9 @@ public:
 		from = from.norm();
 		to = to.norm();
 		vec3<T> 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;