計算幾何-求線段交點算法和代碼(C++語言)

問題描述

問題描述:已知兩條線段P1P2和Q1Q2,判斷P1P2和Q1Q2是否相交,若相交,求出交點。html

兩條線段的位置關係能夠分爲三類:有重合部分、無重合部分但有交點、無交點。算法

clipboard.png

算法核心

算法的步驟以下:ui

1.快速排斥實驗。

設以線段P1P2爲對角線的矩形爲R,設以線段Q1Q2爲對角線的矩形爲T,若是R和T不相交,則兩線段不相交。spa

因此P1P2和Q1Q2相交的必要條件是以他們爲對角線的矩形相交,即:3d

min(p1.x,p2.x) <= max(q1.x,q2.x) &&
min(q1.x,q2.x) <= max(p1.x,p2.x) &&
min(p1.y,p2.y) <= max(q1.y,q2.y) &&
min(q1.y,q2.y) <= max(p1.y,p2.y);

2.跨立實驗。

若是兩線段相交,則兩線段必然相互跨立對方。
線段的跨立到底是什麼意思?向量的跨立是什麼意思?
clipboard.pngcode

a、若P1P2跨立Q1Q2,則矢量(P1-Q1)和(P2-Q1)位於矢量(Q2-Q1)的兩側,即( P1 - Q1 ) × ( Q2 - Q1 ) * ( P2 - Q1 ) × ( Q2 - Q1 ) < 0。htm

等價於
(Q1.x-P1.x,Q1.y-P1.y) × ( Q1.x-Q2.x,Q1.y-Q2.y ) * ( Q1.x-P2.x,Q1.y-P2.y ) × ( Q1.x-Q2.x,Q1.y-Q2.y ) < 0
等價於
((Q1.x-P1.x)*(Q1.y-Q2.y)-(Q1.y-P1.y)*( Q1.x-Q2.x)) * ((Q1.x-P2.x)*(Q1.y-Q2.y)-(Q1.y-P2.y)*(Q1.x-Q2.x)) < 0blog

b、若Q1Q2跨立P1P2,則矢量(Q1-P1)和(Q2-P1)位於矢量(P2-P1)的兩側,即( Q1 - P1 ) × ( P2 - P1 ) * ( Q2 - P1 ) × ( P2 - P1 ) < 0。
等價於
(P1.x-Q1.x,P1.y-Q1.y) × ( P1.x-P2.x,P1.y-P2.y ) * ( P1.x-Q2.x,P1.y-Q2.y ) × ( P1.x-P2.x,P1.y-P2.y ) < 0
等價於
((P1.x-Q1.x)*(P1.y-P2.y)-(P1.y-Q1.y)*(P1.x-P2.x)) * ((P1.x-Q2.x)*(P1.y-P2.y)-(P1.y-Q2.y)*( P1.x-P2.x)) < 0ip

a和b兩個不等式只要知足一個便可。get

排斥實驗和跨立實驗的示例以下圖所示。
clipboard.png

代碼

//跨立判斷
bool isLineSegmentCross(const Point &P1,const Point &P2,const Point &Q1,const Point &Q2)
{
    if(
        ((Q1.x-P1.x)*(Q1.y-Q2.y)-(Q1.y-P1.y)*( Q1.x-Q2.x)) * ((Q1.x-P2.x)*(Q1.y-Q2.y)-(Q1.y-P2.y)*(Q1.x-Q2.x)) < 0 ||
        ((P1.x-Q1.x)*(P1.y-P2.y)-(P1.y-Q1.y)*(P1.x-P2.x)) * ((P1.x-Q2.x)*(P1.y-P2.y)-(P1.y-Q2.y)*( P1.x-P2.x)) < 0
    ) 
        return true;
    else
       return false;
}

3.計算交點。

當斷定兩條線段相交後,能夠進行交點的求解,求交點能夠用平面幾何方法,列點斜式方程來完成。但因爲點斜式方程難以處理斜率爲0的特殊狀況,不方便求解。於是,參用向量法求解交點。

設交點爲(x0,y0),則下列方程組成立:

01162126-6005ab3dc2564e79b76d8a41f920dd24.png

根據以上方程組,消除參數k1和k2,獲得以下方程:

01162126-6005ab3dc2564e79b76d8a41f920dd24.png

而後求解(x0,y0),結果以下所示:

01162126-6005ab3dc2564e79b76d8a41f920dd24.png

#include<stdio.h>
#define N 10002

/**
算法適用於整形點,不適用於浮點型 
**/

typedef struct Point
{
    int x;
    int y;
}Point;

double min(int x, int y)
{
    return x<y?x:y;
}

double max(int x, int y)
{
    return x>y?x:y;
}

//排斥實驗
bool IsRectCross(const Point &p1,const Point &p2,const Point &q1,const Point &q2)
{
    bool ret = min(p1.x,p2.x) <= max(q1.x,q2.x)    &&
                min(q1.x,q2.x) <= max(p1.x,p2.x) &&
                min(p1.y,p2.y) <= max(q1.y,q2.y) &&
                min(q1.y,q2.y) <= max(p1.y,p2.y);
    return ret;
}

