2. 程式人生 > >Delaunay 三角剖分2D（原理 + 原始碼）

Delaunay 三角剖分2D（原理 + 原始碼）

阿新 • • 發佈：2018-12-16

Delaunay 三角剖分

程式碼實現< C++版本 >

程式碼框架

  • 定義點類 vector2.h
  • 定義邊類 edge.h
  • 定義三角形類 triangle.h
  • 定義三角剖分類 delaunay.h
  • 定義比較類(用於三角退化) numeric.h
  • 定義main函式 main.cpp

用到的外部庫

具體實現

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);
            }), 
            
  

           

              
              

            

            
相關推薦

			   
            
            
            
 

    


    
    
Delaunay 三角剖分2D（原理 + 原始碼）

     
 
							
							
							
Delaunay 三角剖分


程式碼實現< C++版本 >

程式碼框架

定義點類  vector2.h
定義邊類    edge.h
定義三角形類    triangle.h
定義三角剖分類    delaunay.h
定義比較類(用於三角退 



  
 

    


    
    
cv2 實現Delaunay三角剖分和Voronoi圖（Delaunay Triangulation and Voronoi Diagram）

     
  
 
 
 本文主要介紹使用cv2模組實現Delaunay三角剖分和Voronoi圖。 
 測試圖片： 
  
 示例程式碼： 
 # 匯入模組
import cv2
import numpy as np
import random

# 檢查點是否在長方形區域內
def rect_contains(rec 



  
 

    


    
    
VC+OpenGL實現空間三維Delaunay三角剖分

     
  
 
  三維建模和等值面的繪製過程中，需要經常使用三角形網格對資料體進行構面。而三角形的生成基於Delaunay三角剖分的演算法實現的。前段時間一直在考慮資料體的任意剖面切割該怎麼做，但是一直被兩個問題所困擾，一個就是交點問題，然後就是對所求交點進行繪製問題（三角形網格面構造）。終於在半個月後有 



  
 

    


    
    
OpenCascade中的Delaunay三角剖分

     
 Delaunay Triangulation in OpenCascade
摘要：本文簡要介紹了Delaunay三角剖分的基礎理論，並使用OpenCascade的三角剖分演算法將邊界BRep表示的幾何體進行三角離散化後在OpenSceneGraph中顯示。 
關鍵字：Delaunay Triangulat 



  
 

    


    
    
OpenCV的Delaunay三角剖分和Voronoi圖的實現

     
 給定平面中的一組點，三角測量指的是將平面細分為三角形，將點作為頂點。在圖1中，我們在左影象上看到一組界標，以及在中間影象中的三角測量。一組點可以有許多可能的三角剖分，但Delaunay三角剖分出眾，因為它有一些不錯的屬性。在Delaunay三角剖分中，選擇三角形使得沒有點在任何三角形的外接圓內。圖2示出了4點 



  
 

    


    
    
Captain Dialog 2009-09-18 VC+OpenGL實現空間三維Delaunay三角剖分

     
 
                
Captain Dialog 2009-09-18
    三維建模和等值面的繪製過程中，需要經常使用三角形網格對資料體進行構面。而三角形的生成基於Delaunay三角剖分的演算法實現的。前段時間一直在考慮資料體的任意剖面切割該怎麼做，但是一直被兩個問題所困擾，一個就是交點 



  
 

    


    
    
mysql 使用 limit 實現底層分頁（附原始碼）

     
  
 
 原理解析： 
 <select id="queryProductList" resultType="com.pojo.Product">
	SELECT *
	FROM
	tb_product
	ORDER BY
	priority DESC
	LIMIT #{rowIndex},#{p 



  
 

    


    
    
自然語言處理之：c++中文分詞（附原始碼）

     
 
								
								            
							
							
							githup地址：https://github.com/jbymy 
一、簡介 
        中文分詞是地然語言處理中的最基礎的環節，到目前為止已經有不少優秀的分詞工具的出現，如“中科院分詞”，“結 



  
 

    


    
    
凸多邊形最優三角剖分（演算法設計：動態規劃）

     
  
 
 一、動態規劃 
       和分治法類似，把原問題劃分成若干個子問題，不同的是，分治法（子問題間互相獨立），動態規劃（子問題不獨立） 
       動態規劃： 
 （1）找出最優解的性質，刻畫其結構特徵 
 （2）遞迴地定義最優值 
  



  
 

    


    
    
opencv平面三角剖分Delaunay

     
  
 
 opencv使用： 
  
  Subdiv2D  實現三角剖分。 
  使用參考：https://blog.csdn.net/czl389/article/details/62264960?fps=1&locationNum=5  
  
 關於三角剖分的基本知識參考： 



  
 

    


    
    
2018.09.30【POJ3348】Cows（凸包）（三角剖分）

     
 
							
							
							傳送門

解析：
讀優沒有寫負數又被卡了半個小時。。。
這裡採用JarrisJarrisJarris步進法求凸包。。主要講一講怎麼求多邊形面積。
思路：
滿足題意的顯然是這些點的凸包，而我們要做的就是求出凸包面積。
那麼怎麼求多邊形面積？
考慮三角剖分，我們將多 



  
 

    


    
    
BZOJ3589 動態樹（樹鏈剖分+容斥原理）

     
 std   class   down   ring   print   color   動態   inf   while   　　顯然容斥後轉化為求樹鏈的交。這個題非常良心的保證了查詢的路徑都是到祖先的，求交就很休閑了。

#include<iostream> 
#include<cstdi 



  
 

    


    
    
8.osg中使用Tesselator分格化（三角剖分）

     
 
                
在球的基礎上進行操作。
效果為與主視窗和子視窗新增一個邊框
        osg::ref_ptr<osg::Vec3Array>vertice2 = new osg::Vec3Array(8);
        //外邊界逆針
        (*verti 



  
 

    


    
    
凸多邊形最優三角剖分（動態規劃）

     
 #include<stdio.h>
#include<stdlib.h>
#include<iostream>
using namespace std;
const int N=7;
int Weight(int **w,int a,int b,int c)
{
    r 



  
 

    


    
    
BZOJ2243 [SDOI2011]染色（樹鏈剖分+線段樹合並）

     
 ech   sca   註意   get   printf   truct   bre   sum   lca   題目鏈接 BZOJ2243
樹鏈剖分+線段樹合並
線段樹合並的一些細節需要註意一下

#include <bits/stdc++.h>

using namespace std;

 



  
 

    


    
    
最優三角剖分

     
 return   min   std   esp   div   技術分享   include   %d   logs   同樣是紫書上的題。
紫書上並沒有給出每一個三角形所貢獻的的權值的計算方法，我這裏就擅作主張，定義成點權的乘積和好了。
那麽做法是DP，這裏註意設狀態的方式（我這麽設是為了使需要求解的問題 



  
 

    


    
    
樹鏈剖分[模板]（洛谷 P3384）

     
 www.   lca   葉子   class   logs   如果   葉子節點   ref   它的   洛谷·[模板]樹鏈剖分
寫在前面
首先，在學樹鏈剖分之前最好先把 LCA、樹形DP、DFS序 這三個知識點學了
如果這三個知識點沒掌握好的話，樹鏈剖分難以理解也是當然的。

樹鏈剖分
樹鏈剖分 就是 



  
 

    


    
    
BZOJ4448 SCOI2015情報傳遞（離線+樹鏈剖分+樹狀數組）

     
 fin   tree   amp   tdi   開始   情報   路徑   stream   from   　　即滋磁單點修改，詢問路徑上小於某數的值有多少個。暴力樹剖套個主席樹即可，然而不太優美。
　　開始覺得可以cdq，然而就變成log^3了。冷靜一下感覺簡直是個弱智，修改本身就是單調的，只要對詢問離 



  
 

    


    
    
bzoj4449: [Neerc2015]Distance on 點分治 三角剖分 平面圖與對偶圖

     
  
 
  
  
 bzoj4449: [Neerc2015]Distance on Triangulation 
 Description 
 給定一個凸n邊形，以及它的三角剖分。再給定q個詢問，每個詢問是一對凸多邊行上的頂點(a,b)，問點a最少經過多少條邊(可以是多邊形上的邊，也可以是剖分上的邊)可以 



  
 

    


    
    
廣義圓方樹+樹鏈剖分+set（Codeforces Round #278 (Div. 1): E. Tourists）

     
  
 
 
   
 前置：雙聯通分量、圓方樹、樹鏈剖分 
 什是是廣義圓方樹 
 圓方樹是針對於仙人掌建樹，而廣義圓方樹是針對無向圖建樹，對於一個無向圖 
 
  無向圖中的所有點 → 廣義圓方樹中的所有圓點  
  無向圖中的一個雙聯通分量 → 廣義圓方樹中的其中一個方點，這個方點向當