path: root/Runtime/Math/AnimationCurve.cpp

#include "UnityPrefix.h"
#include "AnimationCurve.h"
#include "Runtime/Utilities/LogAssert.h"
using namespace std;

#define kOneThird (1.0F / 3.0F)
#define kMaxTan 5729577.9485111479F

int ToInternalInfinity (int pre);
int FromInternalInfinity (int pre);

int ToInternalInfinity (int pre)
{
	if (pre == kRepeat)
		return AnimationCurve::kInternalRepeat;
	else if (pre == kPingPong)
		return AnimationCurve::kInternalPingPong;
	else
		return AnimationCurve::kInternalClamp;
}

int FromInternalInfinity (int pre)
{
	if (pre == AnimationCurve::kInternalRepeat)
		return kRepeat;
	else if (pre == AnimationCurve::kInternalPingPong)
		return kPingPong;
	else
		return kClampForever;
}

template<class T>
KeyframeTpl<T>::KeyframeTpl (float t, const T& v)
{
	time = t;
	value = v;
	inSlope = Zero<T>();
	outSlope = Zero<T>();

#if UNITY_EDITOR
	tangentMode = 0;
#endif
}


template<class T>
void AnimationCurveTpl<T>::InvalidateCache ()
{
	m_Cache.time = std::numeric_limits<float>::infinity ();
	m_Cache.index = 0;
	m_ClampCache.time = std::numeric_limits<float>::infinity ();
	m_ClampCache.index = 0;
}

template<class T>
pair<float, float> AnimationCurveTpl<T>::GetRange () const
{
	if (!m_Curve.empty ())
		return make_pair (m_Curve[0].time, m_Curve.back ().time);
	else
		return make_pair (std::numeric_limits<float>::infinity (), -std::numeric_limits<float>::infinity ());
}



///@TODO: Handle step curves correctly
template<class T>
void AnimationCurveTpl<T>::EvaluateWithoutCache (float curveT, T& output)const
{
	DebugAssertIf (!IsValid ());
	curveT = WrapTime (curveT);

	int lhsIndex, rhsIndex;
	FindIndexForSampling (m_Cache, curveT, lhsIndex, rhsIndex);
	const Keyframe& lhs = m_Curve[lhsIndex];
	const Keyframe& rhs = m_Curve[rhsIndex];

	float dx = rhs.time - lhs.time;
	T m1;
	T m2;
	float t;
	if (dx != 0.0F)
	{
		t = (curveT - lhs.time) / dx;
		m1 = lhs.outSlope * dx;
		m2 = rhs.inSlope * dx;
	}
	else
	{
		t = 0.0F;
		m1 = Zero<T>();
		m2 = Zero<T>();
	}

	output = HermiteInterpolate (t, lhs.value, m1, m2, rhs.value);
	HandleSteppedCurve(lhs, rhs, output);
	DebugAssertIf(!IsFinite(output));
}

template<class T>
inline void EvaluateCache (const typename AnimationCurveTpl<T>::Cache& cache, float curveT, T& output)
{
//	DebugAssertIf (curveT < cache.time - kCurveTimeEpsilon || curveT > cache.timeEnd + kCurveTimeEpsilon);
	float t = curveT - cache.time;
	output = (t * (t * (t * cache.coeff[0] + cache.coeff[1]) + cache.coeff[2])) + cache.coeff[3];
	DebugAssertIf (!IsFinite(output));
}

void SetupStepped (float* coeff, const KeyframeTpl<float>& lhs, const KeyframeTpl<float>& rhs)
{
	// If either of the tangents in the segment are set to stepped, make the constant value equal the value of the left key
	if (lhs.outSlope == std::numeric_limits<float>::infinity() || rhs.inSlope == std::numeric_limits<float>::infinity())
	{
		coeff[0] = 0.0F;
		coeff[1] = 0.0F;
		coeff[2] = 0.0F;
		coeff[3] = lhs.value;
	}
}

void HandleSteppedCurve (const KeyframeTpl<float>& lhs, const KeyframeTpl<float>& rhs, float& value)
{
	if (lhs.outSlope == std::numeric_limits<float>::infinity() || rhs.inSlope == std::numeric_limits<float>::infinity())
		value = lhs.value;
}

void HandleSteppedTangent (const KeyframeTpl<float>& lhs, const KeyframeTpl<float>& rhs, float& tangent)
{
	if (lhs.outSlope == std::numeric_limits<float>::infinity() || rhs.inSlope == std::numeric_limits<float>::infinity())
		tangent = std::numeric_limits<float>::infinity();
}


