Please explain this matrices interpolation code

Question

int rate,frac;
sint16 *frmptr[2];

frac = GetFrames(item, frmptr, &rate);

int GetFrames( ITEM_INFO *item, sint16 *frmptr[], int *rate )
{
    int first;
    int interp;

    // frm - current frame (for example 30 fps, its means one frame from 30)
    // rat - ratio for current animation for example 4

    //get pointer to current frame (start of position and rotation data)
    frmptr[0] = frmptr[1] = anim->frame_ptr;

    first =  frm/rat;
    interp = frm%rat;

    //take this frame
    frmptr[0] += first * frame_size;
    //take next frame
    frmptr[1] = frmptr[0] + frame_size;

    *rate = rat;

    //increase frame counter
    frm++;

    return interp;
    
}

Get position x1, y1, z1 and rotation xr1, yr1, zr1 from frmptr[0];

Get position x2, y2, z2 and rotation xr2, yr2, zr2 from frmptr[1];

Compose two matrices m1 and m2 correspond to position 1,2 and rotation 1,2.

Then interpolate two matrices. Matrices are colmajor.

sint32  *phd_mxptr; /* matrix pointer....*/
sint32  *IMptr;

enum msoff { M00,M01,M02,M03,
     M10,M11,M12,M13,
     M20,M21,M22,M23 };


void    InterpolateMatrix( void )
{
    sint32  *mptr,*iptr;
    sint32  frac,rate;

    //first matrix for interpolation m1
    mptr = phd_mxptr;

    //second matrix for interpolation m2
    iptr = IMptr;

    if ( rate==2 || (frac == 2 && rate == 4))                       // If interpolating Odd Frames
    {                                           // then do simple case
        *(mptr+M00) = (*(iptr+M00) + *(mptr+M00)) >> 1;
        *(mptr+M01) = (*(iptr+M01) + *(mptr+M01)) >> 1;
        *(mptr+M02) = (*(iptr+M02) + *(mptr+M02)) >> 1;
        *(mptr+M03) = (*(iptr+M03) + *(mptr+M03)) >> 1;

        *(mptr+M10) = (*(iptr+M10) + *(mptr+M10)) >> 1;
        *(mptr+M11) = (*(iptr+M11) + *(mptr+M11)) >> 1;
        *(mptr+M12) = (*(iptr+M12) + *(mptr+M12)) >> 1;
        *(mptr+M13) = (*(iptr+M13) + *(mptr+M13)) >> 1;

        *(mptr+M20) = (*(iptr+M20) + *(mptr+M20)) >> 1;
        *(mptr+M21) = (*(iptr+M21) + *(mptr+M21)) >> 1;
        *(mptr+M22) = (*(iptr+M22) + *(mptr+M22)) >> 1;
        *(mptr+M23) = (*(iptr+M23) + *(mptr+M23)) >> 1;
    }
    else if (frac == 1)
    {
        *(mptr+M00) += (*(iptr+M00) - *(mptr+M00)) >> 2;
        *(mptr+M01) += (*(iptr+M01) - *(mptr+M01)) >> 2;
        *(mptr+M02) += (*(iptr+M02) - *(mptr+M02)) >> 2;
        *(mptr+M03) += (*(iptr+M03) - *(mptr+M03)) >> 2;

        *(mptr+M10) += (*(iptr+M10) - *(mptr+M10)) >> 2;
        *(mptr+M11) += (*(iptr+M11) - *(mptr+M11)) >> 2;
        *(mptr+M12) += (*(iptr+M12) - *(mptr+M12)) >> 2;
        *(mptr+M13) += (*(iptr+M13) - *(mptr+M13)) >> 2;

        *(mptr+M20) += (*(iptr+M20) - *(mptr+M20)) >> 2;
        *(mptr+M21) += (*(iptr+M21) - *(mptr+M21)) >> 2;
        *(mptr+M22) += (*(iptr+M22) - *(mptr+M22)) >> 2;
        *(mptr+M23) += (*(iptr+M23) - *(mptr+M23)) >> 2;
    }
    else
    {
        *(mptr+M00) = *(iptr+M00) - ((*(iptr+M00) - *(mptr+M00)) >> 2);
        *(mptr+M01) = *(iptr+M01) - ((*(iptr+M01) - *(mptr+M01)) >> 2);
        *(mptr+M02) = *(iptr+M02) - ((*(iptr+M02) - *(mptr+M02)) >> 2);
        *(mptr+M03) = *(iptr+M03) - ((*(iptr+M03) - *(mptr+M03)) >> 2);

        *(mptr+M10) = *(iptr+M10) - ((*(iptr+M10) - *(mptr+M10)) >> 2);
        *(mptr+M11) = *(iptr+M11) - ((*(iptr+M11) - *(mptr+M11)) >> 2);
        *(mptr+M12) = *(iptr+M12) - ((*(iptr+M12) - *(mptr+M12)) >> 2);
        *(mptr+M13) = *(iptr+M13) - ((*(iptr+M13) - *(mptr+M13)) >> 2);

        *(mptr+M20) = *(iptr+M20) - ((*(iptr+M20) - *(mptr+M20)) >> 2);
        *(mptr+M21) = *(iptr+M21) - ((*(iptr+M21) - *(mptr+M21)) >> 2);
        *(mptr+M22) = *(iptr+M22) - ((*(iptr+M22) - *(mptr+M22)) >> 2);
        *(mptr+M23) = *(iptr+M23) - ((*(iptr+M23) - *(mptr+M23)) >> 2);
    }
}

