文章目录

  • Delaunay 三角剖分
  • 代码实现< C++版本 >
  • 代码实现< Python 版本 >
  • 代码实现< Matlab 版本 >
  • 彩蛋：一个很牛掰的算法学习网站

Delaunay 三角剖分

  • 什么是三角剖分(推荐看这个作者的系列文章，写的相当好)
  • 什么是 Delaunay 三角剖分
    • 空圆特性
    • 最大化最小角特性
  • Delaunay 三角剖分的特点
    • 最接近
    • 唯一性
    • 最优性
    • 最规则
    • 区域性
    • 具有凸多边形的外壳
  • 有什么算法
    • Bowyer-Watson算法 — 逐点插入法
    • Lawson算法
  • 相关博客
    • [图形学]Delaunay三角剖分算法附C++实现
    • [图形学]凸包生成算法附C++实现
    • 笔记：Delaunay三角剖分(Delaunay Triangulation)相关知识
    • Deluanay三角剖分与Voronoi图
    • 强烈推荐这篇：OpenCV——Delaunay三角剖分（C++实现）
    • …

代码实现< C++版本 >

代码框架

  • 定义点类 vector2.h
  • 定义边类 edge.h
  • 定义三角形类 triangle.h
  • 定义三角剖分类 delaunay.h
  • 定义比较类(用于三角退化) numeric.h
  • 定义main函数 main.cpp

用到的外部库

  • SFML库

具体实现

vector2.h

#ifndef H_VECTOR2
#define H_VECTOR2#include "numeric.h"#include <iostream>
#include <cmath>template <typename T>
class Vector2
{public:// Constructors  构造函数Vector2():x(0), y(0){}Vector2(T _x, T _y): x(_x), y(_y){}Vector2(const Vector2 &v): x(v.x), y(v.y){}// Operations  // 计算距离T dist2(const Vector2 &v) const{T dx = x - v.x;T dy = y - v.y;return dx * dx + dy * dy;}T dist(const Vector2 &v) const{return sqrt(dist2(v));}//计算平方和，此函数在 delaunay.h 中判断一点是否在三角形内部T norm2() const{return x * x + y * y;}T x;T y;
};//全特化
template <>
float Vector2<float>::dist(const Vector2<float> &v) const { return hypotf(x - v.x, y - v.y);}template <>
double Vector2<double>::dist(const Vector2<double> &v) const { return hypotf(x - v.x, y - v.y);}template<typename T>
std::ostream &operator << (std::ostream &str, Vector2<T> const &point)
{return str << "Point x: " << point.x << " y: " << point.y;
}template<typename T>
bool operator == (const Vector2<T>& v1, const Vector2<T>& v2)
{return (v1.x == v2.x) && (v1.y == v2.y);
}template<class T>
typename std::enable_if<!std::numeric_limits<T>::is_integer, bool>::typealmost_equal(const Vector2<T>& v1, const Vector2<T>& v2, int ulp=2)
{return almost_equal(v1.x, v2.x, ulp) && almost_equal(v1.y, v2.y, ulp);
}#endif

edge.h

#ifndef H_EDGE
#define H_EDGE#include "vector2.h"template <class T>
class Edge
{public:using VertexType = Vector2<T>;  //相当于 typedef//构造函数Edge(const VertexType &p1, const VertexType &p2) : p1(p1), p2(p2), isBad(false) {}Edge(const Edge &e) : p1(e.p1), p2(e.p2), isBad(false) {}VertexType p1;VertexType p2;bool isBad;
};template <class T>
inline std::ostream &operator << (std::ostream &str, Edge<T> const &e)
{return str << "Edge " << e.p1 << ", " << e.p2;
}template <class T>
inline bool operator == (const Edge<T> & e1, const Edge<T> & e2)
{return     (e1.p1 == e2.p1 && e1.p2 == e2.p2) ||(e1.p1 == e2.p2 && e1.p2 == e2.p1);
}template <class T>
inline bool almost_equal (const Edge<T> & e1, const Edge<T> & e2)
{return (almost_equal(e1.p1, e2.p1) && almost_equal(e1.p2, e2.p2)) ||(almost_equal(e1.p1, e2.p2) && almost_equal(e1.p2, e2.p1));
}#endif

triangle.h

