ACM经典算法之计算几何

转自：http://blog.sina.com.cn/s/blog_93d2ceba010145ci.html

一、（叉乘法求任意多边形面积）

语法：result=polygonarea(Point *polygon,int N);
参数：
*polygon：	多变形顶点数组
N：	多边形顶点数目
返回值：	多边形面积
注意：
	支持任意多边形，凹、凸皆可
	多边形顶点输入时按顺时针顺序排列
源程序：
	typedef struct { double x,y; } Point; double polygonarea(Point *polygon,int N) { int i,j; double area = 0; for (i=0;i<N;i++) { j = (i + 1) % N; area += polygon[i].x * polygon[j].y; area -= polygon[i].y * polygon[j].x; } area /= 2; return(area < 0 ? -area : area); }

二、（求三角形面积）

语法：result=area3(float x1,float y1,float x2,float y2,float x3,float y3);
参数：
x1～3：	三角形3个顶点x坐标
y1～3：	三角形3个顶点y坐标
返回值：	三角形面积
注意：
	需要 math.h
源程序：
	float area3(float x1,float y1,float x2,float y2,float x3,float y3) { float a,b,c,p,s; a=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); b=sqrt((x1-x3)*(x1-x3)+(y1-y3)*(y1-y3)); c=sqrt((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)); p=(a+b+c)/2; s=sqrt(p*(p-a)*(p-b)*(p-c)); return s; }

三、（两矢量间角度）

语法：result=angle(double x1, double y1, double x2, double y2);
参数：
x/y1～2：	两矢量的坐标
返回值：	两的角度矢量
注意：
	返回角度为弧度制，并且以逆时针方向为正方向
	需要 math.h
源程序：
	#define PI 3.1415926 double angle(double x1, double y1, double x2, double y2) { double dtheta,theta1,theta2; theta1 = atan2(y1,x1); theta2 = atan2(y2,x2); dtheta = theta2 - theta1; while (dtheta > PI) dtheta -= PI*2; while (dtheta < -PI) dtheta += PI*2; return(dtheta); }

四、（两点距离（2D、3D））

语法：result=distance_2d(float x1,float x2,float y1,float y2);
参数：
x/y/z1～2：	各点的x、y、z坐标
返回值：	两点之间的距离
注意：
	需要 math.h
源程序：
	float distance_2d(float x1,float x2,float y1,float y2) { return(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))); } float distance_3d(float x1,float x2,float y1,float y2,float z1,float z2) { return(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2))); }

五、（射向法判断点是否在多边形内部）

语法：result=insidepolygon(Point *polygon,int N,Point p);
参数：
*polygon：	多边形顶点数组
N：	多边形顶点个数
p：	被判断点
返回值：	0：点在多边形内部；1：点在多边形外部
注意：
	若p点在多边形顶点或者边上，返回值不确定，需另行判断
	需要 math.h
源程序：
	#define MIN(x,y) (x < y ? x : y) #define MAX(x,y) (x > y ? x : y) typedef struct { double x,y; } Point; int insidepolygon(Point *polygon,int N,Point p) { int counter = 0; int i; double xinters; Point p1,p2; p1 = polygon[0]; for (i=1;i<=N;i++) { p2 = polygon[i % N]; if (p.y > MIN(p1.y,p2.y)) { if (p.y <= MAX(p1.y,p2.y)) { if (p.x <= MAX(p1.x,p2.x)) { if (p1.y != p2.y) { xinters = (p.y-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x; if (p1.x == p2.x || p.x <= xinters) counter++; } } } } p1 = p2; } if (counter % 2 == 0) return(OUTSIDE); else return(INSIDE); }

六、（判断点是否在线段上）

语法：result=Pointonline(Point p1,Point p2,Point p);
参数：
p1、p2：	线段的两个端点
p：	被判断点
返回值：	0：点在不在线段上；1：点在线段上
注意：
	若p线段端点上返回1
	需要 math.h
