Delaunay 三角剖分2D（原理 + 源码）

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')

效果图：

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

演算法筆記