If I understood correctly, TCB spline interpolation is supported in Collada, however, I can find nothing in the specification about how it’s implemented. It also could be that it’s a 3d studio max specific “upgrade” to the collada standard. However, to support loading 3DSM models correctly, i’m forced to include support for TCB splines.
Here’s the interesting part in the collada file:
<source id="Bone14-node-rz_Bone14-node_RotZ.ANGLE-tcbs">
<float_array id="Bone14-node-rz_Bone14-node_RotZ.ANGLE-tcbs-array" count="24">0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5</float_array>
<technique_common>
<accessor source="#Bone14-node-rz_Bone14-node_RotZ.ANGLE-tcbs-array" count="8" stride="3">
<param type="float"/>
<param type="float"/>
<param type="float"/>
</accessor>
</technique_common>
</source>
<source id="Bone14-node-rz_Bone14-node_RotZ.ANGLE-eases">
<float_array id="Bone14-node-rz_Bone14-node_RotZ.ANGLE-eases-array" count="16">0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5</float_array>
<technique_common>
<accessor source="#Bone14-node-rz_Bone14-node_RotZ.ANGLE-eases-array" count="8" stride="2">
<param type="float"/>
<param type="float"/>
</accessor>
</technique_common>
</source>
<sampler id="Bone14-node-rz_Bone14-node_RotZ.ANGLE-sampler">
<input semantic="INPUT" source="#Bone14-node-rz_Bone14-node_RotZ.ANGLE-input"/>
<input semantic="OUTPUT" source="#Bone14-node-rz_Bone14-node_RotZ.ANGLE-output"/>
<input semantic="INTERPOLATION" source="#Bone14-node-rz_Bone14-node_RotZ.ANGLE-interpolations"/>
<input semantic="TCB" source="#Bone14-node-rz_Bone14-node_RotZ.ANGLE-tcbs"/>
<input semantic="EASE_IN_OUT" source="#Bone14-node-rz_Bone14-node_RotZ.ANGLE-eases"/>
</sampler>
I presume the three floats in TCB are Tension, Continuity, and Bias respectively, and the EASE_IN_OUT floats ease_in and ease_out. Now the only problem, what’s the algorithm/mathematics used to get usable results from these values?
My current implementation is as following, but it’s not even close to acting in the way it’s supposed to be (i used http://www.cubic.org/docs/hermite.htm as my reference) :
inline float TCBTangentEquationIncoming(float& pm1, float& p0, float& p1, float& t, float& c, float& b)
{
return ((1-t)*(1-c)*(1+b))/2*(p0-pm1) + ((1-t)*(1+c)*(1-b))/2*(p1-p0);
}
inline float TCBTangentEquationOutgoing(float& pm1, float& p0, float& p1, float& t, float& c, float& b)
{
return ((1-t)*(1+c)*(1+b))/2*(p0-pm1) + ((1-t)*(1-c)*(1-b))/2*(p1-p0);
}
float* CMatrix::ConcatFloat4Array(float* farr)
{
float* V=new float[4];
V[0] = farr[0]*Matrix[0][0] + farr[1]*Matrix[1][0] + farr[2]*Matrix[2][0] + farr[3]*Matrix[3][0];
V[1] = farr[0]*Matrix[0][1] + farr[1]*Matrix[1][1] + farr[2]*Matrix[2][1] + farr[3]*Matrix[3][1];
V[2] = farr[0]*Matrix[0][2] + farr[1]*Matrix[1][2] + farr[2]*Matrix[2][2] + farr[3]*Matrix[3][2];
V[3] = farr[0]*Matrix[0][3] + farr[1]*Matrix[1][3] + farr[2]*Matrix[2][3] + farr[3]*Matrix[3][3];
return V;
}
CMatrix& GetHermiteMatrix()
{
static bool matrixInited=false;
static CMatrix hm;
if(matrixInited)
return hm;
hm.Matrix[0][0] = 2; hm.Matrix[1][0] = -2; hm.Matrix[2][0] = 1; hm.Matrix[3][0] = 1;
hm.Matrix[0][1] = -3; hm.Matrix[1][1] = 3; hm.Matrix[2][1] = -2; hm.Matrix[3][1] = -1;
hm.Matrix[0][2] = 0; hm.Matrix[1][2] = 0; hm.Matrix[2][2] = 1; hm.Matrix[3][2] = 0;
hm.Matrix[0][3] = 1; hm.Matrix[1][3] = 0; hm.Matrix[2][3] = 0; hm.Matrix[3][3] = 0;
matrixInited = true;
return hm;
}
//! pm1 (y coordinate of the previous point)
//! p0 (y coordinate of the current point)
//! p1 (y coordinate of the next point)
//! p2 (y coordinate of the next+1 point)
//! s (interpolation value, use ApproximateCubicBezierParameter?...)
//! to0 (tangent out of (p0))
//! ti1 (tangent in of (p1))
float InterpolateTCB(float pm1, float p0, float p1, float p2, float s, float to0, float ti1)
{
float S[4], C[4];
S[3] = 1;
S[2] = s; // s
S[1] = s*s; // s^2
S[0] = S[1]*s; // s^3
C[0] = p0;
C[1] = p1;
C[2] = to0;
C[3] = ti1;
// according to http://www.cubic.org/docs/hermite.htm
// resulting vector = S * MHermite * C;
// [b](Vector * Matrix? Usually you do Matrix * Vector right?, this is probably where I went wrong in the TCB interpolation code...)[/b]
float* res=GetHermiteMatrix().ConcatFloat4Array(C); // MHermite * C
res[0] *= S[0]; res[1] *= S[1]; res[2] *= S[2]; res[3] *= S[3]; // S * PreviousResult
float y=res[1];
delete res;
return y;
}
float CFormula::Value(float X)
{
unsigned int ps=PointCount();
if(ps == 0)
return 0.0f;
else if(X < points[0].x)
return points[0].y;
else if(X >= points[ps - 1].x)
return points[ps - 1].y;
int segment=this->PointSegmentLocation(X);
if(InterpolationType == INTERPOLATE_CUBICBEZIER)
{
InterpolatePoint& p1=points[segment];
InterpolatePoint& p2=points[segment+1];
float s = ApproximateCubicBezierParameter( X, p1.x, p1.tangentOutX, p2.tangentInX, p2.x );
float s1 = (1-s);
float s2 = s1*s1;
float ss2 = s*s;
return p1.y*s2*s1+3*p1.tangentOutY*s*s2+3*p2.tangentInY*ss2*s1+p2.y*ss2*s;
}
else if (InterpolationType == INTERPOLATE_TCB)
{
// [b]the first and last points are duplicated when it's out of bounds. No idea if Collada/3DSM intended to do the same?![/b]
InterpolatePoint& pm1=points[segment-1<0?segment:segment-1];
InterpolatePoint& p0=points[segment];
InterpolatePoint& p1=points[segment+1];
InterpolatePoint& p2=points[segment+2>=PointCount()?segment+1:segment+2];
// reusing some point parameters here.., you can also see that there's [b]no ease-in and ease-out[/b] yet, as i had no clue how to implement those.
float& t= p0.tangentInX;
float& c= p0.tangentInY;
float& b= p0.tangentOutX;
// ease variables.
//float& easein= p0.tangentOutY;
//float& easeout= p0.easeout;
float to0x = TCBTangentEquationOutgoing(pm1.x, p0.x, p1.x, t, c, b);
float to0y = TCBTangentEquationOutgoing(pm1.y, p0.y, p1.y, t, c, b);
float ti1x = TCBTangentEquationIncoming(p0.x, p1.x, p2.x, t, c, b);
float ti1y = TCBTangentEquationIncoming(p0.y, p1.y, p2.y, t, c, b);
// [b]this is probably very wrong, TCB != CubicBezier...[/b]
float s = ApproximateCubicBezierParameter( X, p0.x, to0x, ti1x, p1.x );
return InterpolateTCB(pm1.y, p0.y, p1.y, p2.y, s, to0y, ti1y);
}
float Factor=(X-points[segment].x) / (points[segment+1].x-points[segment].x);
switch(InterpolationType)
{
case INTERPOLATE_LINEAR:
return InterpolateLinear(points[segment].y, points[segment+1].y, Factor);
case INTERPOLATE_COSINE:
return InterpolateCosine(points[segment].y, points[segment+1].y, Factor);
case INTERPOLATE_CUBIC:
{
float y0 = points[segment-1<0?segment:segment-1].y, y3=points[segment+2>=PointCount()?segment+1:segment+2].y;
return InterpolateCubic(y0, points[segment].y, points[segment+1].y, y3, Factor);
}
case INTERPOLATE_CATMULLROM:
{
float y0 = points[segment-1<0?segment:segment-1].y, y3=points[segment+2>=PointCount()?segment+1:segment+2].y;
return InterpolateCatmullRom(y0, points[segment].y, points[segment+1].y, y3, Factor);
}
case INTERPOLATE_HERMITE:
{
float y0 = points[segment-1<0?segment:segment-1].y, y3=points[segment+2>=PointCount()?segment+1:segment+2].y;
return InterpolateHermite(y0, points[segment].y, points[segment+1].y, y3, HermiteTension, HermiteBias, Factor);
}
}
}
//simply clamps a value between 0 .. 1
float ClampToZeroOne(float value) {
if (value < .0f)
return .0f;
else if (value > 1.0f)
return 1.0f;
else
return value;
}
// the following function i got from the collada specification, but it's intended for the "BEZIER" type of interpolation, which works fine by the way..
/**
* Returns the approximated parameter of a parametric curve for the value X
* @param atX At which value should the parameter be evaluated
* @param P0_X The first interpolation point of a curve segment
* @param C0_X The first control point of a curve segment
* @param C1_X The second control point of a curve segment
* @param P1_x The second interpolation point of a curve segment
* @return The parametric argument that is used to retrieve atX using the parametric function representation of this curve
*/
float ApproximateCubicBezierParameter (
float atX, float P0_X, float C0_X, float C1_X, float P1_X ) {
if (atX - P0_X < VERYSMALL)
return 0.0;
if (P1_X - atX < VERYSMALL)
return 1.0;
long iterationStep = 0;
float u = 0.0f; float v = 1.0f;
//iteratively apply subdivision to approach value atX
while (iterationStep < MAXIMUM_ITERATIONS) {
// de Casteljau Subdivision.
double a = (P0_X + C0_X)*0.5f;
double b = (C0_X + C1_X)*0.5f;
double c = (C1_X + P1_X)*0.5f;
double d = (a + b)*0.5f;
double e = (b + c)*0.5f;
double f = (d + e)*0.5f; //this one is on the curve!
//The curve point is close enough to our wanted atX
if (fabs(f - atX) < APPROXIMATION_EPSILON)
return ClampToZeroOne((u + v)*0.5f);
//dichotomy
if (f < atX) {
P0_X = f;
C0_X = e;
C1_X = c;
u = (u + v)*0.5f;
} else {
C0_X = a; C1_X = d; P1_X = f; v = (u + v)*0.5f;
}
iterationStep++;
}
return ClampToZeroOne((u + v)*0.5f);
}
So, is there any hope for me to get this thing to work?
