zoukankan      html  css  js  c++  java
  • VTK 图形基本操作进阶_点云配准技术(迭代最近点ICP算法)

    1.Iterative Closest Points算法

    点云数据配准最经典的方法是迭代最近点算法(Iterative Closest Points,ICP)。ICP算法是一个迭代的过程,每次迭代中对于源数据点P找到目标点集Q中的最近点,然后给予最小二乘原理求解当前的变换矩阵T。通过不断迭代迭代直至收敛,即完成了点集的配准。

    1.1 基本原理

    ICP算法是大多数点云配准算法的心, 它是一个点对刚性算法。基本思想是: 假设两个点集 P和 X近似对齐, 对 P上的每个点,假设 X上的最近点与之对齐。采用最近点搜索, 在 X上找出 P上各个点对应的最近点, 构成集合 Y, 然后计算一个新的 P到 Y的刚体变换。重复上述过程直到配准收敛。

    1.2 算法步骤(四元数配准)

    假设给两个三维点集 X1 和 X2,ICP方法的配准步骤如下:
    第一步,计算X2中的每一个点在X1 点集中的对应近点;
    第二步,求得使上述对应点对平均距离最小的刚体变换,求得平移参数和旋转参数;
    第三步,对X2使用上一步求得的平移和旋转参数,得到新的变换点集;
    第四步, 如果新的变换点集与参考点集满足两点集的平均距离小于某一给定阈值,则停止迭代计算,否则新的变换点集作为新的X2继续迭代,直到达到目标函数的要求。

    1.3 最近点对查找

    对应点的计算是整个配准过程中耗费时间最长的步骤,查找最近点,利用 k-d tree提高查找速度 K-d tree 法建立点的拓扑关系是基于二叉树的坐标轴分割,构造 k-d tree 的过程就是按照二叉树法则生成,首先按 X 轴寻找分割线,即计算所有点的x值的平均值,以最接近这个平均值的点的x值将空间分成两部分,然后在分成的子空间中按 Y 轴寻找分割线,将其各分成两部分,分割好的子空间在按X轴分割……依此类推,最后直到分割的区域内只有一个点。这样的分割过程就对应于一个二叉树,二叉树的分节点就对应一条分割线,而二叉树的每个叶子节点就对应一个点。这样点的拓扑关系就建立了。

    2.VTK中实现ICP算法实验

    vtk中已经封装了最基本的ICP算法,由类vtkIterativeClosestPointTransform负责。
    具体的示例代码如下所示:
      1 #include <vtkAutoInit.h>
      2 VTK_MODULE_INIT(vtkRenderingOpenGL);
      3 VTK_MODULE_INIT(vtkRenderingFreeType);
      4 VTK_MODULE_INIT(vtkInteractionStyle);
      5  
      6 #include <vtkSmartPointer.h>
      7 #include <vtkPolyDataReader.h>
      8 #include <vtkPolyData.h>
      9 #include <vtkTransform.h>
     10 #include <vtkTransformPolyDataFilter.h>
     11 #include <vtkVertexGlyphFilter.h>
     12 #include <vtkIterativeClosestPointTransform.h>
     13 #include <vtkLandmarkTransform.h>
     14 #include <vtkTransformPolyDataFilter.h>
     15 #include <vtkPolyDataMapper.h>
     16 #include <vtkActor.h>
     17 #include <vtkProperty.h>
     18 #include <vtkAxesActor.h>
     19 #include <vtkRenderer.h>
     20 #include <vtkRenderWindow.h>
     21 #include <vtkRenderWindowInteractor.h>
     22 #include <vtkOrientationMarkerWidget.h> //坐标系交互
     23 int main()
     24 {
     25     vtkSmartPointer <vtkPolyDataReader> reader =
     26         vtkSmartPointer<vtkPolyDataReader>::New();
     27     reader->SetFileName("fran_cut.vtk");
     28     reader->Update();
     29  
     30     //构造浮动数据点集
     31     vtkSmartPointer<vtkPolyData> orig = reader->GetOutput();
     32     vtkSmartPointer<vtkTransform> trans =
     33         vtkSmartPointer<vtkTransform>::New();
     34     trans->Translate(0.2, 0.1, 0.1);
     35     trans->RotateX(10);
     36  
     37     vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter1 =
     38         vtkSmartPointer<vtkTransformPolyDataFilter>::New();
     39     transformFilter1->SetInputData(reader->GetOutput());
     40     transformFilter1->SetTransform(trans);
     41     transformFilter1->Update();
     42     /*********************************************************/
     43     //源数据 与 目标数据
     44     vtkSmartPointer<vtkPolyData> source =
     45         vtkSmartPointer<vtkPolyData>::New();
     46     source->SetPoints(orig->GetPoints());
     47  
     48     vtkSmartPointer<vtkPolyData> target =
     49         vtkSmartPointer<vtkPolyData>::New();
     50     target->SetPoints(transformFilter1->GetOutput()->GetPoints());
     51  
     52     vtkSmartPointer<vtkVertexGlyphFilter>  sourceGlyph =
     53         vtkSmartPointer<vtkVertexGlyphFilter>::New();
     54     sourceGlyph->SetInputData(source);
     55     sourceGlyph->Update();
     56  
     57     vtkSmartPointer<vtkVertexGlyphFilter>  targetGlyph =
     58         vtkSmartPointer<vtkVertexGlyphFilter>::New();
     59     targetGlyph->SetInputData(target);
     60     targetGlyph->Update();
     61  
     62     //进行ICP配准求变换矩阵
     63     vtkSmartPointer<vtkIterativeClosestPointTransform> icptrans =
     64         vtkSmartPointer<vtkIterativeClosestPointTransform>::New();
     65     icptrans->SetSource(sourceGlyph->GetOutput());
     66     icptrans->SetTarget(targetGlyph->GetOutput());
     67     icptrans->GetLandmarkTransform()->SetModeToRigidBody();
     68     icptrans->SetMaximumNumberOfIterations(50);
     69     icptrans->StartByMatchingCentroidsOn();
     70     icptrans->Modified();
     71     icptrans->Update();
     72  
     73     //配准矩阵调整源数据
     74     vtkSmartPointer<vtkTransformPolyDataFilter> solution =
     75         vtkSmartPointer<vtkTransformPolyDataFilter>::New();
     76     solution->SetInputData(sourceGlyph->GetOutput());
     77     solution->SetTransform(icptrans);
     78     solution->Update();
     79     /
     80     vtkSmartPointer<vtkPolyDataMapper> sourceMapper =
     81         vtkSmartPointer<vtkPolyDataMapper>::New();
     82     sourceMapper->SetInputConnection(sourceGlyph->GetOutputPort());
     83     vtkSmartPointer<vtkActor> sourceActor =
     84         vtkSmartPointer<vtkActor>::New();
     85     sourceActor->SetMapper(sourceMapper);
     86     sourceActor->GetProperty()->SetColor(1, 1, 0);
     87     sourceActor->GetProperty()->SetPointSize(2);
     88  
     89     vtkSmartPointer<vtkPolyDataMapper> targetMapper =
     90         vtkSmartPointer<vtkPolyDataMapper>::New();
     91     targetMapper->SetInputConnection(targetGlyph->GetOutputPort());
     92     vtkSmartPointer<vtkActor> targetActor =
     93         vtkSmartPointer<vtkActor>::New();
     94     targetActor->SetMapper(targetMapper);
     95     targetActor->GetProperty()->SetColor(0, 1, 0);
     96     targetActor->GetProperty()->SetPointSize(3);
     97  
     98     vtkSmartPointer<vtkPolyDataMapper> soluMapper =
     99         vtkSmartPointer<vtkPolyDataMapper>::New();
    100     soluMapper->SetInputConnection(solution->GetOutputPort());
    101     vtkSmartPointer<vtkActor> soluActor =
    102         vtkSmartPointer<vtkActor>::New();
    103     soluActor->SetMapper(soluMapper);
    104     soluActor->GetProperty()->SetColor(1, 0, 0);
    105     soluActor->GetProperty()->SetPointSize(2);
    106     //设置坐标系
    107     vtkSmartPointer<vtkAxesActor> axes =
    108         vtkSmartPointer<vtkAxesActor>::New();
    109     ///
    110     vtkSmartPointer<vtkRenderer> render =
    111         vtkSmartPointer<vtkRenderer>::New();
    112     render->AddActor(sourceActor);
    113     render->AddActor(targetActor);
    114     render->AddActor(soluActor);
    115     render->SetBackground(0, 0, 0);
    116  
    117     vtkSmartPointer<vtkRenderWindow> rw =
    118         vtkSmartPointer<vtkRenderWindow>::New();
    119     rw->AddRenderer(render);
    120     rw->SetSize(480, 320);
    121     rw->SetWindowName("Regisration by ICP");
    122  
    123     vtkSmartPointer<vtkRenderWindowInteractor> rwi =
    124         vtkSmartPointer<vtkRenderWindowInteractor>::New();
    125     rwi->SetRenderWindow(rw);
    126     /****************************************************************/
    127     //谨记:顺序不可以颠倒!!!
    128     vtkSmartPointer<vtkOrientationMarkerWidget> widget =
    129         vtkSmartPointer<vtkOrientationMarkerWidget>::New();
    130     widget->SetOutlineColor(1, 1, 1);
    131     widget->SetOrientationMarker(axes);
    132     widget->SetInteractor(rwi); //加入鼠标交互  
    133     widget->SetViewport(0.0, 0.0, 0.3, 0.3);  //设置显示位置  
    134     widget->SetEnabled(1);
    135     widget->InteractiveOn();//开启鼠标交互  
    136     /****************************************************************/
    137     render->ResetCamera();
    138     rw->Render();
    139     rwi->Initialize();
    140     rwi->Start();
    141  
    142     return 0;
    143 }
    这个例子首先读取一个人脸模型。为了方便测试效果,这里对原始模型做了一个平移和旋转变换。
    这里面有个细节应该正视:
    vtkIterativeClosestPointTransform类中设置源点集和目标点集的函数为SetSource()和SetTarget(),其输入数据类型为VTKDataSet,所以集合点集必须进过一定的处理!这里使用vtkVertexGlyphFilter将读入模型和变换后的点集转换为相应的vtkPolyData数据。
    1 vtkSmartPointer<vtkVertexGlyphFilter>  sourceGlyph =
    2         vtkSmartPointer<vtkVertexGlyphFilter>::New();
    3     sourceGlyph->SetInputData(source);
    4     sourceGlyph->Update();
    5  
    6     vtkSmartPointer<vtkVertexGlyphFilter>  targetGlyph =
    7         vtkSmartPointer<vtkVertexGlyphFilter>::New();
    8     targetGlyph->SetInputData(target);
    9     targetGlyph->Update();

    vtkLandmarkTransform类有点不同,输入数据类型就是vtkPoints类型。

     1 //进行ICP配准求变换矩阵
     2     vtkSmartPointer<vtkIterativeClosestPointTransform> icptrans =
     3         vtkSmartPointer<vtkIterativeClosestPointTransform>::New();
     4     icptrans->SetSource(sourceGlyph->GetOutput());
     5     icptrans->SetTarget(targetGlyph->GetOutput());
     6     icptrans->GetLandmarkTransform()->SetModeToRigidBody();
     7     icptrans->SetMaximumNumberOfIterations(50);
     8     icptrans->StartByMatchingCentroidsOn();
     9     icptrans->Modified();
    10     icptrans->Update();
    StartByMatchingCentroidsOn()该函数就是去偏移(中心归一/重心归一)。通过分别计算重心,然后平移,使得两点集中心重合。
    配准后的结果如下图:
    黄色点集是浮动点集;绿色点集是金标准;红色点集是经过配准矩阵调整后的点集。

    2.vtkOrientationMarkerWidget类设置坐标系心得

    1     //谨记:顺序不可以颠倒!!!
    2     vtkSmartPointer<vtkOrientationMarkerWidget> widget =
    3         vtkSmartPointer<vtkOrientationMarkerWidget>::New();
    4     widget->SetOutlineColor(1, 1, 1);
    5     widget->SetOrientationMarker(axes);
    6     widget->SetInteractor(rwi); //加入鼠标交互  
    7     widget->SetViewport(0.0, 0.0, 0.3, 0.3);  //设置显示位置  
    8     widget->SetEnabled(1);
    9     widget->InteractiveOn();//开启鼠标交互
     
  • 相关阅读:
    Java Character类
    JavaMath&Number类
    oracle中创建数据库实例
    Jenkins安装和配置
    crontab计划 -计划任务和日志管理
    centos修改默认字符编码为GBK(服务器端java GBK字符乱码)
    使用vsphere client 克隆虚拟机
    建索引让SQL飞起来
    oracle中对LONG列进行查询
    优化sql,返回行数少情况下,NL比hash快好多
  • 原文地址:https://www.cnblogs.com/ybqjymy/p/14241995.html
Copyright © 2011-2022 走看看