#ifndef H_TRIANGLE
#define H_TRIANGLE#include "vector2.h"
#include "edge.h"
#include "numeric.h"template <class T>
class Triangle
{public:using EdgeType = Edge<T>;using VertexType = Vector2<T>;//构造函数Triangle(const VertexType &_p1, const VertexType &_p2, const VertexType &_p3):   p1(_p1), p2(_p2), p3(_p3),e1(_p1, _p2), e2(_p2, _p3), e3(_p3, _p1), isBad(false){}//判断一点是否在三角形内部bool containsVertex(const VertexType &v) const{// return p1 == v || p2 == v || p3 == v;return almost_equal(p1, v) || almost_equal(p2, v) || almost_equal(p3, v);}//判断一点是否在外接圆内（包括在外接圆上）bool circumCircleContains(const VertexType &v) const{const T ab = p1.norm2();const T cd = p2.norm2();const T ef = p3.norm2();const T circum_x = (ab * (p3.y - p2.y) + cd * (p1.y - p3.y) + ef * (p2.y - p1.y)) / (p1.x * (p3.y - p2.y) + p2.x * (p1.y - p3.y) + p3.x * (p2.y - p1.y));const T circum_y = (ab * (p3.x - p2.x) + cd * (p1.x - p3.x) + ef * (p2.x - p1.x)) / (p1.y * (p3.x - p2.x) + p2.y * (p1.x - p3.x) + p3.y * (p2.x - p1.x));const VertexType circum(half(circum_x), half(circum_y));const T circum_radius = p1.dist2(circum);const T dist = v.dist2(circum);return dist <= circum_radius;}VertexType p1;VertexType p2;VertexType p3;EdgeType e1;EdgeType e2;EdgeType e3;bool isBad;
};template <class T>
inline std::ostream &operator << (std::ostream &str, const Triangle<T> & t)
{return str << "Triangle:" << std::endl << "\t" << t.p1 << std::endl << "\t" << t.p2 << std::endl << "\t" << t.p3 << std::endl << "\t" << t.e1 << std::endl << "\t" << t.e2 << std::endl << "\t" << t.e3 << std::endl;}template <class T>
inline bool operator == (const Triangle<T> &t1, const Triangle<T> &t2)
{return (t1.p1 == t2.p1 || t1.p1 == t2.p2 || t1.p1 == t2.p3) &&(t1.p2 == t2.p1 || t1.p2 == t2.p2 || t1.p2 == t2.p3) &&(t1.p3 == t2.p1 || t1.p3 == t2.p2 || t1.p3 == t2.p3);
}//注意这里的 almost 函数，该函数在 numeric.h 中定义
template <class T>
inline bool almost_equal(const Triangle<T> &t1, const Triangle<T> &t2)
{return (almost_equal(t1.p1 , t2.p1) || almost_equal(t1.p1 , t2.p2) || almost_equal(t1.p1 , t2.p3)) &&(almost_equal(t1.p2 , t2.p1) || almost_equal(t1.p2 , t2.p2) || almost_equal(t1.p2 , t2.p3)) &&(almost_equal(t1.p3 , t2.p1) || almost_equal(t1.p3 , t2.p2) || almost_equal(t1.p3 , t2.p3));
}#endif

delaunay.h