void SetupStepped (Vector3f* coeff, const KeyframeTpl<Vector3f>& lhs, const KeyframeTpl<Vector3f>& rhs)
{
	for (int i=0;i<3;i++)
	{
		// If either of the tangents in the segment are set to stepped, make the constant value equal the value of the left key
		if (lhs.outSlope[i] == std::numeric_limits<float>::infinity() || rhs.inSlope[i] == std::numeric_limits<float>::infinity())
		{
			coeff[0][i] = 0.0F;
			coeff[1][i] = 0.0F;
			coeff[2][i] = 0.0F;
			coeff[3][i] = lhs.value[i];
		}
	}
}

void HandleSteppedCurve (const KeyframeTpl<Vector3f>& lhs, const KeyframeTpl<Vector3f>& rhs, Vector3f& value)
{
	for (int i=0;i<3;i++)
	{
		if (lhs.outSlope[i] == std::numeric_limits<float>::infinity() || rhs.inSlope[i] == std::numeric_limits<float>::infinity())
			value[i] = lhs.value[i];
	}
}

void HandleSteppedTangent (const KeyframeTpl<Vector3f>& lhs, const KeyframeTpl<Vector3f>& rhs, Vector3f& value)
{
	for (int i=0;i<3;i++)
	{
		if (lhs.outSlope[i] == std::numeric_limits<float>::infinity() || rhs.inSlope[i] == std::numeric_limits<float>::infinity())
			value[i] = std::numeric_limits<float>::infinity();
	}
}

void SetupStepped (Quaternionf* coeff, const KeyframeTpl<Quaternionf>& lhs, const KeyframeTpl<Quaternionf>& rhs)
{
	// If either of the tangents in the segment are set to stepped, make the constant value equal the value of the left key
	if (lhs.outSlope[0] == std::numeric_limits<float>::infinity() || rhs.inSlope[0] == std::numeric_limits<float>::infinity() ||
		lhs.outSlope[1] == std::numeric_limits<float>::infinity() || rhs.inSlope[1] == std::numeric_limits<float>::infinity() ||
		lhs.outSlope[2] == std::numeric_limits<float>::infinity() || rhs.inSlope[2] == std::numeric_limits<float>::infinity() ||
		lhs.outSlope[3] == std::numeric_limits<float>::infinity() || rhs.inSlope[3] == std::numeric_limits<float>::infinity() )
	{
		for (int i=0;i<4;i++)
		{
			coeff[0][i] = 0.0F;
			coeff[1][i] = 0.0F;
			coeff[2][i] = 0.0F;
			coeff[3][i] = lhs.value[i];
		}
	}
}

void HandleSteppedCurve (const KeyframeTpl<Quaternionf>& lhs, const KeyframeTpl<Quaternionf>& rhs, Quaternionf& value)
{
	if (lhs.outSlope[0] == std::numeric_limits<float>::infinity() || rhs.inSlope[0] == std::numeric_limits<float>::infinity() ||
		lhs.outSlope[1] == std::numeric_limits<float>::infinity() || rhs.inSlope[1] == std::numeric_limits<float>::infinity() ||
		lhs.outSlope[2] == std::numeric_limits<float>::infinity() || rhs.inSlope[2] == std::numeric_limits<float>::infinity() ||
		lhs.outSlope[3] == std::numeric_limits<float>::infinity() || rhs.inSlope[3] == std::numeric_limits<float>::infinity() )
	{
		value = lhs.value;
	}
}

void HandleSteppedTangent (const KeyframeTpl<Quaternionf>& lhs, const KeyframeTpl<Quaternionf>& rhs, Quaternionf& tangent)
{
	for (int i=0;i<4;i++)
	{
		if (lhs.outSlope[i] == std::numeric_limits<float>::infinity() || rhs.inSlope[i] == std::numeric_limits<float>::infinity())
			tangent[i] = std::numeric_limits<float>::infinity();
	}
}

