POJ 1584：A Round Peg in a Ground Hole _判断点是否在多边形内（5种方法）

Description

The DIY Furniture company specializes in assemble-it-yourself furniture kits. Typically, the pieces of wood are attached to one another using a wooden peg that fits into pre-cut holes in each piece to be attached. The pegs have a circular cross-section and so are intended to fit inside a round hole.
A recent factory run of computer desks were flawed when an automatic grinding machine was mis-programmed. The result is an irregularly shaped hole in one piece that, instead of the expected circular shape, is actually an irregular polygon. You need to figure out whether the desks need to be scrapped or if they can be salvaged by filling a part of the hole with a mixture of wood shavings and glue.
There are two concerns. First, if the hole contains any protrusions (i.e., if there exist any two interior points in the hole that, if connected by a line segment, that segment would cross one or more edges of the hole), then the filled-in-hole would not be structurally sound enough to support the peg under normal stress as the furniture is used. Second, assuming the hole is appropriately shaped, it must be big enough to allow insertion of the peg. Since the hole in this piece of wood must match up with a corresponding hole in other pieces, the precise location where the peg must fit is known.
Write a program to accept descriptions of pegs and polygonal holes and determine if the hole is ill-formed and, if not, whether the peg will fit at the desired location. Each hole is described as a polygon with vertices (x1, y1), (x2, y2), . . . , (xn, yn). The edges of the polygon are (xi, yi) to (x _i+1, y _i+1) for i = 1 . . . n − 1 and (xn, yn) to (x1, y1).

Input

Input consists of a series of piece descriptions. Each piece description consists of the following data:
Line 1 < nVertices > < pegRadius > < pegX > < pegY >
number of vertices in polygon, n (integer)
radius of peg (real)
X and Y position of peg (real)
n Lines < vertexX > < vertexY >
On a line for each vertex, listed in order, the X and Y position of vertex The end of input is indicated by a number of polygon vertices less than 3.

Output

For each piece description, print a single line containing the string:
HOLE IS ILL-FORMED if the hole contains protrusions
PEG WILL FIT if the hole contains no protrusions and the peg fits in the hole at the indicated position
PEG WILL NOT FIT if the hole contains no protrusions but the peg will not fit in the hole at the indicated position

Sample Input

5 1.5 1.5 2.0
1.0 1.0
2.0 2.0
1.75 2.0
1.0 3.0
0.0 2.0
5 1.5 1.5 2.0
1.0 1.0
2.0 2.0
1.75 2.5
1.0 3.0
0.0 2.0
1

Sample Output

HOLE IS ILL-FORMED
PEG WILL NOT FIT

Source

Mid-Atlantic 2003

//题意：给出一些点组成的多边形（顺逆情况不知），判断是否为凸多边形！ 并且在多边形内有一个钉子（圆），判断这个圆是否完全在多边形内。 本题判断点是否在多边形内是重点，下面给出判断点在凸包和任意多边形内的5种方法。

//对于初学计算几何的人来说如何选择好的计算方法和精度控制是件头疼的事。本题wrong了好久，也是栽在了精度控制上！

//本题主要分三个步骤：1、多边形为凸多边形。2、点在凸多边形内（或边上）；3、圆在凸多边形内。 

其中在判断点在凸包内（本题规定了只有是凸多边形才需判定点是否在其内部）的方法各异，而如果并没有规定是凸包，只是判断点是否在多边形内，则更复杂些！ 用了一天的时间研究归纳了判断点是否在凸包内和点是否在一般多边形内的5种方法。

//法1、2：叉积判定、面积法判定（适用于凸多边形）。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<cstdlib>
#include<queue>
#include<stack>
#include<map>
#include<vector>
#include<algorithm>
using namespace std;
#define maxn 10005
#define eps 1e-8
#define max(x,y) (x>y?x:y)
#define min(x,y) (x<y?x:y)

