轉自：https://kuangbin.github.io/2019/04/28/20190428/#more

二維幾何

// `計算幾何模板` const double eps = 1e-8; const double pi = acos(-1.0); const int maxp = 1010; //`Compares a double to zero` int sgn(double x){if(fabs(x) < eps)return 0;if(x < 0)return -1;else return 1; } //square of a double inline double sqr(double x){return x*x;} /** Point* Point() - Empty constructor* Point(double _x,double _y) - constructor* input() - double input* output() - %.2f output* operator == - compares x and y* operator < - compares first by x, then by y* operator - - return new Point after subtracting curresponging x and y* operator ^ - cross product of 2d points* operator * - dot product* len() - gives length from origin* len2() - gives square of length from origin* distance(Point p) - gives distance from p* operator + Point b - returns new Point after adding curresponging x and y* operator * double k - returns new Point after multiplieing x and y by k* operator / double k - returns new Point after divideing x and y by k* rad(Point a,Point b)- returns the angle of Point a and Point b from this Point* trunc(double r) - return Point that if truncated the distance from center to r* rotleft() - returns 90 degree ccw rotated point* rotright() - returns 90 degree cw rotated point* rotate(Point p,double angle) - returns Point after rotateing the Point centering at p by angle radian ccw*/ struct Point{double x,y;Point(){}Point(double _x,double _y){x = _x;y = _y;}void input(){scanf("%lf%lf",&x,&y);}void output(){printf("%.2f %.2f\n",x,y);}bool operator == (Point b)const{return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;}bool operator < (Point b)const{return sgn(x-b.x)== 0?sgn(y-b.y)<0:x<b.x;}Point operator -(const Point &b)const{return Point(x-b.x,y-b.y);}//叉積double operator ^(const Point &b)const{return x*b.y - y*b.x;}//點積double operator *(const Point &b)const{return x*b.x + y*b.y;}//返回長度double len(){return hypot(x,y);//庫函數}//返回長度的平方double len2(){return x*x + y*y;}//返回兩點的距離double distance(Point p){return hypot(x-p.x,y-p.y);}Point operator +(const Point &b)const{return Point(x+b.x,y+b.y);}Point operator *(const double &k)const{return Point(x*k,y*k);}Point operator /(const double &k)const{return Point(x/k,y/k);}//`計算pa 和 pb 的夾角`//`就是求這個點看a,b 所成的夾角`//`測試 LightOJ1203`double rad(Point a,Point b){Point p = *this;return fabs(atan2( fabs((a-p)^(b-p)),(a-p)*(b-p) ));}//`化為長度為r的向量`Point trunc(double r){double l = len();if(!sgn(l))return *this;r /= l;return Point(x*r,y*r);}//`逆時針旋轉90度`Point rotleft(){return Point(-y,x);}//`順時針旋轉90度`Point rotright(){return Point(y,-x);}//`繞著p點逆時針旋轉angle`Point rotate(Point p,double angle){Point v = (*this) - p;double c = cos(angle), s = sin(angle);return Point(p.x + v.x*c - v.y*s,p.y + v.x*s + v.y*c);} }; /** Stores two points* Line() - Empty constructor* Line(Point _s,Point _e) - Line through _s and _e* operator == - checks if two points are same* Line(Point p,double angle) - one end p , another end at angle degree* Line(double a,double b,double c) - Line of equation ax + by + c = 0* input() - inputs s and e* adjust() - orders in such a way that s < e* length() - distance of se* angle() - return 0 <= angle < pi* relation(Point p) - 3 if point is on line* 1 if point on the left of line* 2 if point on the right of line* pointonseg(double p) - return true if point on segment* parallel(Line v) - return true if they are parallel* segcrossseg(Line v) - returns 0 if does not intersect* returns 1 if non-standard intersection* returns 2 if intersects* linecrossseg(Line v) - line and seg* linecrossline(Line v) - 0 if parallel* 1 if coincides* 2 if intersects* crosspoint(Line v) - returns intersection point* dispointtoline(Point p) - distance from point p to the line* dispointtoseg(Point p) - distance from p to the segment* dissegtoseg(Line v) - distance of two segment* lineprog(Point p) - returns projected point p on se line* symmetrypoint(Point p) - returns reflection point of p over se**/ struct Line{Point s,e;Line(){}Line(Point _s,Point _e){s = _s;e = _e;}bool operator ==(Line v){return (s == v.s)&&(e == v.e);}//`根據一個點和傾斜角angle確定直線,0<=angle<pi`Line(Point p,double angle){s = p;if(sgn(angle-pi/2) == 0){e = (s + Point(0,1));}else{e = (s + Point(1,tan(angle)));}}//ax+by+c=0Line(double a,double b,double c){if(sgn(a) == 0){s = Point(0,-c/b);e = Point(1,-c/b);}else if(sgn(b) == 0){s = Point(-c/a,0);e = Point(-c/a,1);}else{s = Point(0,-c/b);e = Point(1,(-c-a)/b);}}void input(){s.input();e.input();}void adjust(){if(e < s)swap(s,e);}//求線段長度double length(){return s.distance(e);}//`返回直線傾斜角 0<=angle<pi`double angle(){double k = atan2(e.y-s.y,e.x-s.x);if(sgn(k) < 0)k += pi;if(sgn(k-pi) == 0)k -= pi;return k;}//`點和直線關系`//`1 在左側`//`2 在右側`//`3 在直線上`int relation(Point p){int c = sgn((p-s)^(e-s));if(c < 0)return 1;else if(c > 0)return 2;else return 3;}// 點在線段上的判斷bool pointonseg(Point p){return sgn((p-s)^(e-s)) == 0 && sgn((p-s)*(p-e)) <= 0;}//`兩向量平行(對應直線平行或重合)`bool parallel(Line v){return sgn((e-s)^(v.e-v.s)) == 0;}//`兩線段相交判斷`//`2 規范相交`//`1 非規范相交`//`0 不相交`int segcrossseg(Line v){int d1 = sgn((e-s)^(v.s-s));int d2 = sgn((e-s)^(v.e-s));int d3 = sgn((v.e-v.s)^(s-v.s));int d4 = sgn((v.e-v.s)^(e-v.s));if( (d1^d2)==-2 && (d3^d4)==-2 )return 2;return (d1==0 && sgn((v.s-s)*(v.s-e))<=0) ||(d2==0 && sgn((v.e-s)*(v.e-e))<=0) ||(d3==0 && sgn((s-v.s)*(s-v.e))<=0) ||(d4==0 && sgn((e-v.s)*(e-v.e))<=0);}//`直線和線段相交判斷`//`-*this line -v seg`//`2 規范相交`//`1 非規范相交`//`0 不相交`int linecrossseg(Line v){int d1 = sgn((e-s)^(v.s-s));int d2 = sgn((e-s)^(v.e-s));if((d1^d2)==-2) return 2;return (d1==0||d2==0);}//`兩直線關系`//`0 平行`//`1 重合`//`2 相交`int linecrossline(Line v){if((*this).parallel(v))return v.relation(s)==3;return 2;}//`求兩直線的交點`//`要保證兩直線不平行或重合`Point crosspoint(Line v){double a1 = (v.e-v.s)^(s-v.s);double a2 = (v.e-v.s)^(e-v.s);return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));}//點到直線的距離double dispointtoline(Point p){return fabs((p-s)^(e-s))/length();}//點到線段的距離double dispointtoseg(Point p){if(sgn((p-s)*(e-s))<0 || sgn((p-e)*(s-e))<0)return min(p.distance(s),p.distance(e));return dispointtoline(p);}//`返回線段到線段的距離`//`前提是兩線段不相交，相交距離就是0了`double dissegtoseg(Line v){return min(min(dispointtoseg(v.s),dispointtoseg(v.e)),min(v.dispointtoseg(s),v.dispointtoseg(e)));}//`返回點p在直線上的投影`Point lineprog(Point p){return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );}//`返回點p關于直線的對稱點`Point symmetrypoint(Point p){Point q = lineprog(p);return Point(2*q.x-p.x,2*q.y-p.y);} }; //圓 struct circle{Point p;//圓心double r;//半徑circle(){}circle(Point _p,double _r){p = _p;r = _r;}circle(double x,double y,double _r){p = Point(x,y);r = _r;}//`三角形的外接圓`//`需要Point的+ / rotate() 以及Line的crosspoint()`//`利用兩條邊的中垂線得到圓心`//`測試：UVA12304`circle(Point a,Point b,Point c){Line u = Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));Line v = Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));p = u.crosspoint(v);r = p.distance(a);}//`三角形的內切圓`//`參數bool t沒有作用，只是為了和上面外接圓函數區別`//`測試：UVA12304`circle(Point a,Point b,Point c,bool t){Line u,v;double m = atan2(b.y-a.y,b.x-a.x), n = atan2(c.y-a.y,c.x-a.x);u.s = a;u.e = u.s + Point(cos((n+m)/2),sin((n+m)/2));v.s = b;m = atan2(a.y-b.y,a.x-b.x) , n = atan2(c.y-b.y,c.x-b.x);v.e = v.s + Point(cos((n+m)/2),sin((n+m)/2));p = u.crosspoint(v);r = Line(a,b).dispointtoseg(p);}//輸入void input(){p.input();scanf("%lf",&r);}//輸出void output(){printf("%.2lf %.2lf %.2lf\n",p.x,p.y,r);}bool operator == (circle v){return (p==v.p) && sgn(r-v.r)==0;}bool operator < (circle v)const{return ((p<v.p)||((p==v.p)&&sgn(r-v.r)<0));}//面積double area(){return pi*r*r;}//周長double circumference(){return 2*pi*r;}//`點和圓的關系`//`0 圓外`//`1 圓上`//`2 圓內`int relation(Point b){double dst = b.distance(p);if(sgn(dst-r) < 0)return 2;else if(sgn(dst-r)==0)return 1;return 0;}//`線段和圓的關系`//`比較的是圓心到線段的距離和半徑的關系`int relationseg(Line v){double dst = v.dispointtoseg(p);if(sgn(dst-r) < 0)return 2;else if(sgn(dst-r) == 0)return 1;return 0;}//`直線和圓的關系`//`比較的是圓心到直線的距離和半徑的關系`int relationline(Line v){double dst = v.dispointtoline(p);if(sgn(dst-r) < 0)return 2;else if(sgn(dst-r) == 0)return 1;return 0;}//`兩圓的關系`//`5 相離`//`4 外切`//`3 相交`//`2 內切`//`1 內含`//`需要Point的distance`//`測試：UVA12304`int relationcircle(circle v){double d = p.distance(v.p);if(sgn(d-r-v.r) > 0)return 5;if(sgn(d-r-v.r) == 0)return 4;double l = fabs(r-v.r);if(sgn(d-r-v.r)<0 && sgn(d-l)>0)return 3;if(sgn(d-l)==0)return 2;if(sgn(d-l)<0)return 1;}//`求兩個圓的交點，返回0表示沒有交點，返回1是一個交點，2是兩個交點`//`需要relationcircle`//`測試：UVA12304`int pointcrosscircle(circle v,Point &p1,Point &p2){int rel = relationcircle(v);if(rel == 1 || rel == 5)return 0;double d = p.distance(v.p);double l = (d*d+r*r-v.r*v.r)/(2*d);double h = sqrt(r*r-l*l);Point tmp = p + (v.p-p).trunc(l);p1 = tmp + ((v.p-p).rotleft().trunc(h));p2 = tmp + ((v.p-p).rotright().trunc(h));if(rel == 2 || rel == 4)return 1;return 2;}//`求直線和圓的交點，返回交點個數`int pointcrossline(Line v,Point &p1,Point &p2){if(!(*this).relationline(v))return 0;Point a = v.lineprog(p);double d = v.dispointtoline(p);d = sqrt(r*r-d*d);if(sgn(d) == 0){p1 = a;p2 = a;return 1;}p1 = a + (v.e-v.s).trunc(d);p2 = a - (v.e-v.s).trunc(d);return 2;}//`得到過a,b兩點，半徑為r1的兩個圓`int gercircle(Point a,Point b,double r1,circle &c1,circle &c2){circle x(a,r1),y(b,r1);int t = x.pointcrosscircle(y,c1.p,c2.p);if(!t)return 0;c1.r = c2.r = r1;return t;}//`得到與直線u相切，過點q,半徑為r1的圓`//`測試：UVA12304`int getcircle(Line u,Point q,double r1,circle &c1,circle &c2){double dis = u.dispointtoline(q);if(sgn(dis-r1*2)>0)return 0;if(sgn(dis) == 0){c1.p = q + ((u.e-u.s).rotleft().trunc(r1));c2.p = q + ((u.e-u.s).rotright().trunc(r1));c1.r = c2.r = r1;return 2;}Line u1 = Line((u.s + (u.e-u.s).rotleft().trunc(r1)),(u.e + (u.e-u.s).rotleft().trunc(r1)));Line u2 = Line((u.s + (u.e-u.s).rotright().trunc(r1)),(u.e + (u.e-u.s).rotright().trunc(r1)));circle cc = circle(q,r1);Point p1,p2;if(!cc.pointcrossline(u1,p1,p2))cc.pointcrossline(u2,p1,p2);c1 = circle(p1,r1);if(p1 == p2){c2 = c1;return 1;}c2 = circle(p2,r1);return 2;}//`同時與直線u,v相切，半徑為r1的圓`//`測試：UVA12304`int getcircle(Line u,Line v,double r1,circle &c1,circle &c2,circle &c3,circle &c4){if(u.parallel(v))return 0;//兩直線平行Line u1 = Line(u.s + (u.e-u.s).rotleft().trunc(r1),u.e + (u.e-u.s).rotleft().trunc(r1));Line u2 = Line(u.s + (u.e-u.s).rotright().trunc(r1),u.e + (u.e-u.s).rotright().trunc(r1));Line v1 = Line(v.s + (v.e-v.s).rotleft().trunc(r1),v.e + (v.e-v.s).rotleft().trunc(r1));Line v2 = Line(v.s + (v.e-v.s).rotright().trunc(r1),v.e + (v.e-v.s).rotright().trunc(r1));c1.r = c2.r = c3.r = c4.r = r1;c1.p = u1.crosspoint(v1);c2.p = u1.crosspoint(v2);c3.p = u2.crosspoint(v1);c4.p = u2.crosspoint(v2);return 4;}//`同時與不相交圓cx,cy相切，半徑為r1的圓`//`測試：UVA12304`int getcircle(circle cx,circle cy,double r1,circle &c1,circle &c2){circle x(cx.p,r1+cx.r),y(cy.p,r1+cy.r);int t = x.pointcrosscircle(y,c1.p,c2.p);if(!t)return 0;c1.r = c2.r = r1;return t;}//`過一點作圓的切線(先判斷點和圓的關系)`//`測試：UVA12304`int tangentline(Point q,Line &u,Line &v){int x = relation(q);if(x == 2)return 0;if(x == 1){u = Line(q,q + (q-p).rotleft());v = u;return 1;}double d = p.distance(q);double l = r*r/d;double h = sqrt(r*r-l*l);u = Line(q,p + ((q-p).trunc(l) + (q-p).rotleft().trunc(h)));v = Line(q,p + ((q-p).trunc(l) + (q-p).rotright().trunc(h)));return 2;}//`求兩圓相交的面積`double areacircle(circle v){int rel = relationcircle(v);if(rel >= 4)return 0.0;if(rel <= 2)return min(area(),v.area());double d = p.distance(v.p);double hf = (r+v.r+d)/2.0;double ss = 2*sqrt(hf*(hf-r)*(hf-v.r)*(hf-d));double a1 = acos((r*r+d*d-v.r*v.r)/(2.0*r*d));a1 = a1*r*r;double a2 = acos((v.r*v.r+d*d-r*r)/(2.0*v.r*d));a2 = a2*v.r*v.r;return a1+a2-ss;}//`求兩圓相交的面積(精度更高)(需要long double)`double areacircle2(circle v){double a=hypot(p.x-v.p.x,p.y-v.p.y),b=r,c=v.r;double s1=pi*r*r,s2=pi*v.r*v.r;if(sgn(a-b-c)>=0)return 0;if(sgn(a+min(b,c)-max(b,c))<=0)return min(s1,s2);else{double cta1=2*acos((a*a+b*b-c*c)/(2*a*b));double cta2=2*acos((a*a+c*c-b*b)/(2*a*c));return cta1/(2*pi)*s1-0.5*sin(cta1)*b*b+cta2/(2*pi)*s2-0.5*sin(cta2)*c*c;}}//`求圓和三角形pab的相交面積`//`測試：POJ3675 HDU3982 HDU2892`double areatriangle(Point a,Point b){if(sgn((p-a)^(p-b)) == 0)return 0.0;Point q[5];int len = 0;q[len++] = a;Line l(a,b);Point p1,p2;if(pointcrossline(l,q[1],q[2])==2){if(sgn((a-q[1])*(b-q[1]))<0)q[len++] = q[1];if(sgn((a-q[2])*(b-q[2]))<0)q[len++] = q[2];}q[len++] = b;if(len == 4 && sgn((q[0]-q[1])*(q[2]-q[1]))>0)swap(q[1],q[2]);double res = 0;for(int i = 0;i < len-1;i++){if(relation(q[i])==0||relation(q[i+1])==0){double arg = p.rad(q[i],q[i+1]);res += r*r*arg/2.0;}else{res += fabs((q[i]-p)^(q[i+1]-p))/2.0;}}return res;} };/** n,p Line l for each side* input(int _n) - inputs _n size polygon* add(Point q) - adds a point at end of the list* getline() - populates line array* cmp - comparision in convex_hull order* norm() - sorting in convex_hull order* getconvex(polygon &convex) - returns convex hull in convex* Graham(polygon &convex) - returns convex hull in convex* isconvex() - checks if convex* relationpoint(Point q) - returns 3 if q is a vertex* 2 if on a side* 1 if inside* 0 if outside* convexcut(Line u,polygon &po) - left side of u in po* gercircumference() - returns side length* getarea() - returns area* getdir() - returns 0 for cw, 1 for ccw* getbarycentre() - returns barycenter**/ struct polygon{int n;Point p[maxp];Line l[maxp];void input(int _n){n = _n;for(int i = 0;i < n;i++)p[i].input();}void add(Point q){p[n++] = q;}void getline(){for(int i = 0;i < n;i++){l[i] = Line(p[i],p[(i+1)%n]);}}struct cmp{Point p;cmp(const Point &p0){p = p0;}bool operator()(const Point &aa,const Point &bb){Point a = aa, b = bb;int d = sgn((a-p)^(b-p));if(d == 0){return sgn(a.distance(p)-b.distance(p)) < 0;}return d > 0;}};//`進行極角排序`//`首先需要找到最左下角的點`//`需要重載號好Point的 < 操作符(min函數要用) `void norm(){Point mi = p[0];for(int i = 1;i < n;i++)mi = min(mi,p[i]);sort(p,p+n,cmp(mi));}//`得到凸包`//`得到的凸包里面的點編號是0$\sim$n-1的`//`兩種凸包的方法`//`注意如果有影響，要特判下所有點共點，或者共線的特殊情況`//`測試 LightOJ1203 LightOJ1239`void getconvex(polygon &convex){sort(p,p+n);convex.n = n;for(int i = 0;i < min(n,2);i++){convex.p[i] = p[i];}if(convex.n == 2 && (convex.p[0] == convex.p[1]))convex.n--;//特判if(n <= 2)return;int &top = convex.n;top = 1;for(int i = 2;i < n;i++){while(top && sgn((convex.p[top]-p[i])^(convex.p[top-1]-p[i])) <= 0)top--;convex.p[++top] = p[i];}int temp = top;convex.p[++top] = p[n-2];for(int i = n-3;i >= 0;i--){while(top != temp && sgn((convex.p[top]-p[i])^(convex.p[top-1]-p[i])) <= 0)top--;convex.p[++top] = p[i];}if(convex.n == 2 && (convex.p[0] == convex.p[1]))convex.n--;//特判convex.norm();//`原來得到的是順時針的點，排序后逆時針`}//`得到凸包的另外一種方法`//`測試 LightOJ1203 LightOJ1239`void Graham(polygon &convex){norm();int &top = convex.n;top = 0;if(n == 1){top = 1;convex.p[0] = p[0];return;}if(n == 2){top = 2;convex.p[0] = p[0];convex.p[1] = p[1];if(convex.p[0] == convex.p[1])top--;return;}convex.p[0] = p[0];convex.p[1] = p[1];top = 2;for(int i = 2;i < n;i++){while( top > 1 && sgn((convex.p[top-1]-convex.p[top-2])^(p[i]-convex.p[top-2])) <= 0 )top--;convex.p[top++] = p[i];}if(convex.n == 2 && (convex.p[0] == convex.p[1]))convex.n--;//特判}//`判斷是不是凸的`bool isconvex(){bool s[3];memset(s,false,sizeof(s));for(int i = 0;i < n;i++){int j = (i+1)%n;int k = (j+1)%n;s[sgn((p[j]-p[i])^(p[k]-p[i]))+1] = true;if(s[0] && s[2])return false;}return true;}//`判斷點和任意多邊形的關系`//` 3 點上`//` 2 邊上`//` 1 內部`//` 0 外部`int relationpoint(Point q){for(int i = 0;i < n;i++){if(p[i] == q)return 3;}getline();for(int i = 0;i < n;i++){if(l[i].pointonseg(q))return 2;}int cnt = 0;for(int i = 0;i < n;i++){int j = (i+1)%n;int k = sgn((q-p[j])^(p[i]-p[j]));int u = sgn(p[i].y-q.y);int v = sgn(p[j].y-q.y);if(k > 0 && u < 0 && v >= 0)cnt++;if(k < 0 && v < 0 && u >= 0)cnt--;}return cnt != 0;}//`直線u切割凸多邊形左側`//`注意直線方向`//`測試：HDU3982`void convexcut(Line u,polygon &po){int &top = po.n;//注意引用top = 0;for(int i = 0;i < n;i++){int d1 = sgn((u.e-u.s)^(p[i]-u.s));int d2 = sgn((u.e-u.s)^(p[(i+1)%n]-u.s));if(d1 >= 0)po.p[top++] = p[i];if(d1*d2 < 0)po.p[top++] = u.crosspoint(Line(p[i],p[(i+1)%n]));}}//`得到周長`//`測試 LightOJ1239`double getcircumference(){double sum = 0;for(int i = 0;i < n;i++){sum += p[i].distance(p[(i+1)%n]);}return sum;}//`得到面積`double getarea(){double sum = 0;for(int i = 0;i < n;i++){sum += (p[i]^p[(i+1)%n]);}return fabs(sum)/2;}//`得到方向`//` 1 表示逆時針，0表示順時針`bool getdir(){double sum = 0;for(int i = 0;i < n;i++)sum += (p[i]^p[(i+1)%n]);if(sgn(sum) > 0)return 1;return 0;}//`得到重心`Point getbarycentre(){Point ret(0,0);double area = 0;for(int i = 1;i < n-1;i++){double tmp = (p[i]-p[0])^(p[i+1]-p[0]);if(sgn(tmp) == 0)continue;area += tmp;ret.x += (p[0].x+p[i].x+p[i+1].x)/3*tmp;ret.y += (p[0].y+p[i].y+p[i+1].y)/3*tmp;}if(sgn(area)) ret = ret/area;return ret;}//`多邊形和圓交的面積`//`測試：POJ3675 HDU3982 HDU2892`double areacircle(circle c){double ans = 0;for(int i = 0;i < n;i++){int j = (i+1)%n;if(sgn( (p[j]-c.p)^(p[i]-c.p) ) >= 0)ans += c.areatriangle(p[i],p[j]);else ans -= c.areatriangle(p[i],p[j]);}return fabs(ans);}//`多邊形和圓關系`//` 2 圓完全在多邊形內`//` 1 圓在多邊形里面，碰到了多邊形邊界`//` 0 其它`int relationcircle(circle c){getline();int x = 2;if(relationpoint(c.p) != 1)return 0;//圓心不在內部for(int i = 0;i < n;i++){if(c.relationseg(l[i])==2)return 0;if(c.relationseg(l[i])==1)x = 1;}return x;} }; //`AB X AC` double cross(Point A,Point B,Point C){return (B-A)^(C-A); } //`AB*AC` double dot(Point A,Point B,Point C){return (B-A)*(C-A); } //`最小矩形面積覆蓋` //` A 必須是凸包(而且是逆時針順序)` //` 測試 UVA 10173` double minRectangleCover(polygon A){//`要特判A.n < 3的情況`if(A.n < 3)return 0.0;A.p[A.n] = A.p[0];double ans = -1;int r = 1, p = 1, q;for(int i = 0;i < A.n;i++){//`卡出離邊A.p[i] - A.p[i+1]最遠的點`while( sgn( cross(A.p[i],A.p[i+1],A.p[r+1]) - cross(A.p[i],A.p[i+1],A.p[r]) ) >= 0 )r = (r+1)%A.n;//`卡出A.p[i] - A.p[i+1]方向上正向n最遠的點`while(sgn( dot(A.p[i],A.p[i+1],A.p[p+1]) - dot(A.p[i],A.p[i+1],A.p[p]) ) >= 0 )p = (p+1)%A.n;if(i == 0)q = p;//`卡出A.p[i] - A.p[i+1]方向上負向最遠的點`while(sgn(dot(A.p[i],A.p[i+1],A.p[q+1]) - dot(A.p[i],A.p[i+1],A.p[q])) <= 0)q = (q+1)%A.n;double d = (A.p[i] - A.p[i+1]).len2();double tmp = cross(A.p[i],A.p[i+1],A.p[r]) *(dot(A.p[i],A.p[i+1],A.p[p]) - dot(A.p[i],A.p[i+1],A.p[q]))/d;if(ans < 0 || ans > tmp)ans = tmp;}return ans; }//`直線切凸多邊形` //`多邊形是逆時針的，在q1q2的左側` //`測試:HDU3982` vector<Point> convexCut(const vector<Point> &ps,Point q1,Point q2){vector<Point>qs;int n = ps.size();for(int i = 0;i < n;i++){Point p1 = ps[i], p2 = ps[(i+1)%n];int d1 = sgn((q2-q1)^(p1-q1)), d2 = sgn((q2-q1)^(p2-q1));if(d1 >= 0)qs.push_back(p1);if(d1 * d2 < 0)qs.push_back(Line(p1,p2).crosspoint(Line(q1,q2)));}return qs; } //`半平面交` //`測試 POJ3335 POJ1474 POJ1279` //*************************** struct halfplane:public Line{double angle;halfplane(){}//`表示向量s->e逆時針(左側)的半平面`halfplane(Point _s,Point _e){s = _s;e = _e;}halfplane(Line v){s = v.s;e = v.e;}void calcangle(){angle = atan2(e.y-s.y,e.x-s.x);}bool operator <(const halfplane &b)const{return angle < b.angle;} }; struct halfplanes{int n;halfplane hp[maxp];Point p[maxp];int que[maxp];int st,ed;void push(halfplane tmp){hp[n++] = tmp;}//去重void unique(){int m = 1;for(int i = 1;i < n;i++){if(sgn(hp[i].angle-hp[i-1].angle) != 0)hp[m++] = hp[i];else if(sgn( (hp[m-1].e-hp[m-1].s)^(hp[i].s-hp[m-1].s) ) > 0)hp[m-1] = hp[i];}n = m;}bool halfplaneinsert(){for(int i = 0;i < n;i++)hp[i].calcangle();sort(hp,hp+n);unique();que[st=0] = 0;que[ed=1] = 1;p[1] = hp[0].crosspoint(hp[1]);for(int i = 2;i < n;i++){while(st<ed && sgn((hp[i].e-hp[i].s)^(p[ed]-hp[i].s))<0)ed--;while(st<ed && sgn((hp[i].e-hp[i].s)^(p[st+1]-hp[i].s))<0)st++;que[++ed] = i;if(hp[i].parallel(hp[que[ed-1]]))return false;p[ed]=hp[i].crosspoint(hp[que[ed-1]]);}while(st<ed && sgn((hp[que[st]].e-hp[que[st]].s)^(p[ed]-hp[que[st]].s))<0)ed--;while(st<ed && sgn((hp[que[ed]].e-hp[que[ed]].s)^(p[st+1]-hp[que[ed]].s))<0)st++;if(st+1>=ed)return false;return true;}//`得到最后半平面交得到的凸多邊形`//`需要先調用halfplaneinsert() 且返回true`void getconvex(polygon &con){p[st] = hp[que[st]].crosspoint(hp[que[ed]]);con.n = ed-st+1;for(int j = st,i = 0;j <= ed;i++,j++)con.p[i] = p[j];} }; //***************************const int maxn = 1010; struct circles{circle c[maxn];double ans[maxn];//`ans[i]表示被覆蓋了i次的面積`double pre[maxn];int n;circles(){}void add(circle cc){c[n++] = cc;}//`x包含在y中`bool inner(circle x,circle y){if(x.relationcircle(y) != 1)return 0;return sgn(x.r-y.r)<=0?1:0;}//圓的面積并去掉內含的圓void init_or(){bool mark[maxn] = {0};int i,j,k=0;for(i = 0;i < n;i++){for(j = 0;j < n;j++)if(i != j && !mark[j]){if( (c[i]==c[j])||inner(c[i],c[j]) )break;}if(j < n)mark[i] = 1;}for(i = 0;i < n;i++)if(!mark[i])c[k++] = c[i];n = k;}//`圓的面積交去掉內含的圓`void init_add(){int i,j,k;bool mark[maxn] = {0};for(i = 0;i < n;i++){for(j = 0;j < n;j++)if(i != j && !mark[j]){if( (c[i]==c[j])||inner(c[j],c[i]) )break;}if(j < n)mark[i] = 1;}for(i = 0;i < n;i++)if(!mark[i])c[k++] = c[i];n = k;}//`半徑為r的圓，弧度為th對應的弓形的面積`double areaarc(double th,double r){return 0.5*r*r*(th-sin(th));}//`測試SPOJVCIRCLES SPOJCIRUT`//`SPOJVCIRCLES求n個圓并的面積，需要加上init\_or()去掉重復圓（否則WA）`//`SPOJCIRUT 是求被覆蓋k次的面積，不能加init\_or()`//`對于求覆蓋多少次面積的問題，不能解決相同圓，而且不能init\_or()`//`求多圓面積并，需要init\_or,其中一個目的就是去掉相同圓`void getarea(){memset(ans,0,sizeof(ans));vector<pair<double,int> >v;for(int i = 0;i < n;i++){v.clear();v.push_back(make_pair(-pi,1));v.push_back(make_pair(pi,-1));for(int j = 0;j < n;j++)if(i != j){Point q = (c[j].p - c[i].p);double ab = q.len(),ac = c[i].r, bc = c[j].r;if(sgn(ab+ac-bc)<=0){v.push_back(make_pair(-pi,1));v.push_back(make_pair(pi,-1));continue;}if(sgn(ab+bc-ac)<=0)continue;if(sgn(ab-ac-bc)>0)continue;double th = atan2(q.y,q.x), fai = acos((ac*ac+ab*ab-bc*bc)/(2.0*ac*ab));double a0 = th-fai;if(sgn(a0+pi)<0)a0+=2*pi;double a1 = th+fai;if(sgn(a1-pi)>0)a1-=2*pi;if(sgn(a0-a1)>0){v.push_back(make_pair(a0,1));v.push_back(make_pair(pi,-1));v.push_back(make_pair(-pi,1));v.push_back(make_pair(a1,-1));}else{v.push_back(make_pair(a0,1));v.push_back(make_pair(a1,-1));}}sort(v.begin(),v.end());int cur = 0;for(int j = 0;j < v.size();j++){if(cur && sgn(v[j].first-pre[cur])){ans[cur] += areaarc(v[j].first-pre[cur],c[i].r);ans[cur] += 0.5*(Point(c[i].p.x+c[i].r*cos(pre[cur]),c[i].p.y+c[i].r*sin(pre[cur]))^Point(c[i].p.x+c[i].r*cos(v[j].first),c[i].p.y+c[i].r*sin(v[j].first)));}cur += v[j].second;pre[cur] = v[j].first;}}for(int i = 1;i < n;i++)ans[i] -= ans[i+1];} };

