Blending matrices

[User137]

I'm not sure if this is called "Interpolating" or something, situation is that we have matrices 1 and 2 (position and rotation of 3D models), then we give these 2 matrix and value between 0..1 to a function that gives third matrix that is something between these 2. Very handy tool for animating stuff.

This is what i have done so far, and it works well only not looking good when 2 objects are more than 45 degrees in different angle.

procedure MatrixBlend(dest,m1,m2: PMatrix; s: single);
var i,j: integer;
begin
for i:=0 to 3 do
for j:=0 to 3 do
dest[i,j]:=m1[i,j]+s*(m2[i,j]-m1[i,j]);
for i:=0 to 2 do
Normalize(dest[i,0],dest[i,1],dest[i,2]);
end;

Is there another kind of way to get steps between 2 animation frames?

[dmantione]

Well, in a lot of situations you would calculate the difference between 2 steps instead of calculating an offset from the start, so you dont need to multiply for each intermediate step.

[User137]

Can you explain a bit more about it? I don't find anything ineffective in current method, the result fron the function can be used straight away by opengl and is very fast. Only lacks the quality...

There may be some way with angles (slow perhaps) or i heard of some quaternions (i know nothing about)... or something else.

Here's some reference:
I have identified my function as linear interpolation
http://www.gamedev.net/reference/articles/article1497.asp

[dmantione]

Well, suppose you have s=0.1, so you have the start + 10 steps, in total 11 positions. Then you would call your code with 0.0, 0.1, 0.2, ...

Now look at the code:


procedure MatrixBlend(dest,m1,m2: PMatrix; s: single);
var i,j: integer;
begin
for i:=0 to 3 do
for j:=0 to 3 do
dest[i,j]:=m1[i,j]+s*(m2[i,j]-m1[i,j]);
for i:=0 to 2 do
Normalize(dest[i,0],dest[i,1],dest[i,2]);
end;


The first time you call it, you increase m1 with 0.1*(m2-m1), the second time with 0.2*(m2-m1), etc.

So, the 11 iterations would be:

output1=m1+0.0*(m2-m1)
output2=m1+0.1*(m2-m1)
output3=m1+0.2*(m2-m1)
output4=m1+0.3*(m2-m1)
output5=m1+0.4*(m2-m1)
output6=m1+0.5*(m2-m1)
output7=m1+0.6*(m2-m1)
output8=m1+0.7*(m2-m1)
output9=m1+0.8*(m2-m1)
output10=m1+0.9*(m2-m1)
output11=m1+1.0*(m2-m1)

That can be simplified:
step=0.1*(m2-m1)
output1=m1
output2=output1+step
output3=output2+step
output4=output3+step
output5=output4+step
output6=output5+step
output7=output6+step
output8=output7+step
output9=output8+step
output10=output9+step
output11=output10+step

In other words, you can save a lot of calculations. The only problem could be the normalize step; some mathematical formulas might require you to use normalized vectors. In most situations this is not a problem.

[LP]

There may be some way with angles (slow perhaps) or i heard of some quaternions (i know nothing about)... or something else.
Indeed you won't get the right result by interpolating matrix values directly. You need to interpolate smoothly translation and rotation, which requires of separating your matrix into position and euler angles and then interpolating them separately. Look more on this here (http://skal.planet-d.net/demo/matrixfaq.htm#Q44). You can also use and interpolate Quaternions (http://mathworld.wolfram.com/Quaternion.html) instead.

By the way, here are two ]Linear Interpolation[/url] and Slerp (http://en.wikipedia.org/wiki/Slerp).

[JSoftware]

i would recommend reading up on quaternions as they are really easy to use(not understand but whtt the heck) :wink: they have the slerp ability which is a very neat thing

[Nitrogen]

Yes, Quaternions are the standard method to use when it comes to animating rotations. You can get hold of a Matrix to Quaternion conversion method on the web and go from there...

[User137]

I've looked on some quaternion links + all the above and must say this will be hard.. very hard, if i ever get it done. That's why i was mainly looking for ready code.

But thanks anyway.

[LP]

[quote="User137"]I've looked on some quaternion ]
Straight from AsphyreMath.pas (this was written by Soulhab, so comments are a bit weird ;)):



//--------------------------------------------------------------------------
const
AsphyreEpsilon = 0.0000001;

