两个简单多边形求交(CGAL)

两个多边形求交的实现需要几个模块
(cgal中有insect函数,但是必须要求使用CGAL::Exact_predicates_exact_constructions_kernel)
本人出于其他想法,没有把基于Exact_predicates_inexact_constructions_kernel核的Polygon_2转换为
Exact_predicates_exact_constructions_kernel的核,而是自己写了一个求交函数自用.

1.线线求交

1.1 cgal版

//cgal实现版本
bool  GetInsertionOfLineAndLine(double bsx,double bsy,double bex,double bey, double sx, double sy, double ex, double ey, double&rx, double &ry) {

		CGAL::Exact_predicates_exact_constructions_kernel::Segment_2 lin(CGAL::Exact_predicates_exact_constructions_kernel::Point_2(bsx, bsy), CGAL::Exact_predicates_exact_constructions_kernel::Point_2(bex, bey));
		CGAL::Exact_predicates_exact_constructions_kernel::Segment_2 seg(CGAL::Exact_predicates_exact_constructions_kernel::Point_2(sx, sy), CGAL::Exact_predicates_exact_constructions_kernel::Point_2(ex, ey));


		CGAL::cpp11::result_of<CGAL::Exact_predicates_exact_constructions_kernel::Intersect_2(CGAL::Exact_predicates_exact_constructions_kernel::Segment_2, CGAL::Exact_predicates_exact_constructions_kernel::Segment_2)>::type
		result = intersection(lin, seg);

		if (result) {
			if (const CGAL::Exact_predicates_exact_constructions_kernel::Segment_2* s = boost::get<CGAL::Exact_predicates_exact_constructions_kernel::Segment_2>(&*result)) {

			}
			else {
				const CGAL::Exact_predicates_exact_constructions_kernel::Point_2* p = boost::get<CGAL::Exact_predicates_exact_constructions_kernel::Point_2 >(&*result);
				rx = to_double((*p).x());
				ry = to_double((*p).y());
				//std::cout << jx << ";" << jy << endl;
				return true;
			}
		}
	return false;
}

2.判断点在多边形内

2.1通用版(判断内角和是否等于360°)

double  GetAngle(Point_2 p0, Point_2 p1, Point_2 p2)
{
	double Pi = 3.1415926;
	double EPS = 0.000001;

	Point_2 a(p0.x() - p1.x(), p0.y() - p1.y());
	Point_2 b(p2.x() - p1.x(), p2.y() - p1.y());

	double angle1 = atan2(a.y(), a.x()) * 180 / Pi;
	double angle2 = atan2(b.y(), b.x()) * 180 / Pi;
	if (angle1*angle2 > -EPS)
		return abs(angle1 - angle2)*Pi / 180;
	else {
		double tmp = abs(angle1) + abs(angle2);
		if (tmp - 180 > EPS)
			tmp = 360 - tmp;
		return tmp * Pi / 180;
	}
}


bool  IsPointInPolygon(const Point_2 p, const Polygon_2 poly)
{
	double EPS = 0.000001;
	double Pi = 3.1415926;

	int length = poly.size();
	double angle = 0;
	for (int i = 0; i < length; ++i)
	{
		angle += GetAngle(poly.vertex(i), p, poly.vertex((i + 1) % length));
	}
	return abs(angle - 2 * Pi) < EPS;
}

bool IsPolygonInPolygon(const Polygon_2 p, const Polygon_2 poly)
{
	int len = p.size();
	for (int i = 0; i < len; ++i)
	{
		if (!IsPointInPolygon(p.vertex(i), poly))
			return false;
	}
	return true;
}

3.多边形求交

3.1 通用版(Polygon逆时针存储vertex)

a.通过线线求交找到两个交点 p1,p2 以及对应的相交边的起始顶点 i1,j1 和 i2,j2
b.通过判断 i1+1是否在polyb中,获取相交部分的逆时针索引vector