/**這段代碼不能獲得正確答案,故註釋
//跨立判斷
bool IsLineSegmentCross(const Point &pFirst1,const Point &pFirst2,const Point &pSecond1,const Point &pSecond2)
{
    long line1,line2;
    line1 = pFirst1.x * (pSecond1.y - pFirst2.y) +
        pFirst2.x * (pFirst1.y - pSecond1.y) +
        pSecond1.x * (pFirst2.y - pFirst1.y);
    line2 = pFirst1.x * (pSecond2.y - pFirst2.y) +
        pFirst2.x * (pFirst1.y - pSecond2.y) + 
        pSecond2.x * (pFirst2.y - pFirst1.y);
    if (((line1 ^ line2) >= 0) && !(line1 == 0 && line2 == 0))
        return false;

    line1 = pSecond1.x * (pFirst1.y - pSecond2.y) +
        pSecond2.x * (pSecond1.y - pFirst1.y) +
        pFirst1.x * (pSecond2.y - pSecond1.y);
    line2 = pSecond1.x * (pFirst2.y - pSecond2.y) + 
        pSecond2.x * (pSecond1.y - pFirst2.y) +
        pFirst2.x * (pSecond2.y - pSecond1.y);
    if (((line1 ^ line2) >= 0) && !(line1 == 0 && line2 == 0))
        return false;
    return true;
}
**/

//跨立判斷
bool IsLineSegmentCross(const Point &P1,const Point &P2,const Point &Q1,const Point &Q2)
{
    if(
        ((Q1.x-P1.x)*(Q1.y-Q2.y)-(Q1.y-P1.y)*( Q1.x-Q2.x)) * ((Q1.x-P2.x)*(Q1.y-Q2.y)-(Q1.y-P2.y)*(Q1.x-Q2.x)) < 0 ||
        ((P1.x-Q1.x)*(P1.y-P2.y)-(P1.y-Q1.y)*(P1.x-P2.x)) * ((P1.x-Q2.x)*(P1.y-P2.y)-(P1.y-Q2.y)*( P1.x-P2.x)) < 0
    ) 
        return true;
    else
       return false;
}

/**
求線段P1P2與Q1Q2的交點。
先進行快速排斥實驗和跨立實驗肯定有交點再進行計算。
交點(x,y)使用引用返回。
沒有驗證過
**/
bool GetCrossPoint(const Point &p1,const Point &p2,const Point &q1,const Point &q2,long &x,long &y)
{
    if(IsRectCross(p1,p2,q1,q2))
    {
        if (IsLineSegmentCross(p1,p2,q1,q2))
        {
            //求交點
            long tmpLeft,tmpRight;
            tmpLeft = (q2.x - q1.x) * (p1.y - p2.y) - (p2.x - p1.x) * (q1.y - q2.y);
            tmpRight = (p1.y - q1.y) * (p2.x - p1.x) * (q2.x - q1.x) + q1.x * (q2.y - q1.y) * (p2.x - p1.x) - p1.x * (p2.y - p1.y) * (q2.x - q1.x);

            x = (int)((double)tmpRight/(double)tmpLeft);

            tmpLeft = (p1.x - p2.x) * (q2.y - q1.y) - (p2.y - p1.y) * (q1.x - q2.x);
            tmpRight = p2.y * (p1.x - p2.x) * (q2.y - q1.y) + (q2.x- p2.x) * (q2.y - q1.y) * (p1.y - p2.y) - q2.y * (q1.x - q2.x) * (p2.y - p1.y); 
            y = (int)((double)tmpRight/(double)tmpLeft);
            return true;
        }
    }
    return false;
}

代碼2-線段求交模板
http://www.cppblog.com/wicbnu/archive/2009/08/24/94225.html

#include <stdio.h>
#include <math.h>
const int N = 100010;
int mark[N];
struct Point
{
    double x,y;
};
struct stline
{
    Point a,b;
} line1,line2, p[N];

int dblcmp(double a,double b)
{
    if (fabs(a-b)<=1E-6) return 0;
    if (a>b) return 1;
    else return -1;
}
//***************點積判點是否在線段上***************
double dot(double x1,double y1,double x2,double y2) //點積
{
    return x1*x2+y1*y2;
}

int point_on_line(Point a,Point b,Point c) //求a點是否是在線段bc上,>0不在,=0與端點重合,<0在。
{
    return dblcmp(dot(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y),0);
}
//**************************************************
double cross(double x1,double y1,double x2,double y2)
{
    return x1*y2-x2*y1;
}
double ab_cross_ac(Point a,Point b,Point c) //ab與ac的叉積
{
    return cross(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);
}
int ab_cross_cd (Point a,Point b,Point c,Point d) //求ab是否與cd相交,交點爲p。1規範相交,0交點是一線段的端點,-1不相交。
{
    double s1,s2,s3,s4;
    int d1,d2,d3,d4;
    Point p;
    d1=dblcmp(s1=ab_cross_ac(a,b,c),0);
    d2=dblcmp(s2=ab_cross_ac(a,b,d),0);
    d3=dblcmp(s3=ab_cross_ac(c,d,a),0);
    d4=dblcmp(s4=ab_cross_ac(c,d,b),0);

//若是規範相交則求交點
    if ((d1^d2)==-2 && (d3^d4)==-2)
    {
        p.x=(c.x*s2-d.x*s1)/(s2-s1);
        p.y=(c.y*s2-d.y*s1)/(s2-s1);
        return 1;
    }

//若是不規範相交
    if (d1==0 && point_on_line(c,a,b)<=0)
    {
        p=c;
        return 0;
    }
    if (d2==0 && point_on_line(d,a,b)<=0)
    {
        p=d;
        return 0;
    }
    if (d3==0 && point_on_line(a,c,d)<=0)
    {
        p=a;
        return 0;
    }
    if (d4==0 && point_on_line(b,c,d)<=0)
    {
        p=b;
        return 0;
    }
//若是不相交
    return -1;
}

Circuit Board

http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=1648
http://www.xuebuyuan.com/552687.html

進階算法 Bentley-Ottmann 時間複雜度O(n^2)

相關文章
相關標籤/搜索