int Fabs(double d) //重点:精度控制，如果d精度很高，如-0.00000000001即使是小于0，但也当做是0，关系到后面数据处理
{
	if(fabs(d)<eps)	return 0;
	else return d>0?1:-1;
}
struct point 
{	
	double x,y;
	bool operator == (const point& p)
	{
		return Fabs(x-p.x)==0&&Fabs(y-p.y)==0;
	}
}p[maxn];
int n;
double pegx,pegy,pegr,max_x,max_y;

double x_multi(point p1,point p2,point p3)  
{
	return (p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y);
}

bool point_is_inside() //叉积判断点在凸包内部！ 只针对于凸多边形。圆心连接每一条边的端点得到的叉积必须同向。 以此可以延伸出面积法判定点是否在凸包内部。这两种方法都局限于在凸多边形
{
	point p1;
	p1.x=pegx,p1.y=pegy;
	int i,flag=1;
	double tmp1=0.0,tmp2;
	for(i=0;i<n;i++)
	{
		tmp2=Fabs(x_multi(p1,p[i],p[(i+1)%n]));
		if(tmp1*tmp2<-eps)
		{
			flag=0;
			break;
		}	
		tmp1=tmp2;
	}
	return flag;
}

double Len_ab(point p1,point p2)
{
	return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}

bool circle_is_inside()  //判断圆是否在凸包内
{
	if(pegr==0.0)
		return true;
	int i;
	double ans;
	point peg,t;
	peg.x=pegx,peg.y=pegy;
	for(i=0;i<n;i++) //判断圆心到多边形的每一条边的最短距离是否不小于半径
	{
		t=peg;
		t.x+=p[i].y-p[(i+1)%n].y;
		t.y+=p[(i+1)%n].x-p[i].x;
		if(Fabs(x_multi(peg,t,p[i])*x_multi(peg,t,p[(i+1)%n]))==1) //如果垂足不在线段上则选择到线段两个端点距离较小的
			ans=Fabs(Len_ab(peg,p[i])-Len_ab(peg,p[(i+1)%n]))==-1?Len_ab(peg,p[i]):Len_ab(peg,p[(i+1)%n]);
		else
			ans=fabs(x_multi(peg,p[(i+1)%n],p[i]))/Len_ab(p[i],p[(i+1)%n]); // 否则利用面积/底边得到最小距离
		if(ans-pegr<-eps)
			return false;
	}
	return true;
}

int main()
{
	int i,j;
	while(scanf("%d",&n)&&n>=3)
	{
		scanf("%lf%lf%lf",&pegr,&pegx,&pegy);
		for(i=0;i<n;i++)
			scanf("%lf%lf",&p[i].x,&p[i].y);
		double tmp1=0.0,tmp2;
		bool flag=true;
		for(i=0;i<n;i++)  //判断是否为凸包，即叉积一直是顺时针或一直是逆时针 
		{
			tmp2=Fabs(x_multi(p[i],p[(i+1)%n],p[(i+2)%n])); //精度控制，否则一直wrong
			if(tmp1*tmp2<-eps)
			{
				flag=false;
				break;
			}
			tmp1=tmp2;	
		}
		if(!flag)
		{
			puts("HOLE IS ILL-FORMED");
			continue;
		}
		
		if(!point_is_inside()||!circle_is_inside()) //判断圆是否在凸多边形内 
		{
			puts("PEG WILL NOT FIT");
			continue;
		}
		puts("PEG WILL FIT");
	}
	return 0;
}

//法3：射线法判定点是否在多边形内部（适用于任意多边形）：做一条水平射线计算与多边形的交点个数num，如果num&1则表示在多边形内，否则在多边形外。 其中射线正好交与多边形端点或者与多边形的边平行需要特判。  适用于任意多边形。

注：以下的处理方法以及代码精确性有待考究，自己分析写的，虽然AC了，但AC也是相对的。个人觉得没有绝对的模板，看是否遇上刁钻的数据，所以更倾向于自己尝试着写。期待有人指出程序不足！

bool Onsegment(point p1,point p2,point p3)  	
{
	double min_x=min(p1.x,p2.x);
	double min_y=min(p1.y,p2.y);
	double max_x=max(p1.x,p2.x);
	double max_y=max(p1.y,p2.y);
	if(p3.x>=min_x&&p3.x<=max_x&&p3.y>=min_y&&p3.y<=max_y)
		return true;
	return false;
}