三維幾何

const double eps = 1e-8; int sgn(double x){if(fabs(x) < eps)return 0;if(x < 0)return -1;else return 1; } struct Point3{double x,y,z;Point3(double _x = 0,double _y = 0,double _z = 0){x = _x;y = _y;z = _z;}void input(){scanf("%lf%lf%lf",&x,&y,&z);}void output(){printf("%.2lf %.2lf %.2lf\n",x,y,z);}bool operator ==(const Point3 &b)const{return sgn(x-b.x) == 0 && sgn(y-b.y) == 0 && sgn(z-b.z) == 0;}bool operator <(const Point3 &b)const{return sgn(x-b.x)==0?(sgn(y-b.y)==0?sgn(z-b.z)<0:y<b.y):x<b.x;}double len(){return sqrt(x*x+y*y+z*z);}double len2(){return x*x+y*y+z*z;}double distance(const Point3 &b)const{return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y)+(z-b.z)*(z-b.z));}Point3 operator -(const Point3 &b)const{return Point3(x-b.x,y-b.y,z-b.z);}Point3 operator +(const Point3 &b)const{return Point3(x+b.x,y+b.y,z+b.z);}Point3 operator *(const double &k)const{return Point3(x*k,y*k,z*k);}Point3 operator /(const double &k)const{return Point3(x/k,y/k,z/k);}//點乘double operator *(const Point3 &b)const{return x*b.x+y*b.y+z*b.z;}//叉乘Point3 operator ^(const Point3 &b)const{return Point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}double rad(Point3 a,Point3 b){Point3 p = (*this);return acos( ( (a-p)*(b-p) )/ (a.distance(p)*b.distance(p)) );}//變換長度Point3 trunc(double r){double l = len();if(!sgn(l))return *this;r /= l;return Point3(x*r,y*r,z*r);} }; struct Line3 {Point3 s,e;Line3(){}Line3(Point3 _s,Point3 _e){s = _s;e = _e;}bool operator ==(const Line3 v){return (s==v.s)&&(e==v.e);}void input(){s.input();e.input();}double length(){return s.distance(e);}//點到直線距離double dispointtoline(Point3 p){return ((e-s)^(p-s)).len()/s.distance(e);}//點到線段距離double dispointtoseg(Point3 p){if(sgn((p-s)*(e-s)) < 0 || sgn((p-e)*(s-e)) < 0)return min(p.distance(s),e.distance(p));return dispointtoline(p);}//`返回點p在直線上的投影`Point3 lineprog(Point3 p){return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );}//`p繞此向量逆時針arg角度`Point3 rotate(Point3 p,double ang){if(sgn(((s-p)^(e-p)).len()) == 0)return p;Point3 f1 = (e-s)^(p-s);Point3 f2 = (e-s)^(f1);double len = ((s-p)^(e-p)).len()/s.distance(e);f1 = f1.trunc(len); f2 = f2.trunc(len);Point3 h = p+f2;Point3 pp = h+f1;return h + ((p-h)*cos(ang)) + ((pp-h)*sin(ang));}//`點在直線上`bool pointonseg(Point3 p){return sgn( ((s-p)^(e-p)).len() ) == 0 && sgn((s-p)*(e-p)) == 0;} }; struct Plane {Point3 a,b,c,o;//`平面上的三個點，以及法向量`Plane(){}Plane(Point3 _a,Point3 _b,Point3 _c){a = _a;b = _b;c = _c;o = pvec();}Point3 pvec(){return (b-a)^(c-a);}//`ax+by+cz+d = 0`Plane(double _a,double _b,double _c,double _d){o = Point3(_a,_b,_c);if(sgn(_a) != 0)a = Point3((-_d-_c-_b)/_a,1,1);else if(sgn(_b) != 0)a = Point3(1,(-_d-_c-_a)/_b,1);else if(sgn(_c) != 0)a = Point3(1,1,(-_d-_a-_b)/_c);}//`點在平面上的判斷`bool pointonplane(Point3 p){return sgn((p-a)*o) == 0;}//`兩平面夾角`double angleplane(Plane f){return acos(o*f.o)/(o.len()*f.o.len());}//`平面和直線的交點，返回值是交點個數`int crossline(Line3 u,Point3 &p){double x = o*(u.e-a);double y = o*(u.s-a);double d = x-y;if(sgn(d) == 0)return 0;p = ((u.s*x)-(u.e*y))/d;return 1;}//`點到平面最近點(也就是投影)`Point3 pointtoplane(Point3 p){Line3 u = Line3(p,p+o);crossline(u,p);return p;}//`平面和平面的交線`int crossplane(Plane f,Line3 &u){Point3 oo = o^f.o;Point3 v = o^oo;double d = fabs(f.o*v);if(sgn(d) == 0)return 0;Point3 q = a + (v*(f.o*(f.a-a))/d);u = Line3(q,q+oo);return 1;} };

