How to compare two curves (array of points)

Question

I have problem to find method to compare two trajectories (curves). The first original contains points (x,y). The second one can be offset, smaller or larger scale, and with rotation - also array with points (x,y)

My first method that i did is to find smallest distance between two points and repeat this process in every iteration, sum of it and divide by number of points - then my result tell me value the average error per point: http://www.mathopenref.com/coorddist.html

And also i find this method: https://help.scilab.org/docs/6.0.0/en_US/fminsearch.html

But i cant figure out how to use it. I would like compare both trajectories but my results have to include rotation, or at least offset for beginning.

My current result is calculate error per point (distance)

1. get coordinate (x,y) second trajectory.
2. in loop i try to find min_distance between (x,y) from 1. and point from original trajectory.
3. add smallest_distance what i found in 2 step.
4. divide sum of smallest distance by number of points from second trajectory.

My result describe average error(distance) per points if we compare with original trajectory.

But i can not figure how to handle if trajectory is rotated, scaled or is shifted.

Please look at my example trajectories:

1. http://pokazywarka.pl/trajectory/

2. http://pokazywarka.pl/trajectory2/

Answer 1

So you need to compare shape of 2 curves invariant on rotation,translation and scale.

Solution

Let assume 2 sinwaves for testing. Both rotated and scaled but with the same aspect ratio and one with added noise. I generated them in C++ like this:

struct _pnt2D
    {
    double x,y;
    // inline
    _pnt2D()    {}
    _pnt2D(_pnt2D& a)   { *this=a; }
    ~_pnt2D()   {}
    _pnt2D* operator = (const _pnt2D *a) { *this=*a; return this; }
    //_pnt2D* operator = (const _pnt2D &a) { ...copy... return this; }
    };
List<_pnt2D> curve0,curve1;         // curves points
_pnt2D p0,u0,v0,p1,u1,v1;           // curves OBBs

const double deg=M_PI/180.0;
const double rad=180.0/M_PI;
void rotate2D(double alfa,double x0,double y0,double &x,double &y)
    {
    double   a=x-x0,b=y-y0,c,s;
    c=cos(alfa);
    s=sin(alfa);
    x=x0+a*c-b*s;
    y=y0+a*s+b*c;
    }

// this code is the init stuff:
int i;
double x,y,a;
_pnt2D p,*pp;
Randomize();
for (x=0;x<2.0*M_PI;x+=0.01)
    {
    y=sin(x);

    p.x= 50.0+(100.0*x);
    p.y=180.0-( 50.0*y);
    rotate2D(+15.0*deg,200,180,p.x,p.y);
    curve0.add(p);

    p.x=150.0+( 50.0*x);
    p.y=200.0-( 25.0*y)+5.0*Random();
    rotate2D(-25.0*deg,250,100,p.x,p.y);
    curve1.add(p);
    }

1. OBB oriented bounding box

   compute OBB which will find the rotation angle and position of both curves so rotate one of them so they start at the same position and has the same orientation.

   If the OBB sizes are too different then the curves are different.

   For above example it yealds this result:

   

   Each OBB is defined by start point P and basis vectors U,V where |U|>=|V| and z coordinate of U x V is positive. That will ensure the same winding for all OBBs. It can be done in OBBox_compute by adding this to the end:

   // |U|>=|V|
   if ((u.x*u.x)+(u.y*u.y)<(v.x*v.x)+(v.y*v.y)) { _pnt2D p; p=u; u=v; v=p; }
   // (U x V).z > 0
   if ((u.x*v.y)-(u.y*v.x)<0.0)
       {
       p0.x+=v.x;
       p0.y+=v.y;
       v.x=-v.x;
       v.y=-v.y;
       }
   

   So curve0 has p0,u0,v0 and curve1 has p1,u1,v1.

   Now we want to rescale,translate and rotate curve1 to match curve0 It can be done like this:

   // compute OBB
   OBBox_compute(p0,u0,v0,curve0.dat,curve0.num);
   OBBox_compute(p1,u1,v1,curve1.dat,curve1.num);
   // difference angle = - acos((U0.U1)/(|U0|.|U1|))
   a=-acos(((u0.x*u1.x)+(u0.y*u1.y))/(sqrt((u0.x*u0.x)+(u0.y*u0.y))*sqrt((u1.x*u1.x)+(u1.y*u1.y))));
   // rotate curve1
   for (pp=curve1.dat,i=0;i<curve1.num;i++,pp++)
    rotate2D(a,p1.x,p1.y,pp->x,pp->y);
   // rotate OBB1
   rotate2D(a,0.0,0.0,u1.x,u1.y);
   rotate2D(a,0.0,0.0,v1.x,v1.y);
   // translation difference = P0-P1
   x=p0.x-p1.x;
   y=p0.y-p1.y;
   // translate curve1
   for (pp=curve1.dat,i=0;i<curve1.num;i++,pp++)
       {
       pp->x+=x;
       pp->y+=y;
       }
   // translate OBB1
   p1.x+=x;
   p1.y+=y;
   // scale difference = |P0|/|P1|
   x=sqrt((u0.x*u0.x)+(u0.y*u0.y))/sqrt((u1.x*u1.x)+(u1.y*u1.y));
   // scale curve1
   for (pp=curve1.dat,i=0;i<curve1.num;i++,pp++)
       {
       pp->x=((pp->x-p0.x)*x)+p0.x;
       pp->y=((pp->y-p0.y)*x)+p0.y;
       }
   // scale OBB1
   u1.x*=x;
   u1.y*=x;
   v1.x*=x;
   v1.y*=x;
   

   You can use Understanding 4x4 homogenous transform matrices to do all this in one step. Here the result:

   

2. sampling

   in case of non uniform or very different point density between curves or between any parts of it you should re-sample your curves to have common point density. You can use linear or polynomial interpolation for this. You also do not need to store the new sampling in memory but instead you could build function that returns point of each curve parametrized by arc-length from start.

   point curve0(double distance);
   point curve1(double distance);
   

3. comparison

   Now you can substract the 2 curves and sum up the abs of the differences. Then divide it by the curve length and threshold the result.

   for (double sum=0.0,l=0.0;d<=bigger_curve_length;l+=step)
    sum+=fabs(curve0(l)-curve1(l));
   sum/=bigger_curve_length;
   if (sum>threshold) curves are different
    else curves match
   

You should try this even with +180deg rotation as the orientation difference from OBB has only half of the true range.

Here few related QAs:

• compare shapes
• How can i produce multi point linear interpolation?