计算两个多边形的重叠面积

#include <iostream>
#include <cmath>
#include <cstring>

using namespace std;

const int maxn = 300;
const double eps = 1e-6;
//位置标识
int dcmp(double x)
{
    if(x > eps) return 1;
    return x < -eps ? -1 : 0;
}

struct Point
{
    double x, y;
};

double cross(Point a, Point b, Point c)  
{
    return (a.x-c.x)*(b.y-c.y)-(b.x-c.x)*(a.y-c.y);  //叉积公式
}

//计算线段ab和cd的交点坐标
Point intersection(Point a, Point b, Point c, Point d)
{
    Point p = a;
    double t =((a.x-c.x)*(c.y-d.y)-(a.y-c.y)*(c.x-d.x))/((a.x-b.x)*(c.y-d.y)-(a.y-b.y)*(c.x-d.x)); 
    p.x +=(b.x-a.x)*t; 
    p.y +=(b.y-a.y)*t;
    cout << "intersection p.x=" << p.x << ", p.y=" << p.y << endl;
    return p;
}

//计算多边形面积，将多边形拆解成连续三个顶点组合成的多个三角形进行计算，这个循环计算一次其实是计算两次多边形的面积。
double PolygonArea(Point p[], int n)
{
    if(n < 3) return 0.0;
    double s = p[0].y * (p[n - 1].x - p[1].x);
    for(int i = 1; i < n - 1; ++ i) {
        s += p[i].y * (p[i - 1].x - p[i + 1].x);
        // cout << "p[i-1].x =" << p[i-1].x << ", p[i-1].y=" << p[i-1].y << endl;
        // cout << "p[i].x =" << p[i].x << ", p[i].y=" << p[i].y << endl;
        // cout << "p[i+1].x =" << p[i+1].x << ", p[i+1].y=" << p[i+1].y << endl;
    }
    s += p[n - 1].y * (p[n - 2].x - p[0].x);
    cout << "s =" << s << endl;
    return fabs(s * 0.5);
}

double CPIA(Point a[], Point b[], int na, int nb)  //ConvexPolygonIntersectArea
{
    Point p[20], tmp[20];
    int tn, sflag, eflag;
    memcpy(p,b,sizeof(Point)*(nb));
    for(int i = 0; i < na && nb > 2; i++)
    {
    	if (i == na - 1) {
    		sflag = dcmp(cross(a[0], p[0],a[i]));
    	} else {
    		sflag = dcmp(cross(a[i + 1], p[0],a[i]));
    	}
        for(int j = tn = 0; j < nb; j++, sflag = eflag)
        {
            if(sflag>=0) {
            	tmp[tn++] = p[j];
            }
            if (i == na - 1) {
            	if (j == nb -1) {
            		eflag = dcmp(cross(a[0], p[0], a[i]));
				} else {
					eflag = dcmp(cross(a[0], p[j + 1], a[i])); //计算下一个连续点在矢量线段的位置
				}
			} else {
				if (j == nb -1) {
					eflag = dcmp(cross(a[i + 1], p[0], a[i]));
				} else {
					eflag = dcmp(cross(a[i + 1], p[j + 1], a[i]));
				}
			}
            if((sflag ^ eflag) == -2) {  //1和-1的异或为-2，也就是两个点分别在矢量线段的两侧
            	if (i == na - 1) {
            		if (j == nb -1) {
            			tmp[tn++] = intersection(a[i], a[0], p[j], p[0]); //求交点
            		} else {
            			tmp[tn++] = intersection(a[i], a[0], p[j], p[j + 1]);
            		}
				} else {
					if (j == nb -1) {
						tmp[tn++] = intersection(a[i], a[i + 1], p[j], p[0]);
					} else {
						tmp[tn++] = intersection(a[i], a[i + 1], p[j], p[j + 1]);
					}
				}
            }
        }
        memcpy(p, tmp, sizeof(Point) * tn);
        nb = tn, p[nb] = p[0];
    }
    if(nb < 3) return 0.0;
    return PolygonArea(p, nb);
}

