简介
本算法可以用来计算任意给定的两个凸多边形的重叠面积。算法思路写起来比较麻烦,希望大家能自己跟着代码走,把思路捋清楚。初次实现这个算法我是用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();
}