//--------------------------------------------------------------------------
function QuaternionToMatrix3(const q: TVector4): TMatrix3;
begin
Result[0, 0]:= 1.0 - (2.0 * q.y * q.y) - (2.0 * q.z * q.z);
Result[0, 1]:= (2.0 * q.x * q.y) + (2.0 * q.w * q.z);
Result[0, 2]:= (2.0 * q.x * q.z) - (2.0 * q.w * q.y);
Result[1, 0]:= (2.0 * q.x * q.y) - (2.0 * q.w * q.z);
Result[1, 1]:= 1.0 - (2.0 * q.x * q.x) - (2.0 * q.z * q.z);
Result[1, 2]:= (2.0 * q.y * q.z) + (2.0 * q.w * q.x);
Result[2, 0]:= (2.0 * q.x * q.z) + (2.0 * q.w * q.y);
Result[2, 1]:= (2.0 * q.y * q.z) - (2.0 * q.w * q.x);
Result[2, 2]:= 1.0 - (2.0 * q.x * q.x) - (2.0 * q.y * q.y);
end;

//--------------------------------------------------------------------------
function MatrixToQuaternion3(const Matrix: TMatrix3): TVector4;
var
Aux : TVector4;
Max : Single;
Index: Integer;
High : Double;
Mult : Double;
begin
// Determine which of w, x, y, z has the largest absolute value.
Aux.w:= Matrix[0, 0] + Matrix[1, 1] + Matrix[2, 2];
Aux.x:= Matrix[0, 0] - Matrix[1, 1] - Matrix[2, 2];
Aux.y:= Matrix[1, 1] - Matrix[0, 0] - Matrix[2, 2];
Aux.z:= Matrix[2, 2] - Matrix[0, 0] - Matrix[1, 1];

Index:= 0;
Max := Aux.w;
if (Aux.x > Max) then
begin
Max := Aux.x;
Index:= 1;
end;
if (Aux.y > Max) then
begin
Max := Aux.y;
Index:= 2;
end;
if (Aux.z > Max) then
begin
Max := Aux.z;
Index:= 3;
end;

// Perform square root and division.
High:= Sqrt(Max + 1.0) * 0.5;
Mult:= 0.25 / High;

// Apply table to compute quaternion values.
case Index of
0: begin
Result.w:= High;
Result.x:= (Matrix[1, 2] - Matrix[2, 1]) * Mult;
Result.y:= (Matrix[2, 0] - Matrix[0, 2]) * Mult;
Result.z:= (Matrix[0, 1] - Matrix[1, 0]) * Mult;
end;
1: begin
Result.x:= High;
Result.w:= (Matrix[1, 2] - Matrix[2, 1]) * Mult;
Result.z:= (Matrix[2, 0] + Matrix[0, 2]) * Mult;
Result.y:= (Matrix[0, 1] + Matrix[1, 0]) * Mult;
end;
2: begin
Result.y:= High;
Result.z:= (Matrix[1, 2] + Matrix[2, 1]) * Mult;
Result.w:= (Matrix[2, 0] - Matrix[0, 2]) * Mult;
Result.x:= (Matrix[0, 1] + Matrix[1, 0]) * Mult;
end;
else
begin
Result.z:= High;
Result.y:= (Matrix[1, 2] + Matrix[2, 1]) * Mult;
Result.x:= (Matrix[2, 0] + Matrix[0, 2]) * Mult;
Result.w:= (Matrix[0, 1] - Matrix[1, 0]) * Mult;
end;
end;
end;

//--------------------------------------------------------------------------
function QuaternionSlerp(const q1, q2: TVector4; Theta: Single): TVector4;
var
CosOmega, SinOmega: Real;
q: TVector4;
Omega: Real;
Coef0, Coef1: Single;
begin
// Compute the "cosine of the angle" between the quaternions,
// using the dot product.
CosOmega:= q1.x * q2.x + q1.y * q2.y + q1.z * q2.z + q1.w * q2.w;

// If negative dot, negate one of the input
// quaterions to take the shorter 4D "arc".
q:= q2;
if &#40;CosOmega < 0&#41; then
begin
CosOmega&#58;= -CosOmega;
q.x&#58;= -q2.x;
q.y&#58;= -q2.y;
q.z&#58;= -q2.z;
q.w&#58;= -q2.w;
end;

// Check if they are very close together to protect against divide-by-zero.
Coef0&#58;= 1.0 - Theta;
Coef1&#58;= Theta;
if &#40;1.0 - CosOmega > AsphyreEpsilon&#41; then
begin
// Spherical interpolation.
Omega&#58;= ArcCos&#40;CosOmega&#41;;
SinOmega&#58;= Sin&#40;Omega&#41;;
Coef0&#58;= Sin&#40;&#40;1.0 - Theta&#41; * Omega&#41; / SinOmega;
Coef1&#58;= Sin&#40;Theta * Omega&#41; / SinOmega;
end;

