版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/OOFFrankDura/article/details/81665905
综述
注意到cgal中的dual没有点到面的对偶,所以我们需要自己遍历完成。
这里使用Edge_Circulator。需要注意的是,Circulator要求边闭合(无射线边)。所以我们使用一个包围盒进行框定。
环境
cgal 4.12
IDE clion
项目管理 cmake
macos
更新
2018.8.14给出第二种计算方案:
使用cgal本身
finite_edge_iterator
的特性.
说明
这是使用powerdiagram进行展示(正则三角剖分后求对偶),如果您希望得到Voronoi可以使用Delaunay三角剖分后做相同的处理。
效果
其中白色区域就是点——绿色点的包围边。
代码
第一种
//作者:SDU窦志扬
//日期:2018/8/14
//联络:[email protected]
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/Weighted_point_2.h>
#include <CGAL/Polygon_2.h>
#include <GLUT/glut.h>
#include <CGAL/bounding_box.h>
#include <fstream>
#include<iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Regular_triangulation_2<K> Regular_triangulation;
typedef K::Vector_2 Vector;
typedef K::Point_2 Point;
typedef K::Weighted_point_2 Wp;
typedef K::Point_2 Point;
typedef CGAL::Regular_triangulation_2<K> Rt;
int window_size = 600;
using namespace std;
std::vector<Regular_triangulation::Weighted_point> wpoints;
std::vector<Point> points;
const int n=100;
const GLfloat Pi=3.1415926536f;
#define MAX_CHAR 128
void MyDisplay(double x,double y ,double R)
{
glBegin(GL_POLYGON);
glColor3f(1,0,0);
R =sqrt(R);
for(int i=0;i<n;i++){
glVertex2f(R*cos(2*Pi/n*i)+x,R*sin(2*Pi/n*i)+y);
}
glEnd();
}
void add_point(int x, int y, double w){
//绘制圆形
wpoints.push_back( Wp(Point(x,y), w) );
points.push_back(Point(x,y) );
MyDisplay(x,y,w);
glPointSize(7);
glBegin(GL_POINTS);
glColor3f( 0.0, 1.0, 0.0 );
glVertex2f( x,y);
glEnd();
}
void display(void)
{
add_point(100, 400,10000); //左中
add_point(200,200,15000);
add_point(240,310,8000); //目标点
add_point(210, 560,20000);
add_point(340, 80, 12000);
add_point(560,120, 10000);
add_point(600, 500,12002);
add_point(450, 550,8000);
Regular_triangulation rt(wpoints.begin(), wpoints.end());
rt.is_valid();
Regular_triangulation::Edge_iterator eit;
Regular_triangulation::Point_iterator pit;
Regular_triangulation::Finite_faces_iterator fit;
//遍历Delaunay的所有边,绘制Delaunay图的对偶图,即Voronoi图
glEnable( GL_LINE_STIPPLE );//使用点画模式,即使用虚线来绘制Voronoi图
glLineStipple( 1, 0x3333 );
glColor3f( 0.0, 1.0,1.0 );
cout << "====all points in rt====" << endl;
for (auto vit = rt.vertices_begin(); vit != rt.vertices_end(); ++vit){
cout << vit->point() << endl;
}
cout << "========" << endl;
vector<Point> dual_points;
cout << "====all points in powerdiagram====" << endl;
//遍历Regular triangulation的所有面
for( fit = rt.faces_begin(); fit != rt.faces_end(); fit ++)
{
dual_points.push_back(rt.dual(fit));
//在其对偶图中所对应的点
}
for (int j = 0; j <dual_points.size() ; ++j) {
cout << dual_points[j] << endl;
}
cout << "========" << endl;
for( eit = rt.edges_begin(); eit != rt.edges_end(); eit ++)
{
CGAL::Object o = rt.dual(eit);
//在其对偶图中所对应的边
if (CGAL::object_cast<K::Segment_2>(&o)) //如果这条边是线段,则绘制线段
{
glBegin(GL_LINES);
glVertex2i( CGAL::object_cast<K::Segment_2>(&o)->source().hx(), CGAL::object_cast<K::Segment_2>(&o)->source().hy() );
glVertex2i( CGAL::object_cast<K::Segment_2>(&o)->target().hx(), CGAL::object_cast<K::Segment_2>(&o)->target().hy() );
glEnd();
}
else if (CGAL::object_cast<K::Ray_2>(&o))//如果这条边是射线,则绘制射线
{
glBegin(GL_LINES);
glVertex2i( CGAL::object_cast<K::Ray_2>(&o)->source().hx(), CGAL::object_cast<K::Ray_2>(&o)->source().hy() );
glVertex2i( CGAL::object_cast<K::Ray_2>(&o)->point(1).hx(), CGAL::object_cast<K::Ray_2>(&o)->point(1).hy() );
glEnd();
}
}
glDisable( GL_LINE_STIPPLE );//关闭点画模式
auto box = CGAL::bounding_box(points.begin(), points.end());
auto center2 = CGAL::ORIGIN + 0.5 * (box.max() - CGAL::ORIGIN) + 0.5 * (box.min() - CGAL::ORIGIN);
double sideLen = max(box.xmax() - box.xmin(), box.ymax() - box.ymin());
Point center = CGAL::ORIGIN;
double maxW = -100000;
map<Point,int> fromSite2ID;
for (int i = 0; i < wpoints.size(); ++i)
{
//fromSite2ID是一个点Point到int的map
//这里对每个点进行标记"i"
fromSite2ID[wpoints[i].point()] = i;
center = center + (wpoints[i].point() - CGAL::ORIGIN);
if (maxW < wpoints[i].weight())
{//负值越大说明
maxW = wpoints[i].weight();
}
}
center = CGAL::ORIGIN + (center - CGAL::ORIGIN) / (double)wpoints.size();
//center计算距离原点的平均偏移量
vector<Wp> exteriorPoints;
Wp p1(center + sideLen * Vector(1, 0), maxW);
exteriorPoints.push_back(p1);
Wp p2(center + sideLen * Vector(-1, 0), maxW);
exteriorPoints.push_back(p2);
Wp p3(center + sideLen * Vector(0, 1), maxW);
exteriorPoints.push_back(p3);
Wp p4(center + sideLen * Vector(0,-1), maxW);
exteriorPoints.push_back(p4);
//用一个更大的区域包围起来(一定包在里面)
rt.insert(exteriorPoints.begin(), exteriorPoints.end());
for (int i = 0; i < exteriorPoints.size(); ++i)
{
fromSite2ID[exteriorPoints[i].point()] = -1; //means the extra sites
}
//遍历显示
cout << "====all edges bounding the point====" << endl;
int con = 0;
for (Rt::Finite_vertices_iterator siteIt = rt.finite_vertices_begin();
siteIt != rt.finite_vertices_end(); ++siteIt)
{
con++;
Point site = siteIt->point().point();
// cout <<"powerdiagram :" << con++<<"--" <<site << endl;
int curSiteID = fromSite2ID.find(site)->second;//找到点对应的id
if (curSiteID == -1){
cout << "这是边界点" << endl;
continue;
}
Rt::Edge_circulator beg = rt.incident_edges(siteIt);
Rt::Edge_circulator it = beg;
vector<Point> cell;
do
{
K::Segment_2 s;
CGAL::Object o = rt.dual(it);
if (CGAL::assign(s, o))
{
cell.push_back(s.source());
}
} while (++it != beg);
glLineWidth(8);
glBegin(GL_LINE_LOOP);
glLineWidth(20);
glColor3f(1,1,1);
for (int j = 0; j <cell.size() ; ++j) {
if (site == Point(240,310)){
glVertex2f(cell[j].x(), cell[j].y());
cout << cell[j] << " --- " << site << endl;
}
}
glEnd();
}
}
void InitEnvironment()
{
//一些初始化操作
glClearColor(0.0,0.0,0.0,0);
glClear(GL_COLOR_BUFFER_BIT);
glPointSize(7);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluOrtho2D(0,window_size,0,window_size);
}
void init_Display(void)
{
glClear(GL_COLOR_BUFFER_BIT);
glColor3f(0.98f, 0.625f, 0.12f);
display();
glFlush();
}
int main(int argc, char *argv[])
{
glutInit(&argc, argv); //初始化GLUT
glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE);
glutInitWindowPosition(300, 100);
glutInitWindowSize(window_size, window_size);
glutCreateWindow("RT");
InitEnvironment(); //初始化
glutDisplayFunc(&init_Display); //回调函数
glutMainLoop(); //持续显示,当窗口改变会重新绘制图形
return 0;
}
第二种
//作者:SDU窦志扬
//日期:2018/8/14
//联络:[email protected]
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/Weighted_point_2.h>
#include <CGAL/Polygon_2.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_2<K> Regular_triangulation;
typedef K::Vector_2 Vector;
typedef K::Point_2 Point;
typedef K::Weighted_point_2 Wp;
typedef K::Point_2 Point;
typedef CGAL::Regular_triangulation_2<K> Rt;
int window_size = 600;
using namespace std;
std::vector<Regular_triangulation::Weighted_point> wpoints;
std::vector<Point> points;
std::set<Point> Point_set;
const int n=100;
const GLfloat Pi=3.1415926536f;
#define MAX_CHAR 128
void MyDisplay(double x,double y ,double R)
{
glBegin(GL_POLYGON);
glColor3f(1,0,0);
R =sqrt(R);
for(int i=0;i<n;i++){
glVertex2f(R*cos(2*Pi/n*i)+x,R*sin(2*Pi/n*i)+y);
}
glEnd();
}
void add_point(int x, int y, double w){
//绘制圆形
wpoints.push_back( Wp(Point(x,y), w) );
points.push_back(Point(x,y) );
Point_set.insert(Point(x,y));
MyDisplay(x,y,w);
glPointSize(7);
glBegin(GL_POINTS);
glColor3f( 0.0, 1.0, 0.0 );
glVertex2f( x,y);
glEnd();
}
void display(void)
{
add_point(100, 400,10000); //左中
add_point(200,200,15000);
add_point(240,310,8000); //目标点
add_point(210, 560,20000);
add_point(340, 80, 12000);
add_point(560,120, 10000);
add_point(600, 500,12002);
add_point(450, 550,8000);
Regular_triangulation rt(wpoints.begin(), wpoints.end());
rt.is_valid();
Regular_triangulation::Finite_edges_iterator eit;
Regular_triangulation::Point_iterator pit;
Regular_triangulation::Finite_faces_iterator fit;
//遍历Delaunay的所有边,绘制Delaunay图的对偶图,即Voronoi图
glEnable( GL_LINE_STIPPLE );//使用点画模式,即使用虚线来绘制Voronoi图
glLineStipple( 1, 0x3333 );
glColor3f( 0.0, 1.0,1.0 );
cout << "====all points in rt====" << endl;
for (auto vit = rt.vertices_begin(); vit != rt.vertices_end(); ++vit){
cout << vit->point() << endl;
}
cout << "========" << endl;
vector<Point> dual_points;
cout << "====all points in powerdiagram====" << endl;
//遍历Regular triangulation的所有面
for( fit = rt.faces_begin(); fit != rt.faces_end(); fit ++)
{
dual_points.push_back(rt.dual(fit));
//在其对偶图中所对应的点
}
for (int j = 0; j <dual_points.size() ; ++j) {
cout << dual_points[j] << endl;
}
cout << "========" << endl;
for( eit = rt.finite_edges_begin(); eit != rt.finite_edges_end(); eit ++)
{
cout << eit->second << endl;
// cout << "processing seg"<< eit->first->vertex(eit->second)->point().point() << endl;
// cout << "processing seg2"<< eit->first->vertex(0)->point().point() << endl;
CGAL::Object o = rt.dual(eit);
//在其对偶图中所对应的边
set<Point>::iterator it;
Point a = (eit->first->vertex(0)->point().point() );
Point b = (eit->first->vertex(1)->point().point() );
it=Point_set.find(a); //查找键值为5的元素
set<Point>::iterator it2;
it2=Point_set.find(eit->first->vertex(1)->point().point()); //查找键值为5的元素
if(it!=Point_set.end()&&it2!=Point_set.end()&& (a==Point(240,310)||Point(240,310)==b)){
// cout << "a" << a<<endl;
// cout << "b" << b<< endl;
if (CGAL::object_cast<K::Segment_2>(&o)) //如果这条边是线段,则绘制线段
{
glBegin(GL_LINES);
glVertex2i( CGAL::object_cast<K::Segment_2>(&o)->source().hx(), CGAL::object_cast<K::Segment_2>(&o)->source().hy() );
glVertex2i( CGAL::object_cast<K::Segment_2>(&o)->target().hx(), CGAL::object_cast<K::Segment_2>(&o)->target().hy() );
glEnd();
}
}
}
glDisable( GL_LINE_STIPPLE );//关闭点画模式
auto box = CGAL::bounding_box(points.begin(), points.end());
auto center2 = CGAL::ORIGIN + 0.5 * (box.max() - CGAL::ORIGIN) + 0.5 * (box.min() - CGAL::ORIGIN);
double sideLen = max(box.xmax() - box.xmin(), box.ymax() - box.ymin());
Point center = CGAL::ORIGIN;
double maxW = -100000;
map<Point,int> fromSite2ID;
for (int i = 0; i < wpoints.size(); ++i)
{
//fromSite2ID是一个点Point到int的map
//这里对每个点进行标记"i"
fromSite2ID[wpoints[i].point()] = i;
center = center + (wpoints[i].point() - CGAL::ORIGIN);
if (maxW < wpoints[i].weight())
{//负值越大说明
maxW = wpoints[i].weight();
}
}
center = CGAL::ORIGIN + (center - CGAL::ORIGIN) / (double)wpoints.size();
//center计算距离原点的平均偏移量
vector<Wp> exteriorPoints;
Wp p1(center + sideLen * Vector(1, 0), maxW);
exteriorPoints.push_back(p1);
Wp p2(center + sideLen * Vector(-1, 0), maxW);
exteriorPoints.push_back(p2);
Wp p3(center + sideLen * Vector(0, 1), maxW);
exteriorPoints.push_back(p3);
Wp p4(center + sideLen * Vector(0,-1), maxW);
exteriorPoints.push_back(p4);
//用一个更大的区域包围起来(一定包在里面)
rt.insert(exteriorPoints.begin(), exteriorPoints.end());
for (int i = 0; i < exteriorPoints.size(); ++i)
{
fromSite2ID[exteriorPoints[i].point()] = -1; //means the extra sites
}
//遍历显示
cout << "====all edges bounding the point====" << endl;
int con = 0;
for (Rt::Finite_vertices_iterator siteIt = rt.finite_vertices_begin();
siteIt != rt.finite_vertices_end(); ++siteIt)
{
con++;
Point site = siteIt->point().point();
// cout <<"powerdiagram :" << con++<<"--" <<site << endl;
int curSiteID = fromSite2ID.find(site)->second;//找到点对应的id
if (curSiteID == -1){
cout << "这是边界点" << endl;
continue;
}
Rt::Edge_circulator beg = rt.incident_edges(siteIt);
Rt::Edge_circulator it = beg;
vector<Point> cell;
do
{
K::Segment_2 s;
CGAL::Object o = rt.dual(it);
if (CGAL::assign(s, o))
{
cell.push_back(s.source());
}
} while (++it != beg);
glLineWidth(8);
// glBegin(GL_LINE_LOOP);
// glLineWidth(20);
// glColor3f(1,1,1);
// for (int j = 0; j <cell.size() ; ++j) {
// if (site == Point(240,310)){
// glVertex2f(cell[j].x(), cell[j].y());
// cout << cell[j] << " --- " << site << endl;
// }
// }
// glEnd();
}
}
void InitEnvironment()
{
//一些初始化操作
glClearColor(0.0,0.0,0.0,0);
glClear(GL_COLOR_BUFFER_BIT);
glPointSize(7);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluOrtho2D(0,window_size,0,window_size);
}
void init_Display(void)
{
glClear(GL_COLOR_BUFFER_BIT);
glColor3f(0.98f, 0.625f, 0.12f);
display();
glFlush();
}
int main(int argc, char *argv[])
{
glutInit(&argc, argv); //初始化GLUT
glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE);
glutInitWindowPosition(300, 100);
glutInitWindowSize(window_size, window_size);
glutCreateWindow("RT");
InitEnvironment(); //初始化
glutDisplayFunc(&init_Display); //回调函数
glutMainLoop(); //持续显示,当窗口改变会重新绘制图形
return 0;
}
往期实现
CGAL单纯的2D Power diagram样例展示(展示暂时使用openGL):
CGAL 2D Power diagram
CGAL单纯的2D Voronoi diagram样例展示:
CGAL 2F Voronoi diagram
CGAL单纯的3D Voronoi diagram样例展示:
CGAL 3D Voronoi