bool Is_intersected(point p1,point p2,point p3,point p4)  //线段相交
{
	double d1=x_multi(p1,p2,p3);
	double d2=x_multi(p1,p2,p4);
	double d3=x_multi(p3,p4,p1);
	double d4=x_multi(p3,p4,p2);
	if(d1*d2<0.0&&d3*d4<0.0)
		return true;
//	if(d1==0.0&&Onsegment(p1,p2,p3))   //由于前面的特判，低处的交点不作为计算
	//	return true;
	if(d2==0.0&&Onsegment(p1,p2,p4))
		return true;
	if(d3==0.0&&Onsegment(p3,p4,p1))
		return true;
	if(d4==0.0&&Onsegment(p3,p4,p2))
		return true;
	return false;
}

double Dot(point p1,point p2,point p3) //点积
{
	return (p2.x-p1.x)*(p3.x-p1.x)+(p2.y-p1.y)*(p3.y-p1.y);
}

int pointonsegment(point p0,point p1,point p2) //判断点是否在线段上
{
	return Fabs(x_multi(p0,p1,p2))==0&&Fabs(Dot(p0,p1,p2))<=0;
}

bool point_is_inside()
{
	int i,num=0;
	point p1,peg,p2,p3;
	p1.x=999999999.0,p1.y=pegy,peg.x=pegx,peg.y=pegy;  //p1坐为在射线极远处的一个点，可以将射线看做线段
	for(i=0;i<n;i++)
	{
 		if(p[i].y==p[(i+1)%n].y)	//如果和多边形的边平行，则判断起点是否在多边形的该边上，避免了和边重合算作无数多个点
		{
			if(pointonsegment(peg,p[i],p[(i+1)%n]))  //判断点是否在线段上
		 		return true;
		}
		else 
		{
			p2=p[i],p3=p[(i+1)%n];
			if(p2.y>p3.y)  //画图知在与多边形端点相交的时候直接计算交的次数都无法直接判定是否在多边形的内外。一条线段的端点有高低之分，此时规定高点的交点为有效交点
				swap(p2,p3);
			if(Is_intersected(peg,p1,p2,p3))
				num++;
		}
	} 
	return num%2==1;
}

//法4：角度和判定法，适用于任意多边形（对自己的解释和做法已做修改，结合黑书知识写的AC代码，）。如果点在多边形内则点连接多边形每条边的到的角度*叉积和为360. 在边上的话为180.

double x_multi(point p1,point p2,point p3)  
{
	return (p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y);
}

double Len_ab(point p1,point p2)
{
	return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}

double Dot(point p1,point p2,point p3)
{
	return (p2.x-p1.x)*(p3.x-p1.x)+(p2.y-p1.y)*(p3.y-p1.y);
}

int pointonsegment(point p1,point p2,point p3)
{
	return Fabs(x_multi(p1,p2,p3))==0&&Fabs(Dot(p1,p2,p3))<=0;
}

bool point_is_inside()
{
	int i;
	double sum=0.0;
	for(i=0;i<n;i++)
	{
		if(p0==p[i]) //点在多边形端点上
			return true;
		if(p[i]==p[(i+1)%n]) //去重点
			continue;
		if(pointonsegment(p0,p[i],p[(i+1)%n]))  //点在多边形边上
			return true;
		double a=Len_ab(p0,p[i]);
		double b=Len_ab(p0,p[(i+1)%n]);
		double c=Len_ab(p[i],p[(i+1)%n]);
		sum+=Fabs(x_multi(p0,p[i],p[(i+1)%n]))*acos((a*a+b*b-c*c)/(2.0*a*b)); //计算角度和，叉积大于0则加上，小于0则减去
	}
	sum=fabs(sum);
	if(Fabs(sum-2.0*pi)==0)
		return true;
	return false;
}