#ifndef H_DELAUNAY
#define H_DELAUNAY#include "vector2.h"
#include "edge.h"
#include "triangle.h"#include <vector>
#include <algorithm>template <class T>
class Delaunay
{public:using TriangleType = Triangle<T>;using EdgeType = Edge<T>;using VertexType = Vector2<T>;//Deluanay 三角剖分核心算法  ---  逐点插入法const std::vector<TriangleType>& triangulate(std::vector<VertexType> &vertices){// 将点云拷贝一份_vertices = vertices;// 计算超级三角形的一些参数T minX = vertices[0].x;T minY = vertices[0].y;T maxX = minX;T maxY = minY;//计算超级三角形的上下左右边界for(std::size_t i = 0; i < vertices.size(); ++i){if (vertices[i].x < minX) minX = vertices[i].x;if (vertices[i].y < minY) minY = vertices[i].y;if (vertices[i].x > maxX) maxX = vertices[i].x;if (vertices[i].y > maxY) maxY = vertices[i].y;}const T dx = maxX - minX;const T dy = maxY - minY;const T deltaMax = std::max(dx, dy);const T midx = half(minX + maxX);const T midy = half(minY + maxY);//尽可能扩大超级三角形范围，可以取比20更大的数字const VertexType p1(midx - 20 * deltaMax, midy - deltaMax);const VertexType p2(midx, midy + 20 * deltaMax);const VertexType p3(midx + 20 * deltaMax, midy - deltaMax);//std::cout << "Super triangle " << std::endl << Triangle(p1, p2, p3) << std::endl;// Create a list of triangles, and add the supertriangle in it//将超级三角形添加至 _triangles_triangles.push_back(TriangleType(p1, p2, p3));//开始依次遍历每个点for(auto p = begin(vertices); p != end(vertices); p++){//std::cout << "Traitement du point " << *p << std::endl;//std::cout << "_triangles contains " << _triangles.size() << " elements" << std::endl;//构造变量用来存储临时新产生的边std::vector<EdgeType> polygon;//依次遍历 _triangles  中的每个三角形for(auto & t : _triangles){//std::cout << "Processing " << std::endl << *t << std::endl;if(t.circumCircleContains(*p))  //如果包含点 p，那么就要产生新的3条边{//std::cout << "Pushing bad triangle " << *t << std::endl;t.isBad = true;  //flag 发生改变，准备接下来在 _triangles  中将其剔除polygon.push_back(t.e1);polygon.push_back(t.e2);polygon.push_back(t.e3);}else{//std::cout << " does not contains " << *p << " in his circum center" << std::endl;}}//更新 _triangles_triangles.erase(std::remove_if(begin(_triangles), end(_triangles), [](TriangleType &t){return t.isBad;}), end(_triangles));   //std::remove_if只移动不删除，erase将其删除，这里还用到了C++11的 lambda 表达式for(auto e1 = begin(polygon); e1 != end(polygon); ++e1)  //这个用来删除重复的边{for(auto e2 = e1 + 1; e2 != end(polygon); ++e2){if(almost_equal(*e1, *e2)){e1->isBad = true;e2->isBad = true;}}}//更新 polygonpolygon.erase(std::remove_if(begin(polygon), end(polygon), [](EdgeType &e){return e.isBad;}), end(polygon));for(const auto e : polygon)_triangles.push_back(TriangleType(e.p1, e.p2, *p));}//删除超级三角形_triangles.erase(std::remove_if(begin(_triangles), end(_triangles), [p1, p2, p3](TriangleType &t){return t.containsVertex(p1) || t.containsVertex(p2) || t.containsVertex(p3);}), end(_triangles));for(const auto t : _triangles){_edges.push_back(t.e1);_edges.push_back(t.e2);_edges.push_back(t.e3);}return _triangles;}const std::vector<TriangleType>& getTriangles() const { return _triangles; }const std::vector<EdgeType>& getEdges() const { return _edges; }const std::vector<VertexType>& getVertices() const { return _vertices; }private:std::vector<TriangleType> _triangles;std::vector<EdgeType> _edges;std::vector<VertexType> _vertices;
};
#endif

numeric.h

#ifndef H_NUMERIC
#define H_NUMERIC#include <cmath>
#include <limits>/*** @brief use of machine epsilon to compare floating-point values for equality* http://en.cppreference.com/w/cpp/types/numeric_limits/epsilon*///其目的就是用来比较数据是否相等，变相的 float1 - float2  < 0.000001
template<class T>
typename std::enable_if<!std::numeric_limits<T>::is_integer, bool>::typealmost_equal(T x, T y, int ulp=2)
{// the machine epsilon has to be scaled to the magnitude of the values used// and multiplied by the desired precision in ULPs (units in the last place)return std::abs(x-y) <= std::numeric_limits<T>::epsilon() * std::abs(x+y) * ulp// unless the result is subnormal|| std::abs(x-y) < std::numeric_limits<T>::min();
}template<class T>
T half(T x){}template <>
float half(float x){return 0.5f * x;}template <>
double half(double x){return 0.5 * x;}#endif

main.cpp