?平面最近點對

const int MAXN = 100010; const double eps = 1e-8; const double INF = 1e20; struct Point{double x,y;void input(){scanf("%lf%lf",&x,&y);} }; double dist(Point a,Point b){return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)); } Point p[MAXN]; Point tmpt[MAXN]; bool cmpx(Point a,Point b){return a.x < b.x || (a.x == b.x && a.y < b.y); } bool cmpy(Point a,Point b){return a.y < b.y || (a.y == b.y && a.x < b.x); } double Closest_Pair(int left,int right){double d = INF;if(left == right)return d;if(left+1 == right)return dist(p[left],p[right]);int mid = (left+right)/2;double d1 = Closest_Pair(left,mid);double d2 = Closest_Pair(mid+1,right);d = min(d1,d2);int cnt = 0;for(int i = left;i <= right;i++){if(fabs(p[mid].x - p[i].x) <= d)tmpt[cnt++] = p[i];}sort(tmpt,tmpt+cnt,cmpy);for(int i = 0;i < cnt;i++){for(int j = i+1;j < cnt && tmpt[j].y - tmpt[i].y < d;j++)d = min(d,dist(tmpt[i],tmpt[j]));}return d; } int main(){int n;while(scanf("%d",&n) == 1 && n){for(int i = 0;i < n;i++)p[i].input();sort(p,p+n,cmpx);printf("%.2lf\n",Closest_Pair(0,n-1));}return 0; }

