Is it pointcloud or polyline?
Is the sampling,scale the same for both curves?
You can use this How to compare two shapes? and find parts where the sequence of angles is the "same". this is invariant on rotation,scaling and translation. So it should lead to direct match which points from one curve match to which points on the other in order.
In case your polylines do not match exactly you can use correlation coefficient on the angles or threshold the sum of differences of the angles/number of angles.
In case your order of points might be in reverse between the polylines you need to check both combinations...
After that its just a matter of computing:
scale
just by dividing sizes between 2 distant points (forming big line) on both curves
rotation
use atan2 on the same lines from #1 and substract the angles
translation
just substract single point from both curves (but after applying scale and rotation on the point belonging to curve which is being aligned).
Or you can use transform matrix approach instead.
Here small simple C++/GL/VCL example of brute force search:
//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#pragma hdrstop
#include "Unit1.h"
#include "gl_simple.h"
#include "OpenGLrep4d_double.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
#define polysep -1.1e10,-1.1e10
#define polyend +1.1e10,+1.1e10
double pol0[]=
{
3237,253 ,3652,425 ,polysep,2272,258 ,2384,409 ,2560,545 ,2730,637 ,2837,814 ,
2882,942 ,2891,1069,2837,1221,2770,1323,2597,1462,2563,1537,2566,1590,2591,1645,
2686,1721,2778,1740,3013,1707,3102,1709,3276,1887,3648,2127,3722,2212,3764,2295,
3869,2437,3950,2515,3995,2538,4166,2559,4484,2280,4534,2255,4555,2260,4593,2382,
4664,2491,4703,2671,4738,2754,4811,2862,4906,2935,4991,2942,5196,2877,5342,2869,
5545,2782,5715,2808,5869,2810,5992,2832,6105,2891,6369,3082,6583,3116,6526,3203,
6465,3242,6366,3340,6192,3414,5922,3457,5732,3576,5732,3618,5815,3727,5871,3833,
5968,3921,6029,3944,6034,4000,6015,4043,5687,4191,4672,4806,4562,4838,3258,5475,
3166,6035,3135,6398,3152,6425,polysep,6584,3115,6745,3059,6802,3054,6979,3053,
7129,3089,polysep,7344,3139,7299,3158,7269,3150,polysep,6900,3306,6756,3397,
6756,3436,polysep,7378,3441,7345,3474,7237,3513,polysep,7205,3512,7168,3544,
7079,3578,7060,3603,polysep,7269,3580,7216,3601,7180,3701,polysep,7190,3649,
7057,3648,6984,3678,polysep,7040,3578,7019,3607,6915,3653,6899,3679,polysep,
6876,3658,6857,3682,6741,3734,6726,3758,polysep,6293,3755,6186,3833,6042,3901,
polysep,6551,3798,6536,3822,6400,3893,polysep,6539,3862,6442,4793,polysep,
6372,3879,6217,3970,polysep,6200,3948,6056,4046,polysep,6036,4078,6013,4615,
polysep,7294,3867,7245,4069,7237,4346,polyend
};
double pol1[]=
{
235,1853 ,485,2032 ,600,2091 ,682,2183 ,813,2414 ,914,2493 ,1051,2513,1390,2214,
1421,2201,1495,2222,1522,2261,1541,2363,1609,2460,1680,2718,1750,2827,1830,2884,
1917,2881,2091,2825,2223,2819,2399,2737,2489,2732,2895,2782,3080,2878,3193,2973,
3301,3032,3426,3042,3614,3000,3890,2980,4027,3029,4073,3062,4089,3097,4082,3168,
3788,3130,3660,3139,3548,3195,3409,3299,3307,3404,3112,3489,2860,3520,2732,3590,
2720,3611,2728,3641,2796,3772,2907,3871,2903,3900,polysep,2908,3871,2968,3882,
2984,3928,2991,3992,2954,4087,2638,4224,2203,4503,1588,4863,1479,4895,1184,5035,
252,5489 ,235,5518 ,polysep,235,1912 ,257,1964 ,228,2035 ,polysep,258,1962 ,
332,1963 ,polysep,2152,4342,2038,4436,1872,4528,polysep,1795,4616,1818,4568,
1881,4574,polysep,1441,4864,1303,4904,927,5098 ,831,5141 ,767,5152 ,polysep,
597,5017 ,669,4938 ,692,4942 ,polysep,838,4944 ,802,4981 ,855,5057 ,904,5054 ,
polysep,843,5043 ,786,5062 ,795,5084 ,polysep,784,5063 ,732,5034 ,641,5086 ,
637,5132 ,polysep,640,5086 ,581,5084 ,550,5107 ,585,5207 ,polysep,742,5031 ,
760,4997 ,799,4981 ,polysep,1554,4604,1427,4654,polysep,1140,4622,1133,4647,
1287,4669,1428,4634,polyend
};
const int pnts0=(sizeof(pol0)/sizeof(pol0[0]))>>1;
const int pnts1=(sizeof(pol1)/sizeof(pol1[0]))>>1;
double da0[pnts0]; // delta angles of pol0
double da1[pnts1]; // delta angles of pol1
double M[16]; // align pol1 to pol0 matrix so pol0=M*pol1
double V[16]; // just GL view matrix so pol0,pol1 are fully visible
int a0=20,a1=30,b0=20,b1=30; // match
//---------------------------------------------------------------------------
void compute_MV()
{
//--------------------------------------------------------
int i,j,k;
double x0,y0,x1,y1,x,y;
double x2,y2,x3,y3;
//--------------------------------------------------------
/*
// reverse po1[] in case your polylines are in reverse in between pol0/pol1
for (i=0,j=pnts1+pnts1-4;i<j;i+=2,j-=2)
{
x=pol1[i+0]; pol1[i+0]=pol1[j+0]; pol1[j+0]=x;
y=pol1[i+1]; pol1[i+1]=pol1[j+1]; pol1[j+1]=y;
}
// rotate po1[] to better test rotation correction
x=5.0*deg; x0=cos(x); y0=sin(x);
for (i=0;i<pnts1+pnts1;i+=2)
{
x=pol1[i+0];
y=pol1[i+1];
pol1[i+0]=+(x*x0)+(y*y0);
pol1[i+1]=-(x*y0)+(y*x0);
}
*/
//--------------------------------------------------------
// M,V set to unit matrices
for (i=0;i<16;i++) M[i]=0.0; for (i=0;i<16;i+=5) M[i]=1.0;
for (i=0;i<16;i++) V[i]=0.0; for (i=0;i<16;i+=5) V[i]=1.0;
// pol0 AABB
x0=y0=+1e10;
x1=y1=-1e10;
for (i=0;;)
{
x=pol0[i]; i++;
y=pol0[i]; i++;
if (x>+1e10) break;
if (x<-1e10) continue;
if (x0>x) x0=x;
if (x1<x) x1=x;
if (y0>y) y0=y;
if (y1<y) y1=y;
}
// pol0+pol1 AABB
for (i=0;;)
{
x=pol1[i]; i++;
y=pol1[i]; i++;
if (x>+1e10) break;
if (x<-1e10) continue;
if (x0>x) x0=x;
if (x1<x) x1=x;
if (y0>y) y0=y;
if (y1<y) y1=y;
}
// view V scale
V[0] =+1.8/(x1-x0);
V[5] =-1.8/(y1-y0);
V[10]=+1.0;
V[15]= 1.0;
// view V offset
V[12]=-0.5*(x1+x0)*V[0];
V[13]=-0.5*(y1+y0)*V[5];
//--------------------------------------------------------
// compute a0
for (i=0,j=0,k=0;;k++)
{
da0[k]=0.0;
x0=x1; x1=pol0[i]; i++;
y0=y1; y1=pol0[i]; i++;
if (!j){ x0=x1; y0=y1; } j++;
if (x1>+1e10){ da0[k]=1000.0; break; }
if (x1<-1e10){ da0[k]=1000.0; j=0; continue; }
da0[k]=atanxy(x1-x0,y1-y0);
}
// compute da0
for (x1=0.0,i=0;i<pnts0;i++)
{
x0=x1; x1=da0[i]; i++;
if (x1>999.0){ x1=0.0; continue; }
x=x1-x0;
while (x>+pi) x-=pi2;
while (x<-pi) x+=pi2;
da0[i]=x;
}
//--------------------------------------------------------
// compute a1
for (i=0,j=0,k=0;;k++)
{
da1[k]=0.0;
x0=x1; x1=pol1[i]; i++;
y0=y1; y1=pol1[i]; i++;
if (!j){ x0=x1; y0=y1; } j++;
if (x1>+1e10){ da1[k]=1000.0; break; }
if (x1<-1e10){ da1[k]=1000.0; j=0; continue; }
da1[k]=atanxy(x1-x0,y1-y0);
}
// compute da1
for (x1=0.0,i=0;i<pnts1;i++)
{
x0=x1; x1=da1[i]; i++;
if (x1>999.0){ x1=0.0; continue; }
x=x1-x0;
while (x>+pi) x-=pi2;
while (x<-pi) x+=pi2;
da1[i]=x;
}
//--------------------------------------------------------
// brute force window search ignoring separators
int n=15; // search window size
a0=-1; b0=-1; x1=-1.0;
for (i=0;i<pnts0-n;i++)
for (j=0;j<pnts1-n;j++)
{
for (x0=0.0,k=0;k<n;k++)
{
y0=da0[i+k]; if (y0>999.0) break;
y1=da1[j+k]; if (y1>999.0) break;
x0+=fabs(y1-y0);
}
if (k<n) continue; // skip shorter polylines
x0/=n;
if ((a0<0)||(x1>x0)){ x1=x0; a0=i; b0=j; }
}
a1=a0+n; b1=b0+n;
//--------------------------------------------------------
// line0
i=a0+a0;
x0=pol0[i]; i++;
y0=pol0[i]; i++;
i=a1+a1;
x1=pol0[i]; i++;
y1=pol0[i]; i++;
// line1
i=b0+b0;
x2=pol1[i]; i++;
y2=pol1[i]; i++;
i=b1+b1;
x3=pol1[i]; i++;
y3=pol1[i]; i++;
// trasform
x=sqrt(((x1-x0)*(x1-x0))+((y1-y0)*(y1-y0)));
y=sqrt(((x3-x2)*(x3-x2))+((y3-y2)*(y3-y2)));
double scale=x/y;
double ang=atanxy(x1-x0,y1-y0)-atanxy(x3-x2,y3-y2);
M[0]=+scale*cos(ang);
M[1]=+scale*sin(ang);
M[4]=-scale*sin(ang);
M[5]=+scale*cos(ang);
// offset (x,y,?) = M*(x2,y2,0.0,1.0)
x=(M[0]*x2)+(M[4]*y2);
y=(M[1]*x2)+(M[5]*y2);
M[12]=x0-x;
M[13]=y0-y;
}
//---------------------------------------------------------------------------
void gl_draw() // main rendering code
{
int i,j;
double x,y,r;
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glDisable(GL_CULL_FACE);
glDisable(GL_DEPTH_TEST);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glLoadMatrixd(V);
// render pol0 (purple)
glBegin(GL_LINE_STRIP);
for (i=0;;)
{
if (i==a0+a0) glColor3f(0.2,0.4,0.8);
if ((!i)||(i==a1+a1)) glColor3f(0.8,0.0,0.6);
x=pol0[i]; i++;
y=pol0[i]; i++;
if (x>+1e10) break;
if (x<-1e10){ glEnd(); glBegin(GL_LINE_STRIP); continue; }
glVertex2d(x,y);
}
glEnd();
// render a0 start point (yellow)
glBegin(GL_LINES);
glColor3f(1.0,1.0,0.0); r=100.0;
x=pol0[a0+a0 ];
y=pol0[a0+a0+1];
glVertex2d(x-r,y-r);
glVertex2d(x+r,y+r);
glVertex2d(x+r,y-r);
glVertex2d(x-r,y+r);
glEnd();
for (j=0;j<2;j++) // render pol1 twice
{
if (j)
{
// apply M matrix to view on second pass (align)
glMatrixMode(GL_MODELVIEW);
glMultMatrixd(M);
}
// render pol1 (white)
glBegin(GL_LINE_STRIP);
for (i=0;;)
{
if (i==b0+b0) glColor3f(0.2,0.9,0.2);
if ((!i)||(i==b1+b1)) glColor3f(0.7,0.7,0.7);
x=pol1[i]; i++;
y=pol1[i]; i++;
if (x>+1e10) break;
if (x<-1e10){ glEnd(); glBegin(GL_LINE_STRIP); continue; }
glVertex2d(x,y);
}
glEnd();
// render b0 start point (yellow)
glBegin(GL_LINES);
glColor3f(1.0,1.0,0.0); r=100.0;
x=pol1[b0+b0 ];
y=pol1[b0+b0+1];
glVertex2d(x-r,y-r);
glVertex2d(x+r,y+r);
glVertex2d(x+r,y-r);
glVertex2d(x-r,y+r);
glEnd();
}
glFlush();
SwapBuffers(hdc);
}
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
{
// application init
gl_init(Handle);
compute_MV();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
{
// application exit
gl_exit();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
{
// window resize
gl_resize(ClientWidth,ClientHeight);
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
{
// window repaint
gl_draw();
}
//---------------------------------------------------------------------------
Just ignore the VCL stuff... The important thing is compute_MV() function which computes M,V matrices. M is your aligning transform and V is just view transform to <-1,+1> range.
The atanxy(x,y) is my own function but its the same as atan2(y,x) which is common math function however that one needs to handle edge cases when x or y is zero otherwise exception might be thrown (that is what my atanxy does).
Used OpenGL matrix math is described in here and OpenGL stuff is in my complete GL+GLSL+VAO/VBO C++ example
Here preview of result
The pol0[] is in purple and pol1[] is in white twice once in raw state and second is aligned to pol0[]. The code is far from perfect and the search would need much more tweaking but I think its a good start point.
The search itself simply check constant size n window in both polylines and compute the abs difference between their delta angles and remembers the smallest distance positions ... Even if it is O(n.pnts0*pnts1) its pretty fast (less than 0.2 sec on my machine) ...
There is a lot of room for improvement like you do not need to step both polylines by single point its enough that just one loop is stepping by one point and the other might step by half or full window size... You can test for more window sizes... you can tweak the comparison and much more ...
In case your data is not sampled with the same average distance between points you should resample your data to common segment size first!
