一个计算任意凸多边形重叠面积的算法

任意凸多边形重叠面积计算

简介

本算法可以用来计算任意给定的两个凸多边形的重叠面积。算法思路写起来比较麻烦,希望大家能自己跟着代码走,把思路捋清楚。初次实现这个算法我是用fortran语言编写的,但考虑到除了做科学计算的群体以外似乎用fortran的不多,所以这边在原代码的基础上将它改写成了C++,由于基本是从Fortran逐行改过来的,所以一些代码细节可能会看着别扭。另外,算法的思路是适用于任意两个凸多边形重叠面积计算的,不过下面代码中我是假设了多边形类型只包括三角形和四边形(当时我的任务中只涉及到这两个类型),所以很多变量没用动态类型,因为边数最多为4。
好了,不多说了,直接贴代码。。算法的重点在 “重叠形状的重建”那一块。

完整代码

#include <stdio.h>
#include <iostream>
#include <vector>
#include <math.h>

using namespace std;

typedef struct{
	int nVer;
	double coord[4][2];
} t_polygon;

bool PointInPolygon(double* coord, t_polygon* polygon) {
    
	bool res;
	double vec1[2], vec2[2], ref, tmp;
	int i, nVer, point2;

	nVer = polygon->nVer;
	vec1[0] = polygon->coord[1][0] - polygon->coord[0][0];
	vec1[1] = polygon->coord[1][1] - polygon->coord[0][1];
	vec2[0] = coord[0] - polygon->coord[0][0];
	vec2[1] = coord[1] - polygon->coord[0][1];

	ref = vec1[0] * vec2[1] - vec1[1] * vec2[0];

	res = true;

	for (i = 2; i <= nVer; i++){
		point2 = i + 1;
		if (i == nVer) point2 = 1;
		vec1[0] = polygon->coord[point2 - 1][0] - polygon->coord[i - 1][0];
		vec1[1] = polygon->coord[point2 - 1][1] - polygon->coord[i - 1][1];
		vec2[0] = coord[0] - polygon->coord[i-1][0];
		vec2[1] = coord[1] - polygon->coord[i-1][1];
		tmp = vec1[0] * vec2[1] - vec1[1] * vec2[0];
		if (tmp*ref < 0.0) res = false;
	}
	
	return res;
}

double CalArea(int nVer, double** coord){
	double area;
	int i, point1, point2;

	area = 0.0;
	for (i = 1; i <= nVer; i++){
		point1 = i;
		point2 = i + 1;
		if (i == nVer) point2 = 1;

		area += ((coord[point2-1][0] + coord[point1-1][0])*(coord[point2-1][1] - coord[point1-1][1]));
	}

	area = 0.5*fabs(area);
	return area;
}