#include <iostream>
#include <vector>
#include <iterator>
#include <algorithm>
#include <array>#include <SFML/Graphics.hpp>   //注意这是外部库//包含头文件
#include "vector2.h"
#include "triangle.h"
#include "delaunay.h"//随机生成二维点
float RandomFloat(float a, float b) {const float random = ((float) rand()) / (float) RAND_MAX;const float diff = b - a;const float r = random * diff;return a + r;
}int main(int argc, char * argv[])
{int numberPoints = 40;/*if (argc==1){numberPoints = (int) roundf(RandomFloat(4, numberPoints));}else if (argc>1){numberPoints = atoi(argv[1]);}*/srand (time(NULL));std::cout << "Generating " << numberPoints << " random points" << std::endl;std::vector<Vector2<float> > points;for(int i = 0; i < numberPoints; ++i) {points.push_back(Vector2<float>(RandomFloat(0, 800), RandomFloat(0, 600)));}Delaunay<float> triangulation;const std::vector<Triangle<float> > triangles = triangulation.triangulate(points);std::cout << triangles.size() << " triangles generated\n";const std::vector<Edge<float> > edges = triangulation.getEdges();std::cout << " ========= ";std::cout << "\nPoints : " << points.size() << std::endl;for(const auto &p : points)std::cout << p << std::endl;std::cout << "\nTriangles : " << triangles.size() << std::endl;for(const auto &t : triangles)std::cout << t << std::endl;std::cout << "\nEdges : " << edges.size() << std::endl;for(const auto &e : edges)std::cout << e << std::endl;// SFML windowsf::RenderWindow window(sf::VideoMode(800, 600), "Delaunay triangulation");// Transform each points of each vector as a rectanglestd::vector<sf::RectangleShape*> squares;for(const auto p : points) {sf::RectangleShape *c1 = new sf::RectangleShape(sf::Vector2f(4, 4));c1->setPosition(p.x, p.y);squares.push_back(c1);}// Make the linesstd::vector<std::array<sf::Vertex, 2> > lines;for(const auto &e : edges) {lines.push_back({{sf::Vertex(sf::Vector2f(e.p1.x + 2, e.p1.y + 2)),sf::Vertex(sf::Vector2f(e.p2.x + 2, e.p2.y + 2))}});}while (window.isOpen()){sf::Event event;while (window.pollEvent(event)){if (event.type == sf::Event::Closed)window.close();}window.clear();// Draw the squaresfor(const auto &s : squares) {window.draw(*s);}// Draw the linesfor(const auto &l : lines) {window.draw(l.data(), 5, sf::Lines);}window.display();}return 0;
}

效果图：

代码实现< Python 版本 >

原文链接：https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.Delaunay.html

import numpy as np
from scipy.spatial import Delaunay
import matplotlib.pyplot as plt# Triangulation of a set of points:
points = np.array([[0, 0], [0, 1.1], [1, 0.3], [1, 1]]) # 定义三角点
tri = Delaunay(points) # 三角剖分# We can plot it:
plt.triplot(points[:,0], points[:,1], tri.simplices.copy())
plt.plot(points[:,0], points[:,1], 'o')
plt.show()# Point indices and coordinates for the two triangles forming the triangulation:
print(tri.simplices) # 每个三角面对应的点的索引index
print(points[tri.simplices]) # 每个三角面所包含的坐标点# Triangle 0 is the only neighbor of triangle 1, and it’s opposite to vertex 1 of triangle 1:
print(tri.neighbors[1]) # 第一个三角面周围有几个邻居三角形，这里只有 1 个
print(points[tri.simplices[1,0]]) # 第 1 个三角面的 X 坐标
print(points[tri.simplices[1,1]]) # 第 1 个三角面的 Y 坐标
print(points[tri.simplices[1,2]]) # 第 1 个三角面的 Z 坐标# We can find out which triangle points are in:
p = np.array([(0.1, 0.2), (1.5, 0.5)]) # 判断两个点是都在三角网内部
print(tri.find_simplex(p))# We can also compute barycentric(重心) coordinates in triangle 1 for these points:
b = tri.transform[1,:2].dot(p - tri.transform[1,2])
print(tri)
print(np.c_[b, 1 - b.sum(axis=1)])

效果图：

代码实现< Matlab 版本 >

原文链接：https://www.mathworks.com/help/matlab/ref/delaunaytriangulation.html

%%
close all
number = 20;
x = 10*rand(1,number);
y = 10*rand(1,number);tri = delaunay(x,y);figure
hold on
plot(x, y, 'r*')for ii = 1:size(tri, 1)plot( [x(tri(ii,1)) x(tri(ii,2))], [y(tri(ii,1)) y(tri(ii,2))], 'b' )plot( [x(tri(ii,2)) x(tri(ii,3))], [y(tri(ii,2)) y(tri(ii,3))], 'b' )plot( [x(tri(ii,1)) x(tri(ii,3))], [y(tri(ii,1)) y(tri(ii,3))], 'b' )
end
set(gca, 'box', 'on')
print(gcf,'-dpng','delaunary.png')

效果图：

彩蛋：一个很牛掰的算法学习网站

演算法筆記

