2D折线配准对齐算法

如何解决2D折线配准对齐算法

我有两个来源的多条折线（由点集定义）-看到渲染的图像，其中一个来源的折线是白色，而另一个来源的折线是紫色。

这两个部分可以对齐，但是我不知道使用哪种算法。

我尝试使用ICP（迭代最近点），但是问题是线太远了。我尝试创建修改-而不是直接使用最近的点，而是使用相似角度附近的最近的点。但是，要检测相似的角度是行不通的，因为线条不一样-一个角度很锐利，而另一个角度却很平滑。

您知道如何对齐布景吗？数据没有刻度，只有平移（主要是）和轻微旋转（+-20度）

折线的原始数据。有多条折线，每条折线之间用===

隔开

格式：

x1 y1
x2 y2
...

紫色：

白色

解决方法

是点云还是折线？
两条曲线的采样比例是否相同？

您可以使用此How to compare two shapes?来查找角度顺序为“相同”的零件。这在旋转，缩放和平移方面是不变的。因此，它应该导致直接匹配，一条曲线中的哪些点与另一条曲线中的哪些点顺序匹配。

如果折线不完全匹配，则可以在角度上使用correlation coefficient或对角度之和/角度数量进行阈值确定。

如果折线之间的点顺序可能相反，则需要检查两种组合...

之后，这仅是计算问题：

规模

仅需在两条曲线的2个遥远点之间划分大小（形成大线）即可
旋转

在与＃1 相同的行上使用atan2并减去角度
翻译

仅从两条曲线中减去单个点（但在对要对齐的曲线的点上应用比例和旋转之后）。

或者您也可以使用transform matrix approach。

以下是暴力搜索的简单C ++ / GL / VCL小示例：

//---------------------------------------------------------------------------
#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,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,6369,3082,6583,3116,6526,3203,6465,3242,6366,3340,6192,3414,5922,3457,5732,3576,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,6584,3115,6745,3059,6802,3054,6979,3053,7129,3089,7344,3139,7299,3158,7269,3150,6900,3306,6756,3397,3436,7378,3441,7345,3474,7237,3513,7205,3512,7168,3544,7079,3578,7060,3603,3580,7216,3601,7180,3701,7190,3649,7057,6984,3678,7040,7019,3607,6915,3653,6899,3679,6876,3658,6857,3682,6741,3734,6726,3758,6293,3755,6186,6042,3901,6551,3798,6536,3822,6400,3893,6539,3862,6442,4793,6372,3879,6217,3970,6200,3948,6056,4046,6036,4078,6013,4615,7294,3867,7245,4069,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,2825,2223,2819,2399,2737,2489,2732,2895,3080,2878,3193,2973,3301,3032,3426,3042,3614,3000,3890,2980,4027,3029,4073,3062,4089,3097,4082,3168,3788,3130,3660,3548,3195,3409,3299,3307,3404,3112,3489,2860,3520,3590,2720,3611,2728,3641,2796,3772,2907,3871,2903,3900,2908,2968,3882,2984,3928,2991,3992,2954,4087,2638,4224,2203,4503,1588,4863,1479,4895,1184,5035,252,5489,235,5518,1912,257,1964,228,2035,1962,332,1963,2152,4342,2038,4436,1872,4528,1795,4616,1818,4568,1881,4574,1441,4864,1303,4904,927,5098,831,5141,767,5152,597,5017,669,4938,692,4942,838,4944,802,4981,855,5057,904,5054,843,5043,786,5062,795,5084,784,5063,732,5034,641,5086,5132,640,581,550,5107,585,5207,742,5031,760,4997,799,1554,4604,1427,4654,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,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,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,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.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,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();
    }
//---------------------------------------------------------------------------

只需忽略VCL内容...重要的是compute_MV()函数，它计算M,V矩阵。 M是您的对齐变换，V只是视图变换到<-1,+1>范围。

atanxy(x,y)是我自己的函数，但与atan2(y,x)相同，后者是常用的数学函数，但是当x或y为零时，需要处理边缘情况，否则可能会引发异常（即我的atanxy做什么）。

here中描述了使用的OpenGL矩阵数学，complete GL+GLSL+VAO/VBO C++ example中包含了OpenGL内容

此处预览结果

在原始状态下，pol0[]呈紫色，pol1[]呈白色两次，第二次与pol0[]对齐。代码远非完美，搜索需要更多调整，但我认为这是一个很好的起点。

搜索本身只需检查两条折线中恒定大小的n窗口并计算其delta angles之间的绝对差，并记住最小的距离位置...即使它是O(n.pnts0*pnts1)相当快（我的机器上不到0.2秒）...

还有很多改进的余地，就像您不需要将两条折线一步一步地移动一样，只需要一个循环就一步一步地移动，另一个循环就可以将窗口大小的一半或全开。测试更多的窗口大小...您可以调整比较以及更多...

如果未以相同的点间平均距离采样数据，则应先将数据采样到共同的分段大小！