double OverlappedArea(t_polygon* polygon1, t_polygon* polygon2){
	  
	//calculate the overlapped area of two arbitrary convex polygon.Here, only triangel and quadrangle are considered.
	double olp_area;
	t_polygon* polygon_tmp;
	double (*coord1)[2], (*coord2)[2], (*coord_tmp)[2];
	double** coord_new = nullptr;
	double coord_intpoin[8][2];      //coordinate of the intersection point.最多有8个交点。
	int connections[4][4] = { { -1, -1, -1, -1 }, { -1, -1, -1, -1 }, { -1, -1, -1, -1 }, { -1, -1, -1, -1 } },  //used to indicate which intersection point it is.
		intpoin_on_edge1[4][2] = { { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 } },      //used to indicate which edge intersect with the current edge.For polygon1.
		intpoin_on_edge2[4][2] = { { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 } };		 //used to indicate which edge intersect with the current edge.For polygon2.
	int (*intpoin_on_edge)[2];		 // a temp pointer
	int tag[16] = {0};				 // 4 vertices + 4 vertices + 8 intpoins
	vector<int> vec_vtx;			 // store vertices
	bool inPolygon;
	int i, j, iVer, i_next, i_tmp, iEdge,
		idxA, idxB, idxP, idxQ;		 // edge(A, B), and edge(P, Q)
	int nVer, nVer1, nVer2, counter;
	double PA[2], QA[2], norm[2], coord_m[2], temp, t, d, e, t_old;
	double eps = 1e-14;
	int error;

	coord1 = polygon1->coord;
	coord2 = polygon2->coord;
	nVer1 = polygon1->nVer;
	nVer2 = polygon2->nVer;
	
	counter = 0;

	//determine intpoin_on_edge1
	for (i = 1; i <= nVer1; i++){
		idxP = i;
		idxQ = i + 1;
		if (idxP == nVer1) idxQ = 1;
		for (j = 1; j <= nVer2; j++){
			idxA = j;
			idxB = j + 1;
			if (idxA == nVer2) idxB = 1;
			norm[0] = coord2[idxA - 1][1] - coord2[idxB - 1][1]; //AB = (x,y), then norm = (-y,x)
			norm[1] = coord2[idxB - 1][0] - coord2[idxA - 1][0];
			norm[0] = norm[0] / sqrt(norm[0] * norm[0] + norm[1] * norm[1]);
			norm[1] = norm[1] / sqrt(norm[0] * norm[0] + norm[1] * norm[1]);
			
			PA[0] = coord2[idxA-1][0] - coord1[idxP-1][0];
			PA[1] = coord2[idxA - 1][1] - coord1[idxP - 1][1];

			QA[0] = coord2[idxA-1][0] - coord1[idxQ-1][0];
			QA[1] = coord2[idxA - 1][1] - coord1[idxQ - 1][1];
			d = PA[0] * norm[0] + PA[1] * norm[1];
			e = QA[0] * norm[0] + QA[1] * norm[1];

			if (d*e >= 0.0) continue;
			
			t = d / (d - e);
			coord_m[0] = (1.0 - t)*coord1[idxP-1][0] + t*coord1[idxQ-1][0];
			coord_m[1] = (1.0 - t)*coord1[idxP - 1][1] + t*coord1[idxQ - 1][1];
			if ((coord_m[0] - coord2[idxA-1][0])*(coord_m[0] - coord2[idxB-1][0]) +
				(coord_m[1] - coord2[idxA-1][1])*(coord_m[1] - coord2[idxB-1][1]) <= 0.0) {

				counter += 1;
				coord_intpoin[counter - 1][0] = coord_m[0];
				coord_intpoin[counter - 1][1] = coord_m[1];
				connections[j-1][i-1] = counter;
			
				if (intpoin_on_edge1[i-1][0]==-1){
					intpoin_on_edge1[i-1][0] = j;
					t_old = t;
				}
				else{
					if (t < t_old){
						intpoin_on_edge1[i-1][1] = intpoin_on_edge1[i-1][0];
						intpoin_on_edge1[i-1][0] = j;
					}else
						intpoin_on_edge1[i-1][1] = j;
					break;
				}
			}
		}
	}

	/*for (i = 0; i < 4; i++){
		for (j = 0; j < 4; j++){
			printf("%d  ", connections[i][j]);
		}
		printf("\n");
	}*/

	//determine intpoin_on_edge2
	for (i = 1; i <= nVer2; i++){
		idxP = i;
		idxQ = i + 1;
		if (idxP == nVer2) idxQ = 1;

		for (j = 1; j <= nVer1; j++){
			idxA = j;
			idxB = j + 1;
			if (idxA == nVer1) idxB = 1;

			norm[0] = coord1[idxA-1][1] - coord1[idxB-1][1];
			norm[1] = coord1[idxB-1][0] - coord1[idxA-1][0];
			norm[0] = norm[0] / sqrt(norm[0] * norm[0] + norm[1] * norm[1]);
			norm[1] = norm[1] / sqrt(norm[0] * norm[0] + norm[1] * norm[1]);

			PA[0] = coord1[idxA - 1][0] - coord2[idxP - 1][0];
			PA[1] = coord1[idxA - 1][1] - coord2[idxP - 1][1];
			QA[0] = coord1[idxA - 1][0] - coord2[idxQ - 1][0];
			QA[1] = coord1[idxA - 1][1] - coord2[idxQ - 1][1];
			d = PA[0] * norm[0] + PA[1] * norm[1];
			e = QA[0] * norm[0] + QA[1] * norm[1];

			if (d*e >= 0.0) continue;

			t = d / (d - e);
			coord_m[0] = (1.0 - t)*coord2[idxP-1][0] + t*coord2[idxQ-1][0];
			coord_m[1] = (1.0 - t)*coord2[idxP - 1][1] + t*coord2[idxQ - 1][1];
			if ((coord_m[0] - coord1[idxA - 1][0])*(coord_m[0] - coord1[idxB - 1][0]) +
				(coord_m[1] - coord1[idxA - 1][1])*(coord_m[1] - coord1[idxB - 1][1]) <= 0.0){

				if (connections[i-1][j-1] == -1){
					counter += 1;
					coord_intpoin[counter-1][0] = coord_m[0];
					coord_intpoin[counter - 1][1] = coord_m[1];
					connections[i-1][j-1] = counter;
				}
				if (intpoin_on_edge2[i-1][0]==-1){
					intpoin_on_edge2[i-1][0] = j;
					t_old = t;
				}
				else{
					if (t < t_old){
						intpoin_on_edge2[i - 1][1] = intpoin_on_edge2[i - 1][0];
						intpoin_on_edge2[i - 1][0] = j;
					}else
						intpoin_on_edge2[i-1][1] = j;
					break;
				}
			}
		}
	}
	
	// if there is no intersection point
	if (counter == 0){
		coord_new = new double*[4];
		for (i = 0; i < 4; i++){
			coord_new[i] = new double[2];
		}
		for (i = 0; i < nVer1; i++){
			coord_new[i][0] = coord1[i][0];
			coord_new[i][1] = coord1[i][1];
		}
		temp = CalArea(nVer1,coord_new);
		for (i = 0; i < nVer2; i++){
			coord_new[i][0] = coord2[i][0];
			coord_new[i][1] = coord2[i][1];
		}
		olp_area = CalArea(nVer2, coord_new);
		if (temp < olp_area) olp_area = temp;

		if (coord_new != nullptr){
			for (i = 0; i < 4; i++){
				delete[] coord_new[i];
			}
			delete[] coord_new;
		}

		return olp_area;
	}

	//constructe the overlapped shape
	i = 1;
	iEdge = 1;
	inPolygon = PointInPolygon(coord1[0], polygon2);
	do {
		if (i <= 4){
			iVer = i;
			//iEdge = i;
			if (i != 4)
				i_next = i + 1;
			else
				i_next = 1;
			coord_tmp = coord1;
			polygon_tmp = polygon2;
			intpoin_on_edge = intpoin_on_edge1;
		}
		else if (i <= 8){
			iVer = i - 4;
			//iEdge = i;
			if (i != 8)
				i_next = i + 1;
			else
				i_next = 5;
			coord_tmp = coord2;
			polygon_tmp = polygon1;
			intpoin_on_edge = intpoin_on_edge2;
		}
		else {//i>8
			switch(iEdge){
			case 1:
			case 2:
			case 3:
				iVer = iEdge + 1;
				coord_tmp = coord1;
				polygon_tmp = polygon2;
				break;
			case 4:
				iVer = 1;
				coord_tmp = coord1;
				polygon_tmp = polygon2;
				break;
			case 5:
			case 6:
			case 7:
				iVer = iEdge + 1 - 4;
				coord_tmp = coord2;
				polygon_tmp = polygon1;
				break;
			case 8:
				iVer = 1;
				coord_tmp = coord2;
				polygon_tmp = polygon1;
				break;
			}
		}

		if (i <= 8){
			if (inPolygon){
				vec_vtx.push_back(i);
				tag[i-1] = 1;
				if (intpoin_on_edge[iVer-1][0]==-1){
					i = i_next;
					iEdge = i_next;
					inPolygon = true;
				}
				else{
					if (i <= 4){
						i = 8 + connections[intpoin_on_edge[iVer - 1][0]-1][iVer - 1];
						iEdge = 4 + intpoin_on_edge[iVer - 1][0];
					}
					else{
						i = 8 + connections[iVer-1][intpoin_on_edge[iVer-1][0]-1];
						iEdge = intpoin_on_edge[iVer-1][0];
					}
				}
			}
			else{
				if (intpoin_on_edge[iVer - 1][0] == -1){
					i = i_next;
					iEdge = i_next;
					inPolygon = false;
				}
				else{
					if (i <= 4)
						//iEdge = i,
						i = 8 + connections[intpoin_on_edge[iVer - 1][0]-1][iVer - 1];
					else
						//iEdge = i,
						i = 8 + connections[iVer - 1][intpoin_on_edge[iVer - 1][0]-1];
				}
			}
		}
		else{ //i>8
			vec_vtx.push_back(i);
			tag[i - 1] = 1;
			if (PointInPolygon(coord_tmp[iVer-1], polygon_tmp)){
				if (iEdge <= 4)
					i = iVer;
				else
					i = iVer + 4;
				inPolygon = true;
			}
			else{
				if (iEdge <= 4){
					i = 8 + connections[intpoin_on_edge1[iEdge - 1][1]-1][iEdge - 1];
					iEdge = 4 + intpoin_on_edge1[iEdge - 1][1];
				}
				else{
					i = 8 + connections[iEdge-4-1][intpoin_on_edge2[iEdge - 4 - 1][1]-1];
					iEdge = intpoin_on_edge2[iEdge-4-1][1];
				}
			}
		}
	} while (tag[i - 1] == 0);

	//calculate the polygon area according to its vertices 
	nVer = vec_vtx.size();
	coord_new = new double* [nVer];
	for (i = 0; i < nVer; i++){
		coord_new[i] = new double [2];
	}

	for (i = 1; i <= nVer; i++){
		j = vec_vtx[i-1];
		if (j <= 4){
			coord_new[i - 1][0] = coord1[j - 1][0];
			coord_new[i - 1][1] = coord1[j - 1][1];
		}
		else if(j<=8){
			coord_new[i - 1][0] = coord2[j - 4 - 1][0];
			coord_new[i - 1][1] = coord2[j - 4 - 1][1];
		}
		else{
			coord_new[i - 1][0] = coord_intpoin[j - 8 - 1][0];
			coord_new[i - 1][1] = coord_intpoin[j - 8 - 1][1];
		}
	}

	olp_area = CalArea(nVer, coord_new);

	if (coord_new != nullptr){
		for (i = 0; i < nVer; i++){
			delete[] coord_new[i];
		}
		delete[] coord_new;
	}

	return olp_area;
}
	