三維凸包：HDU4273

const double eps = 1e-8; const int MAXN = 550; int sgn(double x){if(fabs(x) < eps)return 0;if(x < 0)return -1;else return 1; } struct Point3{double x,y,z;Point3(double _x = 0, double _y = 0, double _z = 0){x = _x;y = _y;z = _z;}void input(){scanf("%lf%lf%lf",&x,&y,&z);}bool operator ==(const Point3 &b)const{return sgn(x-b.x) == 0 && sgn(y-b.y) == 0 && sgn(z-b.z) == 0;}double len(){return sqrt(x*x+y*y+z*z);}double len2(){return x*x+y*y+z*z;}double distance(const Point3 &b)const{return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y)+(z-b.z)*(z-b.z));}Point3 operator -(const Point3 &b)const{return Point3(x-b.x,y-b.y,z-b.z);}Point3 operator +(const Point3 &b)const{return Point3(x+b.x,y+b.y,z+b.z);}Point3 operator *(const double &k)const{return Point3(x*k,y*k,z*k);}Point3 operator /(const double &k)const{return Point3(x/k,y/k,z/k);}//點乘double operator *(const Point3 &b)const{return x*b.x + y*b.y + z*b.z;}//叉乘Point3 operator ^(const Point3 &b)const{return Point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);} }; struct CH3D{struct face{//表示凸包一個面上的三個點的編號int a,b,c;//表示該面是否屬于最終的凸包上的面bool ok;};//初始頂點數int n;Point3 P[MAXN];//凸包表面的三角形數int num;//凸包表面的三角形face F[8*MAXN];int g[MAXN][MAXN];//叉乘Point3 cross(const Point3 &a,const Point3 &b,const Point3 &c){return (b-a)^(c-a);}//`三角形面積*2`double area(Point3 a,Point3 b,Point3 c){return ((b-a)^(c-a)).len();}//`四面體有向面積*6`double volume(Point3 a,Point3 b,Point3 c,Point3 d){return ((b-a)^(c-a))*(d-a);}//`正：點在面同向`double dblcmp(Point3 &p,face &f){Point3 p1 = P[f.b] - P[f.a];Point3 p2 = P[f.c] - P[f.a];Point3 p3 = p - P[f.a];return (p1^p2)*p3;}void deal(int p,int a,int b){int f = g[a][b];face add;if(F[f].ok){if(dblcmp(P[p],F[f]) > eps)dfs(p,f);else {add.a = b;add.b = a;add.c = p;add.ok = true;g[p][b] = g[a][p] = g[b][a] = num;F[num++] = add;}}}//遞歸搜索所有應該從凸包內刪除的面void dfs(int p,int now){F[now].ok = false;deal(p,F[now].b,F[now].a);deal(p,F[now].c,F[now].b);deal(p,F[now].a,F[now].c);}bool same(int s,int t){Point3 &a = P[F[s].a];Point3 &b = P[F[s].b];Point3 &c = P[F[s].c];return fabs(volume(a,b,c,P[F[t].a])) < eps &&fabs(volume(a,b,c,P[F[t].b])) < eps &&fabs(volume(a,b,c,P[F[t].c])) < eps;}//構建三維凸包void create(){num = 0;face add;//***********************************//此段是為了保證前四個點不共面bool flag = true;for(int i = 1;i < n;i++){if(!(P[0] == P[i])){swap(P[1],P[i]);flag = false;break;}}if(flag)return;flag = true;for(int i = 2;i < n;i++){if( ((P[1]-P[0])^(P[i]-P[0])).len() > eps ){swap(P[2],P[i]);flag = false;break;}}if(flag)return;flag = true;for(int i = 3;i < n;i++){if(fabs( ((P[1]-P[0])^(P[2]-P[0]))*(P[i]-P[0]) ) > eps){swap(P[3],P[i]);flag = false;break;}}if(flag)return;//**********************************for(int i = 0;i < 4;i++){add.a = (i+1)%4;add.b = (i+2)%4;add.c = (i+3)%4;add.ok = true;if(dblcmp(P[i],add) > 0)swap(add.b,add.c);g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = num;F[num++] = add;}for(int i = 4;i < n;i++)for(int j = 0;j < num;j++)if(F[j].ok && dblcmp(P[i],F[j]) > eps){dfs(i,j);break;}int tmp = num;num = 0;for(int i = 0;i < tmp;i++)if(F[i].ok)F[num++] = F[i];}//表面積//`測試：HDU3528`double area(){double res = 0;if(n == 3){Point3 p = cross(P[0],P[1],P[2]);return p.len()/2;}for(int i = 0;i < num;i++)res += area(P[F[i].a],P[F[i].b],P[F[i].c]);return res/2.0;}double volume(){double res = 0;Point3 tmp = Point3(0,0,0);for(int i = 0;i < num;i++)res += volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);return fabs(res/6);}//表面三角形個數int triangle(){return num;}//表面多邊形個數//`測試：HDU3662`int polygon(){int res = 0;for(int i = 0;i < num;i++){bool flag = true;for(int j = 0;j < i;j++)if(same(i,j)){flag = 0;break;}res += flag;}return res;}//重心//`測試：HDU4273`Point3 barycenter(){Point3 ans = Point3(0,0,0);Point3 o = Point3(0,0,0);double all = 0;for(int i = 0;i < num;i++){double vol = volume(o,P[F[i].a],P[F[i].b],P[F[i].c]);ans = ans + (((o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0)*vol);all += vol;}ans = ans/all;return ans;}//點到面的距離//`測試：HDU4273`double ptoface(Point3 p,int i){double tmp1 = fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p));double tmp2 = ((P[F[i].b]-P[F[i].a])^(P[F[i].c]-P[F[i].a])).len();return tmp1/tmp2;} }; CH3D hull; int main() {while(scanf("%d",&hull.n) == 1){for(int i = 0;i < hull.n;i++)hull.P[i].input();hull.create();Point3 p = hull.barycenter();double ans = 1e20;for(int i = 0;i < hull.num;i++)ans = min(ans,hull.ptoface(p,i));printf("%.3lf\n",ans);}return 0; }

叉乘性質：

總結

以上是生活随笔為你收集整理的(转)计算几何模板 - kuangbin的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。