imathquat.h

对gif

    // Given a set of quaternion keys: q0, q1, q2, q3,
    // this routine does the interpolation between
    // q1 and q2 by constructing two intermediate
    // quaternions: qa and qb. The qa and qb are 
    // computed by the intermediate function to 
    // guarantee the continuity of tangents across
    // adjacent cubic segments. The qa represents in-tangent
    // for q1 and the qb represents the out-tangent for q2.
    // 
    // The q1 q2 is the cubic segment being interpolated. 
    // The q0 is from the previous adjacent segment and q3 is 
    // from the next adjacent segment. The q0 and q3 are used
    // in computing qa and qb.
    // 

    Quat<T> qa = intermediate (q0, q1, q2);
    Quat<T> qb = intermediate (q1, q2, q3);
    Quat<T> result = squad(q1, qa, qb, q2, t);

    return result;
}

template<class T>
Quat<T> squad(const Quat<T> &q1, const Quat<T> &qa,
	      const Quat<T> &qb, const Quat<T> &q2,
	      T t)
{
    // Spherical Quadrangle Interpolation -
    // from Advanced Animation and Rendering
    // Techniques by Watt and Watt, Page 366:
    // It constructs a spherical cubic interpolation as 
    // a series of three spherical linear interpolations 
    // of a quadrangle of unit quaternions. 
    //     
  
    Quat<T> r1 = slerp(q1, q2, t);
    Quat<T> r2 = slerp(qa, qb, t);
    Quat<T> result = slerp(r1, r2, 2*t*(1-t));

    return result;
}

template<class T>
Quat<T> intermediate(const Quat<T> &q0, const Quat<T> &q1, const Quat<T> &q2)
{
    // From advanced Animation and Rendering
    // Techniques by Watt and Watt, Page 366:
    // computing the inner quadrangle 
    // points (qa and qb) to guarantee tangent
    // continuity.
    // 
    Quat<T> q1inv = q1.inverse();
    Quat<T> c1 = q1inv*q2;
    Quat<T> c2 = q1inv*q0;
    Quat<T> c3 = (T) (-0.25) * (c2.log() + c1.log());
    Quat<T> qa = q1 * c3.exp();
    qa.normalize();
    return qa;
}

template <class T>
inline Quat<T> Quat<T>::log() const
{
    //
    // For unit quaternion, from Advanced Animation and 
    // Rendering Techniques by Watt and Watt, Page 366:
    //

    T theta = Math<T>::acos (std::min (r, (T) 1.0));

    if (theta == 0)
	return Quat<T> (0, v);
    
    T sintheta = Math<T>::sin (theta);
    
    T k;
    if (abs (sintheta) < 1 && abs (theta) >= limits<T>::max() * abs (sintheta))
	k = 1;
    else
	k = theta / sintheta;

    return Quat<T> ((T) 0, v.x * k, v.y * k, v.z * k);
}

template <class T>
inline Quat<T> Quat<T>::exp() const
{
    //
    // For pure quaternion (zero scalar part):
    // from Advanced Animation and Rendering
    // Techniques by Watt and Watt, Page 366:
    //

    T theta = v.length();
    T sintheta = Math<T>::sin (theta);
    
    T k;
    if (abs (theta) < 1 && abs (sintheta) >= limits<T>::max() * abs (theta))
	k = 1;
    else
	k = sintheta / theta;

    T costheta = Math<T>::cos (theta);

    return Quat<T> (costheta, v.x * k, v.y * k, v.z * k);
}

template <class T>
inline T Quat<T>::angle() const
{
    return 2.0*Math<T>::acos(r);
}

template <class T>
inline Vec3<T> Quat<T>::axis() const
{
    return v.normalized();
}

template <class T>
inline Quat<T>& Quat<T>::setAxisAngle(const Vec3<T>& axis, T radians)
{
    r = Math<T>::cos(radians/2);
    v = axis.normalized() * Math<T>::sin(radians/2);
    return *this;
}


template <class T>
Quat<T>&
Quat<T>::setRotation(const Vec3<T>& from, const Vec3<T>& to)
{
    //
    // Create a quaternion that rotates vector from into vector to,
    // such that the rotation is around an axis that is the cross
    // product of from and to.
    //
    // This function calls function setRotationInternal(), which is
    // numerically accurate only for rotation angles that are not much
    // greater than pi/2.  In order to achieve good accuracy for angles
    // greater than pi/2, we split large angles in half, and rotate in
    // two steps.
    //

    //
    // Normalize from and to, yielding f0 and t0.
    //

    Vec3<T> f0 = from.normalized();
    Vec3<T> t0 = to.normalized();

    if ((f0 ^ t0) >= 0)
    {
	//
	// The rotation angle is less than or equal to pi/2.
	//

	setRotationInternal (f0, t0, *this);
    }
    else
    {
	//
	// The angle is greater than pi/2.  After computing h0,
	// which is halfway between f0 and t0, we rotate first
	// from f0 to h0, then from h0 to t0.
	//

	Vec3<T> h0 = (f0 + t0).normalized();

	if ((h0 ^ h0) != 0)
	{
	    setRotationInternal (f0, h0, *this);

	    Quat<T> q;
	    setRotationInternal (h0, t0, q);

	    *this *= q;
	}
	else
	{
	    //
	    // f0 and t0 point in exactly opposite directions.
	    // Pick an arbitrary axis that is orthogonal to f0,
	    // and rotate by pi.
	    //

	    r = T (0);

	    Vec3<T> f02 = f0 * f0;

	    if (f02.x <= f02.y && f02.x <= f02.z)
		v = (f0 % Vec3<T> (1, 0, 0)).normalized();
	    else if (f02.y <= f02.z)
		v = (f0 % Vec3<T> (0, 1, 0)).normalized();
	    else
		v = (f0 % Vec3<T> (0, 0, 1)).normalized();
	}
    }

    return *this;
}