int PolygonRcs(t_polygon* polygon, int nVer, double coord[4][2]){
	int i, j;
	polygon->nVer = nVer;
	for (i = 0; i < nVer; i++){
		for (j = 0; j < 2;j++)
			polygon->coord[i][j] = coord[i][j];
	}
	return 0;
}


void main(){

	t_polygon *polygon1, *polygon2;
	double olp_area;
	//double coord1[4][2] = { { 0.0, 0.0 }, { 0.5, 0.8660254 }, { 1.0, 0.0 }, { 0.0, 0.0 } };
	//double coord2[4][2] = { { 0.5, 0.5 }, { 1.5, 0.5 }, { 1.5, -0.5 }, { 0.5, -0.5 } };
	double coord1[4][2] = { { 0.0, 0.0 }, { 0.5, 0.8660254 }, { 1.0, 0.0 }, { 0.0, 0.0 } };
	double coord2[4][2] = { { 0.1, 0.5 }, { 0.9, 0.5 }, { 0.9, -0.5 }, { 0.1, -0.5 } };

	polygon1 = new t_polygon;
	polygon2 = new t_polygon;

	PolygonRcs(polygon1, 3, coord1);
	PolygonRcs(polygon2, 4, coord2);

	olp_area = OverlappedArea(polygon1, polygon2);

	printf("overlapped area = %f\n", olp_area);
	
	getchar();
}
  • 0
    点赞
  • 8
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值