template<class T>
void AnimationCurveTpl<T>::CalculateCacheData (Cache& cache, int lhsIndex, int rhsIndex, float timeOffset) const
{
	const Keyframe& lhs = m_Curve[lhsIndex];
	const Keyframe& rhs = m_Curve[rhsIndex];
//	DebugAssertIf (timeOffset < -0.001F || timeOffset - 0.001F > rhs.time - lhs.time);
	cache.index = lhsIndex;
	cache.time = lhs.time + timeOffset;
	cache.timeEnd = rhs.time + timeOffset;
	cache.index = lhsIndex;

	float dx, length;
	T dy;
	T m1, m2, d1, d2;

	dx = rhs.time - lhs.time;
	dx = max(dx, 0.0001F);
	dy = rhs.value - lhs.value;
	length = 1.0F / (dx * dx);

	m1 = lhs.outSlope;
	m2 = rhs.inSlope;
	d1 = m1 * dx;
	d2 = m2 * dx;

	cache.coeff[0] = (d1 + d2 - dy - dy) * length / dx;
	cache.coeff[1] = (dy + dy + dy - d1 - d1 - d2) * length;
	cache.coeff[2] = m1;
	cache.coeff[3] = lhs.value;
	SetupStepped(cache.coeff, lhs, rhs);

	DebugAssertIf(!IsFinite(cache.coeff[0]));
	DebugAssertIf(!IsFinite(cache.coeff[1]));
	DebugAssertIf(!IsFinite(cache.coeff[2]));
	DebugAssertIf(!IsFinite(cache.coeff[3]));
}

// When we look for the next index, how many keyframes do we just loop ahead instead of binary searching?
#define SEARCH_AHEAD 3
///@TODO: Cleanup old code to completely get rid of this
template<class T>
int AnimationCurveTpl<T>::FindIndex (const Cache& cache, float curveT) const
{
	#if SEARCH_AHEAD >= 0
	int cacheIndex = cache.index;
	if (cacheIndex != -1)
	{
		// We can not use the cache time or time end since that is in unwrapped time space!
		float time = m_Curve[cacheIndex].time;

		if (curveT > time)
		{
			if (cacheIndex + SEARCH_AHEAD < static_cast<int>(m_Curve.size()))
			{
				for (int i=0;i<SEARCH_AHEAD;i++)
				{
					if (curveT < m_Curve[cacheIndex + i + 1].time)
						return cacheIndex + i;
				}
			}
		}
		else
		{
			if (cacheIndex - SEARCH_AHEAD >= 0)
			{
				for (int i=0;i<SEARCH_AHEAD;i++)
				{
					if (curveT > m_Curve[cacheIndex - i - 1].time)
						return cacheIndex - i - 1;
				}
			}
		}
	}

	#endif

	///@ use cache to index into next if not possible use binary search
	const_iterator i = std::lower_bound (m_Curve.begin (), m_Curve.end (), curveT, KeyframeCompare());
	int index = distance (m_Curve.begin (), i);
	index--;
	index = min<int> (m_Curve.size () - 2, index);
	index = max<int> (0, index);

	return index;
}

///@TODO: Cleanup old code to completely get rid of this
template<class T>
int AnimationCurveTpl<T>::FindIndex (float curveT) const
{
	pair<float, float> range = GetRange ();
	if (curveT <= range.first || curveT >= range.second)
		return -1;

	const_iterator i = std::lower_bound (m_Curve.begin (), m_Curve.end (), curveT, KeyframeCompare());
	AssertIf (i == m_Curve.end ());
	int index = distance (m_Curve.begin (), i);
	index--;
	index = min<int> (m_Curve.size () - 2, index);
	index = max<int> (0, index);

	AssertIf (curveT < m_Curve[index].time || curveT > m_Curve[index+1].time);
	return index;
}