源程序：
	#define MIN(x,y) (x < y ? x : y) #define MAX(x,y) (x > y ? x : y) typedef struct { double x,y; } Point; int FC(double x1,double x2) { if (x1-x2<0.000002&&x1-x2>-0.000002) return 1; else return 0; } int Pointonline(Point p1,Point p2,Point p) { double x1,y1,x2,y2; x1=p.x-p1.x; x2=p2.x-p1.x; y1=p.y-p1.y; y2=p2.y-p1.y; if (FC(x1*y2-x2*y1,0)==0) return 0; if ((MIN(p1.x,p2.x)<=p.x&&p.x<=MAX(p1.x,p2.x))&& (MIN(p1.y,p2.y)<=p.y&&p.y<=MAX(p1.y,p2.y))) return 1; else return 0; }

七、（判断两线段是否相交）

语法：result=sectintersect(Point p1,Point p2,Point p3,Point p4);
参数：
p1～4：	两条线段的四个端点
返回值：	0：两线段不相交；1：两线段相交；2两线段首尾相接
注意：
	p1!=p2;p3!=p4;
源程序：
	#define MIN(x,y) (x < y ? x : y) #define MAX(x,y) (x > y ? x : y) typedef struct { double x,y; } Point; int lineintersect(Point p1,Point p2,Point p3,Point p4) { Point tp1,tp2,tp3; if ((p1.x==p3.x&&p1.y==p3.y)||(p1.x==p4.x&&p1.y==p4.y)||(p2.x==p3.x&&p2.y==p3.y)||(p2.x==p4.x&&p2.y==p4.y)) return 2; //快速排斥试验 if ((MIN(p1.x,p2.x)<p3.x&&p3.x<MAX(p1.x,p2.x)&&MIN(p1.y,p2.y)<p3.y<MAX(p1.y,p2.y))|| (MIN(p1.x,p2.x)<p4.x&&p3.x<MAX(p1.x,p2.x)&&MIN(p1.y,p2.y)<p3.y<MAX(p1.y,p2.y))) ;else return 0; //跨立试验 tp1.x=p1.x-p3.x; tp1.y=p1.y-p3.y; tp2.x=p4.x-p3.x; tp2.y=p4.y-p3.y; tp3.x=p2.x-p3.x; tp3.y=p2.y-p3.y; if ((tp1.x*tp2.y-tp1.y*tp2.x)*(tp2.x*tp3.y-tp2.y*tp3.x)>=0) return 1; else return 0; }

八、（判断线段与直线是否相交）

语法：result=lineintersect(Point p1,Point p2,Point p3,Point p4);
参数：
p1、p2：	线段的两个端点
p3、p4：	直线上的两个点
返回值：	0：线段直线不相交；1：线段和直线相交
注意：
	如线段在直线上，返回 1
源程序：
	typedef struct { double x,y; } Point; int lineintersect(Point p1,Point p2,Point p3,Point p4) { Point tp1,tp2,tp3; tp1.x=p1.x-p3.x; tp1.y=p1.y-p3.y; tp2.x=p4.x-p3.x; tp2.y=p4.y-p3.y; tp3.x=p2.x-p3.x; tp3.y=p2.y-p3.y; if ((tp1.x*tp2.y-tp1.y*tp2.x)*(tp2.x*tp3.y-tp2.y*tp3.x)>=0) return 1; else return 0; }

九、（点到线段最短距离）

语法：result=mindistance(Point p1,Point p2,Point q);
参数：
p1、p2：	线段的两个端点
q：	判断点
返回值：	点q到线段p1p2的距离
注意：
	需要 math.h
源程序：
	#define MIN(x,y) (x < y ? x : y) #define MAX(x,y) (x > y ? x : y) typedef struct { double x,y; } Point; double mindistance(Point p1,Point p2,Point q) { int flag=1; double k; Point s; if (p1.x==p2.x) {s.x=p1.x;s.y=q.y;flag=0;} if (p1.y==p2.y) {s.x=q.x;s.y=p1.y;flag=0;} if (flag) { k=(p2.y-p1.y)/(p2.x-p1.x); s.x=(k*k*p1.x+k*(q.y-p1.y)+q.x)/(k*k+1); s.y=k*(s.x-p1.x)+p1.y; } if (MIN(p1.x,p2.x)<=s.x&&s.x<=MAX(p1.x,p2.x)) return sqrt((q.x-s.x)*(q.x-s.x)+(q.y-s.y)*(q.y-s.y)); else return MIN(sqrt((q.x-p1.x)*(q.x-p1.x)+(q.y-p1.y)*(q.y-p1.y)),sqrt((q.x-p2.x)*(q.x-p2.x)+(q.y-p2.y)*(q.y-p2.y))); }