template <class T>
void
Quat<T>::setRotationInternal (const Vec3<T>& f0, const Vec3<T>& t0, Quat<T> &q)
{
    //
    // The following is equivalent to setAxisAngle(n,2*phi),
    // where the rotation axis, is orthogonal to the f0 and
    // t0 vectors, and 2*phi is the angle between f0 and t0.
    //
    // This function is called by setRotation(), above; it assumes
    // that f0 and t0 are normalized and that the angle between
    // them is not much greater than pi/2.  This function becomes
    // numerically inaccurate if f0 and t0 point into nearly
    // opposite directions.
    //

    //
    // Find a normalized vector, h0, that is half way between f0 and t0.
    // The angle between f0 and h0 is phi.
    //

    Vec3<T> h0 = (f0 + t0).normalized();

    //
    // Store the rotation axis and rotation angle.
    //

    q.r = f0 ^ h0;	//  f0 ^ h0 == cos (phi)
    q.v = f0 % h0;	// (f0 % h0).length() == sin (phi)
}


template<class T>
Matrix33<T> Quat<T>::toMatrix33() const
{
    return Matrix33<T>(1. - 2.0 * (v.y * v.y + v.z * v.z),
			    2.0 * (v.x * v.y + v.z * r),
			    2.0 * (v.z * v.x - v.y * r),

			    2.0 * (v.x * v.y - v.z * r),
		       1. - 2.0 * (v.z * v.z + v.x * v.x),
			    2.0 * (v.y * v.z + v.x * r),

			    2.0 * (v.z * v.x + v.y * r),
			    2.0 * (v.y * v.z - v.x * r),
		       1. - 2.0 * (v.y * v.y + v.x * v.x));
}

template<class T>
Matrix44<T> Quat<T>::toMatrix44() const
{
    return Matrix44<T>(1. - 2.0 * (v.y * v.y + v.z * v.z),
			    2.0 * (v.x * v.y + v.z * r),
			    2.0 * (v.z * v.x - v.y * r),
			    0.,
			    2.0 * (v.x * v.y - v.z * r),
		       1. - 2.0 * (v.z * v.z + v.x * v.x),
			    2.0 * (v.y * v.z + v.x * r),
			    0.,
			    2.0 * (v.z * v.x + v.y * r),
			    2.0 * (v.y * v.z - v.x * r),
		       1. - 2.0 * (v.y * v.y + v.x * v.x),
			    0.,
			    0.,
			    0.,
			    0.,
			    1.0 );
}


template<class T>
inline Matrix33<T> operator* (const Matrix33<T> &M, const Quat<T> &q)
{
    return M * q.toMatrix33();
}

template<class T>
inline Matrix33<T> operator* (const Quat<T> &q, const Matrix33<T> &M)
{
    return q.toMatrix33() * M;
}

template<class T>
std::ostream& operator<< (std::ostream &o, const Quat<T> &q)
{
    return o << "(" << q.r
	     << " " << q.v.x
	     << " " << q.v.y
	     << " " << q.v.z
	     << ")";

}

template<class T>
inline Quat<T> operator* (const Quat<T>& q1, const Quat<T>& q2)
{
    // (S1+V1) (S2+V2) = S1 S2 - V1.V2 + S1 V2 + V1 S2 + V1 x V2
    return Quat<T>( q1.r * q2.r - (q1.v ^ q2.v),
		    q1.r * q2.v + q1.v * q2.r + q1.v % q2.v );
}

template<class T>
inline Quat<T> operator/ (const Quat<T>& q1, const Quat<T>& q2)
{
    return q1 * q2.inverse();
}

template<class T>
inline Quat<T> operator/ (const Quat<T>& q,T t)
{
    return Quat<T>(q.r/t,q.v/t);
}

template<class T>
inline Quat<T> operator* (const Quat<T>& q,T t)
{
    return Quat<T>(q.r*t,q.v*t);
}

template<class T>
inline Quat<T> operator* (T t, const Quat<T>& q)
{
    return Quat<T>(q.r*t,q.v*t);
}

template<class T>
inline Quat<T> operator+ (const Quat<T>& q1, const Quat<T>& q2)
{
    return Quat<T>( q1.r + q2.r, q1.v + q2.v );
}

template<class T>
inline Quat<T> operator- (const Quat<T>& q1, const Quat<T>& q2)
{
    return Quat<T>( q1.r - q2.r, q1.v - q2.v );
}

template<class T>
inline Quat<T> operator~ (const Quat<T>& q)
{
    return Quat<T>( q.r, -q.v );	// conjugate: (S+V)* = S-V
}

template<class T>
inline Quat<T> operator- (const Quat<T>& q)
{
    return Quat<T>( -q.r, -q.v );
}

template<class T>
inline Vec3<T> operator* (const Vec3<T>& v, const Quat<T>& q)
{
    Vec3<T> a = q.v % v;
    Vec3<T> b = q.v % a;
    return v + T (2) * (q.r * a + b);
}

#if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
#pragma warning(default:4244)
#endif

} // namespace Imath

#endif