template<class T>
void AnimationCurveTpl<T>::FindIndexForSampling (const Cache& cache, float curveT, int& lhs, int& rhs) const
{
	AssertIf (curveT < GetRange ().first || curveT > GetRange ().second);
	int actualSize = m_Curve.size();
	const Keyframe* frames = &m_Curve[0];

	// Reference implementation:
	// (index is the last value that is equal to or smaller than curveT)
	#if 0
	int foundIndex = 0;
	for (int i=0;i<actualSize;i++)
	{
		if (frames[i].time <= curveT)
			foundIndex = i;
	}

	lhs = foundIndex;
	rhs = min<int>(lhs + 1, actualSize - 1);
	AssertIf (curveT < m_Curve[lhs].time || curveT > m_Curve[rhs].time);
	AssertIf(frames[rhs].time == curveT && frames[lhs].time != curveT)
	return;
	#endif


	#if SEARCH_AHEAD > 0
	int cacheIndex = cache.index;
	if (cacheIndex != -1)
	{
		// We can not use the cache time or time end since that is in unwrapped time space!
		float time = m_Curve[cacheIndex].time;

		if (curveT > time)
		{
			for (int i=0;i<SEARCH_AHEAD;i++)
			{
				int index = cacheIndex + i;
				if (index + 1 < actualSize && frames[index + 1].time > curveT)
				{
					lhs = index;

					rhs = min<int>(lhs + 1, actualSize - 1);
					AssertIf (curveT < frames[lhs].time || curveT > frames[rhs].time);
					AssertIf(frames[rhs].time == curveT && frames[lhs].time != curveT);
					return;
				}
			}
		}
		else
		{
			for (int i=0;i<SEARCH_AHEAD;i++)
			{
				int index = cacheIndex - i;
				if (index >= 0 && curveT >= frames[index].time)
				{
					lhs = index;
					rhs = min<int>(lhs + 1, actualSize - 1);
					AssertIf (curveT < frames[lhs].time || curveT > m_Curve[rhs].time);
					AssertIf(frames[rhs].time == curveT && frames[lhs].time != curveT);
					return;
				}
			}
		}
	}

	#endif

	// Fall back to using binary search
	// upper bound (first value larger than curveT)
	int __len = actualSize;
	int __half;
	int __middle;
	int __first = 0;
	while (__len > 0)
	{
		__half = __len >> 1;
		__middle = __first + __half;

		if (curveT < frames[__middle].time)
			__len = __half;
		else
	    {
			__first = __middle;
			++__first;
			__len = __len - __half - 1;
	    }
	}

	// If not within range, we pick the last element twice
	lhs = __first - 1;
	rhs = min(actualSize - 1, __first);

	AssertIf(lhs < 0 || lhs >= actualSize);
	AssertIf(rhs < 0 || rhs >= actualSize);

	AssertIf (curveT < m_Curve[lhs].time || curveT > m_Curve[rhs].time);
	AssertIf(frames[rhs].time == curveT && frames[lhs].time != curveT);
}

template<class T>
void AnimationCurveTpl<T>::SetPreInfinity (int pre)
{
	m_PreInfinity = ToInternalInfinity(pre);
	InvalidateCache ();
}

template<class T>
void AnimationCurveTpl<T>::SetPostInfinity (int post)
{
	m_PostInfinity = ToInternalInfinity(post);
	InvalidateCache ();
}

template<class T>
int AnimationCurveTpl<T>::GetPreInfinity () const
{
	return FromInternalInfinity(m_PreInfinity);
}

template<class T>
int AnimationCurveTpl<T>::GetPostInfinity () const
{
	return FromInternalInfinity(m_PostInfinity);
}

template<class T>
T AnimationCurveTpl<T>::EvaluateClamp (float curveT) const
{
	T output;
	if (curveT >= m_ClampCache.time && curveT < m_ClampCache.timeEnd)
	{
//		AssertIf (!CompareApproximately (EvaluateCache (m_Cache, curveT), EvaluateWithoutCache (curveT), 0.001F));
		EvaluateCache<T> (m_ClampCache, curveT, output);
		return output;
	}
	else
	{
		DebugAssertIf (!IsValid ());

		float begTime = m_Curve[0].time;
		float endTime = m_Curve.back().time;

		if (curveT > endTime)
		{
			m_ClampCache.time = endTime;
			m_ClampCache.timeEnd = std::numeric_limits<float>::infinity ();
			m_ClampCache.coeff[0] = m_ClampCache.coeff[1] = m_ClampCache.coeff[2] = Zero<T>();
			m_ClampCache.coeff[3] = m_Curve[m_Curve.size()-1].value;
		}
		else if (curveT < begTime)
		{
			m_ClampCache.time = curveT - 1000.0F;
			m_ClampCache.timeEnd = begTime;
			m_ClampCache.coeff[0] = m_ClampCache.coeff[1] = m_ClampCache.coeff[2] = Zero<T>();
			m_ClampCache.coeff[3] = m_Curve[0].value;
		}
		else
		{
			int lhs, rhs;
			FindIndexForSampling (m_ClampCache, curveT, lhs, rhs);
			CalculateCacheData (m_ClampCache, lhs, rhs, 0.0F);
		}

//		AssertIf (!CompareApproximately (EvaluateCache (m_Cache, curveT), EvaluateWithoutCache (curveT), 0.001F));
		EvaluateCache<T> (m_ClampCache, curveT, output);
		return output;
	}
}

