Piecewise linear integer curve interpolation in C#/Unity3D

Solution 1:

I would use this interpolation cubic:

x=a0+a1*t+a2*t*t+a3*t*t*t
y=b0+b1*t+b2*t*t+b3*t*t*t

where a0..a3 are computed like this:

d1=0.5*(p2.x-p0.x);
d2=0.5*(p3.x-p1.x);
a0=p1.x;
a1=d1;
a2=(3.0*(p2.x-p1.x))-(2.0*d1)-d2;
a3=d1+d2+(2.0*(-p2.x+p1.x));

b0 .. b3 are computed in same way but use y coordinates of course
p0..p3 are control points for cubic interpolation curve
t = < 0.0 , 1.0 > is curve parameter from p1 to p2

This ensures that position and first derivation is continuous (c1). If you want to do this on integer math then just scale ai,bi ant t accordingly. You can also add as many dimensions as you need in the same manner

Now you need some parameter to go through your interpolation points for example u = <0 , N-1>

p(0..N-1) are your control points list
u = 0 means start point p(0)
u = N-1 means end point p(N-1)
P0..P3 are control points used for interpolation

So you need to compute t and select which points to use for interpolation

    double t=u-floor(u); // fractional part between control points
    int i=floor(u);       // integer part points to starting control point used
         if (i<1)     { P0=p(  0),P1=p(  0),P2=p(  1),P3=p(  2); }               // handle start edge case
    else if (i==N-1) { P0=p(N-2),P1=p(N-1),P2=p(N-1),P3=p(N-1); }  // handle end edge case
    else if (i>=N-2) { P0=p(N-3),P1=p(N-2),P2=p(N-1),P3=p(N-1); }  // handle end edge case
    else              { P0=p(i-1),P1=p(i  ),P2=p(i+1),P3=p(i+2); }

    (x,y) = interpolation (P0,P1,P2,P3,t);

If you want to do this on integer math then just scale u,t accordingly. If N<3 then use linear interpolation ... or duplicate end points until N>=3

[edit1] linear interpolation approach

struct pnt { int x,y; };

pnt interpolate (pnt *p,int N,int x)
    {
    int i,j;
    pnt p;
    for (j=1,i=N-1;j<i;j<<=1); j>>=1; if (!j) j=1; // this just determine max mask for binary search ... can do it on p[] size change
    for (i=0;j;j>>=1) // binary search by x coordinate output is i as point index with  p[i].x<=x
        {
        i|=j;
        if (i>=N) { i-=j; continue; }
        if (p[i].x==x) break;
        if (p[i].x> x) i-=j;
        }
    p.x=x;
    p.y=p[i].y+((p[i+1].y-p[i].y)*(x-p[i].x)/(p[i+1].x-p[i].x))
    return p;
    }

add edge cases handling like x is out of points bound or point list is too small