fixed a bug in the quaternion class that caused an ambiguity for rotation angles…

… greater than 90-degrees

fixed a bug in the quaternion class that caused an ambiguity for rotation angles…
… greater than 90-degrees
David Mayerich
1 parent 8718a1a0
Showing 3 changed files with 55 additions and 13 deletions Show diff stats
stim/image/image.h
stim/math/quaternion.h
stim/ui/progressbar.h
@@ -248,7 +248,7 @@ public:
 			outfile.close();
 		}
-		cimg_library::CImg<T> cimg(R[1], R[2], 1, R[0]);
+		cimg_library::CImg<T> cimg((unsigned int)R[1], (unsigned int)R[2], 1, (unsigned int)R[0]);
 		get_noninterleaved(cimg.data());
 		cimg.save(filename.c_str());
 	}
@@ -24,6 +24,11 @@ public:
 		z=z/length;
 	}
+	//calculate the quaternion length (norm)
+	CUDA_CALLABLE T norm() {
+		return sqrt(w*w + x * x + y * y + z * z);
+	}
+
 	CUDA_CALLABLE void CreateRotation(T theta, T ux, T uy, T uz){
 		vec3<T> u(ux, uy, uz);
@@ -46,10 +51,12 @@ public:
 		from = from.norm();
 		to = to.norm();
 		vec3<T> r = from.cross(to);			//compute the rotation vector
-		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
+		//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)
+		T cos_theta = from.dot(to);			//calculate the cosine between the two vectors
+		T theta = acos(cos_theta);			//calculate the angle between the two vectors
 		if(theta == (T)0){
 			return;
 		}
@@ -117,12 +124,34 @@ public:
 	}
 	CUDA_CALLABLE matrix_sq<T, 4> toMatrix4(){
-
-		matrix_sq<T, 4> result;
-	    T wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
+		matrix_sq<T, 4> R;
+	    T s, wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
+		s = sqrt(norm());
+		xx = x * x;		xy = x * y;		xz = x * z;
+		yy = y * y;		yz = y * z;
+		zz = z * z;
+		wx = w * x;		wy = w * y;		wz = w * z;
+
+		R(0, 0) = 1 - 2 * s * (yy + zz);
+		R(0, 1) = 2 * s * (xy - wz);
+		R(0, 2) = 2 * s * (xz + wy);
+		R(1, 0) = 2 * s * (xy + wz);
+		R(1, 1) = 1 - 2 * s * (xx + zz);
+		R(1, 2) = 2 * s * (yz - wx);
+		R(2, 0) = 2 * s * (xz - wy);
+		R(2, 1) = 2 * s * (yz + wx);
+		R(2, 2) = 1 - 2 * s * (xx + yy);
+
+		R(0, 3) = 0;
+		R(1, 3) = 0;
+		R(2, 3) = 0;
+		R(3, 0) = 0;
+		R(3, 1) = 0;
+		R(3, 2) = 0;
+		R(3, 3) = 1;
 	    // calculate coefficients
-	    x2 = x + x; y2 = y + y;
+	    /*x2 = x + x; y2 = y + y;
 	    z2 = z + z;
 	    xx = x * x2; xy = x * y2; xz = x * z2;
 	    yy = y * y2; yz = y * z2; zz = z * z2;
@@ -143,9 +172,9 @@ public:
 		result(2, 2) = 1 - (xx + yy);
-		result(3, 3) = 1;
+		result(3, 3) = 1;*/
-		return result;
+		return R;
 	}
@@ -157,6 +186,20 @@ public:
 		w=c;  x=i;  y=j;  z=k;
 	}
+	// create a pure quaternion from a vector
+	CUDA_CALLABLE quaternion(vec3<T> v){
+		w = 0; x = v[0]; y = v[1]; z = v[2];
+	}
+
+	CUDA_CALLABLE quaternion<T> conj(){
+		quaternion<T> c;
+		c.w = w;
+		c.x = -x;
+		c.y = -y;
+		c.z = -z;
+		return c;
+	}
+
 };
 }	//end rts namespace
@@ -7,7 +7,6 @@
 using namespace std;
 void stimPrintTime(unsigned long long s) {
-	char buffer[26];
 	std::cout << std::fixed << std::showpoint << std::setprecision(1);
 	if (s >= 60) {											//if there are more than 60 seconds
 		unsigned long long m = s / 60;						//get the number of minutes
@@ -49,12 +48,12 @@ static void rtsProgressBar(unsigned int percent, double elapsed_s = 0)
 	cout << "\r"; // carriage return back to beginning of line
 	cout << bar.str() << " " << slash[x] << " " << percent << " %"; // print the bars and percentage
 	if (elapsed_s > 0) {
-		stimPrintTime(elapsed_s);
+		stimPrintTime((unsigned long long)elapsed_s);
 		if (percent > 0) {
 			std::cout << " / ";
 			double s_per_percent = elapsed_s / percent;
 			double total = s_per_percent * 100;
-			stimPrintTime(total);
+			stimPrintTime((unsigned long long)total);
 		}
 	}
 	cout.flush();