CGAL-三维Delaunay/Voronoi图/Powerdiagram

综述

3D Vornoi和Powerdiagram的例子比较少。这里给出一些方法。

说明

这里给出两个版本
版本一
使用一般的RT完成
（注意在add_point中最后一个参数带权重则为powerdiagram，如果为0就为Voronoidiagram）
支持Delaunay三角的顶点输出
支持显示所有Voronoi顶点输出
（一般使用推荐使用版本一）

版本二
使用LCC map完成，数据查询更为方便，但是数据结构较难组织
支持Delaunay三角的顶点输出
支持显示所有Voronoi顶点输出
支持显示所有有限空间划分的顶点输出

更新

2018.8.14
2018.8.16 提供了一种更为普遍的方法

有网友和我交流问“3D powerdiagram”是否可以这样组织（通过三角正则化求对偶）？
事实上是不能的。
因为LCC没有Weighted_Point的类型，也就是说不支持。
所以还是要使用：
CGAL::Exact_predicates_inexact_constructions_kernel
下的
Weighted_Point;

环境

clion
maxos

代码

版本一代码

//作者：SDU窦志扬
//日期：2018/8/16
//联络：sdudzy@163.com
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Weighted_point_3.h>
#include <GLUT/glut.h>
#include <CGAL/bounding_box.h>
#include <fstream>
#include<iostream>
#include<set>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Regular_triangulation_3<K> Regular_triangulation;
typedef K::Vector_3                                     Vector;
typedef K::Point_3 Point;
typedef K::Weighted_point_3 Wp;
typedef K::Point_3                                     Point;
typedef CGAL::Regular_triangulation_3<K>               Rt;
using namespace std;
vector<Regular_triangulation::Weighted_point> wpoints;
vector<Point> points;
vector<double> X, Y, Z;
set<Point> Voronoi_vert;
void add_point(double x, double y, double z , double w){wpoints.push_back( Wp(Point(x,y,z), w) );points.push_back(Point(x,y,z) );X.push_back(x);Y.push_back(y);Z.push_back(z);}void display(void)
{add_point(100+0, 100+0,100+0,0); //左中add_point(100+100,100+0,100+0,0);add_point(100+0,100+0,100+400,0);   //目标点add_point(100+0, 100+200,100+0,0);add_point(300+0, 70+200,100+0,0);add_point(20+0, 19+200,80+0,0);add_point(100+0, 100+200,100+0,0);add_point(100+10, 100+20,100+20,000);Regular_triangulation rt(wpoints.begin(), wpoints.end());rt.is_valid();Regular_triangulation::Edge_iterator eit;Regular_triangulation::Point_iterator pit;//遍历Delaunay的所有边，绘制Delaunay图的对偶图，即Voronoi图cout << "====all points in rt====" << endl;for (auto vit = rt.vertices_begin(); vit != rt.vertices_end(); ++vit){cout << vit->point() << endl;}cout << "==== Voronoi vert ====" << endl;Rt::Finite_edges_iterator fei;for (fei = rt.finite_edges_begin(); fei != rt.finite_edges_end(); fei++){//遍历每个边auto face = rt.incident_facets(*fei);auto face_start = face;//找到与边相邻的面（更准确的，遍历每个邻接面）do {CGAL::Object  dualedge = rt.dual(*face);//找到面对应的边if (CGAL::object_cast<Rt::Segment>(&dualedge)){Point source = CGAL::object_cast<Rt::Segment>(&dualedge)->source();Point target = CGAL::object_cast<Rt::Segment>(&dualedge)->target();Voronoi_vert.insert(source);Voronoi_vert.insert(target);} else if (CGAL::object_cast<Rt::Ray>(&dualedge))//如果这条边是射线，则绘制射线{Point target = CGAL::object_cast<Rt::Ray>(&dualedge)->source();Voronoi_vert.insert(target);}++face;} while (face_start != face);}}
int main(int argc, char *argv[]) {display();set<Point>::iterator iter=Voronoi_vert.begin();while(iter!=Voronoi_vert.end()){cout<<*iter<<endl;++iter;}return 0;
}

版本二代码

/** 作者:窦志扬* 日期:2018年8月13日* 联络:sdudzy@163.com* 院校:山东大学* */
#include <CGAL/Linear_cell_complex_for_combinatorial_map.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_3_to_lcc.h>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
typedef CGAL::Linear_cell_complex_for_combinatorial_map<3> LCC_3;
typedef LCC_3::Dart_handle                                 Dart_handle;
typedef LCC_3::Point                                          Point;
typedef CGAL::Delaunay_triangulation_3<LCC_3::Traits> Triangulation;
using namespace std;
void remove_the_infinite_volume(LCC_3 &alcc, Dart_handle adart)
{/** 该方法参考cgal 方法* */std::stack<Dart_handle> toremove;LCC_3::size_type mark_toremove=alcc.get_new_mark();toremove.push(adart);CGAL::mark_cell<LCC_3,3>(alcc, adart, mark_toremove);//首先记录需要去除的无限区域for (LCC_3::Dart_of_cell_range<3>::iteratorit=alcc.darts_of_cell<3>(adart).begin(),itend=alcc.darts_of_cell<3>(adart).end(); it!=itend; ++it){if ( !alcc.is_marked(alcc.beta(it,3), mark_toremove) ){CGAL::mark_cell<LCC_3,3>(alcc, alcc.beta(it,3), mark_toremove);toremove.push(alcc.beta(it,3));}}while( !toremove.empty() ){Dart_handle tem = toremove.top();alcc.remove_cell<3>(tem);toremove.pop();}CGAL_assertion(alcc.is_without_boundary(1) && alcc.is_without_boundary(2));std::cout<<"Voronoi subdvision, only finite volumes:"<<std::endl<<"  ";alcc.display_characteristics(std::cout) << ", valid="<< alcc.is_valid()<< std::endl; std::cout<<"Voronoi剖分中所有的有限空间点： "<< endl;for (LCC_3::Vertex_attribute_const_range::iteratorv=alcc.vertex_attributes().begin(),vend=alcc.vertex_attributes().end();v!=vend; ++v)std::cout <<"v "<< alcc.point_of_vertex_attribute(v) << endl;std::cout<<std::endl;
}template<typename LCC, typename TR>
void transform_dart_to_their_dual(LCC& alcc, LCC& adual,std::map<typename TR::Cell_handle,typename LCC::Dart_handle>& assoc)
{typename LCC::Dart_range::iterator it1=alcc.darts().begin();typename LCC::Dart_range::iterator it2=adual.darts().begin();std::map<typename LCC::Dart_handle, typename LCC::Dart_handle> dual;for ( ; it1!=alcc.darts().end(); ++it1, ++it2 ){dual[it1]=it2;}for ( typename std::map<typename TR::Cell_handle, typename LCC::Dart_handle>::iterator it=assoc.begin(), itend=assoc.end(); it!=itend; ++it){assoc[it->first]=dual[it->second];}
}template<typename LCC, typename TR>
void set_geometry_of_dual(LCC& alcc, TR& tr,std::map<typename TR::Cell_handle,typename LCC::Dart_handle>& assoc)
{for ( typename std::map<typename TR::Cell_handle, typename LCC::Dart_handle>::iterator it=assoc.begin(), itend=assoc.end(); it!=itend; ++it){if ( !tr.is_infinite(it->first) )alcc.set_vertex_attribute(it->second,alcc.create_vertex_attribute(tr.dual(it->first)));elsealcc.set_vertex_attribute(it->second,alcc.create_vertex_attribute());}
}int main()
{char*input1 =  "/Users/frankdura/Desktop/222.txt";vector<Point> boundaryPoints;  //装载边界点ifstream in_boundary(input1);Point p;while (in_boundary >> p){   //读入边界点cout << p << endl;boundaryPoints.push_back(p);}cout << "read in vertexs: " << boundaryPoints.size() << endl;in_boundary.close();//请输入您的剖分点Triangulation T;//用于三角剖分T.insert(boundaryPoints.begin(), boundaryPoints.end());LCC_3 lcc;std::map<Triangulation::Cell_handle,LCC_3::Dart_handle > vol_to_dart;Dart_handle dh=CGAL::import_from_triangulation_3<LCC_3, Triangulation>(lcc, T, &vol_to_dart);std::cout<<"Delaunay triangulation :"<<std::endl<<"  ";lcc.display_characteristics(std::cout) << ", valid="<< lcc.is_valid() << std::endl;CGAL::set_ascii_mode(std::cout);std::cout<<"狄洛尼三角剖分顶点: "<< endl;for (LCC_3::Vertex_attribute_const_range::iteratorv=lcc.vertex_attributes().begin(),vend=lcc.vertex_attributes().end();v!=vend; ++v)std::cout << lcc.point_of_vertex_attribute(v) << "\n\r";std::cout<<std::endl;cout <<"end"<< endl;LCC_3 dual_lcc;Dart_handle ddh=lcc.dual(dual_lcc, dh);// Here, dual_lcc is the 3D Voronoi diagram.CGAL_assertion(dual_lcc.is_without_boundary());transform_dart_to_their_dual<LCC_3,Triangulation>(lcc, dual_lcc, vol_to_dart);set_geometry_of_dual<LCC_3,Triangulation>(dual_lcc, T, vol_to_dart);// 5) Display the dual_lcc characteristics.std::cout<<"Voronoi剖分信息 :"<<std::endl<<"  ";dual_lcc.display_characteristics(std::cout) << ", valid="<< dual_lcc.is_valid()<< std::endl;std::cout<<"Voronoi剖分中所有点: "<< endl;for (LCC_3::Vertex_attribute_const_range::iteratorv=dual_lcc.vertex_attributes().begin(),vend=dual_lcc.vertex_attributes().end();v!=vend; ++v)std::cout <<"v " <<dual_lcc.point_of_vertex_attribute(v) << endl;remove_the_infinite_volume(dual_lcc, ddh);return EXIT_SUCCESS;
}

CGAL-三维Delaunay/Voronoi图/Powerdiagram相关推荐

Voronoi图 | 泰森多边形
概念 Voronoi图,又叫泰森多边形.沃罗诺伊图或Dirichlet图,是一个关于空间划分的基础数据结构. 它是由一组由连接两邻点直线的垂直平分线组成的连续多边形组成. N个在平面上有区别的点,按照 ...
TIN的构建、Delaunay三角网、Voronoi图
一.TIN的三角剖分准则 (1)空外接圆准则过每个三角形的外接圆均不包含点集的其余任何点. (2)最大最小角准则两三角形中的最小内角>交换z凸四边形对角线后三角形的最小角. (3)最短距离和 ...
如何使用OpenCV进行Delaunay三角剖分和Voronoi图
图1.左:使用dlib检测到具有标志性建筑的奥巴马总统图像.中心:地标的Delaunay三角剖分.右:对应的Voronoi图. 俄国数学家鲍里斯·尼古拉耶维奇·德劳内(Boris Nikolaevic ...
基于C++(MFC)的二维Delaunay三角剖分与Voronoi图的算法及代码
一. Delaunay三角网 Delaunay三角网的特性: (1)空圆性,任一三角形外接圆内部不包含其他点. (2)最接近:以最近临的三点形成三角形,且各线段(三角形的边)皆不相交. (3)唯一性: ...
数字图像处理 Delaunay三角剖分和Voronoi图
一.什么是Delaunay三角剖分? 给定平面中的一组点,三角测量是指将平面细分为三角形,以这些点为顶点.在下图1中,我们在左图像中看到一组地标,在中间图像中看到三角剖分.一组点可以有许多可能的三角剖 ...
Voronoi图（泰森多边形）和Delaunay三角形
荷兰气候学家A.H.Thiessen提出了一个根据离散分布的气象站的降雨量来计算平均降雨量的方法,即将所有的相邻气象站连成三角形,做这些三角形的各边的垂直平分线,于是每个气象站的周围的 ...
加权voronoi图 matlab,加权voronoi图matlab
广义Voronoi图的快速生成算法_电力/水利_工程科技_专业资料.27 卷第 ... Voronoi图理论与应用新成... 3页免费 Voronoi图理论与应用新成... 3页免费 Vorono ...
泰森多边形（Voronoi图）
二维Delaunay(德洛内)三角网剖分的matlab实现 https://blog.csdn.net/weixin_42943114/article/details/82262122 泰森多边形(V ...
泰森多边形(Voronoi图)生成算法
一.文档目的本文描述了在geomodel模块中,生成泰森多边形所使用的算法. 二.概述 GIS和地理分析中经常采用泰森多边形进行快速插值,和分析地理实体的影响区域,是解决邻接度问题的又一常用工具. ...

CGAL-三维Delaunay/Voronoi图/Powerdiagram

综述

说明

更新

环境

代码

版本一代码

版本二代码

CGAL-三维Delaunay/Voronoi图/Powerdiagram相关推荐

最新文章

热门文章