VTK 三角剖分 Delaunay2D (二)

标签: 计算机图形学  GIS  vtk  Delaunay2D  三角剖分  限制边界环  最短对角线

在这里插入图片描述

Description

随机生成上下两个三角网,由两个三角网的边界生成侧面三角网,组合成封闭的三角网。

在这里插入图片描述
在这里插入图片描述

Code

#include <vtkSmartPointer.h>

#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>

#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkNamedColors.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>

#include <vtkDelaunay2D.h>
#include <vtkFeatureEdges.h>
#include <vtkStripper.h>
#include <vtkCleanPolyData.h>
#include <vtkAppendPolyData.h>
#include <vtkVertexGlyphFilter.h>

int main(int, char *[])
{
  /*
  * 模拟生成第一个三角网
  * 1. 随机生成点集
  * 2. 获取生成三角网的边界环点集(有序)
  */
  vtkSmartPointer<vtkPoints> ps1 =vtkSmartPointer<vtkPoints>::New();
 
  unsigned int GridSize = 10;

  for(unsigned int x = 0; x < GridSize; x++)
  {
    for(unsigned int y = 0; y < GridSize; y++)
    {
		ps1->InsertNextPoint(x, y, vtkMath::Random(-.25, .25));
	}
  }

  vtkSmartPointer<vtkPolyData> polydata1 =vtkSmartPointer<vtkPolyData>::New();
  polydata1->SetPoints(ps1);

  vtkSmartPointer<vtkDelaunay2D> delaunay1 =vtkSmartPointer<vtkDelaunay2D>::New();
  delaunay1->SetInputData(polydata1);
  delaunay1->Update();

  //获取三角网边界点集
  vtkSmartPointer<vtkFeatureEdges> featureEdges1 =vtkSmartPointer<vtkFeatureEdges>::New();
  featureEdges1->SetInputConnection(delaunay1->GetOutputPort());
  featureEdges1->BoundaryEdgesOn();
  featureEdges1->FeatureEdgesOff();
  featureEdges1->ManifoldEdgesOff();
  featureEdges1->NonManifoldEdgesOff();
  featureEdges1->Update();
  
  //使点集有序
  vtkSmartPointer<vtkStripper> stripper1 =vtkSmartPointer<vtkStripper>::New();
  stripper1->SetInputConnection(featureEdges1->GetOutputPort());

  vtkSmartPointer<vtkCleanPolyData> cleanPolyData1 =
	  vtkSmartPointer<vtkCleanPolyData>::New();
  cleanPolyData1->SetInputConnection(stripper1->GetOutputPort());
  cleanPolyData1->Update();

  vtkPoints *pt = nullptr;
  double* p;
  unsigned int num1, num2;

  pt = cleanPolyData1->GetOutput()->GetPoints();
  num1 = pt->GetNumberOfPoints();

  //拷贝到边界点集中
  vtkSmartPointer<vtkPoints> boundary1 = vtkSmartPointer<vtkPoints>::New();
  //boundary1->SetNumberOfPoints(num1);
 // boundary1->DeepCopy(pt);
  //DeepCopy 拷贝节点会无序
  for (unsigned int i = 0; i < num1; ++i) {
	  p = pt->GetPoint(i);
	  boundary1->InsertNextPoint(p);
  }

  /*
  * 模拟生成第一个三角网
  * 1. 随机生成点集
  * 2. 获取生成三角网的边界环点集(有序)
  */
  vtkSmartPointer<vtkPoints> pts2 = vtkSmartPointer<vtkPoints>::New();
  GridSize = 9;
  for (unsigned int x =1; x < GridSize; x++)
  {
	  for (unsigned int y =1; y < GridSize; y++)
	  {
		  pts2->InsertNextPoint(x+ vtkMath::Random(-.25, 1.25), y+ vtkMath::Random(-.25, 1.25), vtkMath::Random(1.5, 1.75));
		 
	  }
  }

  vtkSmartPointer<vtkPolyData> polydata2 = vtkSmartPointer<vtkPolyData>::New();
  polydata2->SetPoints(pts2);
  vtkSmartPointer<vtkDelaunay2D> delaunay2 = vtkSmartPointer<vtkDelaunay2D>::New();
  delaunay2->SetInputData(polydata2);
  delaunay2->Update();


  vtkSmartPointer<vtkFeatureEdges> featureEdges2 =vtkSmartPointer<vtkFeatureEdges>::New();
  featureEdges2->SetInputConnection(delaunay2->GetOutputPort());
  featureEdges2->BoundaryEdgesOn();
  featureEdges2->FeatureEdgesOff();
  featureEdges2->ManifoldEdgesOff();
  featureEdges2->NonManifoldEdgesOff();
  featureEdges2->Update();

  vtkSmartPointer<vtkStripper> stripper2 =vtkSmartPointer<vtkStripper>::New();
  stripper2->SetInputConnection(featureEdges2->GetOutputPort());

  vtkSmartPointer<vtkCleanPolyData> cleanPolyData2 =vtkSmartPointer<vtkCleanPolyData>::New();
  cleanPolyData2->SetInputConnection(stripper2->GetOutputPort());
  cleanPolyData2->Update();
  pt = cleanPolyData2->GetOutput()->GetPoints();
  num2 = pt->GetNumberOfPoints();
  //将第二个边界点集拷贝到边界点集中
  vtkSmartPointer<vtkPoints> boundary2 = vtkSmartPointer<vtkPoints>::New();
  //boundary2->SetNumberOfPoints(num2);
  //boundary2->DeepCopy(pt);
  for (unsigned int i = 0; i < num1; ++i) {
	  p = pt->GetPoint(i);
	  boundary2->InsertNextPoint(p);

  }

   /*
  * 由上下两个边界环生成三角网
  * 原则:按照最短对角线原则
  * 方法:上下两个边界环,依次组成四边形,比较对角线的大小,生成三角形
  */
  vtkPoints *boundary = vtkPoints::New();
 
  for (unsigned int i = 0; i < num1; ++i) {
	  p = boundary1->GetPoint(i);
	  boundary->InsertNextPoint(p);
  }
  
  for (unsigned int i = 0; i < num2; ++i) {
	  p = boundary2->GetPoint(i);
	  boundary->InsertNextPoint(p);
  }
  
  vtkIdType pts[3];               //保存三角形三个顶点的索引
  vtkPolyData* meshPloyData = vtkPolyData::New();
  vtkCellArray* triangles = vtkCellArray::New();
  triangles->Allocate(triangles->EstimateSize(2 * (num1+num2), 3));
  meshPloyData->SetPoints(boundary);
  meshPloyData->SetPolys(triangles);
  
  double d0 = 0.0,d1=0.0;         //保存四边形两条对角线的长度
  double p1[3],p2[3],p3[3],p4[3]; //四边形的四个顶点
  unsigned int n1=0, n2=0;        //上下边界环当前遍历的位置
  unsigned int start1=0, start2=0;//上下边界环遍历的起始位置
  
  start1 = 0;// 设置第一条边界环遍历的起始位置 start1=0
  boundary1->GetPoint(0, p1);
  //寻找第二条边的起始位置 start2 (距离最近的点)
  d1 = VTK_DOUBLE_MAX;
  for (unsigned int i = 0; i < num2; ++i) {
	  boundary2->GetPoint(i, p2);
	  d0= vtkMath::Distance2BetweenPoints(p1, p2);
	  if (d0 < d1) {
		  d1 = d0;
		  start2 = i;
	  }
  }

  n2 = start2;
  bool flag = false;
 do {
	 // 组成四边形的四个点
	 // p1,p2 一条边界环的点
	 // p3,p4 另一条边界环的点

	 boundary1->GetPoint(n1, p1);
	 boundary1->GetPoint((n1+1)%num1, p2);

	 boundary2->GetPoint(n2, p3);
	 boundary2->GetPoint((n2+1)%num2, p4);

	 d0 = vtkMath::Distance2BetweenPoints(p1, p4);
	 d1 = vtkMath::Distance2BetweenPoints(p2, p3);

	 pts[0] = n1%num1;
	 pts[1] = n2%num2 + num1; 

	 if (d0<=d1) {
		  pts[2] = (n2 + 1) % num2 + num1;
		  n2 = (n2 + 1)% num2;

	  }else{
		  pts[2] = (n1 + 1) % num1; 
		  n1 = (n1 + 1) % num1;
	  }
	  triangles->InsertNextCell(3, pts);

	  flag = !((n1%num1 == start1||(n1 + 1) % num1 == start1) && (n2%num2 == start2||(n2 + 1) % num2 == start2));

  }while (flag);

 //最后两个三角形的情况
 if ((n1+1)%num1 == start1 && ( n2+1)%num2 == start2) {
	
	 boundary1->GetPoint(n1, p1);
	 boundary1->GetPoint((n1 + 1) % num1, p2);
	 boundary2->GetPoint(n2, p3);
	 boundary2->GetPoint((n2 + 1) % num2, p4);

	 d0 = vtkMath::Distance2BetweenPoints(p1, p4);
	 d1 = vtkMath::Distance2BetweenPoints(p2, p3);

	 if (d0 <= d1) {

		 pts[0] = n1%num1;
		 pts[1] = n2%num2 + num1; 
		 pts[2] = (n2 + 1) % num2 + num1; 
		 triangles->InsertNextCell(3, pts);

		 pts[0] = n1%num1;
		 pts[1] = (n2 + 1) % num2 + num1;
		 pts[2] = (n1 + 1) % num1;
		 triangles->InsertNextCell(3, pts);
	 }
	 else {

		 pts[0] = n1%num1;
		 pts[1] = n2%num2 + num1;
		 pts[2] = (n1 + 1) % num1;
		 triangles->InsertNextCell(3, pts);

		 pts[0] = (n1 + 1) % num1;
		 pts[1] = (n2 + 1) % num2 + num1;
		 pts[2] = n2%num2 + num1;
		 triangles->InsertNextCell(3, pts);
	 }
	 
 }
 else {  // 最后一个三角形的情况

	 if (n1 == 0) {

		 pts[0] = n1%num1;
		 pts[1] = n2%num2 + num1;
		 pts[2] = (n2 + 1) % num2 + num1;

	 }
	 else {
		 pts[0] = n1%num1;
		 pts[1] = n2%num2 + num1;
		 pts[2] = (n1 + 1) % num1;
	 }
	 triangles->InsertNextCell(3, pts);
 }

   meshPloyData->BuildLinks();
  // combine two poly data
  // 上下三角网和侧面三角网组合,生成封闭的体
   vtkSmartPointer<vtkAppendPolyData> appendFilter =vtkSmartPointer<vtkAppendPolyData>::New();
   appendFilter->AddInputData(delaunay1->GetOutput());
   appendFilter->AddInputData(delaunay2->GetOutput());
   appendFilter->AddInputData(meshPloyData);
   appendFilter->Update();

  /*vtkSmartPointer<vtkCleanPolyData> cleanPolyData = vtkSmartPointer<vtkCleanPolyData>::New();
  cleanPolyData->SetInputConnection(appendFilter->GetOutputPort());
  cleanPolyData->Update();*/
  // Visualize
  vtkSmartPointer<vtkNamedColors> colors =vtkSmartPointer<vtkNamedColors>::New();

  vtkSmartPointer<vtkPolyDataMapper> meshMapper =vtkSmartPointer<vtkPolyDataMapper>::New();
  meshMapper->SetInputConnection(appendFilter->GetOutputPort());
  /*meshMapper->SetInputData(meshPloyData);
  meshMapper->Update();*/

  vtkSmartPointer<vtkActor> meshActor =vtkSmartPointer<vtkActor>::New();
  meshActor->SetMapper(meshMapper);
  meshActor->GetProperty()->SetColor(colors->GetColor3d("Banana").GetData());
  meshActor->GetProperty()->EdgeVisibilityOn();

  vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =vtkSmartPointer<vtkVertexGlyphFilter>::New();
  glyphFilter->SetInputData(appendFilter->GetOutput());

  vtkSmartPointer<vtkPolyDataMapper> pointMapper =vtkSmartPointer<vtkPolyDataMapper>::New();
  pointMapper->SetInputConnection(glyphFilter->GetOutputPort());

  vtkSmartPointer<vtkActor> pointActor =vtkSmartPointer<vtkActor>::New();
  pointActor->GetProperty()->SetColor(colors->GetColor3d("Tomato").GetData());
  pointActor->GetProperty()->SetPointSize(5);
  pointActor->SetMapper(pointMapper);

  vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
  vtkSmartPointer<vtkRenderWindow> renderWindow =vtkSmartPointer<vtkRenderWindow>::New();
  renderWindow->AddRenderer(renderer);
  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =vtkSmartPointer<vtkRenderWindowInteractor>::New();
  renderWindowInteractor->SetRenderWindow(renderWindow);

  renderer->AddActor(meshActor);
  renderer->AddActor(pointActor);

  renderer->SetBackground(colors->GetColor3d("Mint").GetData());

  renderWindow->Render();
  renderWindowInteractor->Start();

  return EXIT_SUCCESS;
}

算法示意图

在这里插入图片描述

参考

VTK 三角剖分 Delaunay2D (一)

版权声明:本文为mrbaolong原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/mrbaolong/article/details/104876411