Delaunay三角网

三角网生长法（课本上的方法按原意实现）

//控制台演示
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <ctime>
#define PI 3.14159

struct POINT
{	
	double x;
	double y;
};
typedef POINT VECTOR;

struct SEGMENT
{
	POINT * p_begin;
	POINT * p_end;
	bool bExpandable;//记录边可否扩展
};

struct TRIANGLE
{
	POINT * vertex1;
	POINT * vertex2;
	POINT * vertex3;
};

//方便调试的功能
std::ostream & operator<<(std::ostream & os, const POINT & pt)
{
	os<<"("<<pt.x<<", "<<pt.y<<")";
	return os;
}
std::ostream & operator<<(std::ostream & os, const SEGMENT & seg)
{
	os<<"--: "<<*seg.p_begin<<" -> "<<*seg.p_end;
	return os;
}
std::ostream & operator<<(std::ostream & os, const TRIANGLE & tri)
{
	os<<"△: "<<*tri.vertex1<<" , "<<*tri.vertex2<<" , "<<*tri.vertex3;
	return os;
}

int main()
{
	//数据来源于《空间分析》课本P205
	clock_t begin = clock();
	std::cout<<"数据来源于《空间分析》课本P205 采用内核生长算法"<<std::endl;
	const int size = 14;
	POINT pts[size] = {{5.5, 2}, {13.5, 4.5}, {22.5, 2.5}, {22.6, 7.5}, {15.3, 8.7}, {2.1, 12}, {30.8, 10.6}, {6, 14.5}, {10.5, 18.8}, {22.4, 19}, {25.8, 21.5}, {14, 22.5}, {3.4, 24}, {14.2, 27}};

	std::vector<SEGMENT> segment_list;//边表
	std::vector<TRIANGLE> triangle_list;//三角形表

	//最近两点
	auto GetDis = [] (const POINT & p1, const POINT & p2) ->double {return sqrt(pow(p1.x - p2.x, 2) + pow(p1.y - p2.y, 2));};
	double min_dis = GetDis(pts[0], pts[1]);
	int min_i = 0;
	int min_j = 1;
	for(int i = 0; i < size - 1; i++)
	{
		for(int j = i + 1; j < size; j++)
		{
			double dis = GetDis(pts[i], pts[j]);
			if(dis < min_dis)
			{
				min_dis = dis;
				min_i = i;
				min_j = j;
			}
		}
	}//至此，找到p1, p2
	SEGMENT seg_first = {pts + min_i, pts + min_j};
	std::cout<<"起始边："<<seg_first<<std::endl;

	//原理上的缺陷警告：若第一条边就在边界会出错

	//检测边是否为边界
	auto IsEdge = [] (const SEGMENT & seg/*线段结构*/, const POINT * pts/*点数组*/, int size/*数组大小*/) ->bool
	{
		bool pt_on_left = false;//是否有点落在seg_first的左边
		bool pt_on_right = false;//是否有点落在seg_first的右边
		for(int i = 0; i < size; i++)
		{
			VECTOR right_90 = {seg.p_end->y - seg.p_begin->y, -(seg.p_end->x - seg.p_begin->x)};
			VECTOR ap = {pts[i].x - seg.p_begin->x, pts[i].y - seg.p_begin->y};
			if(right_90.x * ap.x + right_90.y * ap.y < 0)
			{
				pt_on_left = true;
			}
			else if(right_90.x * ap.x + right_90.y * ap.y > 0)
			{
				pt_on_right = true;
			}
			if(pt_on_right && pt_on_left)//如果已经发现两边都有点了，那么一定不是边界线，结束循环
				break;
		}
		return !(pt_on_left && pt_on_right);
	};//{SEGMENT s = {pts, pts + 1};std::cout<<std::boolalpha<<IsEdge(s, pts, size)<<std::endl;}

	//检测第一条边是否为边界，若为边界，不能用该生长算法
	if(IsEdge(seg_first, pts, size))
	{
		std::cout<<"EXIT_FAILURE: seg_first is on the edge"<<std::endl;
		exit(EXIT_FAILURE);
	}

	auto IsLeft = [] (const POINT & pt_begin, const POINT & pt_end, const POINT & pt) ->bool
	{
		//VECTOR base = {pt_end.x - pt_begin.x, pt_end.y - pt_begin.y};
		//VECTOR right_90 = {base.y, -base.x};
		VECTOR right_90 = {pt_end.y - pt_begin.y, -(pt_end.x - pt_begin.x)};
		VECTOR begin_to_pt = {pt.x - pt_begin.x, pt.y - pt_begin.y};
		return right_90.x * begin_to_pt.x + right_90.y * begin_to_pt.y < 0;
	};

	//计算pt到向量首位连线的夹角余弦值
	auto CalcCos = [] (const POINT & pt_begin, const POINT & pt_end, const POINT & pt) ->double
	{
		if(pt.x == pt_begin.x && pt.y == pt_begin.y || pt.x == pt_end.x && pt.y == pt_end.y)
		{
			return 1.0;//如果比较点是起始点赋个最大值
		}
		VECTOR v1 = {pt_begin.x - pt.x, pt_begin.y - pt.y};
		VECTOR v2 = {pt_end.x - pt.x, pt_end.y - pt.y};
		double cos = (v1.x * v2.x + v1.y * v2.y) / (sqrt(pow(v1.x, 2) + pow(v1.y, 2)) * sqrt(pow(v2.x, 2) + pow(v2.y, 2)));
		return cos;
	};/*{POINT p0 = {0,0}, p1 = {1, 0}, p2 = {1, 0};std::cout<<CalcCos(p0, p1, p2)<<std::endl;}*/

	double left_min_cos = 1.0;
	double right_min_cos = 1.0;//最大角cos值，cos最小值
	int left_i;
	int right_i;
	for(int i = 0; i < size; i++)
	{
		double cos = CalcCos(*seg_first.p_begin, *seg_first.p_end, pts[i]);
		if(IsLeft(*seg_first.p_begin, *seg_first.p_end, pts[i]))
		{
			if(cos < left_min_cos)//找左边夹角最大
			{
				left_min_cos = cos;//更新左边cos最小值
				left_i = i;
			}
		}
		else
		{
			if(cos < right_min_cos)//找右边夹角最大
			{
				right_min_cos = cos;
				right_i = i;
			}
		}
	}//至此，找到p3, p4
	std::cout<<"左边P3："<<pts[left_i]<<std::endl;
	std::cout<<"右边P4："<<pts[right_i]<<std::endl;
	
	segment_list.resize(4);
	for(int i = 0; i < 4; i++)
	{
		segment_list[i].bExpandable = true;//内核边界标记为可扩展边
	}
	segment_list[0].p_begin = seg_first.p_begin;//三角网内核
	segment_list[0].p_end = pts + left_i;
	segment_list[1].p_begin = pts + left_i;
	segment_list[1].p_end = seg_first.p_end;
	segment_list[2].p_begin = seg_first.p_end; 
	segment_list[2].p_end = pts + right_i;
	segment_list[3].p_begin = pts + right_i;
	segment_list[3].p_end = seg_first.p_begin;

	TRIANGLE tri1 = {seg_first.p_begin, seg_first.p_end, pts + left_i};
	TRIANGLE tri2 = {seg_first.p_begin, seg_first.p_end, pts + right_i};
	triangle_list.push_back(tri1);
	triangle_list.push_back(tri2);

	std::for_each(segment_list.begin(), segment_list.end(), [] (const SEGMENT & s) {std::cout<<s<<std::endl;});
	std::for_each(triangle_list.begin(), triangle_list.end(), [] (const TRIANGLE & t) {std::cout<<t<<std::endl;});

	auto GetLeftExpandPointPtr/*得到点结构指针*/ = [&IsLeft, &CalcCos] (const SEGMENT s, POINT * pts, int size) ->POINT *
	{
		int index;
		double min_cos = 1;
		for(int i = 0; i < size; i++)
		{
			if(IsLeft(*s.p_begin, *s.p_end, pts[i]))
			{
				double cos = CalcCos(*s.p_begin, *s.p_end, pts[i]);
				if(cos < min_cos)
				{
					min_cos = cos;
					index = i;
				}
			}
		}
		return pts + index;
	};

	for(int i = 0; i < segment_list.size(); i++)
	{
		if(!segment_list[i].bExpandable)//跳过本身被标记为不可扩展的边
		{	
			continue;//跳过本身被标记为不可扩展的边
		}
		else if(IsEdge(segment_list[i], pts, size))//如果该边是边界
		{
			segment_list[i].bExpandable = false;//将这条边标记为不可扩展边并跳过
			continue;//并跳过
		}

		POINT * tpt = GetLeftExpandPointPtr(segment_list[i], pts, size);
		TRIANGLE t_triangle = {segment_list[i].p_begin, segment_list[i].p_end, tpt};
		triangle_list.push_back(t_triangle);//加入三角形表

		SEGMENT ts1 = {segment_list[i].p_begin, tpt, true};
		SEGMENT ts2 = {tpt, segment_list[i].p_end, true};
		//寻找之前是否已经包括这两条条边的反向边，有则标记为不可扩展边
		std::vector<SEGMENT>::iterator its1 = std::find_if(segment_list.begin(), segment_list.end(), [&ts1] (const SEGMENT & s) ->bool 
		{
			if(ts1.p_begin == s.p_end && ts1.p_end == s.p_begin)
			{
				return true;
			}
			else
				return false;
		});
		if(its1 != segment_list.end())//如果找得到反向边
		{
			its1->bExpandable = false;//把找到的那条标记为不可扩展
			ts1.bExpandable = false;//把自己也标记为不可扩展
		}
		std::vector<SEGMENT>::iterator its2 = std::find_if(segment_list.begin(), segment_list.end(), [&ts2] (const SEGMENT & s) ->bool
		{
			if(ts2.p_begin == s.p_end && ts2.p_end == s.p_begin)
			{
				return true;
			}
			else
				return false;
		});
		if(its2 != segment_list.end())
		{
			its2->bExpandable = false;
			ts2.bExpandable = false;
		}
		
		//如果可以扩展就添加到边表
		if(ts1.bExpandable)
		{
			segment_list.push_back(ts1);
		}
		if(ts2.bExpandable)
		{
			segment_list.push_back(ts2);
		}

		//扩展之后要把自己设置为不可扩展
		segment_list[i].bExpandable = false;
	}

	//显示结果
	std::cout<<std::endl;
	std::cout<<"显示结果："<<std::endl;
	int count = 0;
	std::for_each(triangle_list.begin(), triangle_list.end(), [&count] (const TRIANGLE & t) {std::cout<<t<<std::endl; count++;});
	std::cout<<"剖分三角形共"<<count<<"个"<<std::endl;

	std::cout<<"用时(包括打印时间)："<<clock() - begin<<"ms"<<std::endl;
	return 0;
}

为了更加直观，我用MFC做了下图形

此外还可以胡乱加些点