//法5：改进弧长法——权威算法！精度很高！ 以下对于该算法的解析摘自 http://hi.baidu.com/zhuangxie1013/item/7b7e1443abde1a9e833ae1ae

该算法只需O(1)的附加空间,时间复杂度为O(n),系数很小;最大的优点是具有很高的精度,只需做乘法和减法,若针对整数坐标则完全没有精度问题.而且实现起来也非常简单,比转角法和射线法都要好写且不易出错.
   有关"弧长法"的介绍:"弧长法要求多边形是有向多边形,一般规定沿多边形的正向,边的左侧为多边形的内侧域.以被测点为圆心作单位圆,将全部有向边向单位圆作径向投影,并计算其中单位圆上弧长的代数和,若代数和为0,则点在多边形外部;若代数和为2π,则点在多边形内部;若代数和为π,则点在多边形上."
   根据上面的介绍,其实弧长法就是转角法,但它的改进方法比较比较厉害:将坐标原点平移到被测点P,这个新坐标系将平面划分为4个象限,对每个多边形顶点P,只考虑其所在的象限，然后按邻接顺序访问多边形的各个顶点P,分析P[i]和P[i+1],有下列三种情况:
      (1) P[i+1]在P[i]的下一象限,此时弧长和加π/2;
     （2）P[i+1]在P[i]的上一象限,此时弧长和减π/2;
     （3）P[i+1]在P[i]的相对象限,首先计算f=p[i+1].y*p[i].x-p[i+1].x*p[i].y(叉积),若f=0,则点在多边形上;若f<0,弧长和减π;若f>0,弧长和加π.最后对算出的代数和和上述的情况一样判断即可.
实现的时候还有两点要注意:
      1> 若P的某个坐标为0时,一律当正号处理;
      2> 若被测点和多边形的顶点重合时要特殊处理.
   还有一个问题那就是当多边形的某条边在坐标轴上而且两个顶点分别在原点的两侧时会出错,如边(3,0)……>(-3,0),按以上的处理,象限分别是第一和第二,这样会使代数和加π/2,有可能导致最后结果是被测点在多边形外,而实际上被测点是在多边形上(该边穿过该点).
   对于这点,处理办法是:每次计算P[i]和P[i+1]时,就计算叉积和点积,判断该点是否在该边上,是则判断结束,否则继续上述过程,这样牺牲了时间,但保证了正确性.
   具体实现的时候,由于只需知道当前点和上一点的象限位置,所以附加空间只需O(1).实现的时候可以把上述的"π/2"改成1,"π"改成2,这样便可以完全使用整数进行计算,不必考虑顶点的顺序,逆时针和顺时针都可以处理,只是最后的代数和符号不同而已

int get_tmp(point p0)
{
	return p0.x>=0?(p0.y>=0?0:3):(p0.y>=0?1:2);
}
bool point_is_inside()
{
	int tmp1,tmp2,sum=0,i;
	point p0,p1;
	p0.x=pegx,p0.y=pegy;
	p1.x=p[0].x-p0.x,p1.y=p[0].y-p0.y;
	tmp1=get_tmp(p1);
	
	for(i=0;i<n;i++)
	{
		if(p[i]==p0)
			break;
		int t0=Fabs(x_multi(p0,p[i],p[(i+1)%n]));
		int t1=Fabs((p[i].x-p0.x)*(p[(i+1)%n].x-p0.x));
		int t2=Fabs((p[i].y-p0.y)*(p[(i+1)%n].y-p0.y));
		
		if(!t0&&t1<=0&&t2<=0)  //被测点在多边形边上
			break;
		p1.x=p[(i+1)%n].x-p0.x,p1.y=p[(i+1)%n].y-p0.y;
		tmp2=get_tmp(p1);  //计算象限
		switch((tmp2-tmp1+4)%4)
		{
			case 1:{ sum++; break; }
			case 2:
			{
				if(t0>0) sum+=2;
				else sum-=2;
				break;
			}
			case 3: { sum--; break; }
		}
		tmp1=tmp2;
	}
	if(i<n||sum) //被测点在多边形边上或者在多边形内部
		return true;
	return false;
}