geometry - 二维折线配准（对齐）算法

问题描述

我有来自两个来源的几条折线（由一组点定义） - 请参阅渲染图像，其中一个来源的折线是白色和另一个来源的紫色。

这两个部分可以对齐在一起，但是我不知道要使用什么算法。

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

你知道如何对齐集合吗？数据没有比例，只有平移（主要）和轻微旋转（+- 20 度）

折线的原始数据。有多条折线，每条折线由===

格式：

x1 y1
x2 y2
...

紫色的：

白色的

标签： geometrypattern-matching2d

解决方案

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

您可以使用这个如何比较两个形状？并找到角度序列“相同”的部分。这在旋转、缩放和平移上是不变的。因此，它应该导致直接匹配一条曲线上的哪些点按顺序匹配另一条曲线上的哪些点。

如果您的折线不完全匹配，您可以在角度上使用相关系数或阈值角度/角度差的总和。

如果您的点顺序可能在多段线之间颠倒，您需要检查这两种组合......

之后，它只是一个计算问题：

规模

只需在两条曲线上的两个远点（形成大线）之间划分大小
回转

在#1atan2的同一行上使用并减去角度
翻译

只需从两条曲线中减去一个点（但在对属于正在对齐的曲线的点应用缩放和旋转之后）。

或者您可以改用变换矩阵方法。

这里是蛮力搜索的简单 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,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();
    }
//---------------------------------------------------------------------------

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

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

此处描述了使用的 OpenGL 矩阵数学，OpenGL 的内容在我的完整 GL+GLSL+VAO/VBO C++ 示例中

这里预览结果

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

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

有很大的改进空间，例如您不需要将两条多段线按单点步进，它足以使一个循环步进一个点，而另一个循环可能步进一半或整个窗口大小...您可以测试更多窗口大小...您可以调整比较等等...

如果您的数据没有以相同的点之间的平均距离进行采样，您应该首先将数据重新采样为常见的段大小！

geometry - 二维折线配准（对齐）算法

问题描述

解决方案

推荐阅读