Delaunay 三角剖分2D（原理 + 源码）相关推荐

  1. Delaunay 三角剖分3D（原理 + 源码）

    文章目录 Delaunay 三角剖分 二维与三维的区别 代码实现< Python版本 > 代码实现< Matlab 版本 > Delaunay 三角剖分 原理部分请参照我的上一 ...

  2. k线顶分型 python_顶底分型－（K线分类及顶底分型的一种数学原理 源码 贴图）...

    好股票软件下载网(www.goodgupiao.com)提示:您正在下载的是:顶底分型-(K线分类及顶底分型的一种数学原理 源码 贴图) 参考缠论,研究了很多天终于将顶底分型进行了具体的数学量化,涵盖 ...

  3. Unity Fog 原理 源码分析 案例

    Unity Fog 原理 源码分析 案例 效果图 简述 背景知识 clip空间坐标的范围 d3d (near,0), unity remapping to (0,far) 越靠近相机z值越小 open ...

  4. 仿猎豹垃圾清理 实现原理+源码

    仿猎豹垃圾清理(实现原理+源码) 转载请注明出处: 仿猎豹垃圾清理(实现原理+源码) 前几天无意打开猎豹内存大师, 发现它的垃圾清理很强大, 效果也不错, 闲着就研究了下. 不过.. 结果貌似和我想象 ...

  5. 仿猎豹垃圾清理(实现原理+源码)

    仿猎豹垃圾清理(实现原理+源码) 转载请注明出处: 仿猎豹垃圾清理(实现原理+源码) 前几天无意打开猎豹内存大师, 发现它的垃圾清理很强大, 效果也不错, 闲着就研究了下. 不过.. 结果貌似和我想象 ...

  6. 仿iOS猎豹垃圾清理(实现原理+源码)

    转载请注明出处: 仿猎豹垃圾清理(实现原理+源码) 前几天无意打开猎豹内存大师, 发现它的垃圾清理很强大, 效果也不错, 闲着就研究了下. 不过.. 结果貌似和我想象的不太一样.怎么说呢, 听我下文一 ...

  7. 【原理+源码详细解读】从Transformer到ViT

    文章目录 参考文献 简介 Transformer架构 Position Encoding Self-attention Multi-head Self-attention Masked Multi-H ...

  8. SpringBoot（2.4.0）自动配置原理(源码)

    一.从@SpringBootApplication讲起 源码 @Target({ElementType.TYPE}) @Retention(RetentionPolicy.RUNTIME) @Docu ...

  9. Mybatis Interceptor 拦截器原理 源码分析

    Mybatis采用责任链模式,通过动态代理组织多个拦截器(插件),通过这些拦截器可以改变Mybatis的默认行为(诸如SQL重写之类的),由于插件会深入到Mybatis的核心,因此在编写自己的插件前最 ...

最新文章

  1. HTML5 学习之地理定位
  2. 怎么用python做表格-零基础小白怎么用Python做表格？
  3. 文件创建与文件格式的修改
  4. MainWindow 简介
  5. Wpf拖动按钮实现（二）
  6. 【POJ - 3169】 Layout（差分约束+spfa）（当板子记？）
  7. python 绘制柱状图
  8. 分享]人生忠告——七天改变人生影响世界
  9. HP server ILO
  10. 一周总结汇总_2016-09-25
  11. 超神能力：云库局面分析
  12. oracle存货管理系统,存货管理系统的功能模块有哪些?
  13. 数字信号处理重要学习资源
  14. 如何在linux上运行asp网站,linux上搭建asp网站
  15. 如何测试计算机的运行速度,如何查看cpu运行速度
  16. Operator基础：2: Operator SDK安装
  17. 开源软件的许可证（License）
  18. CCNA之EIGRP(IGRP)
  19. 如何打jar包和运行jar包
  20. 自然语言处理及计算语言学相关术语中英对译表

热门文章

  1. SAP案例教程FI财务后台配置
  2. 腾讯人口密度热力图_从腾讯位置大数据，看中国的超级城市
  3. 微信群管理助手哪里弄的?
  4. 【办公协作软件】万彩办公大师教程丨图片OCR工具的应用
  5. Android面试题解答(结尾有彩蛋)
  6. vue 二级路由嵌套和二级路由高亮问题
  7. 如何开启红米手机4X的ROOT超级权限
  8. 硬盘的读写速度如何计算
  9. k8s 应用更新策略：灰度发布和蓝绿发布
  10. 交换机 tagged 与 untagged 的关系