二维几何相关（未完全测试）
基础+凸包+旋转卡壳+半平面交
1 /************************************************************************* 2 > File Name: board.cpp 3 > Author: yijiull 4 > Mail: 1147161372@qq.com 5 > Created Time: 2017年09月22日 星期五 08时32分54秒 6 ************************************************************************/ 7 #include <iostream> 8 #include <cstring> 9 #include <cstdio> 10 #include <bits/stdc++.h> 11 using namespace std; 12 const double eps = 1e-12; 13 const int inf = 0x3f3f3f3f; 14 struct Point { 15 double x,y; 16 Point (double x = 0, double y = 0) : x(x), y(y) {} 17 }; 18 typedef Point Vector; 19 Vector operator + (Vector a, Vector b) { 20 return Vector (a.x + b.x, a.y + b.y); 21 } 22 Vector operator * (Vector a, double s) { 23 return Vector (a.x * s, a.y * s); 24 } 25 Vector operator / (Vector a, double p) { 26 return Vector (a.x / p, a.y / p); 27 } 28 Vector operator - (Point a, Point b) { 29 return Vector (a.x - b.x, a.y - b.y); 30 } 31 bool operator < (Point a, Point b) { 32 return a.x < b.x || (a.x == b.x && a.y < b.y); 33 } 34 int dcmp (double x) { 35 if(fabs(x) < eps) return 0; 36 return x < 0 ? -1 : 1; 37 } 38 bool operator == (const Point &a, const Point &b) { 39 return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0; 40 } 41 double Angel (Vector a) { 42 return atan2(a.y, a.x); 43 } 44 double Dot(Vector a, Vector b) { 45 return a.x * b.x + a.y * b.y; 46 } 47 double Length (Vector a) { 48 return sqrt(Dot(a, a)); 49 } 50 double Angle (Vector a, Vector b) { 51 return acos(Dot(a, b) / Length(a) / Length(b)); 52 } 53 double Cross (Vector a, Vector b) { 54 return a.x * b.y - a.y * b.x; 55 } 56 double Area2 (Point a, Point b, Point c) { 57 return Cross(b - a, c - a); 58 } 59 Vector Rotate (Vector a, double rad) { //右手逆时针旋转rad(弧度) 60 return Vector (a.x * cos(rad) - a.y * sin(rad), a.x * sin(rad) + a.y * cos(rad)); 61 } 62 Vector Normal (Vector a) { 63 double L = Length(a); 64 return Vector(-a.y / L, a.x / L); 65 } 66 //两直线交点 67 Point GetLineIntersection (Point p, Vector v, Point q, Vector w) { 68 Vector u = p - q; 69 double t1 = Cross(w, u) / Cross(v, w); 70 double t2 = Cross(v, u) / Cross(v, w); 71 return p + v * t1; // return q + w * t2; 72 } 73 //点到直线的dis 74 double DistanceToLine(Point p, Point a, Point b) { 75 Vector v1 = b - a, v2 = p - a; 76 return fabs(Cross(v1, v2)) / Length(v1); 77 } 78 //点到线段的dis 79 double DistanceToSegment(Point p, Point a, Point b) { 80 if(a == b) return Length(a - p); 81 Vector v1 = b - a, v2 = p - a, v3 = p - b; 82 if(dcmp(Dot(v1, v2)) < 0) return Length(v2); 83 else if(dcmp(Dot(v1, v3)) > 0) return Length(v3); 84 else return fabs(Cross(v1, v2)) / Length(v1); 85 } 86 //点在直线的投影 87 Point GetLineProjection(Point p, Point a, Point b) { 88 Vector v = b - a; 89 return a + v * (Dot(v, p-a) / Dot(v, v)); 90 } 91 //判断两线段是否规范相交 92 bool SegmentProperIntersection(Point a1, Point b1, Point a2, Point b2) { 93 double c1 = Cross(b1 - a1, a2 - a1), c2 = Cross(b1 - a1, b2 - a1), 94 c3 = Cross(b2 - a2, a1 - a2), c4 = Cross(b2 - a2, b1 - a2); 95 return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4) < 0; 96 } 97 //判断点是否在线段上(即∠APB等于π) 98 bool OnSegment(Point p, Point a, Point b) { 99 return dcmp(Cross(a - p, b - p)) == 0 && dcmp(Dot(a - p, b - p)) < 0; 100 } 101 //多边形面积 102 double ConvexPolygonArea(Point *p, int n) { 103 double area = 0; 104 for(int i = 1; i < n-1; i++) { 105 area += Cross(p[i] - p[0], p[i+1] - p[0]); 106 } 107 return area / 2; 108 } 109 //判断点是否在多边形内or外or上 110 bool isPointInPolygon(Point p,Point *poly, int n) { 111 int cnt = 0; 112 poly[n] = poly[0]; 113 for(int i = 0; i < n; i++) { 114 if(OnSegment(p, poly[i], poly[i+1])) return 0; //在边界上 115 int k = dcmp(Cross(poly[i+1]-poly[i], p - poly[i])); 116 int d1 = dcmp(poly[i].y - p.y); 117 int d2 = dcmp(poly[i+1].y - p.y); 118 if(k > 0 && d1 <= 0 && d2 > 0) cnt++; 119 if(k < 0 && d2 <= 0 && d1 >= 0) cnt--; 120 } 121 if(cnt) return -1; //内部 122 return 1; //外部 123 } 124 //凸包 125 int ConvexHull(Point *p, int n, Point *ch) { 126 int m = 0; 127 sort(p, p+n); 128 for(int i = 0; i < n; i++) { 129 while(m > 1 && Cross(ch[m-1] - ch[m-2], p[i] - ch[m-2]) <= 0) m--; //去掉在边上的点 130 ch[m++] = p[i]; 131 } 132 int k = m; 133 for(int i = n-2; i >= 0; i--) { 134 while(m > k && Cross(ch[m-1] - ch[m-2], p[i] - ch[m-2]) <= 0) m--; // 135 ch[m++] = p[i]; 136 } 137 if(m > 1) m--; 138 return m; 139 } 140 //判断是否是稳定凸包 141 bool IsStableConvexHull(Point *ch, int n) { 142 if(n == 0) return 0; 143 ch[n] = ch[0]; 144 ch[n+1] = ch[1]; 145 for(int i = 0; i < n; i++) { 146 if(Cross(ch[i] - ch[(i-1+n)%n], ch[i+1] - ch[i]) != 0 && Cross(ch[i+1] - ch[i] , ch[i+2] - ch[i]) != 0) 147 return 0; 148 } 149 return 1; 150 } 151 //旋转卡壳求凸包直径 152 double RotateCalipers(Point *ch, int n) { 153 double ans = - inf; 154 ch[n] = ch[0]; 155 int q = 1; 156 for(int i = 0; i < n; i++) { 157 while(Cross(ch[i+1] - ch[i], ch[q] - ch[i]) < Cross(ch[i+1] - ch[i], ch[q+1] -ch[i])) q = (q+1)%n; 158 ans = max(ans, max(Length(ch[i] - ch[q]), Length(ch[i+1] -ch[q+1]))); 159 } 160 return ans; 161 } 162 163 //半平面交 164 struct Line{ 165 Point p; 166 Vector v; 167 double rad; 168 Line () {} 169 Line (Point p, Vector v) : p(p), v(v) { 170 rad = atan2(v.y,v.x); 171 } 172 bool operator < (const Line &L) const { 173 return rad < L.rad; 174 } 175 }; 176 bool OnLeft(Line L, Point p) { 177 return Cross(L.v, p - L.p) > 0; 178 } 179 Point GetLineIntersection (Line a, Line b) { 180 Vector u = a.p - b.p; 181 double t = Cross(b.v, u) / Cross(a.v, b.v); 182 return a.p + a.v*t; 183 } 184 185 int HalfplaneIntersection(Line *L, int n,Point *poly) { 186 sort(L, L+n); 187 int first,last; 188 Point *p = new Point[n]; 189 Line *q = new Line[n]; //双端队列 190 q[first = last = 0] = L[0]; 191 for(int i = 1; i < n; i++) { 192 while(first < last && !OnLeft(L[i], p[last-1])) last--; //去尾 193 while(first < last && !OnLeft(L[i], p[first])) first++; 194 q[++last] = L[i]; 195 if(dcmp(Cross(q[last].v, q[last-1].v)) == 0) { 196 last--; 197 if(OnLeft(q[last], L[i].p)) q[last] = L[i]; 198 } 199 if(first < last) p[last-1] = GetLineIntersection (q[last-1],q[last]); 200 } 201 while(first < last && !OnLeft(q[first], p[last-1])) last--; //删除无用平面 202 if(last - first <= 1) return 0; //空集 203 p[last] = GetLineIntersection (q[last], q[first]); 204 205 int m = 0; 206 for(int i = first; i <= last; i++) poly[m++] = p[i]; 207 return m; 208 } 209 210 int main(){ 211 212 }
1 //线段AB与圆的交点 2 void LineIntersectionCircle(Point A, Point B, Circle cr, Point *p, int &num){ 3 double x0 = cr.o.x, y0 = cr.o.y; 4 double x1 = A.x, y1 = A.y, x2 = B.x, y2 = B.y; 5 double dx = x2 - x1, dy = y2 - y1; 6 double a = dx * dx + dy * dy; 7 double b = 2 * dx *(x1 - x0) + 2 * dy * (y1 -y0); 8 double c = (x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0) - cr.r *cr.r; 9 double delta = b*b - 4*a*c; 10 num = 0; 11 if(dcmp(delta) >= 0){ 12 double t1 = (-b - sqrt(delta)) / (2*a); 13 double t2 = (-b + sqrt(delta)) / (2*a); 14 if(dcmp(t1-1) <= 0 && dcmp(t1) >= 0) { 15 p[num++] = Point(x1 + t1*dx, y1 + t1*dy); 16 } 17 if(dcmp(t2-1) <= 0 && dcmp(t2) >= 0) { 18 p[num++] = Point(x1 + t2*dx, y1 + t2*dy); 19 } 20 } 21 } 22 //扇形面积 23 double SectorArea(Point a, Point b, Circle cr){ 24 double theta = Angle(a - cr.o) - Angle(b - cr.o); 25 while(theta <= 0) theta += 2*pi; 26 while(theta > 2*pi) theta -= 2*pi; 27 theta = min(theta, 2*pi - theta); 28 return cr.r * cr.r * theta / 2; 29 } 30 31 double cal(Point a, Point b, Circle cr){ 32 Point tp[3]; 33 int num = 0; 34 Vector ta = a - cr.o; 35 Vector tb = b - cr.o; 36 bool ina = (Length(ta) - cr.r) < 0; 37 bool inb = (Length(tb) - cr.r) < 0; 38 if(ina){ 39 if(inb){ 40 return fabs(Cross(ta, tb))/2; 41 } else{ 42 LineIntersectionCircle(a, b, cr, tp, num); 43 return SectorArea(b, tp[0], cr) + fabs(Cross(ta, tp[0] - cr.o))/2; 44 } 45 } else{ 46 if(inb){ 47 LineIntersectionCircle(a, b, cr, tp, num); 48 return SectorArea(a, tp[0], cr) + fabs(Cross(tb, tp[0] - cr.o))/2; 49 } else { 50 LineIntersectionCircle(a, b, cr, tp, num); 51 if(num == 2) { 52 return SectorArea(a, tp[0], cr) + SectorArea(b, tp[1], cr) + fabs(Cross(tp[0]-cr.o, tp[1]-cr.o))/2; 53 } else { 54 return SectorArea(a, b, cr); 55 } 56 } 57 } 58 } 59 //圆与多边形的面积并 60 double CirclePolyArea(Point *p, int n, Circle cr){ 61 double res = 0; 62 p[n] = p[0]; 63 for(int i = 0; i < n; i++){ 64 int sgn = dcmp(Cross(p[i] - cr.o, p[i+1] - cr.o)); 65 if(sgn){ 66 res += sgn * cal(p[i], p[i+1], cr); 67 } 68 } 69 return res; 70 }
自适应辛普森积分
1 double F(double x) { 2 //Simpson公式用到的函数 3 } 4 double simpson(double a, double b) { //三点Simpson法，这里要求F是一个全局函数 5 double c = a + (b - a) / 2; 6 return (F(a) + 4 * F(c) + F(b))*(b - a) / 6; 7 } 8 double asr(double a, double b, double eps, double A) { //自适应Simpson公式（递归过程）。已知整个区间[a,b]上的三点Simpson值A 9 double c = a + (b - a) / 2; 10 double L = simpson(a, c), R = simpson(c, b); 11 if (fabs(L + R - A) <= 15 * eps)return L + R + (L + R - A) / 15.0; 12 return asr(a, c, eps / 2, L) + asr(c, b, eps / 2, R); 13 } 14 double asr(double a, double b, double eps) { //自适应Simpson公式（主过程） 15 return asr(a, b, eps, simpson(a, b)); 16 }
多边形重心
1 /* 2 * 求多边形重心 3 * INIT: pnt[]已按顺时针(或逆时针)排好序; | CALL: res = bcenter(pnt, n); 4 */ 5 Point bcenter(Point pnt[], int n) { 6 Point p, s; 7 double tp, area = 0, tpx = 0, tpy = 0; 8 p.x = pnt[0].x; 9 p.y = pnt[0].y; 10 for (int i = 1; i <= n; ++i) { 11 // Point:0 ~ n - 1 12 s.x = pnt[(i == n) ? 0 : i].x; 13 s.y = pnt[(i == n) ? 0 : i].y; 14 tp = (p.x * s.y - s.x * p.y); 15 area += tp / 2; 16 tpx += (p.x + s.x) * tp; 17 tpy += (p.y + s.y) * tp; 18 p.x = s.x; 19 p.y = s.y; 20 } 21 s.x = tpx / (6 * area); 22 s.y = tpy / (6 * area); 23 return s; 24 }
相关公式（转自f_zyj）
1 已知圆锥表面积S求最大体积V 2 3 V = S * sqrt(S / (72 * Pi)) 4 5 三角形重点 6 7 设三角形的三条边为a, b, c, 且不妨假设a <= b <= c. 8 面积 9 10 三角形面积可以根据海伦公式求得： 11 12 s = sqrt(p * (p - a) * (p - b) * (p - c)); 13 p = (a + b + c) / 2; 14 15 关键点与A, B, C三顶点距离之和 16 费马点 17 18 该点到三角形三个顶点的距离之和最小。 19 有个有趣的结论: 20 若三角形的三个内角均小于120度,那么该点连接三个顶点形成的三个角均为120度;若三角形存在一个内角大于120度,则该顶点就是费马点。 21 计算公式如下: 22 若有一个内角大于120度(这里假设为角C),则距离为a + b;若三个内角均小于120度,则距离为sqrt((a * a + b * b + c * c + 4 * sqrt(3.0) * s) / 2)。 23 内心 24 25 角平分线的交点。 26 令x = (a + b - c) / 2, y = (a - b + c) / 2, z = (-a + b + c) / 2, h = s / p. 27 计算公式为sqrt(x * x + h * h) + sqrt(y * y + h * h) + sqrt(z * z + h * h)。 28 重心 29 30 中线的交点。 31 计算公式如下: 32 2.0 / 3 * (sqrt((2 * (a * a + b * b) - c * c) / 4) 33 + sqrt((2 * (a * a + c * c) - b * b) / 4) + sqrt((2 * (b * b + c * c) - a * a) / 4))。 34 垂心 35 36 垂线的交点。 37 计算公式如下: 38 3 * (c / 2 / sqrt(1 - cosC * cosC))。 39 外心 40 41 三点求圆心坐标。 42 43 Point waixin(Point a, Point b, Point c) 44 { 45 double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1 * a1 + b1 * b1) / 2; 46 double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2 * a2 + b2 * b2) / 2; 47 double d = a1 * b2 - a2 * b1; 48 return Point(a.x + (c1 * b2 - c2 * b1) / d, a.y + (a1 * c2 -a2 * c1) / d); 49 }
圆的面积交
1 double area[maxn]; 2 #define sqr(x) (x)*(x) 3 int dcmp(double x) { 4 if (x < -eps) return -1; else return x > eps; 5 } 6 struct cp { 7 double x, y, r, angle; 8 int d; 9 cp(){} 10 cp(double xx, double yy, double ang = 0, int t = 0) { 11 x = xx; y = yy; angle = ang; d = t; 12 } 13 void get() { 14 scanf("%lf%lf%lf", &x, &y, &r); 15 d = 1; 16 } 17 }cir[maxn], tp[maxn * 2]; 18 double dis(cp a, cp b) { 19 return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y)); 20 } 21 double cross(cp p0, cp p1, cp p2) { 22 return (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x); 23 } 24 int CirCrossCir(cp p1, double r1, cp p2, double r2, cp &cp1, cp &cp2) { 25 double mx = p2.x - p1.x, sx = p2.x + p1.x, mx2 = mx * mx; 26 double my = p2.y - p1.y, sy = p2.y + p1.y, my2 = my * my; 27 double sq = mx2 + my2, d = -(sq - sqr(r1 - r2)) * (sq - sqr(r1 + r2)); 28 if (d + eps < 0) return 0; if (d < eps) d = 0; else d = sqrt(d); 29 double x = mx * ((r1 + r2) * (r1 - r2) + mx * sx) + sx * my2; 30 double y = my * ((r1 + r2) * (r1 - r2) + my * sy) + sy * mx2; 31 double dx = mx * d, dy = my * d; sq *= 2; 32 cp1.x = (x - dy) / sq; cp1.y = (y + dx) / sq; 33 cp2.x = (x + dy) / sq; cp2.y = (y - dx) / sq; 34 if (d > eps) return 2; else return 1; 35 } 36 bool circmp(const cp& u, const cp& v) { 37 return dcmp(u.r - v.r) < 0; 38 } 39 bool cmp(const cp& u, const cp& v) { 40 if (dcmp(u.angle - v.angle)) return u.angle < v.angle; 41 return u.d > v.d; 42 } 43 double calc(cp cir, cp cp1, cp cp2) { 44 double ans = (cp2.angle - cp1.angle) * sqr(cir.r) 45 - cross(cir, cp1, cp2) + cross(cp(0, 0), cp1, cp2); 46 return ans / 2; 47 } 48 void CirUnion(cp cir[], int n) { 49 cp cp1, cp2; 50 sort(cir, cir + n, circmp); 51 for (int i = 0; i < n; ++i) 52 for (int j = i + 1; j < n; ++j) 53 if (dcmp(dis(cir[i], cir[j]) + cir[i].r - cir[j].r) <= 0) 54 cir[i].d++; 55 for (int i = 0; i < n; ++i) { 56 int tn = 0, cnt = 0; 57 for (int j = 0; j < n; ++j) { 58 if (i == j) continue; 59 if (CirCrossCir(cir[i], cir[i].r, cir[j], cir[j].r, 60 cp2, cp1) < 2) continue; 61 cp1.angle = atan2(cp1.y - cir[i].y, cp1.x - cir[i].x); 62 cp2.angle = atan2(cp2.y - cir[i].y, cp2.x - cir[i].x); 63 cp1.d = 1; tp[tn++] = cp1; 64 cp2.d = -1; tp[tn++] = cp2; 65 if (dcmp(cp1.angle - cp2.angle) > 0) cnt++; 66 } 67 tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, pi, -cnt); 68 tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, -pi, cnt); 69 sort(tp, tp + tn, cmp); 70 int p, s = cir[i].d + tp[0].d; 71 for (int j = 1; j < tn; ++j) { 72 p = s; s += tp[j].d; 73 area[p] += calc(cir[i], tp[j - 1], tp[j]); 74 } 75 } 76 }
三维几何相关
1 /************************************************************************* 2 > File Name: board.cpp 3 > Author: yijiull 4 > Mail: 1147161372@qq.com 5 > Created Time: 2017年09月22日 星期五 16时51分54秒 6 ************************************************************************/ 7 #include <iostream> 8 #include <cstring> 9 #include <cstdio> 10 #include <bits/stdc++.h> 11 using namespace std; 12 const double eps = 1e-12; 13 const int inf = 0x3f3f3f3f; 14 struct Point3{ 15 double x, y, z; 16 Point3(double x=0, double y=0, double z=0) : x(x), y(y), z(z) {} 17 }; 18 int dcmp(double x) { 19 if(fabs(x) < eps) return 0; 20 return x < 0 ? -1 : 1; 21 } 22 typedef Point3 Vector3; 23 Vector3 operator + (Vector3 a, Vector3 b) { 24 return Vector3(a.x + b.x, a.y + b.y, a.z + b.z); 25 } 26 Vector3 operator - (Point3 a, Point3 b) { 27 return Vector3(a.x - b.x, a.y - b.y, a.z - b.z); 28 } 29 Vector3 operator * (Vector3 a, double p) { 30 return Vector3(a.x * 3, a.y * 3, a.z * 3); 31 } 32 Vector3 operator / (Vector3 a, double p) { 33 return Vector3(a.x / p, a.y / p, a.z / p); 34 } 35 double Dot(Vector3 a, Vector3 b) { 36 return a.x * b.x + a.y * b.y + a.z * b.z; 37 } 38 bool operator == (Point3 a, Point3 b) { 39 return a.x==b.x && a.y==b.y && a.z==b.z; 40 } 41 double Length(Vector3 a) { 42 return sqrt(Dot(a,a)); 43 } 44 double Angle(Vector3 a, Vector3 b) { 45 return acos(Dot(a,b) / Length(a) / Length(b)) ; 46 } 47 Vector3 Normal(Vector3 a) { 48 return a / Length(a); 49 } 50 //点p到平面p0-n的距离,n是单位向量 51 double DistanceToPlane(const Point3 & p, const Point3 &p0, const Vector3 &n) { 52 return fabs(Dot(p-p0, n)); //如果不取绝对值，得到的是有向距离 53 } 54 55 //点P在平面p0-n上的投影，ｎ必须为单位向量 56 Point3 GetPlaneProjection(const Point3 p, const Point3 &p0, const Vector3 &n) { 57 return p - n * Dot(p - p0, n); 58 } 59 //直线p1-p2到平面p0-n的交点，假定交点唯一存在 60 Point3 LinePlaneIntersection(Point3 p1, Point3 p2, Point3 p0, Vector3 n) { 61 Vector3 v = p2 - p1; 62 double t = (Dot(n, p0 -p1) / Dot(n, p2 -p1)); //判断分母是否为０ 63 return p1 + v * t; //如果是线段，判断ｔ是不是在０和１之间 64 } 65 66 Vector3 Cross(Vector3 a, Vector3 b) { 67 return Vector3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x); 68 } 69 double Area2(Point3 a, Point3 b, Point3 c) { 70 return Length(Cross(b-a, c-a)); 71 } 72 //点p在三角形p0,p1,p2中 73 bool PointInTri(Point3 p, Point3 p0, Point3 p1, Point3 p2) { 74 double area1 = Area2(p, p0, p1); 75 double area2 = Area2(p, p1, p2); 76 double area3 = Area2(p, p2, p0); 77 return dcmp(area1 + area2 + area3 - Area2(p0, p1, p2)) == 0; 78 } 79 //三角形p0,p1,p2是否和线段ab相交 80 bool TriSegIntersection(Point3 p0, Point3 p1, Point3 p2, Point3 a, Point3 b, Point3 &p) { 81 Vector3 n = Cross(p1 - p0, p2 - p0); 82 if(dcmp(Dot(n, b - a)) == 0) return false; //无交点 83 double t = Dot(n, p0 - a) / Dot(n, b - a); 84 if(dcmp(t) < 0 || dcmp(t-1) > 0) return false; //交点不在直线上 85 p = a + (b - a) * t; 86 return PointInTri(p, p0, p1, p2); //是否在三角形内 87 } 88 //点p到直线ab的距离 89 double DistanceToLine(Point3 p, Point3 a, Point3 b) { 90 Vector3 v1 = b - a, v2 = p - a; 91 return Length(Cross(v1, v2)) / Length(v1); 92 } 93 //点p到线段ab的距离 94 double DisstanceToSegment(Point3 p, Point3 a, Point3 b) { 95 if(a == b) return Length(p - a); 96 Vector3 v1 = b - a, v2 = p - a, v3 = p - b; 97 if(dcmp(Dot(v1, v2)) < 0) return Length(v2); 98 else if(dcmp(Dot(v1, v3)) > 0) return Length(v3); 99 else return Length(Cross(v1, v2)) / Length(v1); 100 } 101 //四面体abcd的体积 102 double Volume6(Point3 a, Point3 b, Point3 c, Point3 d) { 103 return Dot(d-a, Cross(b-a, c-a)); 104 } 105 106 //三维凸包 107 int vis[200][200]; 108 struct Face{ 109 int v[3]; 110 Vector3 normal(Point3 *p) const { 111 return Cross(p[v[1]] - p[v[0]], p[v[2]] - p[v[0]]); 112 } 113 int cansee(Point3 *p, int i) const { 114 return Dot(p[i] - p[v[0]], normal(p)) > 0 ? 1 : 0; 115 } 116 }; 117 //倍增法求三维凸包 118 //没有考虑特殊情况(如四点共面),实践中,请在调用之前对输入点进行微小扰动 119 vector<Face> CH3D(Point3 *p, int n) { 120 vector<Face> cur; 121 memset(vis, 0, sizeof(vis)); 122 //由于已经扰动,前三个点不共线 123 cur.push_back((Face){ {0, 1, 2} } ); 124 cur.push_back((Face){ {2, 1, 0} } ); 125 for(int i = 3; i < n; i++) { 126 vector<Face> next; 127 for(int j = 0; j < cur.size(); j++) { 128 Face &f = cur[j]; 129 int res = f.cansee(p, i); 130 if(!res) next.push_back(f); 131 for(int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res; 132 } 133 for(int j = 0; j < cur.size(); j++) { 134 for(int k = 0; k < 3; k++) { 135 int a = cur[j].v[k], b = cur[j].v[(k+1)%3]; 136 if(vis[a][b] != vis[b][a] && vis[a][b]) //(a, b) 是分界线, 左边对(a, b) 可见 137 next.push_back((Face) { {a, b, i} } ); 138 } 139 } 140 cur = next; 141 } 142 return cur; 143 } 144 //扰动!!! 145 double rand01() { 146 return rand() / (double)RAND_MAX; 147 } 148 double randeps() { 149 return (rand01() - 0.5) * eps; 150 } 151 Point3 add_noise(Point3 p) { 152 return Point3(p.x + randeps(), p.y + randeps(), p.z + randeps()); 153 } 154 155 int main(){ 156 157 }