VTK编程指南＜十六＞：VTK表面重建之Delaunay三角剖分

1、无边界限制

  三角剖分是一种应用非常广泛的重建技术。三角剖分将一些散乱的点云数据划分为一系列的三角形网格。最常用的三角剖分技术是Delaunay三角剖分。Delaunay三角剖分具有许多优良的性质，即最大化最小角特性，即所有可能的三角形剖分中，其生成的三角形最小角的角度最大。所以，Delaunay三角剖分无论在哪个区域开始构建，最终生成的三角网格是唯一的。

  VTK的vtkDelaunay2D类实现了二维三角剖分。该类的输入数据集为一个vtkPointSet或其它类表示的三维空间点集，其中输出为一个三角网格vtkPolyData数据。虽然输入的三维数据， 但是计算仅使用XY平面数据进行平面三角剖分，而忽略Z方向数据。 当然，也可以为vtkDelaunay2D设置一个投影转换从而在投影平面上进行三角剖分。需要注意的是，在不加任何限制的情况下，该类生成的平面网格为一个凸包。下面通过示例演示vtkDelaunay2D， 将其生成的数据信息用于模拟地形数据（程序执行结果如图所示）：

#include <vtkAutoInit.h>
VTK_MODULE_INIT(vtkRenderingOpenGL2);// VTK was built with vtkRenderingOpenGL2
VTK_MODULE_INIT(vtkInteractionStyle);#include <vtkSmartPointer.h>
#include <vtkProperty.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkPointData.h>
#include <vtkDelaunay2D.h>
#include <vtkMath.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkVertexGlyphFilter.h>int main(int, char *[])
{unsigned int gridSize = 10;vtkSmartPointer<vtkPoints> points =vtkSmartPointer<vtkPoints>::New();for(unsigned int x = 0; x < gridSize; x++){for(unsigned int y = 0; y < gridSize; y++){points->InsertNextPoint(x, y, vtkMath::Random(0.0, 3.0));}}vtkSmartPointer<vtkPolyData> polydata =vtkSmartPointer<vtkPolyData>::New();polydata->SetPoints(points);vtkSmartPointer<vtkDelaunay2D> delaunay =vtkSmartPointer<vtkDelaunay2D>::New();delaunay->SetInputData(polydata);delaunay->Update();vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =vtkSmartPointer<vtkVertexGlyphFilter>::New();glyphFilter->SetInputData(polydata);glyphFilter->Update();vtkSmartPointer<vtkPolyDataMapper> pointsMapper =vtkSmartPointer<vtkPolyDataMapper>::New();pointsMapper->SetInputData(glyphFilter->GetOutput());vtkSmartPointer<vtkActor> pointsActor =vtkSmartPointer<vtkActor>::New();pointsActor->SetMapper(pointsMapper);pointsActor->GetProperty()->SetPointSize(3);pointsActor->GetProperty()->SetColor(1,0,0);vtkSmartPointer<vtkPolyDataMapper> triangulatedMapper =vtkSmartPointer<vtkPolyDataMapper>::New();triangulatedMapper->SetInputData(delaunay->GetOutput());vtkSmartPointer<vtkActor> triangulatedActor =vtkSmartPointer<vtkActor>::New();triangulatedActor->SetMapper(triangulatedMapper);vtkSmartPointer<vtkRenderer> renderer =vtkSmartPointer<vtkRenderer>::New();renderer->AddActor(pointsActor);renderer->AddActor(triangulatedActor);renderer->SetBackground(1.0, 1.0, 1.0);vtkSmartPointer<vtkRenderWindow> renderWindow =vtkSmartPointer<vtkRenderWindow>::New();renderWindow->AddRenderer(renderer);vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =vtkSmartPointer<vtkRenderWindowInteractor>::New();renderWindowInteractor->SetRenderWindow(renderWindow);renderWindow->SetSize(640, 320);renderWindow->Render();renderWindow->SetWindowName("PolyDataDelaunay2D");renderWindow->Render();renderWindowInteractor->Start();return EXIT_SUCCESS;
}

2、边界限制

  vtkDelaunay2D 还支持加入边界限制。用户需要设置另外一个vtkPolyData数据，其中部分的线段、闭合或非闭合的线段集合作为边界条件控制三角剖分的过程。过程中组成这些边界的点的索引必须与原始点集数据一致。加入边界条件后，最后的剖分结果可能不再满足Delaunay准则。在上述例子的基础上，加入一个多边形边界来限制三角剖分，这里定义了一个vtkPolyData类型的数据boundary，其点数据与上例中的points一致，其单元数据为一个多边形。通过vtkDelaunay2D的SetSource()函数设置边界数据，运行结果如图6-22a所示，该边界多边形内部的数据并未进行三角剖分。对于边界多边形数据的内部或者外部，与多边形点的顺序有关。这里采用右手坐标系，从Z轴向下看去，如果多边形的点顺序为逆时针，则对多边形内部数据进行剖分；而如果顺时针方向，则对多边形外部数据进行剖分，此时该边界多边形可以看作一个孔洞。将上例中的点数据反向，则剖分结果如图所示。