template<class T>
T AnimationCurveTpl<T>::Evaluate (float curveT) const
{
	int lhs, rhs;
	T output;
	if (curveT >= m_Cache.time && curveT < m_Cache.timeEnd)
	{
//		AssertIf (!CompareApproximately (EvaluateCache (m_Cache, curveT), EvaluateWithoutCache (curveT), 0.001F));
		EvaluateCache<T> (m_Cache, curveT, output);
		return output;
	}
	// @TODO: Optimize IsValid () away if by making the non-valid case always use the m_Cache codepath
	else if (IsValid ())
	{
		float begTime = m_Curve[0].time;
		float endTime = m_Curve.back().time;
		float wrappedTime;

		if (curveT >= endTime)
		{
			if (m_PostInfinity == kInternalClamp)
			{
				m_Cache.time = endTime;
				m_Cache.timeEnd = std::numeric_limits<float>::infinity ();
				m_Cache.coeff[0] = m_Cache.coeff[1] = m_Cache.coeff[2] = Zero<T>();
				m_Cache.coeff[3] = m_Curve[m_Curve.size()-1].value;
			}
			else if (m_PostInfinity == kInternalRepeat)
			{
				wrappedTime = Repeat (curveT, begTime, endTime);

				FindIndexForSampling (m_Cache, wrappedTime, lhs, rhs);
				CalculateCacheData (m_Cache, lhs, rhs, curveT - wrappedTime);
			}
			///@todo optimize pingpong by making it generate a cache too
			else
			{
				EvaluateWithoutCache (curveT, output);
				return output;
			}
		}
		else if (curveT < begTime)
		{
			if (m_PreInfinity == kInternalClamp)
			{
				m_Cache.time = curveT - 1000.0F;
				m_Cache.timeEnd = begTime;
				m_Cache.coeff[0] = m_Cache.coeff[1] = m_Cache.coeff[2] = Zero<T>();
				m_Cache.coeff[3] = m_Curve[0].value;
			}
			else if (m_PreInfinity == kInternalRepeat)
			{
				wrappedTime = Repeat (curveT, begTime, endTime);
				FindIndexForSampling (m_Cache, wrappedTime, lhs, rhs);
				CalculateCacheData (m_Cache, lhs, rhs, curveT - wrappedTime);
			}
			///@todo optimize pingpong by making it generate a cache too
			else
			{
				EvaluateWithoutCache (curveT, output);
				return output;
			}
		}
		else
		{
			FindIndexForSampling (m_Cache, curveT, lhs, rhs);
			CalculateCacheData (m_Cache, lhs, rhs, 0.0F);
		}

		//		AssertIf (!CompareApproximately (EvaluateCache (m_Cache, curveT), EvaluateWithoutCache (curveT), 0.001F));
		EvaluateCache<T> (m_Cache, curveT, output);
		return output;
	}
	else
	{
		if (m_Curve.size () == 1)
			return m_Curve.begin()->value;
		else
			return Zero<T> ();
	}
}

template<class T>
float AnimationCurveTpl<T>::WrapTime (float curveT) const
{
	DebugAssertIf (!IsValid ());

	float begTime = m_Curve[0].time;
	float endTime = m_Curve.back().time;

	if (curveT < begTime)
	{
		if (m_PreInfinity == kInternalClamp)
			curveT = begTime;
		else if (m_PreInfinity == kInternalPingPong)
			curveT = PingPong (curveT, begTime, endTime);
		else
			curveT = Repeat (curveT, begTime, endTime);
	}
	else if (curveT > endTime)
	{
		if (m_PostInfinity == kInternalClamp)
			curveT = endTime;
		else if (m_PostInfinity == kInternalPingPong)
			curveT = PingPong (curveT, begTime, endTime);
		else
			curveT = Repeat (curveT, begTime, endTime);
	}
	return curveT;
}

template<class T>
int AnimationCurveTpl<T>::AddKey (const Keyframe& key)
{
	InvalidateCache ();

	iterator i = std::lower_bound (m_Curve.begin (), m_Curve.end (), key);

	// is not included in container and value is not a duplicate
	if (i == end () || key < *i)
	{
		iterator ii = m_Curve.insert (i, key);
		return std::distance (m_Curve.begin (), ii);
	}
	else
		return -1;
}

template<class T>
void AnimationCurveTpl<T>::RemoveKeys (iterator begin, iterator end)
{
	InvalidateCache ();
	m_Curve.erase (begin, end);
}

void ScaleCurveValue (AnimationCurve& curve, float scale)
{
	for (int i=0;i<curve.GetKeyCount ();i++)
	{
		curve.GetKey (i).value *= scale;
		curve.GetKey (i).inSlope *= scale;
		curve.GetKey (i).outSlope *= scale;
	}
	
	curve.InvalidateCache();
}

void OffsetCurveValue (AnimationCurve& curve, float offset)
{
	for (int i=0;i<curve.GetKeyCount ();i++)
		curve.GetKey (i).value += offset;
	
	curve.InvalidateCache();
}