double SPIA(Point a[], Point b[], int na, int nb)    //SimplePolygonIntersectArea 调用此函数
{
    int i, j;
    Point t1[na], t2[nb];
    double res = 0, num1, num2;
    t1[0] = a[0], t2[0] = b[0];
    for(i = 2; i < na; i++)
    {
        t1[1] = a[i-1], t1[2] = a[i];
        num1 = dcmp(cross(t1[1], t1[2],t1[0]));  //根据差积公式来计算t1[2]在矢量线段（t1[0], t1[1]）的左侧还是右侧，
                                                 //值为负数在矢量线段左侧，值为正数在矢量线段右侧
        if(num1 < 0) swap(t1[1], t1[2]);  // 按逆时针进行排序
        for(j = 2; j < nb; j++)
        {
            t2[1] = b[j - 1], t2[2] = b[j];
            num2 = dcmp(cross(t2[1], t2[2],t2[0]));
            if(num2 < 0) swap(t2[1], t2[2]);  
            res += CPIA(t1, t2, 3, 3) * num1 * num2; 
        }
    }
    cout << "Sum::res=" <<res << endl;
    return res;
}

Point p1[maxn], p2[maxn];
int n1, n2;
int main()
{

    while(cin>>n1>>n2)
    {
        for(int i = 0; i < n1; i++) scanf("%lf%lf", &p1[i].x, &p1[i].y);
        for(int i = 0; i < n2; i++) scanf("%lf%lf", &p2[i].x, &p2[i].y);
        double Area = SPIA(p1, p2, n1, n2);
        cout << "Area=" << Area << endl;
        double A1 = PolygonArea(p1, n1);
        double A2 = PolygonArea(p2, n2);
        cout << "A1 =" << A1 << ", A2=" << A2 << endl;
    }

    return 0;
}

其中一些接口方法说明：
calculateOverlappedArea(ipc_point_t *a, ipc_point_t *b, int na, int nb)方法用于计算重叠面积。
参数含义：
ipc_point_t *a：检测区域的点的集合
ipc_point_t *b：目标物点的集合
int na：检测区域的点的数量
int nb：目标物点的数量
这个接口方法的设计思想：
将检测区域和目标的边框分别拆解成多个三角形，分别进行计算两个区域内不同三角形的重叠面积的总和，最后算出的结果也就是两个区域的重叠面积。

PolygonIntersectArea(ipc_point_t *a, ipc_point_t *b, int na, int nb)：计算两个三角形重叠部分的面积。
ipc_point_t *a：取自检测区域三个点组成的一个三角形点的数组
ipc_point_t *b：取自检测目标三个点组成的一个三角形点的数组
int na：数组a的数量，也就是三角形点的数量
int nb：数组b的数量，也就是三角形点的数量
这个接口方法的设计思想：
1） 方法形参中的数组a和b实则上是可以理解为组成的三角形A和B，并设置一个临时数组p用来存储两个三角形重叠部分的顶点坐标，数组p的初始值为数组b的拷贝。
2） 从三角形 A中取出一条边（也就是数组a里面的两个点），将A中的一条边与数组p的其中一条边进行计算，计算两条边是否有交点，如果有交点就将该交点拷贝放到数组p中，如果没有交点并且p取中的点在A中取出两个点组成的矢量线段的右侧，也将该点拷贝到数组p中，这个时候的数组p可能包含多个交点，组成的是一个多边形，这样第一轮下来就是三角形A中的其中一条边与数组p中的每一条边（即三角形B中的三条边）都会进行一次计算是否有交点。
3） 循环依次从三角形A中其他边，与新的数组p依次根据第 2）步骤进行计算。
4） 三角形A的每个边都与每次得到的新数组p比较结束之后，如果新数组p内的点的个数大于等于3，计算这个新数组p内所有点组成的多边形的面积。

intersection(ipc_point_t a, ipc_point_t b, ipc_point_t c, ipc_point_t d)：计算两个线段的交点
两个线段的顶点分别为(a, b), (c, d)

PolygonArea(ipc_point_t *p, int n)：计算多边形的面积，将多边形拆解成连续三个顶点组合成的多个三角形进行计算，这个循环计算一次其实是计算两次多边形的面积。
ipc_point_t *p：多边形的顶点集合
int n：多边形的顶点数