1. IT分享知识网首页

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

老牧童 • • 未分类 • 阅读 122

Delaunay 三角剖分2D（原理 + 源码）Delaunay三角剖分什么是三角剖分(推荐看这个作者的系列文章，写的相当好)什么是Delaunay三角剖分空圆特性最大化最小角特性Delaunay三角剖分的特点最接近唯一性最优性最规则区域性具有凸多边形的外壳有什么方法Bowyer-Watson算法—逐点插入法Lawson算法Lawson算法相关博客[图形学]Delaunay三…

大家好，欢迎来到IT知识分享网。

文章目录

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>::type
    almost_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;
                        }
                    }
                }

				//更新 polygon
                polygon.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>::type
    almost_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 window
    sf::RenderWindow window(sf::VideoMode(800, 600), "Delaunay triangulation");

    // Transform each points of each vector as a rectangle
    std::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 lines
    std::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 squares
        for(const auto &s : squares) { 
   
            window.draw(*s);
        }

        // Draw the lines
        for(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')

效果图：

在这里插入图片描述

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

演算法筆記

免责声明：本站所有文章内容,图片，视频等均是来源于用户投稿和互联网及文摘转载整编而成，不代表本站观点，不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益，请在线联系站长,一经查实,本站将立刻删除。 本文来自网络,若有侵权，请联系删除，如若转载，请注明出处：https://yundeesoft.com/12068.html

赞 (0)
老牧童老牧童
0

相关推荐

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

关注微信