计算几何札记-线段交

叉积求面积和直线交点坐标模板：

int n;
struct Point
{
    double x;
    double y;
};
Point a[MAX], b[MAX], c[MAX], d[MAX], p[MAX][MAX];
double area[MAX][MAX], ans;
Point GetPoint(Point p1, Point p2, Point p3, Point p4)//已知四个点求直线相交的点的坐标
//（这种求法包括了直线斜率不存在或为0的情况下）
{
    Point p0;
    double A1 = p2.y - p1.y;
    double B1 = p1.x - p2.x;
    double C1 = p1.y * (-B1) - p1.x * A1;
    double A2 = p4.y - p3.y;
    double B2 = p3.x - p4.x;
    double C2 = p3.y *(-B2) - p3.x * A2;
    p0.x = (C2 * B1 - C1 * B2) / (A1 * B2 - A2 * B1);
    p0.y = (C2 * A1 - C1 * A2) / (B1 * A2 - B2 * A1);
    return p0;
}
void GetPoints()
{
    int i, j;
    //矩形四个角的坐标
    p[1][1].x = 0.0;
    p[1][1].y = 0.0;
    p[1][n + 2].x = 0.0;
    p[1][n + 2].y = 1.0;
    p[n + 2][1].x = 1.0;
    p[n + 2][1].y = 0.0;
    p[n + 2][n + 2].x = 1.0;
    p[n + 2][n + 2].y = 1.0;
    //矩形四条边的上点的坐标
    for (i = 2; i < n + 2; ++i)
    {
        p[i][1] = a[i - 1];//y为0的边（下）
        p[i][n + 2] = b[i - 1];//x为1的边（右）
        p[1][i] = c[i - 1];//x为0的边（左）
        p[n + 2][i] = d[i -1];//y为1的边（上）
    }
    for (i = 2; i < n + 2; ++i)//获取矩阵中间的点的坐标
    {
        for (j = 2; j < n + 2; ++j)
        {
            p[i][j] = GetPoint(p[i][1],p[i][n+2],p[1][j],p[n+2][j]);
        }
    }
}
double GetArea(Point p1, Point p2, Point p3, Point p4)//已知四个点的坐标求四边形的面积
{
    int i, j;
    double ar = 0;
    Point ac[4];
    ac[0] = p1;
    ac[1] = p2;
    ac[2] = p3;
    ac[3] = p4;
    for (i = 0; i < 4; ++i)//叉积求四边形面积
    {
        j = (i + 1) % 4;
        ar += ac[i].x * ac[j].y - ac[i].y * ac[j].x;
    }
    ar = ar / 2.0;//要除以2.0
    if (ar < 0)//面积是正数
    {
        ar = - ar;
    }
    return ar;
}

void GetAreas()
{
    memset(area, 0, sizeof(area));
    int i, j;
    for (i = 1; i <= n + 1; ++i)
    {
        for (j = 1; j <= n + 1; ++j)//枚举四个点组成的四边形的面积
        {
            area[i][j] = GetArea(p[i][j], p[i + 1][j], p[i + 1][j + 1], p[i][j + 1]);
            if (ans < area[i][j])//找出最大的网格面积
            {
                ans = area[i][j];//更新
            }
        }
    }
}

 

计算几何的线段交分为两个步骤：快速排斥和跨立

[以下摘抄]

（1）       快速排斥试验

    设以线段 P1P2 为对角线的矩形为R1， 设以线段 Q1Q2 为对角线的矩形为R2，如果R1和R2不相交，则两线段不会有交点；

 

（2）       跨立试验。

如果两线段相交，则两线段必然相互跨立对方，所谓跨立，指的是一条线段的两端点分别位于另一线段所在直线的两边。判断是否跨立，还是要用到矢量叉积的几何意义。以图3为例，若P1P2跨立Q1Q2 ，则矢量 ( P1 - Q1 ) 和( P2 - Q1 )位于矢量( Q2 - Q1 ) 的两侧，即：

 

( P1 - Q1 ) × ( Q2 - Q1 ) * ( P2 - Q1 ) × ( Q2 - Q1 ) < 0

 

上式可改写成：

 

( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) > 0

 

