点、线、面相关的算法

* 需要包含的头文件 */ ＃i nclude <cmath > /* 常用的常量定义 */ const double INF = 1E200    const double EP = 1E-10 const int MAXV = 300 const double PI = 3.14159265 /* 基本几何结构 */ struct POINT { double x; double y; POINT(double a=0, double b=0) { x=a; y=b;} //constructor }; struct LINESEG { POINT s; POINT e; LINESEG(POINT a, POINT b) { s=a; e=b;} LINESEG() { } }; struct LINE           // 直线的解析方程 a*x+b*y+c=0  为统一表示，约定 a >= 0 {    double a;    double b;    double c; LINE(double d1=1, double d2=-1, double d3=0) {a=d1; b=d2; c=d3;} }; /********************/ *                    * *   点的基本运算     * *                    * /********************/ double dist(POINT p1,POINT p2)                // 返回两点之间欧氏距离 { return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) ); } bool equal_point(POINT p1,POINT p2)           // 判断两个点是否重合  { return ( (abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP) ); } /****************************************************************************** r=multiply(sp,ep,op),得到(sp-op)*(ep-op)的叉积 r>0:ep在矢量opsp的逆时针方向； r=0：opspep三点共线； r<0:ep在矢量opsp的顺时针方向 *******************************************************************************/ double multiply(POINT sp,POINT ep,POINT op) { return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y)); } /******************************************************************************* r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积，如果两个矢量都非零矢量 r<0:两矢量夹角为锐角；r=0：两矢量夹角为直角；r>0:两矢量夹角为钝角 *******************************************************************************/ double dotmultiply(POINT p1,POINT p2,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,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)<=0 ) ) ); } // 返回点p以点o为圆心逆时针旋转alpha(单位：弧度)后所在的位置 POINT rotate(POINT o,double alpha,POINT p) { POINT tp; p.x-=o.x; p.y-=o.y; tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x; tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y; return tp; } /* 返回顶角在o点，起始边为os，终止边为oe的夹角(单位：弧度) 角度小于pi，返回正值 角度大于pi，返回负值 可以用于求线段之间的夹角 */ double angle(POINT o,POINT s,POINT e) { double cosfi,fi,norm; double dsx = s.x - o.x; double dsy = s.y - o.y; double dex = e.x - o.x; double dey = e.y - o.y; cosfi=dsx*dex+dsy*dey; norm=(dsx*dsx+dey*dey)*(dex*dex+dey*dey); cosfi /= sqrt( norm ); if (cosfi >=  1.0 ) return 0;   if (cosfi <= -1.0 ) return -3.1415926; fi=acos(cosfi); if (dsx*dey-dsy*dex>0) return fi;      // 说明矢量os 在矢量 oe的顺时针方向   return -fi; }                           /*****************************/   *                             *   *      线段及直线的基本运算   *   *                             *   /*****************************/ /* 判断点与线段的关系,用途很广泛 本函数是根据下面的公式写的，P是点C到线段AB所在直线的垂足                 AC dot AB         r =     ---------                  ||AB||^2              (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)           = -------------------------------                           L^2     r has the following meaning:         r=0      P = A         r=1      P = B         r<0 P is on the backward extension of AB r>1      P is on the forward extension of AB         0<r<1 P is interior to AB */ double relation(POINT p,LINESEG l) { LINESEG tl; tl.s=l.s; tl.e=p; return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e)); } // 求点C到线段AB所在直线的垂足 P POINT perpendicular(POINT p,LINESEG l) { double r=relation(p,l); 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 ptolinesegdist(POINT p,LINESEG l,POINT &np) { double r=relation(p,l); if(r<0) { np=l.s; return dist(p,l.s); } if(r>1) { np=l.e; return dist(p,l.e); } np=perpendicular(p,l); return dist(p,np); } // 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别  double ptoldist(POINT p,LINESEG l) { return abs(multiply(p,l.e,l.s))/dist(l.s,l.e); } /* 计算点到折线集的最近距离,并返回最近点. 注意：调用的是ptolineseg()函数 */ double ptopointset(int vcount,POINT pointset[],POINT p,POINT &q) { int i; double cd=double(INF),td; LINESEG l; POINT tq,cq; for(i=0;i<vcount-1;i++) { l.s=pointset[i]; l.e=pointset[i+1]; td=ptolinesegdist(p,l,tq); if(td<cd) { cd=td; cq=tq; } } q=cq; return cd; } /* 判断圆是否在多边形内.ptolineseg()函数的应用2 */ bool CircleInsidePolygon(int vcount,POINT center,double radius,POINT polygon[]) { POINT q; double d; q.x=0; q.y=0; d=ptopointset(vcount,polygon,center,q); if(d<radius||fabs(d-radius)<EP) return true; else return false; } /* 返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意：如果想从余弦求夹角的话，注意反余弦函数的定义域是从 0到pi */ double cosine(LINESEG l1,LINESEG l2) { return (((l1.e.x-l1.s.x)*(l2.e.x-l2.s.x) + (l1.e.y-l1.s.y)*(l2.e.y-l2.s.y))/(dist(l1.e,l1.s)*dist(l2.e,l2.s))) ); } // 返回线段l1与l2之间的夹角 单位：弧度 范围(-pi，pi) double lsangle(LINESEG l1,LINESEG l2) { POINT o,s,e; o.x=o.y=0; s.x=l1.e.x-l1.s.x; s.y=l1.e.y-l1.s.y; e.x=l2.e.x-l2.s.x; e.y=l2.e.y-l2.s.y; return angle(o,s,e); } // 如果线段u和v相交(包括相交在端点处)时，返回true bool intersect(LINESEG u,LINESEG v) { return( (max(u.s.x,u.e.x)>=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)*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相交)&&(交点不是双方的端点) 时返回true    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相交时返回true；方法：判断线段u是否跨立线段v  bool intersect_l(LINESEG u,LINESEG v) { return multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0; } // 根据已知两点坐标，求过这两点的直线解析方程： a*x+b*y+c = 0  (a >= 0)  LINE makeline(POINT p1,POINT p2) {    LINE tl;    int sign = 1;    tl.a=p2.y-p1.y;    if(tl.a<0) { 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,水平线返回 0,竖直线返回 1e200