+ /// Return a quaternion that rotates the from vector onto the to vector.
+ static SGQuat fromRotateTo(const SGVec3<T>& from, const SGVec3<T>& to)
+ {
+ T nfrom = norm(from);
+ T nto = norm(to);
+ if (nfrom <= SGLimits<T>::min() || nto <= SGLimits<T>::min())
+ return SGQuat::unit();
+
+ return SGQuat::fromRotateToNorm((1/nfrom)*from, (1/nto)*to);
+ }
+
+ /// Return a quaternion that rotates v1 onto the i1-th unit vector
+ /// and v2 into a plane that is spanned by the i2-th and i1-th unit vector.
+ static SGQuat fromRotateTo(const SGVec3<T>& v1, unsigned i1,
+ const SGVec3<T>& v2, unsigned i2)
+ {
+ T nrmv1 = norm(v1);
+ T nrmv2 = norm(v2);
+ if (nrmv1 <= SGLimits<T>::min() || nrmv2 <= SGLimits<T>::min())
+ return SGQuat::unit();
+
+ SGVec3<T> nv1 = (1/nrmv1)*v1;
+ SGVec3<T> nv2 = (1/nrmv2)*v2;
+ T dv1v2 = dot(nv1, nv2);
+ if (fabs(fabs(dv1v2)-1) <= SGLimits<T>::epsilon())
+ return SGQuat::unit();
+
+ // The target vector for the first rotation
+ SGVec3<T> nto1 = SGVec3<T>::zeros();
+ SGVec3<T> nto2 = SGVec3<T>::zeros();
+ nto1[i1] = 1;
+ nto2[i2] = 1;
+
+ // The first rotation can be done with the usual routine.
+ SGQuat q = SGQuat::fromRotateToNorm(nv1, nto1);
+
+ // The rotation axis for the second rotation is the
+ // target for the first one, so the rotation axis is nto1
+ // We need to get the angle.
+
+ // Make nv2 exactly orthogonal to nv1.
+ nv2 = normalize(nv2 - dv1v2*nv1);
+
+ SGVec3<T> tnv2 = q.transform(nv2);
+ T cosang = dot(nto2, tnv2);
+ T cos05ang = T(0.5)+T(0.5)*cosang;
+ if (cos05ang <= 0)
+ cosang = 0;
+ cos05ang = sqrt(cos05ang);
+ T sig = dot(nto1, cross(nto2, tnv2));
+ T sin05ang = T(0.5)-T(0.5)*cosang;
+ if (sin05ang <= 0)
+ sin05ang = 0;
+ sin05ang = copysign(sqrt(sin05ang), sig);
+ q *= SGQuat::fromRealImag(cos05ang, sin05ang*nto1);
+
+ return q;
+ }
+
+
+ // Return a quaternion which rotates the vector given by v
+ // to the vector -v. Other directions are *not* preserved.
+ static SGQuat fromChangeSign(const SGVec3<T>& v)
+ {
+ // The vector from points to the oposite direction than to.
+ // Find a vector perpendicular to the vector to.
+ T absv1 = fabs(v(0));
+ T absv2 = fabs(v(1));
+ T absv3 = fabs(v(2));
+
+ SGVec3<T> axis;
+ if (absv2 < absv1 && absv3 < absv1) {
+ T quot = v(1)/v(0);
+ axis = (1/sqrt(1+quot*quot))*SGVec3<T>(quot, -1, 0);
+ } else if (absv1 < absv2 && absv3 < absv2) {
+ T quot = v(2)/v(1);
+ axis = (1/sqrt(1+quot*quot))*SGVec3<T>(0, quot, -1);
+ } else if (absv1 < absv3 && absv2 < absv3) {
+ T quot = v(0)/v(2);
+ axis = (1/sqrt(1+quot*quot))*SGVec3<T>(-1, 0, quot);
+ } else {
+ // The all zero case.
+ return SGQuat::unit();
+ }
+
+ return SGQuat::fromRealImag(0, axis);
+ }
+