C++几何库.docx

1、C+几何库#include #include Geometry.h#include math.h#include Overlap.h#define max(x,y) (x)(y)?(x):(y)#define min(x,y) (x)(y)?(x):(y)/* 常量定义 */const double INF = 1E200;const double EP = 1E-10;const int MAXV = 300;const double PI = 3.14159265;/浮点误差的处理int dblcmp(double d)if(fabs(d)0) ?1 :-1 ;/ 返回两点之间欧氏距离do

2、uble dist(d_POINT p1,d_POINT p2)return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) );/ 判断两个点是否重合bool equal_point(d_POINT p1,d_POINT p2) return ( (abs(p1.x-p2.x)EP)&(abs(p1.y-p2.y)0:sp在矢量op ep的顺时针方向；r=0：op sp ep三点共线；r0: sp在矢量op ep的逆时针方向 */double multiply(d_POINT sp,d_POINT ep,d_POINT op)

3、return(sp.x-op.x)*(ep.y-op.y) - (ep.x-op.x)*(sp.y-op.y);double amultiply(d_POINT sp,d_POINT ep,d_POINT op)return fabs(sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y);/*矢量(p1-op)和(p2-op)的点积r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积如果两个矢量都非零矢量r 0: 两矢量夹角为钝角 */double dotmultiply(d_POINT p1,d_POINT p2

4、,d_POINT p0) return (p1.x-p0.x)*(p2.x-p0.x) + (p1.y-p0.y)*(p2.y-p0.y);/* 判断点p是否在线段l上条件：(p在线段l所在的直线上)& (点p在以线段l为对角线的矩形内) */bool online(LINESEG l,d_POINT p)return (multiply(l.e, p, l.s)=0)& ( ( (p.x-l.s.x) * (p.x-l.e.x) =0 ) & ( (p.y-l.s.y)*(p.y-l.e.y) = 1.0 ) return 0; if (cosfi 0) return fi;/ 说明矢量os

5、 在矢量 oe的顺时针方向 return -fi;/* 判断点C在线段AB所在的直线l上垂足P的与线段AB的关系本函数是根据下面的公式写的，P是点C到线段AB所在直线的垂足AC dot ABr = -|AB|2(Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)= -L2r has the following meaning:r=0 P = Ar=1 P = Br1 P is on the forward extension of AB0r1 P is interior to AB*/double relation(d_POINT c,LINESEG l)LINESEG tl;tl.

6、s=l.s;tl.e=c;return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e);/ 求点C到线段AB所在直线的垂足 Pd_POINT perpendicular(d_POINT p,LINESEG l)double r=relation(p,l);d_POINT tp;tp.x=l.s.x+r*(l.e.x-l.s.x);tp.y=l.s.y+r*(l.e.y-l.s.y);return tp;/* 求点p到线段l的最短距离返回线段上距该点最近的点np 注意：np是线段l上到点p最近的点，不一定是垂足 */double pt

7、olinesegdist(d_POINT p,LINESEG l,d_POINT &np) double r=relation(p,l); if(r1) np=l.e; return dist(p,l.e); np=perpendicular(p,l); return dist(p,np);/ 求点p到线段l所在直线的距离/请注意本函数与上个函数的区别double ptoldist(d_POINT p,LINESEG l) return abs(multiply(p,l.e,l.s)/dist(l.s,l.e);/* 计算点到折线集的最近距离,并返回最近点.注意：调用的是ptolineseg(

8、)函数 */double ptopointset(int vcount, d_POINT pointset, d_POINT p, d_POINT &q) int i; double cd=double(INF),td; LINESEG l; d_POINT tq,cq; for(i=0;ivcount-1;i+) l.s=pointseti; l.e=pointseti+1; td=ptolinesegdist(p,l,tq); if(tdcd) cd=td; cq=tq; q=cq; return cd;/* 判断圆是否在多边形内*/bool CircleInsidePolygon(int

9、 vcount,d_POINT center,double radius,d_POINT polygon)d_POINT q;double d;q.x=0;q.y=0;d=ptopointset(vcount,polygon,center,q);if(dradius|fabs(d-radius)=min(v.s.x,v.e.x)& /排斥实验(max(v.s.x,v.e.x)=min(u.s.x,u.e.x)&(max(u.s.y,u.e.y)=min(v.s.y,v.e.y)&(max(v.s.y,v.e.y)=min(u.s.y,u.e.y)&(multiply(v.s,u.e,u.s)*

10、multiply(u.e,v.e,u.s)=0)& /跨立实验(multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)=0);/ 判断线段u和v相交（不包括双方的端点）bool intersect_A(LINESEG u,LINESEG v)return (intersect(u,v) &(!online(u,v.s) &(!online(u,v.e) &(!online(v,u.e) &(!online(v,u.s);/ 判断线段v所在直线与线段u相交/方法：判断线段u是否跨立线段vbool intersect_l(LINESEG u,LINESEG v)re

11、turn multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)=0;/ 已知两点坐标p1和p2，求过这两点的直线解析方程： a*x+b*y+c = 0 (a = 0)LINE makeLine(d_POINT p1,d_POINT p2) LINE tl; int sign = 1; tl.a=p2.y-p1.y; if(tl.a0) sign = -1; tl.a=sign*tl.a; tl.b=sign*(p1.x-p2.x); tl.c=sign*(p1.y*p2.x-p1.x*p2.y); return tl;/ 根据直线解析方程返回直线的斜率k,水平

12、线返回 0,竖直线返回 1e200double slope(LINE l)if(abs(l.a) 1e-20)return 0;if(abs(l.b) 1e-20)return INF;return -(l.a/l.b);/ 返回直线的倾斜角alpha ( 0 - pi)/ 注意：atan()返回的是 -PI/2 PI/2double alpha(LINE l)if(abs(l.a) EP)return 0;if(abs(l.b)0)return atan(k);elsereturn PI+atan(k);/ 求点p关于直线l的对称点d_POINT symmetry(LINE l,d_POIN

13、T p)d_POINT tp;tp.x=(l.b*l.b-l.a*l.a)*p.x-2*l.a*l.b*p.y-2*l.a*l.c)/(l.a*l.a+l.b*l.b);tp.y=(l.a*l.a-l.b*l.b)*p.y-2*l.a*l.b*p.x-2*l.b*l.c)/(l.a*l.a+l.b*l.b);return tp;/ 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交，返回true，且返回交点pbool lineInterSect(LINE l1,LINE l2,d_POINT &p) / 是 L1，L2 double d=l1

14、.a*l2.b-l2.a*l1.b; if(abs(d)EP) / 不相交 return false; p.x = (l2.c*l1.b-l1.c*l2.b)/d; p.y = (l2.a*l1.c-l1.a*l2.c)/d; return true;/ 如果线段l1和l2相交，返回true且交点由(inter)返回，否则返回falsebool intersection(LINESEG l1,LINESEG l2,d_POINT &inter) LINE ll1,ll2; ll1=makeLine(l1.s,l1.e); ll2=makeLine(l2.s,l2.e); if(lineInte

15、rSect(ll1,ll2,inter) return online(l1,inter); else return false;/如果无特别说明，输入多边形顶点要求按逆时针排列/ 返回多边形面积(signed)；/ 输入顶点按逆时针排列时，返回正值；否则返回负值double area_of_polygon(int vcount,d_POINT polygon)int i;double s;if (vcount3)return 0;s=polygon0.y*(polygonvcount-1.x-polygon1.x);for (i=1;i0;/*射线法判断点q与多边形polygon的位置关系要求

16、polygon为简单多边形，顶点时针排列如果点在多边形内： 返回0如果点在多边形边上：返回1如果点在多边形外： 返回2 */int insidepolygon(int vcount,d_POINT Polygon, d_POINT q) int c=0,i; LINESEG l1,l2; l1.s=q; l1.e=q;l1.e.x=double(INF); /n=vcount; for (i=0;il2.e.y) c+;/l2的一个端点在l1上且该端点是两端点中纵坐标较大的那个 if(!online(l1,l2.e)& online(l1,l2.s) & l2.e.yl2.e.y) c+;/忽

17、略平行边 if(c%2 = 1) return 0; else return 2;/判断点q在凸多边形polygon内/ 点q是凸多边形polygon内包括边上时，返回true/ 注意：多边形polygon一定要是凸多边形bool InsideConvexPolygon(int vcount,d_POINT polygon,d_POINT q) d_POINT p; LINESEG l; int i; p.x=0; p.y=0; for(i=0;ivcount;i+) / 寻找一个肯定在多边形polygon内的点p：多边形顶点平均值 p.x+=polygoni.x; p.y+=polygoni

18、.y; p.x /= vcount; p.y /= vcount; for(i=0;ivcount;i+) l.s=polygoni; l.e=polygon(i+1)%vcount; if(multiply(p,l.e,l.s)*multiply(q,l.e,l.s)0) /* 点p和点q在边l的两侧，说明点q肯定在多边形外 */ return false; return true;/*寻找凸包的graham 扫描法PointSet为输入的点集；ch为输出的凸包上的点集，按照逆时针方向排列;n为PointSet中的点的数目len为输出的凸包上的点的个数 */void Graham_scan(

19、d_POINT PointSet,d_POINT ch,int n,int &len)int i,j,k=0,top=2;d_POINT tmp;/ 选取PointSet中y坐标最小的点PointSetk，如果这样的点有多个，则取最左边的一个for(i=1;in;i+)if ( PointSeti.yPointSetk.y | (PointSeti.y=PointSetk.y)& (PointSeti.xPointSetk.x) )k=i;tmp=PointSet0;PointSet0=PointSetk;PointSetk=tmp; / 现在PointSet中y坐标最小的点在PointSet

20、0for (i=1;in-1;i+) /* 对顶点按照相对PointSet0的极角从小到大进行排序，极角相同的按照距离PointSet0从近到远进行排序 */k=i;for (j=i+1;j0 | / 极角更小(multiply(PointSetj,PointSetk,PointSet0)=0) & /*极角相等，距离更短 */ dist(PointSet0,PointSetj)dist(PointSet0,PointSetk) )k=j;tmp=PointSeti;PointSeti=PointSetk;PointSetk=tmp;ch0=PointSet0;ch1=PointSet1;ch

21、2=PointSet2;for (i=3;i=0) top-;ch+top=PointSeti;len=top+1;/ 求凸多边形的重心,要求输入多边形按逆时针排序d_POINT gravitycenter(int vcount,d_POINT polygon)d_POINT tp;double x,y,s,x0,y0,cs,k;x=0;y=0;s=0;for(int i=1;ivcount-1;i+)x0=(polygon0.x+polygoni.x+polygoni+1.x)/3;y0=(polygon0.y+polygoni.y+polygoni+1.y)/3; /求当前三角形的重心cs

22、=multiply(polygoni,polygoni+1,polygon0)/2;/三角形面积可以直接利用该公式求解if(abs(s)1e-20)x=x0;y=y0;s+=cs;continue;k=cs/s; /求面积比例x=(x+k*x0)/(1+k);y=(y+k*y0)/(1+k);s += cs;tp.x=x;tp.y=y;return tp;/*所谓凸多边形的直径，即凸多边形任两个顶点的最大距离。下面的算法仅耗时O(n)，是一个优秀的算法。 输入必须是一个凸多边形，且顶点必须按顺序（顺时针、逆时针均可）依次输入。若输入不是凸多边形而是一般点集，则要先求其凸包。 就是先求出所有跖对，然后求出每个跖对的距离，取最大者。点数要多于5个*/void Diameter(d_POINT ch,int n,double &dia) int znum=0,i,j,k=1; int zdMAXV2; double tmp