轉(zhuǎn)角法判斷點(diǎn)和多邊形的關(guān)系大家都知道，原理比較簡(jiǎn)單，在多邊形內(nèi)掃過(guò)的轉(zhuǎn)角一定是360度，在邊界上和外面則不一定。
實(shí)現(xiàn)起來(lái)也比較麻煩，浮點(diǎn)誤差比較大，而且還要考慮些特殊情況。
在網(wǎng)上找到一種叫做改進(jìn)弧長(zhǎng)法的算法，原理和轉(zhuǎn)角法類似，但是做了很多重要的改進(jìn)。比如，計(jì)算轉(zhuǎn)角改成了計(jì)算叉積，根據(jù)叉積決定
旋轉(zhuǎn)方向，還要根據(jù)計(jì)算下一個(gè)點(diǎn)的象限決定偏轉(zhuǎn)多少，每次偏轉(zhuǎn)的都是90度的倍數(shù)。
該算法可以方便判斷出點(diǎn)在多邊形內(nèi)，還是邊界上，還是在多邊形外面。
摘自別人對(duì)該算法的描述如下：
首先從該書(shū)中摘抄一段弧長(zhǎng)法的介紹：“弧長(zhǎng)法要求多邊形是有向多邊形，一般規(guī)定沿多邊形的正向，邊的左側(cè)為多邊形的內(nèi)側(cè)域。
以被測(cè)點(diǎn)為圓心作單位圓，將全部有向邊向單位圓作徑向投影，并計(jì)算其中單位圓上弧長(zhǎng)的代數(shù)和。若代數(shù)和為0，則點(diǎn)在多邊形外部；
若代數(shù)和為2π則點(diǎn)在多邊形內(nèi)部；若代數(shù)和為π，則點(diǎn)在多邊形上。” 按書(shū)上的這個(gè)介紹，其實(shí)弧長(zhǎng)法就是轉(zhuǎn)角法。但它的改進(jìn)方法比較厲害：將坐標(biāo)原點(diǎn)平移到被測(cè)點(diǎn)P，這個(gè)新坐標(biāo)系將平面劃分為4個(gè)
象限，對(duì)每個(gè)多邊形頂點(diǎn)P ，只考慮其所在的象限，然后按鄰接順序訪問(wèn)多邊形的各個(gè)頂點(diǎn)P，分析P和P[i+1]，有下列三種情況：（1）P[i+1]在P的下一象限。此時(shí)弧長(zhǎng)和加π/2；（2）P[i+1]在P的上一象限。此時(shí)弧長(zhǎng)和減π/2；（3）P[i+1]在Pi的相對(duì)象限。首先計(jì)算f=y[i+1]*x-x[i+1]*y（叉積），若f=0，則點(diǎn)在多邊形上；若f<0，弧長(zhǎng)和減π；若f>0，弧長(zhǎng)和加π。 最后對(duì)算出的代數(shù)和和上述的情況一樣判斷即可。 實(shí)現(xiàn)的時(shí)候還有兩點(diǎn)要注意，第一個(gè)是若P的某個(gè)坐標(biāo)為0時(shí)，一律當(dāng)正號(hào)處理；第二點(diǎn)是若被測(cè)點(diǎn)和多邊形的頂點(diǎn)重合時(shí)要特殊處理。
以上就是書(shū)上講解的內(nèi)容，其實(shí)還存在一個(gè)問(wèn)題。那就是當(dāng)多邊形的某條邊在坐標(biāo)軸上而且兩個(gè)頂點(diǎn)分別在原點(diǎn)的兩側(cè)時(shí)會(huì)出錯(cuò)。
如邊(3,0)-(-3,0)，按以上的處理，象限分別是第一和第二，這樣會(huì)使代數(shù)和加π/2，有可能導(dǎo)致最后結(jié)果是被測(cè)點(diǎn)在多邊形外。而實(shí)際上
被測(cè)點(diǎn)是在多邊形上（該邊穿過(guò)該點(diǎn)）。 對(duì)于這點(diǎn)，我的處理辦法是：每次算P和P[i+1]時(shí)，就計(jì)算叉積和點(diǎn)積，判斷該點(diǎn)是否在該邊上，是則判斷結(jié)束，否則繼續(xù)上述過(guò)程。
這樣犧牲了時(shí)間，但保證了正確性。 具體實(shí)現(xiàn)的時(shí)候，由于只需知道當(dāng)前點(diǎn)和上一點(diǎn)的象限位置，所以附加空間只需O(1)。實(shí)現(xiàn)的時(shí)候可以把上述的“π/2”改成1，“π”改成2，
這樣便可以完全使用整數(shù)進(jìn)行計(jì)算。不必考慮頂點(diǎn)的順序，逆時(shí)針和順時(shí)針都可以處理，只是最后的代數(shù)和符號(hào)不同而已。整個(gè)算法編寫(xiě)
起來(lái)非常容易。
代碼如下：
#include <stdio.h>
#include <math.h>
#define MAX (100 + 10)
struct Point
{
double x,y;
};
Point pts[MAX];
const int OUT = 0;
const int IN = 1;
const int EDGE = 2;
const double fPre = 1e-8;
int DblCmp(double fD)
{
if (fabs(fD) < fPre)
{
return 0;
}
else
{
return fD > 0 ? 1 : -1;
}
}
int GetQuadrant(Point p)
{
return DblCmp(p.x) >= 0 ? (DblCmp(p.y) >= 0 ? 0 : 3) :
(DblCmp(p.y) >= 0 ? 1 : 2);
}
double Det(double fX1, double fY1, double fX2, double fY2)
{
return fX1 * fY2 - fX2 * fY1;
}
int PtInPolygon(Point* pts, int nN, Point p)
{
int i, j, k;
for (j = 0; j < nN; ++j)
{
pts[j].x -= p.x;
pts[j].y -= p.y;
}
int nA1, nA2;
int nSum = 0;
nA1 = GetQuadrant(pts[0]);
for (i = 0; i < nN; ++i)
{
k = (i + 1) % nN;
if (DblCmp(pts[k].x) == 0 && DblCmp(pts[k].y) == 0)
{
break;//與頂點(diǎn)重合
}
int nC = DblCmp(Det(pts[i].x, pts[i].y,
pts[k].x, pts[k].y));
if (!nC && DblCmp(pts[i].x * pts[k].x) <= 0
&& DblCmp(pts[i].y * pts[k].y) <= 0)
{
break;//邊上
}
nA2 = GetQuadrant(pts[k]);
if ((nA1 + 1) % 4 == nA2)
{
nSum += 1;
}
else if ((nA1 + 2) % 4 == nA2)
{
if (nC > 0)
{
nSum += 2;
}
else
{
nSum -= 2;
}
}
else if ((nA1 + 3) % 4 == nA2)
{
nSum -= 1;
}
nA1 = nA2;
}
for (j = 0; j < nN; ++j)
{
pts[j].x += p.x;
pts[j].y += p.y;
}
if (i < nN)
{
return EDGE;
}
else if (nSum)//逆時(shí)針nSum == 4, 順時(shí)針nSum == -4
{
return IN;
}
else
{
return OUT;
}
}
int main()
{
int nN, nM;
int nCase = 1;
while (scanf("%d%d", &nN, &nM), nN)
{
if (nCase > 1)
{
printf("\n");
}
for (int i = 0; i < nN; ++i)
{
scanf("%lf%lf", &pts[i].x, &pts[i].y);
}
printf("Problem %d:\n", nCase++);
for (int i = 0; i < nM; ++i)
{
Point p;
scanf("%lf%lf", &p.x, &p.y);
if (PtInPolygon(pts, nN, p))
{
printf("Within\n");
}
else
{
printf("Outside\n");
}
}
}
return 0;
}