十、（求两直线的交点）

语法：result=mindistance(Point p1,Point p2,Point q);
参数：
p1～p4：	直线上不相同的两点
*p：	通过指针返回结果
返回值：	1：两直线相交；2：两直线平行
注意：
	如需要判断两线段交点，检验k和对应k1（注释中）的值是否在0～1之间，用在0～1之间的那个求交点
源程序：
	typedef struct { double x,y; } Point; int linecorss(Point p1,Point p2,Point p3,Point p4,Point *p) { double k; //同一直线 if ((p4.x-p3.x)*(p1.y-p3.y)-(p4.y-p3.y)*(p1.x-p3.x)==0&& (p2.x-p1.x)*(p1.y-p3.y)-(p2.y-p1.y)*(p1.x-p3.x)==0) return 2; //平行，不同一直线 if ((p4.y-p3.y)*(p2.x-p1.x)-(p4.x-p3.x)*(p2.y-p1.y)==0) return 0; k=((p4.x-p3.x)*(p1.y-p3.y)-(p4.y-p3.y)*(p1.x-p3.x))/((p4.y-p3.y)*(p2.x-p1.x)-(p4.x-p3.x)*(p2.y-p1.y)); //k1=((p2.x-p1.x)*(p1.y-p3.y)-(p2.y-p1.y)*(p1.x-p3.x))/((p4.y-p3.y)*(p2.x-p1.x)-(p4.x-p3.x)*(p2.y-p1.y)); (*p).x=p1.x+k*(p2.x-p1.x); (*p).y=p1.y+k*(p2.y-p1.y); return 1;//有交点}

十一、（判断一个封闭图形是凹集还是凸集）

语法：result=convex(Point *p,int n);
参数：
*p：	封闭曲线顶点数组
n：	封闭曲线顶点个数
返回值：	1：凸集；-1：凹集；0：曲线不符合要求无法计算
注意：
	默认曲线为简单曲线：无交叉、无圈
源程序：
	typedef struct { double x,y; } Point; int convex(Point *p,int n) { int i,j,k; int flag = 0; double z; if (n < 3) return(0); for (i=0;i<n;i++) { j = (i + 1) % n; k = (i + 2) % n; z = (p[j].x - p[i].x) * (p[k].y - p[j].y); z -= (p[j].y - p[i].y) * (p[k].x - p[j].x); if (z < 0) flag |= 1; else if (z > 0) flag |= 2; if (flag == 3) return －1; //CONC***E } if (flag != 0) return 1; //CONVEX else return 0; }

十二、（Graham扫描法寻找凸包）

语法：Graham_scan(Point PointSet[],Point ch[],int n,int &len);
参数：
PointSet[]：	输入的点集
ch[]：	输出的凸包上的点集，按照逆时针方向排列
n：	PointSet中的点的数目
len：	输出的凸包上的点的个数
返回值：	null
源程序：
	struct Point{ float x,y; }; float multiply(Point p1,Point p2,Point p0) { return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y)); } float distance(Point p1,Point p2) { return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y))); } void Graham_scan(Point PointSet[],Point ch[],int n,int &len) { int i,j,k=0,top=2; Point tmp; for(i=1;i<n;i++) if ((PointSet[i].y<PointSet[k].y)||((PointSet[i].y==PointSet[k].y)&&(PointSet[i].x<PointSet[k].x))) k=i; tmp=PointSet[0]; PointSet[0]=PointSet[k]; PointSet[k]=tmp; for (i=1;i<n-1;i++) { k=i; for (j=i+1;j<n;j++) if ( (multiply(PointSet[j],PointSet[k],PointSet[0])>0) || ((multiply(PointSet[j],PointSet[k],PointSet[0])==0) &&(distance(PointSet[0],PointSet[j])<distance(PointSet[0],PointSet[k]))) ) k=j; tmp=PointSet[i]; PointSet[i]=PointSet[k]; PointSet[k]=tmp; } ch[0]=PointSet[0]; ch[1]=PointSet[1]; ch[2]=PointSet[2]; for (i=3;i<n;i++) { while (multiply(PointSet[i],ch[top],ch[top-1])>=0) top--; ch[++top]=PointSet[i]; } len=top+1; }