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/art...rticle1497.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. You can also use and interpolate Quaternions instead.

By the way, here are two ]Linear Interpolation[/url] and Slerp.

[JSoftware]

i would recommend reading up on quaternions as they are really easy to use(not understand but whtt the heck) 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(dest,m1,m2: PMatrix; s: single);
var q,q1,q2: TVector4; m: TMatrix; i,j: integer;
begin
  q1:=MatrixToQuaternion3(m1^);
  q2:=MatrixToQuaternion3(m2^);
  q:=QuaternionSlerp(q1,q2,s);
  m:=QuaternionToMatrix3(q);
  dest^:=m;
  for i:=0 to 3 do // Interpolate position
    for j:=0 to 3 do
      if (i=3) or (j=3) then
        dest[i,j]:=m1[i,j]+s*(m2[i,j]-m1[i,j]);
end;

To optimize this even further, quaternions could be pre-calculated for each matrix.