void ScaleCurveTime (AnimationCurve& curve, float scale)
{
	for (int i=0;i<curve.GetKeyCount ();i++)
	{
		curve.GetKey (i).time *= scale;
		curve.GetKey (i).inSlope /= scale;
		curve.GetKey (i).outSlope /= scale;
	}
	curve.InvalidateCache();
}

void OffsetCurveTime (AnimationCurve& curve, float offset)
{
	for (int i=0;i<curve.GetKeyCount ();i++)
		curve.GetKey (i).time += offset;
	curve.InvalidateCache();
}


/*

Calculating tangents from a hermite spline () Realtime rendering page 56



> On this first pass we're stuck with linear keyframing, because we just
> don't have the cycles in game to go to splines. I know this would help a
> lot, but it isn't an option.

In Granny I do successive least-squares approximations do the data.  I
take the array of samples for a given channel (ie., position) which is 2x
oversampled or better from the art tool.  I then start by doing a
least-square solve for a spline of arbitrary degree with knots at either
end.  I compute the error over the spline from the original samples, and
add knots in the areas of highest error.  Repeat and salt to taste.

Since it's fairly easy to write a single solver that solves for any degree
of spline, I use the same solver for linear keyframes as I do for
quadratic keyframes, cubic keyframes, or even "0th order keframes", which
is to say if you don't want to interpolate _at all_, the solver can still
place the "stop-motion" key frames in the best locations.  But hopefully
no one is still doing that kind of animation.

If I were doing this over again (which I probably will at some point), I
would probably use some kind of weird waveletty scheme now.  I decided not
to do that originally, because I had about a month to do the entire
spline/solver/reduction thing the first time, and it had to be highly
optimized.  So I didn't want to have to learn wavelets and spline solvers
at the same time.  Next time I'll spend some time on wavelets and probably
use some sort of hierachical reduction scheme instead, primarily so you
can change how densely your keyframes are placed at run-time, and so you
can get hard edges easier (motions which are intended to have sharp
discontinuities aren't handled well at all by my current scheme).

- Casey



--

> Is anybody looking at the way errors add up?  eg. if you do some
> reduction on the hip, and the thigh, and the ankle bone channels, the
> result is a foot that moves quite a bit differently than it should.

This has been on my list for a long time.  I don't think it's a simple
case of error analysis though.  I think it's more a case for using
discrete skeletal changes or integrated IK, because hey, if what you care
about is having a particular thing stay in one place, then it seems to me
that the best way to compress that is by just saying "this stays in one
place", rather than trying to spend a lot of data on the joint curves
necessary to make that so.

> Also, is anybody doing things like using IK in the reducer, to make sure
> that even in the reduced version the feet stay in exactly the same spot?

That doesn't work with splines, unfortunately, because the splines are
continuous, and you will always have the feet slipping as a result (if not
at the keyframes, then in between for sure).  So you end up with the
problem that the IK reducer would need to shove a metric assload of keys
into the streams, which defeats the compression.

From Ian:

> You said you are doing this on linear data as well. I can see that this
> would work, but do you find you get a fairly minimal result? I can
> envisage cases where you'd get many unneeded keyframes from initial
> 'breaks' at points of large error.

I'm not sure what you mean by this.  The error is controllable, so you say
how much you are willing to accept.  You get a minimal spline for the
error that you ask for, but no more - obviously in the areas of high
error, it has to add more keys, but that's what you want.  The objective
of compression, at least in my opinion, is not to remove detail, but
rather to more efficiently store places where there is less detail.

> I think I may be missing the point. Do you do least squares on the
> already fitted cures (ie, a chi squared test) then curve fit separately,
> or do you use least squares to fit your actual spline to the data?

The latter.  The incremental knot addition sets up the t_n's of an
arbitrary degree spline.  You have a matrix A that has sample-count rows
and knot-count columns.  The vector you're looking for is x, the spline
vector, which has knot-count rows.  You're producing b, the vector of
samples, which has sample-count rows.  So it's Ax = b, but A is
rectangular.  You solve via A^T Ax = A^Tb, a simple least squares problem.
You can solve it any way you like.  I chose a straightforward
implementation, because there's no numerical difficulties with these
things.

The version that comes with Granny is a highly optimized A^T Ax = A^Tb
solver, which constructs A^T A and A^T b directly and sparsely (so it is
O(n) in the number of samples, instead of O(n^3)).  The A^T A is band
diagonal, because of the sparsity pattern of A, so I use a sparse banded
cholesky solver on the back end.  It's _extremely_ fast.  In fact, the
stupid part of Granny's solver is actually the other part (the error
analysis + knot addition), because it adds very few knots per cycle, so on
a really long animation it can call the least-squares Ax = b solver many
hundreds of times for a single animation, and it _still_ never takes more
than 10 seconds or so even on animations with many hundred frames.  So
really, the right place to fix currently is the stupid knot addition
alogirithm, which really should be improved quite a bit.

> I have not used least squared before and from the reading I have done
> (numerical recipes, last night), it seem that the equation of the
> line/spline are inherent in the method and it needs to be reformulated
> to use a different line/spline. (The version I read was for fitting a
> straight line to a data-set, which I know won't work for quaternions
> slerps anyway)

Quaternion lerps are great, and work great with splines.  I use them
throughout all of Granny - there is no slerping anywhere.  Slerping is not
a very useful thing in a run-time engine, in my opinion, unless you are
dealing with orientations that are greater the 90 degrees apart.  It
allows me to use the same solver for position, rotation, AND scale/shear,
with no modifications.

- Casey


-----
> When is it best to use quaternions and when to use normal matrix
> rotations?

In Granny, I use quaternions for everything except the final composite
phase where rotations (and everything else) are multiplied by their
parents to produce the final world-space results.  Since we support
scale/shear, orientation, and position, it tends to be fastest at the
composition stage to convert the orientation from quaternion to matrix,
and then do all the matrix concatenation together.  I'm not sure I would
do the same if I didn't support scale/shear.  It might be a bit more fun
and efficient to go ahead and do the whole pipe in quaternion and only
convert to matrices at the very end (maybe somebody else has played with
that and can comment?)

> And for character animation with moving joints quaternions is usually
> used, but why?

There's lots of nice things about quaternions:

1) You can treat the linearly if you want to (so you can blend animations
quickly and easily)

2) You can treat them non-linearly if you want to (so you can do exact
geodesic interpolation at a constant speed)

3) You can convert them to a matrix with no transcendentals

4) They are compact, and can be easily converted to a 3-element
representation when necessary (ie., any one of the four values can be
easily re-generated from the other 3)

5) They cover 720 degrees of rotation, not 360 (if you do a lot of work
with character animation, you will appreciate this!), and the
neighborhooding operator is extremely simple and fast (inner product and a
negation)

6) They are easy to visualize (they're not much more complicated than
angle/axis) and manipulate (ie., the axis of rotation is obvious, you can
transform the axis of rotation without changing the amount of rotation,
and even just use a regular 3-vector transform, etc.)

7) You can easily transform vectors with them without converting to a
matrix if you want to

8) They are easy to renormalize, and they can never become "unorthogonal"

9) No gimbal lock

10) There is a simple, fast, and relevant distance metric (the inner
product)

> Is it just by trial and error that we decide when we have to use
> quarternions?

Well, I can't speak for everyone, but I have been extremely careful in my
selection of quaternions.  With Granny 1 I did quaternions to get some
practice with them, but when I did 2 I did a lot of research and wrote
code to visualize the actions of various rotational representations and
operations, and I can very confidently say there's no better choice for
general character animation than quaternions, at least among the
representations that I know (euler, angle/axis, matrix, quaternion, exp
map).  There are some other mappings that I've come across since I worked
everything out (like the Rational Map), that I have not experimented with,
because so far quaternions have been working superbly so I haven't had
much cause to go hunting.

> "Quarternions have the ability to have smooth interpollated rotations",
> I've made smooth interpollated rotations with normal matrix rotations
> though.

Well, it's not so much what you can do as how easy it is to do it.
Quaternions can be interpolated directly as a smooth geodesic (via slerp)
or via a nonlinear geodesic (lerp).  The latter is particularly powerful
because it distributes - it's a linear operator.  So you can literally
plug it in to anything that would work, like splines, multi-point blends,
etc., and it will "just work" (well, that's a bit of an exaggeration,
because there is a neighborhooding concern, but it's always easy to do).

> "They take less room, 4 elements versus 9 and some operations
> are cheaper in terms of CPU cycles". I accept this reason except someone
> said that quarternions use more CPU cycles, so I'm not sure who to
> believe.

It depends what you're doing.  Quaternions take more CPU cycles to
transform a vector than a matrix does.  But quaternions are MUCH less
expensive for lots of other operations, like, for example, splining.
This is why it's very useful to use quaternions for everything except the
final stage of your pipe, where matrices are more appropriate.

> "Quarternions are not susceptible to gimbal lock. Gimbal lock shows its
> face when two axes point in the same direction." So if no axes face in
> the same direction then this is not a reason to use quarternions.
> ...
> Am I right in saying that if no axes face the same direction it is
> possible to represent all rotations with normal matrix rotations?

Matrix rotations are not subject to gimbal lock.  That's Euler angles.

- Casey


Least squares is very simple.  Given a set of data points (2x oversampled in caseys case) and a knot vector, compute the coefficients for the control points of the spline that minimize squared error.  This involves solving a linear system where each row corresponds to a data point (di - they are uniformly space so each point has a corresponding ti which is in the domain of the curve) and each column corresponds to a control point for the curve (cj the unkowns).  So matrix element Aij is Bj(ti) where Bj is the basis function for the jth control point.  The right hand side is just a vector of the di and the x's are the unkown control points.  Sine you are compressing things there will be many more rows then columns in this linear system so you don't have to worry about over fitting problems (which you would if there were more DOF...)

-Peter-Pike

	-----Original Message-----
	From: Ian Elsley [mailto:ielsley@kushgames.com]
	Sent: Tue 12/17/2002 10:35 AM
	To: gdalgorithms-list@lists.sourceforge.net
	Cc:
	Subject: RE: [Algorithms] Re: Keyframe reduction



	Casey

	I think I may be missing the point. Do you do least squares on the
	already fitted cures (ie, a chi squared test) then curve fit separately,
	or do you use least squares to fit your actual spline to the data?

	I have not used least squared before and from the reading I have done
	(numerical recipes, last night), it seem that the equation of the
	line/spline are inherent in the method and it needs to be reformulated
	to use a different line/spline. (The version I read was for fitting a
	straight line to a data-set, which I know won't work for quaternions
	slerps anyway)

	I see the elegance of this method and would love to work it out, but
	I've got a few blanks here.

	Any pointers?

	Thanks in advance,

	Ian






         q
  q' = -----
        |q|

For renormalizing, since you are very close to 1, I usually use the
tangent-line approximation, which was suggested by Checker a long long
time ago for 3D vectors and works just peachy for 4D ones as well:

inline void NormalizeCloseToOne4(float *Dest)
{
    float const Sum = (Dest[0] * Dest[0] +
                       Dest[1] * Dest[1] +
                       Dest[2] * Dest[2] +
                       Dest[3] * Dest[3]);

    float const ApproximateOneOverRoot = (3.0f - Sum) * 0.5f;

    Dest[0] *= ApproximateOneOverRoot;
    Dest[1] *= ApproximateOneOverRoot;
    Dest[2] *= ApproximateOneOverRoot;
    Dest[3] *= ApproximateOneOverRoot;
}


*/