当 ( P1 - Q1 ) × ( Q2 - Q1 ) = 0 时，说明线段P1P2和Q1Q2共线（但是不一定有交点）。同理判断Q1Q2跨立P1P2的依据是：

 

( Q1 - P1 ) × ( P2 - P1 ) * ( Q2 - P1 ) × ( P2 - P1 ) < 0

 

具体情况如下图所示：

 

图3 直线段跨立试验示意图

 

        根据矢量叉积的几何意义，跨立试验只能证明线段的两端点位于另一个线段所在直线的两边，但是不能保证是在另一直线段的两端，因此，跨立试验只是证明两条线段有交点的必要条件，必需和快速排斥试验一起才能组成直线段相交的充分必要条件。根据以上分析，两条线段有交点的完整判断依据就是：1)以两条线段为对角线的两个矩形有交集；2)两条线段相互跨立。

        判断直线段跨立用计算叉积算法的CrossProduct()函数即可，还需要一个判断两个矩形是否有交的算法。矩形求交也是最简单的求交算法之一，原理就是根据两个矩形的最大、最小坐标判断。图4展示了两个矩形没有交集的各种情况：

 

图4 矩形没有交集的几种情况

 

图5展示了两个矩形有交集的各种情况：

图5 矩形有交集的几种情况

 

从图4和图5可以分析出来两个矩形有交集的几何坐标规律，就是在x坐标方向和y坐标方向分别满足最大值最小值法则，简单解释这个法则就是每个矩形在每个方向上的坐标最大值都要大于另一个矩形在这个坐标方向上的坐标最小值，否则在这个方向上就不能保证一定有位置重叠。由以上分析，判断两个矩形是否相交的算法就可以如下实现：

186 bool IsRectIntersect(const Rect& rc1,const Rect& rc2) 187 { 188     return ( (std::max(rc1.p1.x, rc1.p2.x)>= std::min(rc2.p1.x, rc2.p2.x)) 189              && (std::max(rc2.p1.x, rc2.p2.x)>= std::min(rc1.p1.x, rc1.p2.x)) 190              && (std::max(rc1.p1.y, rc1.p2.y)>= std::min(rc2.p1.y, rc2.p2.y)) 191              && (std::max(rc2.p1.y, rc2.p2.y)>= std::min(rc1.p1.y, rc1.p2.y))); 192 }

 

        完成了排斥试验和跨立试验的算法，最后判断直线段是否有交点的算法就水到渠成了：

204 bool IsLineSegmentIntersect(const LineSeg& ls1,const LineSeg& ls2) 205 { 206     if(IsLineSegmentExclusive(ls1, ls2))//排斥实验 207     { 208         return false; 209     } 210     //( P1 - Q1 ) ×'a1?( Q2 - Q1 ) 211     double p1xq = CrossProduct(ls1.ps.x- ls2.ps.x, ls1.ps.y- ls2.ps.y, 212                                ls2.pe.x- ls2.ps.x, ls2.pe.y- ls2.ps.y); 213     //( P2 - Q1 ) ×'a1?( Q2 - Q1 ) 214     double p2xq = CrossProduct(ls1.pe.x- ls2.ps.x, ls1.pe.y- ls2.ps.y, 215                                ls2.pe.x- ls2.ps.x, ls2.pe.y- ls2.ps.y); 216  217     //( Q1 - P1 ) ×'a1?( P2 - P1 ) 218     double q1xp = CrossProduct(ls2.ps.x- ls1.ps.x, ls2.ps.y- ls1.ps.y, 219                                ls1.pe.x- ls1.ps.x, ls1.pe.y- ls1.ps.y); 220     //( Q2 - P1 ) ×'a1?( P2 - P1 ) 221     double q2xp = CrossProduct(ls2.pe.x- ls1.ps.x, ls2.pe.y- ls1.ps.y, 222                                ls1.pe.x- ls1.ps.x, ls1.pe.y- ls1.ps.y); 223  224     //跨立实验 225     return ( (p1xq* p2xq <= 0.0) && (q1xp * q2xp <=0.0) ); 226 }

 IsLineSegmentExclusive()函数就是调用IsRectIntersect()函数根据结果做排斥判断，此处不再列出代码。