In this example used fixed point math.

I dont understand process interpolation of matrices - 3 cases from InterpolateMatrix() function. Please explain.

Answer 1

This is more a less a component-wise interpolation of 2 matrices of integers using some binary hacks to optimize divide by 2 and divide by 4. You can rewrite the 3 cases like this:

1. (A + B) >> 1 => (1/2 * A) + (1/2 * B)
2. B + (A - B) >> 2 => (1/4 * A) + (3/4 * B)
3. A - (A - B) >> 2 => (3/4 * A) + (1/4 * B)

So it's splitting the interpolation in 3 cases each 25% along the linear interp of the 2 matrices.

From the looks of it this an attempt to inerpolate position/rotation matrices, and while the positional component will work relatively fine this is not a correct way to interpolate rotation matrices. It might be good enough for an extremely performance-critical application though.

To do it correctly the easiest way is to convert the rotation to a quaternion and then using a slerp operation.

Answer 2

Shift on 1 = division on 2 (i.e. 0.5). Shift on 2 is division on 4 (i.e. 0.25)

Key is this part of code:

    first =  frm/rat;
    interp = frm%rat;

For example we do count the frame each time. And now number frame is 8 (totaly could be 30 fps - frame per second). And Frame Rate = 4.

For example currect frame is 8. i.e. frm = 8. Then:

int first = 8/4 = 2;
int interp = 8%4 = 0;

For example the current frame frm= 9:

int first = 9/4 = 2;
int interp = 9%4 = 1;

And next frame is 10:

int first = 10/4 = 2;
int interp = 10%4 = 2;

And next frame is 11:

int first = 10/4 = 2;
int interp = 10%4 = 3;

And for example next frame is 12:

int first = 12/4 = 3;
int interp = 12%4 = 0;

Four frames first variable stays 2 and four frames interp variable is 0,1,3.

The first case:

if ( rate==2 || (frac == 2 && rate == 4))                

       // If interpolating Odd Frames
    {  

If frac == 2 i.e. first 8 + interp 2 = 10 - it's half of frame (Frame Rate = 4 i.e. half of frame = 2) then we do interpolation:

m1 = (m2 + m1) * 0.5;

Second case:

else if (frac == 1)
{
m1 = m1 + (m2 - m1) * 0.25;

if frac == 1 i.e. 1 of 4 = 0.25, i.e. 1/4 = 0.25:

And third case obviously frac == 3 i.e. 3 (third) part of 4 (four) = 3/4 = 0.75, i.e. from m2 - 0.25 of(m2- m1):

else // if (t == 0.75)
  {
    m1 = m2 - (m2 - m1) * 0.25; // = lerp(m1, m2, 0.75)
  }

As result we can do so:

if(frac == 2) //i.e. 0.5 (2 of 4)
{
    //matrix interpolation
}
else if(frac == 1) //i.e. 0.25 (1 of 4)
{
    //matrix interpolation
}
else // if(frac == 3) i.e. 0.75 (3 of 4)
{
    //matrix interpolation
}

I.e. it is my answer:

if (t == 0.5)
  {
    m1 = (m2 + m1) * 0.5; // = lerp(m1, m2, 0.5)
  }
  else if (t == 0.25)
  {
    m1 = m1 + (m2 - m1) * 0.25; // = lerp(m1, m2, 0.25)
  }
  else // if (t == 0.75)
  {
    m1 = m2 - (m2 - m1) * 0.25; // = lerp(m1, m2, 0.75)
  }

Shorterner case of code above in my question is:

matrix4x4 m1, m2;

    if ( rate==2 || (frac == 2 && rate == 4))
    {
        m1 = (m2 + m1) >> 1;
    }
    else if (frac == 1)
    {
        m1 = m1 + (m2 - m1) >> 2;   
    }
    else
    {
        m1 = m2 - (m2 - m1) >> 2;
    }