#include <vtkAutoInit.h>
VTK_MODULE_INIT(vtkRenderingOpenGL2);// VTK was built with vtkRenderingOpenGL2
VTK_MODULE_INIT(vtkInteractionStyle);#include <vtkSmartPointer.h>
#include <vtkProperty.h>
#include <vtkPolygon.h>
#include <vtkCellArray.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkPointData.h>
#include <vtkDelaunay2D.h>
#include <vtkMath.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkVertexGlyphFilter.h>int main(int, char *[])
{vtkSmartPointer<vtkPoints> points =vtkSmartPointer<vtkPoints>::New();unsigned int gridSize = 10;for(unsigned int x = 0; x < gridSize; x++){for(unsigned int y = 0; y < gridSize; y++){points->InsertNextPoint(x, y, vtkMath::Random(0.0, 3.0));}}vtkSmartPointer<vtkPolyData> polydata =vtkSmartPointer<vtkPolyData>::New();polydata->SetPoints(points);vtkSmartPointer<vtkPolygon> poly =vtkSmartPointer<vtkPolygon>::New();/*poly->GetPointIds()->InsertNextId(22);poly->GetPointIds()->InsertNextId(23);poly->GetPointIds()->InsertNextId(24);poly->GetPointIds()->InsertNextId(25);poly->GetPointIds()->InsertNextId(35);poly->GetPointIds()->InsertNextId(45);poly->GetPointIds()->InsertNextId(44);poly->GetPointIds()->InsertNextId(43);poly->GetPointIds()->InsertNextId(42);poly->GetPointIds()->InsertNextId(32);*/poly->GetPointIds()->InsertNextId(32);poly->GetPointIds()->InsertNextId(42);poly->GetPointIds()->InsertNextId(43);poly->GetPointIds()->InsertNextId(44);poly->GetPointIds()->InsertNextId(45);poly->GetPointIds()->InsertNextId(35);poly->GetPointIds()->InsertNextId(25);poly->GetPointIds()->InsertNextId(24);poly->GetPointIds()->InsertNextId(23);poly->GetPointIds()->InsertNextId(22);vtkSmartPointer<vtkCellArray> cell =vtkSmartPointer<vtkCellArray>::New();cell->InsertNextCell(poly);vtkSmartPointer<vtkPolyData> boundary =vtkSmartPointer<vtkPolyData>::New();boundary->SetPoints(points);boundary->SetPolys(cell);vtkSmartPointer<vtkDelaunay2D> delaunay =vtkSmartPointer<vtkDelaunay2D>::New();delaunay->SetInputData(polydata);delaunay->SetSourceData(boundary);delaunay->Update();vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =vtkSmartPointer<vtkVertexGlyphFilter>::New();glyphFilter->SetInputData(polydata);glyphFilter->Update();vtkSmartPointer<vtkPolyDataMapper> pointsMapper =vtkSmartPointer<vtkPolyDataMapper>::New();pointsMapper->SetInputData(glyphFilter->GetOutput());vtkSmartPointer<vtkActor> pointsActor =vtkSmartPointer<vtkActor>::New();pointsActor->SetMapper(pointsMapper);pointsActor->GetProperty()->SetPointSize(3);pointsActor->GetProperty()->SetColor(1,0,0);vtkSmartPointer<vtkPolyDataMapper> triangulatedMapper =vtkSmartPointer<vtkPolyDataMapper>::New();triangulatedMapper->SetInputData(delaunay->GetOutput());vtkSmartPointer<vtkActor> triangulatedActor =vtkSmartPointer<vtkActor>::New();triangulatedActor->SetMapper(triangulatedMapper);vtkSmartPointer<vtkRenderer> renderer =vtkSmartPointer<vtkRenderer>::New();renderer->AddActor(pointsActor);renderer->AddActor(triangulatedActor);renderer->SetBackground(1.0, 1.0, 1.0);vtkSmartPointer<vtkRenderWindow> renderWindow =vtkSmartPointer<vtkRenderWindow>::New();renderWindow->AddRenderer(renderer);vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =vtkSmartPointer<vtkRenderWindowInteractor>::New();renderWindowInteractor->SetRenderWindow(renderWindow);renderWindow->SetSize(640, 480);renderWindow->Render();renderWindow->SetWindowName("PolyDataConstrainedDelaunay2D");renderWindow->Render();renderWindowInteractor->Start();return EXIT_SUCCESS;
}