bool  Intersect(const Polygon_2 polya, const Polygon_2 polyb, Polygon_2 *intersect)
{
	int lena = polya.size();
	int lenb = polyb.size();

	double *_polyax = new double[lena];
	double *_polyay = new double[lena];

	double *_polybx = new double[lenb];
	double *_polyby = new double[lenb];

	for (int i = 0; i < lena; ++i)
	{
		_polyax[i] = polya.vertex(i).x();
		_polyay[i] = polya.vertex(i).y();
	}
	for (int i = 0; i < lenb; ++i)
	{
		_polybx[i] = polyb.vertex(i).x();
		_polyby[i] = polyb.vertex(i).y();
	}

	bool first = false;
	bool second = false;
	int i1, i2, j1, j2;
	double p1x, p1y, p2x, p2y;

	bool IsInsertion = false;
	
	for (int i = 0; i < lena; ++i)
	{
		for (int j = 0; j < lenb; ++j)
		{
			IsInsertion = false;
			double jx, jy;
			IsInsertion = GetInsertionOfLineAndLine(_polyax[i], _polyay[i], _polyax[(i+1)%lena], _polyay[(i + 1) % lena], _polybx[j], _polyby[j], _polybx[(j + 1) % lenb], _polyby[(j + 1) % lenb], jx, jy);
			if (IsInsertion) {
				if (!first) {
					first = true;
					i1 = i;	j1 = j;
					p1x = jx; p1y = jy;
				}
				else {
					second = true;
					i2 = i;	j2 = j;
					p2x = jx; p2y = jy;
				}
				break;
			}
		}
	}

	if (!first || !second)
		return false;

	Polygon_2 result; int s1, s2; 
	//if vertex(i1+1) is located in poly2,get p1->i1+1->...->i2->p2->j2+1->...->j1
	if (this->IsPointInPolygonByBao(polya.vertex((i1+1)%lena), polyb))	
	{
		s1 = (i1 + 1) % lena;
		s2 = (j2 + 1) % lenb;

		result.push_back(Point_2(p1x, p1y));
		if (s1 <= i2)
		{
			for (int i = s1; i <= i2; ++i)
				result.push_back(polya.vertex(i));
		}
		else {
			for (int i = s1; i < lena; ++i)
				result.push_back(polya.vertex(i));
			for(int i=0;i<=i2;++i)
				result.push_back(polya.vertex(i));
		}
		result.push_back(Point_2(p2x, p2y));
		if(j1 != j2)
			if (s2 <= j1)
			{
				for (int i = s2; i <= j1; ++i)
					result.push_back(polyb.vertex(i));
			}
			else {
				for (int i = s2; i < lenb; ++i)
					result.push_back(polyb.vertex(i));
				for (int i = 0; i <= j1; ++i)
					result.push_back(polyb.vertex(i));
			}

	}
	//if vertex(i1+1) isn't located in poly2,get p1->j1+1->...->j2->p2->i2+1->...->i1
	else {
		s1 = (j1 + 1) % lenb;
		s2 = (i2 + 1) % lena;

		result.push_back(Point_2(p1x, p1y));

		if(j1!=j2)
			if (s1 <= j2)
			{
				for (int i = s1; i <= j2; ++i)
					result.push_back(polyb.vertex(i));
			}
			else {
				for (int i = s1; i < lenb; ++i)
					result.push_back(polyb.vertex(i));
				for (int i = 0; i <= j2; ++i)
					result.push_back(polyb.vertex(i));
			}
		result.push_back(Point_2(p2x, p2y));

	
		if (s2 <= i1)
		{
			for (int i = s2; i <= i1; ++i)
				result.push_back(polya.vertex(i));
		}
		else {
			for (int i = s2; i < lena; ++i)
				result.push_back(polya.vertex(i));
			for (int i = 0; i <= i1; ++i)
				result.push_back(polya.vertex(i));
		}
	}


	*intersect = result;
	return true;
}