// Interpolate.
Result.x&#58;= &#40;Coef0 * q1.x&#41; + &#40;Coef1 * q.x&#41;;
Result.y&#58;= &#40;Coef0 * q1.y&#41; + &#40;Coef1 * q.y&#41;;
Result.z&#58;= &#40;Coef0 * q1.z&#41; + &#40;Coef1 * q.z&#41;;
Result.w&#58;= &#40;Coef0 * q1.w&#41; + &#40;Coef1 * q.w&#41;;
end;

[User137]

Cool! Thank you very much, it worked. Also learned that quaternions cover only rotation. Anyways, here is current one that does not bug visually:


procedure MatrixBlend3&#40;dest,m1,m2&#58; PMatrix; s&#58; single&#41;;
var q,q1,q2&#58; TVector4; m&#58; TMatrix; i,j&#58; integer;
begin
q1&#58;=MatrixToQuaternion3&#40;m1^&#41;;
q2&#58;=MatrixToQuaternion3&#40;m2^&#41;;
q&#58;=QuaternionSlerp&#40;q1,q2,s&#41;;
m&#58;=QuaternionToMatrix3&#40;q&#41;;
dest^&#58;=m;
for i&#58;=0 to 3 do // Interpolate position
for j&#58;=0 to 3 do
if &#40;i=3&#41; or &#40;j=3&#41; then
dest&#91;i,j&#93;&#58;=m1&#91;i,j&#93;+s*&#40;m2&#91;i,j&#93;-m1&#91;i,j&#93;&#41;;
end;
To optimize this even further, quaternions could be pre-calculated for each matrix.

[Paulius]

Only the last row defines position, lineary interpolating w components makes no sense

[JSoftware]

scaling?

[User137]

Well, it is not just interpolating w components but resetting those values to the matrix. QuaternionToMatrix3 only touches first 3x3 part. If i replace 0..3 to 0..2 it doesn't work, either copying only other i or j row doesn't work. Must have the full matrix initialized.

This will get cleaner when i edit those functions using pointers instead of slowly taking and returning full matrixes and vectors.

[LP]

This will get cleaner when i edit those functions using pointers instead of slowly taking and returning full matrixes and vectors.
You can do that, but you lose elegancy *and* usability. You won't be able to do something like that:


// Delphi 2005 or lower code
RawMtx&#58;= MatMul&#40;MatScale&#40;Vector3&#40;0.5, 0.5, 0.5&#41;&#41;, MatMul&#40;MatTransl&#40;Vector3&#40;100.0, 50.0, 0.0&#41;&#41;, MatRotateX&#40;Pi / 2&#41;&#41;&#41;;

// Delphi 2006
RawMtx&#58;= Vector3&#40;0.5, 0.5, 0.5&#41; * &#40;MatTransl&#40;Vector3&#40;100.0, 50.0, 0.0&#41;&#41; * MatRotateX&#40;Pi / 2&#41;&#41;;


The only real optimization you can do is replacing return values by pointers, but as I said, you won't be able to use functions as parameters to other functions. You won't be able to optimize parameters much, since they are already passed as pointers (notice "const" in declaration).

[XProger]

Lifepower,
Delphi 2005-2006 have inline functions :)

[User137]

I don't mind usability as long as it needs short code and is faster.

This may be final version:
procedure MatrixBlend(dest,m1,m2: PMatrix; s: single);
var q,q1,q2: TVector4; i: integer;
begin
// Initialize matrix
for i:=0 to 3 do dest[i,3]:=m1[i,3];
// Interpolate rotation
MatrixToQuat(m1,@q1);
MatrixToQuat(m2,@q2);
QuatSlerp(@q1,@q2,@q,s);
QuatToMatrix(@q,dest);
for i:=0 to 2 do // Interpolate position
dest[3,i]:=m1[3,i]+s*(m2[3,i]-m1[3,i]);
end;

[LP]

Lifepower,
Delphi 2005-2006 have inline functions :)
I know. :)

By the way, I've done some tests with Delphi inline functions and most of the time performance was actually dropping when using inline functions. Also, I'm not quite fond of the code generation when it comes to inline functions (somehow inline code is treated as separate in code optimization). They should have implemented macro support instead.