#define INSTANTIATE(T) \
template EXPORT_COREMODULE std::pair<float, float> AnimationCurveTpl<T>::GetRange() const; \
template EXPORT_COREMODULE T AnimationCurveTpl<T>::Evaluate (float curveT) const; \
template EXPORT_COREMODULE void AnimationCurveTpl<T>::RemoveKeys (iterator begin, iterator end);\
template EXPORT_COREMODULE int AnimationCurveTpl<T>::AddKey (const Keyframe& key);\
template EXPORT_COREMODULE int AnimationCurveTpl<T>::FindIndex (float time) const; \
template EXPORT_COREMODULE int AnimationCurveTpl<T>::FindIndex (const Cache& cache, float time) const; \
template EXPORT_COREMODULE void AnimationCurveTpl<T>::InvalidateCache (); \
template EXPORT_COREMODULE void AnimationCurveTpl<T>::SetPreInfinity (int mode); \
template EXPORT_COREMODULE void AnimationCurveTpl<T>::SetPostInfinity (int mode); \
template EXPORT_COREMODULE int AnimationCurveTpl<T>::GetPreInfinity () const; \
template EXPORT_COREMODULE int AnimationCurveTpl<T>::GetPostInfinity () const; \
template EXPORT_COREMODULE KeyframeTpl<T>::KeyframeTpl (float time, const T& value);

INSTANTIATE(float)
INSTANTIATE(Vector3f)
INSTANTIATE(Quaternionf)

template EXPORT_COREMODULE void AnimationCurveTpl<float>::CalculateCacheData (Cache& cache, int lhs, int rhs, float timeOffset) const;
template EXPORT_COREMODULE float AnimationCurveTpl<float>::EvaluateClamp (float curveT) const;
template EXPORT_COREMODULE Quaternionf AnimationCurveTpl<Quaternionf>::EvaluateClamp (float curveT) const;
template EXPORT_COREMODULE Vector3f AnimationCurveTpl<Vector3f>::EvaluateClamp (float curveT) const;