1、凸包

1 inline bool cmp(const POINT &a, const POINT &b) { 2 if(a.y == b.y) return a.x < b.x; 3 return a.y < b.y; 4 } 5 //turn left 6 inline bool Cross(POINT &sp, POINT &ep, POINT &op) { 7 return (sp.x - op.x) * (ep.y - op.y) - (ep.x - op.x) * (sp.y - op.y) >= 0; 8 } 9 10 inline double dist(POINT &a, POINT &b) { 11 return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); 12 } 13 14 void Graham_scan() { 15 std::sort(p, p + n, cmp); 16 top = 1; 17 stk[0] = 0; stk[1] = 1; 18 for(int i = 2; i < n; ++i) { 19 while(top && Cross(p[i], p[stk[top]], p[stk[top - 1]])) --top; 20 stk[++top] = i; 21 } 22 int len = top; 23 stk[++top] = n - 2; 24 for(int i = n - 3; i >= 0; --i) { 25 while(top != len && Cross(p[i], p[stk[top]], p[stk[top - 1]])) --top; 26 stk[++top] = i; 27 } 28 } View Code

2、旋轉(zhuǎn)卡殼（POJ 3384）

1 #include <cstdio> 2 #include <cstring> 3 #include <cmath> 4 #include <algorithm> 5 using namespace std; 6 7 #define EPS 1e-8 8 #define MAXN 1000 9 10 inline int sgn(double x) { 11 if(fabs(x) < EPS) return 0; 12 return x > 0 ? 1 : -1; 13 } 14 15 struct Point { 16 double x, y; 17 Point(double xx = 0, double yy = 0): x(xx), y(yy) {} 18 bool operator == (const Point &b) const { 19 return sgn(x - b.x) == 0 && sgn(y - b.y) == 0; 20 } 21 }; 22 //cross 23 inline double operator ^ (const Point &a, const Point &b) { 24 return a.x * b.y - a.y * b.x; 25 } 26 27 inline Point operator - (const Point &a, const Point &b) { 28 return Point(a.x - b.x, a.y - b.y); 29 } 30 31 struct Line { 32 Point s, e; 33 double ag; 34 }; 35 36 struct polygon { 37 Point v[MAXN]; 38 int n; 39 } pg, res; 40 41 inline double dist(Point &a, Point &b) { 42 return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); 43 } 44 45 inline double Cross(Point o, Point s, Point e) { 46 return (s - o) ^ (e - o); 47 } 48 //cross_point 49 Point operator * (const Line &a, const Line &b) { 50 Point res; 51 double u = Cross(a.s, a.e, b.s), v = Cross(a.e, a.s, b.e); 52 res.x = (b.s.x * v + b.e.x * u)/(u + v); 53 res.y = (b.s.y * v + b.e.y * u)/(u + v); 54 return res; 55 } 56 57 int parallel(Line a, Line b) { 58 double u = (a.e.x - a.s.x) * (b.e.y - b.s.y) - (a.e.y - a.s.y) * (b.e.x - b.s.x); 59 return sgn(u) == 0; 60 } 61 62 inline void set_vector(double x1, double y1, double x2, double y2, Line &v) { 63 v.s.x = x1; v.s.y = y1; 64 v.e.x = x2; v.e.y = y2; 65 v.ag = atan2(y2 - y1, x2 - x1); 66 } 67 68 Line vct[MAXN], deq[MAXN]; 69 70 bool cmp(const Line &a, const Line &b) { 71 if(sgn(a.ag - b.ag) == 0) 72 return sgn(Cross(b.s, b.e, a.s)) < 0; 73 return a.ag < b.ag; 74 } 75 76 int half_planes_cross(Line *v, int vn) { 77 int i, n; 78 //sort(v, v + vn, cmp); 79 for(i = n = 1; i < vn; ++i) { 80 if(sgn(v[i].ag - v[i-1].ag) == 0) continue; 81 v[n++] = v[i]; 82 } 83 int head = 0, tail = 1; 84 deq[0] = v[0], deq[1] = v[1]; 85 for(i = 2; i < n; ++i) { 86 if(parallel(deq[tail - 1], deq[tail]) || parallel(deq[head], deq[head + 1])) 87 return false; 88 while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[tail - 1] * deq[tail])) > 0) 89 --tail; 90 while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[head] * deq[head + 1])) > 0) 91 ++head; 92 deq[++tail] = v[i]; 93 } 94 while(head < tail && sgn(Cross(deq[head].s, deq[head].e, deq[tail - 1] * deq[tail])) > 0) 95 --tail; 96 while(head < tail && sgn(Cross(deq[tail].s, deq[tail].e, deq[head] * deq[head + 1])) > 0) 97 ++head; 98 if(tail <= head + 1) return false; 99 res.n = 0; 100 for(i = head; i < tail; ++i) 101 res.v[res.n++] = deq[i] * deq[i + 1]; 102 res.v[res.n++] = deq[head] * deq[tail]; 103 res.n = unique(res.v, res.v + res.n) - res.v; 104 res.v[res.n] = res.v[0]; 105 return true; 106 } 107 108 void moving(Line v[], int vn, double r) { 109 for(int i = 0; i < vn; ++i) { 110 double dx = v[i].e.x - v[i].s.x, dy = v[i].e.y - v[i].s.y; 111 dx = dx / dist(v[i].s, v[i].e) * r; 112 dy = dy / dist(v[i].s, v[i].e) * r; 113 v[i].s.x += dy; v[i].e.x += dy; 114 v[i].s.y -= dx; v[i].e.y -= dx; 115 } 116 } 117 118 int ix, jx; 119 120 double dia_roataing_calipers() { 121 double dia = 0; 122 ix = jx = 0; 123 int q = 1; 124 for(int i = 0; i < res.n; ++i) { 125 while(sgn(Cross(res.v[i+1], res.v[i], res.v[q+1]) - Cross(res.v[i+1], res.v[i], res.v[q])) > 0) 126 q = (q + 1) % res.n; 127 if(sgn(dist(res.v[i], res.v[q]) - dia) > 0) { 128 dia = dist(res.v[i], res.v[q]); 129 ix = i; jx = q; 130 } 131 if(sgn(dist(res.v[i+1], res.v[q]) - dia) > 0) { 132 dia = dist(res.v[i+1], res.v[q]); 133 ix = i+1; jx = q; 134 } 135 } 136 return dia; 137 } 138 139 int main() { 140 int n; 141 double r; 142 while(scanf("%d%lf", &n, &r) != EOF) { 143 for(int i = 0; i < n; ++i) scanf("%lf%lf", &pg.v[i].x, &pg.v[i].y); 144 pg.v[n] = pg.v[0]; 145 for(int i = 0; i < n; ++i) 146 set_vector(pg.v[i].x, pg.v[i].y, pg.v[i+1].x, pg.v[i+1].y, vct[i]); 147 moving(vct, n, r); 148 half_planes_cross(vct, n); 149 dia_roataing_calipers(); 150 printf("%.4f %.4f %.4f %.4f\n", res.v[ix].x, res.v[ix].y, res.v[jx].x, res.v[jx].y); 151 } 152 } View Code

3、最小圓覆蓋

1 #include <cstdio> 2 #include <cmath> 3 4 const int MAXN = 1010; 5 6 struct Point { 7 double x, y; 8 Point(double xx = 0, double yy = 0): x(xx), y(yy) {} 9 }; 10 11 double operator ^ (const Point &a, const Point &b) { 12 return a.x * b.y - a.y * b.x; 13 } 14 15 Point operator - (const Point &a, const Point &b) { 16 return Point(a.x - b.x, a.y - b.y); 17 } 18 19 double dist(const Point &a, const Point &b) { 20 double x = a.x - b.x, y = a.y - b.y; 21 return sqrt(x * x + y * y); 22 } 23 24 struct Circle { 25 double r; 26 Point centre; 27 }; 28 29 struct Triangle { 30 Point t[3]; 31 double Area() { 32 return fabs((t[1] - t[0]) ^ (t[2] - t[0]))/2; 33 } 34 }; 35 36 Circle circumcircleOfTriangle(Triangle t) { 37 Circle tmp; 38 double a, b, c, c1, c2; 39 Point A(t.t[0].x, t.t[0].y); 40 Point B(t.t[1].x, t.t[1].y); 41 Point C(t.t[2].x, t.t[2].y); 42 a = dist(B, C); 43 b = dist(A, C); 44 c = dist(A, B); 45 tmp.r = a * b * c / t.Area() / 4; 46 c1 = (A.x * A.x + A.y * A.y - B.x * B.x - B.y * B.y) / 2; 47 c2 = (A.x * A.x + A.y * A.y - C.x * C.x - C.y * C.y) / 2; 48 tmp.centre.x = (c1 * (A.y - C.y) - c2 * (A.y - B.y)) / ((A.x - B.x) * (A.y - C.y) - (A.x - C.x) * (A.y - B.y)); 49 tmp.centre.y = (c1 * (A.x - C.x) - c2 * (A.x - B.x)) / ((A.y - B.y) * (A.x - C.x) - (A.y - C.y) * (A.x - B.x)); 50 return tmp; 51 } 52 53 Circle c; 54 Point a[MAXN]; 55 56 Circle MinCircle2(int tce, Triangle ce) { 57 Circle tmp; 58 if(tce == 0) tmp.r = -2; 59 else if(tce == 1) { 60 tmp.centre = ce.t[0]; 61 tmp.r = 0; 62 } 63 else if(tce == 2) { 64 tmp.r = dist(ce.t[0], ce.t[1]) / 2; 65 tmp.centre.x = (ce.t[0].x + ce.t[1].x) / 2; 66 tmp.centre.y = (ce.t[0].y + ce.t[1].y) / 2; 67 } 68 else if(tce == 3) tmp = circumcircleOfTriangle(ce); 69 return tmp; 70 } 71 72 void MinCircle(int t, int tce, Triangle ce) { 73 Point tmp; 74 c = MinCircle2(tce, ce); 75 if(tce == 3) return; 76 for(int i = 1; i <= t; ++i) { 77 if(dist(a[i], c.centre) > c.r) { 78 ce.t[tce] = a[i]; 79 MinCircle(i - 1, tce + 1, ce); 80 tmp = a[i]; 81 for(int j = i; j >= 2; --j) a[j] = a[j - 1]; 82 a[1] = tmp; 83 } 84 } 85 } 86 87 void run(int n) { 88 Triangle ce; 89 MinCircle(n, 0, ce); 90 printf("%.2f %.2f %.2f\n", c.centre.x, c.centre.y, c.r); 91 } 92 93 int main() { 94 int n; 95 while(scanf("%d", &n) != EOF) { 96 if(n == 0) break; 97 for(int i = 1; i <= n; ++i) 98 scanf("%lf%lf", &a[i].x, &a[i].y); 99 run(n); 100 } 101 } View Code

4、半平面交+最遠(yuǎn)點(diǎn)對(duì)

1 #include <cstdio> 2 #include <cstring> 3 #include <cmath> 4 #include <algorithm> 5 using namespace std; 6 7 #define EPS 1e-8 8 #define MAXN 1000 9 10 inline int sgn(double x) { 11 if(fabs(x) < EPS) return 0; 12 return x > 0 ? 1 : -1; 13 } 14 15 struct Point { 16 double x, y; 17 Point(double xx = 0, double yy = 0): x(xx), y(yy) {} 18 bool operator == (const Point &b) const { 19 return sgn(x - b.x) == 0 && sgn(y - b.y) == 0; 20 } 21 }; 22 //cross 23 inline double operator ^ (const Point &a, const Point &b) { 24 return a.x * b.y - a.y * b.x; 25 } 26 27 inline Point operator - (const Point &a, const Point &b) { 28 return Point(a.x - b.x, a.y - b.y); 29 } 30 31 struct Line { 32 Point s, e; 33 double ag; 34 }; 35 36 struct polygon { 37 Point v[MAXN]; 38 int n; 39 } pg, res; 40 41 inline double dist(Point &a, Point &b) { 42 return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); 43 } 44 45 inline double Cross(Point o, Point s, Point e) { 46 return (s - o) ^ (e - o); 47 } 48 //cross_point 49 Point operator * (const Line &a, const Line &b) { 50 Point res; 51 double u = Cross(a.s, a.e, b.s), v = Cross(a.e, a.s, b.e); 52 res.x = (b.s.x * v + b.e.x * u)/(u + v); 53 res.y = (b.s.y * v + b.e.y * u)/(u + v); 54 return res; 55 } 56 57 int parallel(Line a, Line b) { 58 double u = (a.e.x - a.s.x) * (b.e.y - b.s.y) - (a.e.y - a.s.y) * (b.e.x - b.s.x); 59 return sgn(u) == 0; 60 } 61 62 inline void set_vector(double x1, double y1, double x2, double y2, Line &v) { 63 v.s.x = x1; v.s.y = y1; 64 v.e.x = x2; v.e.y = y2; 65 v.ag = atan2(y2 - y1, x2 - x1); 66 } 67 68 Line vct[MAXN], deq[MAXN]; 69 70 bool cmp(const Line &a, const Line &b) { 71 if(sgn(a.ag - b.ag) == 0) 72 return sgn(Cross(b.s, b.e, a.s)) < 0; 73 return a.ag < b.ag; 74 } 75 76 int half_planes_cross(Line *v, int vn) { 77 int i, n; 78 //sort(v, v + vn, cmp); 79 for(i = n = 1; i < vn; ++i) { 80 if(sgn(v[i].ag - v[i-1].ag) == 0) continue; 81 v[n++] = v[i]; 82 } 83 int head = 0, tail = 1; 84 deq[0] = v[0], deq[1] = v[1]; 85 for(i = 2; i < n; ++i) { 86 if(parallel(deq[tail - 1], deq[tail]) || parallel(deq[head], deq[head + 1])) 87 return false; 88 while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[tail - 1] * deq[tail])) > 0) 89 --tail; 90 while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[head] * deq[head + 1])) > 0) 91 ++head; 92 deq[++tail] = v[i]; 93 } 94 while(head < tail && sgn(Cross(deq[head].s, deq[head].e, deq[tail - 1] * deq[tail])) > 0) 95 --tail; 96 while(head < tail && sgn(Cross(deq[tail].s, deq[tail].e, deq[head] * deq[head + 1])) > 0) 97 ++head; 98 if(tail <= head + 1) return false; 99 res.n = 0; 100 for(i = head; i < tail; ++i) 101 res.v[res.n++] = deq[i] * deq[i + 1]; 102 res.v[res.n++] = deq[head] * deq[tail]; 103 res.n = unique(res.v, res.v + res.n) - res.v; 104 res.v[res.n] = res.v[0]; 105 return true; 106 } 107 108 void moving(Line v[], int vn, double r) { 109 for(int i = 0; i < vn; ++i) { 110 double dx = v[i].e.x - v[i].s.x, dy = v[i].e.y - v[i].s.y; 111 dx = dx / dist(v[i].s, v[i].e) * r; 112 dy = dy / dist(v[i].s, v[i].e) * r; 113 v[i].s.x += dy; v[i].e.x += dy; 114 v[i].s.y -= dx; v[i].e.y -= dx; 115 } 116 } 117 118 int ix, jx; 119 120 double dia_roataing_calipers() { 121 double dia = 0; 122 ix = jx = 0; 123 int q = 1; 124 for(int i = 0; i < res.n; ++i) { 125 while(sgn(Cross(res.v[i+1], res.v[i], res.v[q+1]) - Cross(res.v[i+1], res.v[i], res.v[q])) > 0) 126 q = (q + 1) % res.n; 127 if(sgn(dist(res.v[i], res.v[q]) - dia) > 0) { 128 dia = dist(res.v[i], res.v[q]); 129 ix = i; jx = q; 130 } 131 if(sgn(dist(res.v[i+1], res.v[q]) - dia) > 0) { 132 dia = dist(res.v[i+1], res.v[q]); 133 ix = i+1; jx = q; 134 } 135 } 136 return dia; 137 } 138 139 int main() { 140 int n; 141 double r; 142 while(scanf("%d%lf", &n, &r) != EOF) { 143 for(int i = 0; i < n; ++i) scanf("%lf%lf", &pg.v[i].x, &pg.v[i].y); 144 pg.v[n] = pg.v[0]; 145 for(int i = 0; i < n; ++i) 146 set_vector(pg.v[i].x, pg.v[i].y, pg.v[i+1].x, pg.v[i+1].y, vct[i]); 147 moving(vct, n, r); 148 half_planes_cross(vct, n); 149 dia_roataing_calipers(); 150 printf("%.4f %.4f %.4f %.4f\n", res.v[ix].x, res.v[ix].y, res.v[jx].x, res.v[jx].y); 151 } 152 } View Code

?

其他1：

#include <cmath> #include <cstdio> #include <cstring> #include <algorithm> using namespace std;typedef double TYPE; #define Abs(x) (((x)>0)?(x):(-(x))) #define Sgn(x) (((x)<0)?(-1):(1)) #define Max(a,b) (((a)>(b))?(a):(b)) #define Min(a,b) (((a)<(b))?(a):(b)) #define EPS 1e-8 #define Infinity 1e+10 #define PI acos(-1.0)//3.14159265358979323846 TYPE Deg2Rad(TYPE deg){return (deg * PI / 180.0);} TYPE Rad2Deg(TYPE rad){return (rad * 180.0 / PI);} TYPE Sin(TYPE deg){return sin(Deg2Rad(deg));} TYPE Cos(TYPE deg){return cos(Deg2Rad(deg));} TYPE ArcSin(TYPE val){return Rad2Deg(asin(val));} TYPE ArcCos(TYPE val){return Rad2Deg(acos(val));} TYPE Sqrt(TYPE val){return sqrt(val);}//點(diǎn) struct POINT {TYPE x, y;POINT(TYPE xx = 0, TYPE yy = 0): x(xx), y(yy) {} }; // 兩個(gè)點(diǎn)的距離 TYPE Distance(const POINT &a, const POINT &b) {return Sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); } //線段 struct SEG {POINT a ,b;SEG() {};SEG(POINT aa, POINT bb): a(aa), b(bb) {} }; //直線(兩點(diǎn)式) struct LINE {POINT a, b;LINE() {};LINE(POINT aa, POINT bb): a(aa), b(bb) {} }; //直線(一般式) struct LINE2 {TYPE A,B,C;LINE2() {};LINE2(TYPE AA, TYPE BB, TYPE CC): A(AA), B(BB), C(CC) {} }; //兩點(diǎn)式化一般式 LINE2 Line2line(const LINE &L) {LINE2 L2;L2.A = L.b.y - L.a.y;L2.B = L.a.x - L.b.x;L2.C = L.b.x * L.a.y - L.a.x * L.b.y;return L2; } // 引用返回直線 Ax + By + C =0 的系數(shù) void Coefficient(const LINE &L, TYPE &A, TYPE &B, TYPE &C) {A = L.b.y - L.a.y;B = L.a.x - L.b.x;C = L.b.x * L.a.y - L.a.x * L.b.y; } void Coefficient(const POINT &p,const TYPE &a,TYPE &A,TYPE &B,TYPE &C) {A = Cos(a);B = Sin(a);C = - (p.y * B + p.x * A); } //直線相交的交點(diǎn) POINT Intersection(const LINE &A, const LINE &B) {TYPE A1, B1, C1;TYPE A2, B2, C2;Coefficient(A, A1, B1, C1);Coefficient(B, A2, B2, C2);POINT I;I.x = - (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1);I.y = (A2 * C1 - A1 * C2) / (A1 * B2 - A2 * B1);return I; } //點(diǎn)到直線的距離 TYPE Ptol(const POINT &p,const LINE2 &L) {return fabs(L.A * p.x + L.B * p.y + L.C)/Sqrt(L.A * L.A + L.B * L.B); } //兩線段叉乘 TYPE Cross2(const SEG &p, const SEG &q) {return (p.b.x - p.a.x) * (q.b.y - q.a.y) - (q.b.x - q.a.x) * (p.b.y - p.a.y); } //有公共端點(diǎn)O叉乘,并判拐,若正p0->p1->p2拐向左 TYPE Cross(const POINT &a, const POINT &b, const POINT &o) {return (a.x - o.x) * (b.y - o.y) - (b.x - o.x) * (a.y - o.y); } //判等(值，點(diǎn)，直線) bool IsEqual(const TYPE &a, const TYPE &b) {return (Abs(a - b) < EPS); } bool IsEqual(const POINT & a, const POINT & b) {return IsEqual(a.x, b.x) && IsEqual(a.y, b.y); } bool IsEqual(const LINE & A, const LINE & B) {TYPE A1, B1, C1;TYPE A2, B2, C2;Coefficient(A, A1, B1, C1);Coefficient(B, A2, B2, C2);return IsEqual(A1 * B2, A2 * B1) && IsEqual(A1 * C2, A2 * C1) && IsEqual(B1 * C2, B2 * C1); } // 判斷點(diǎn)是否在線段上 bool IsOnSeg(const SEG &seg, const POINT &p) {return (IsEqual(p, seg.a) || IsEqual(p, seg.b)) ||(((p.x - seg.a.x) * (p.x - seg.b.x) < 0 ||(p.y - seg.a.y) * (p.y - seg.b.y) < 0) &&(IsEqual(Cross(seg.b, p, seg.a), 0))); }//判斷兩條線斷是否相交，端點(diǎn)重合算相交 bool IsIntersect(const SEG &u, const SEG &v) {return (Cross(v.a, u.b, u.a) * Cross(u.b, v.b, u.a) >= 0) &&(Cross(u.a, v.b, v.a) * Cross(v.b, u.b, v.a) >= 0) &&(Max(u.a.x, u.b.x) >= Min(v.a.x, v.b.x)) &&(Max(v.a.x, v.b.x) >= Min(u.a.x, u.b.x)) &&(Max(u.a.y, u.b.y) >= Min(v.a.y, v.b.y)) &&(Max(v.a.y, v.b.y) >= Min(u.a.y, u.b.y)); } //判斷線段P和直線Q是否相交，端點(diǎn)重合算相交 bool Isonline(const SEG &p, const LINE &q) {SEG p1,p2,q0;p1.a = q.a; p1.b = p.a;p2.a = q.a; p2.b = p.b;q0.a = q.a; q0.b = q.b;return (Cross2(p1,q0) * Cross2(q0,p2) >= 0); } //判斷兩條線斷是否平行 bool IsParallel(const LINE &A, const LINE &B) {TYPE A1, B1, C1;TYPE A2, B2, C2;Coefficient(A, A1, B1, C1);Coefficient(B, A2, B2, C2);return (A1 * B2 == A2 * B1); } //判斷兩條直線是否相交 bool IsIntersect(const LINE & A, const LINE & B) {return !IsParallel(A, B); } // 矩形 struct RECT {POINT a; // 左下點(diǎn)POINT b; // 右上點(diǎn)RECT() {}RECT(const POINT &aa, const POINT &bb): a(aa), b(bb) {} }; //矩形化標(biāo)準(zhǔn) RECT Stdrect(const RECT &q) {RECT p = q;if(p.a.x > p.b.x) swap(p.a.x , p.b.x);if(p.a.y > p.b.y) swap(p.a.y , p.b.y);return p; } //根據(jù)下標(biāo)返回矩形的邊 SEG Edge(const RECT &rect, int idx) {SEG edge;while (idx < 0) idx += 4;switch (idx % 4) {case 0: //下邊edge.a = rect.a;edge.b = POINT(rect.b.x, rect.a.y);break;case 1: //右邊edge.a = POINT(rect.b.x, rect.a.y);edge.b = rect.b;break;case 2: //上邊edge.a = rect.b;edge.b = POINT(rect.a.x, rect.b.y);break;case 3: //左邊edge.a = POINT(rect.a.x, rect.b.y);edge.b = rect.a;break;default:break;}return edge; } //矩形的面積 TYPE Area(const RECT &rect) {return (rect.b.x - rect.a.x) * (rect.b.y - rect.a.y); } //兩個(gè)矩形的公共面積 TYPE CommonArea(const RECT &A, const RECT &B) {TYPE area = 0.0;POINT LL(Max(A.a.x, B.a.x), Max(A.a.y, B.a.y));POINT UR(Min(A.b.x, B.b.x), Min(A.b.y, B.b.y));if((LL.x <= UR.x) && (LL.y <= UR.y)) {area = Area(RECT(LL, UR));}return area; } // 圓 struct CIRCLE {TYPE x, y, r;CIRCLE() {}CIRCLE(TYPE xx, TYPE yy, TYPE zz): x(xx), y(yy), r(zz) {} }; //圓心 POINT Center(const CIRCLE &circle) {return POINT(circle.x, circle.y); } //圓的面積 TYPE Area(const CIRCLE &circle) {return PI * circle.r * circle.r; } //兩個(gè)圓的公共面積 TYPE CommonArea(const CIRCLE &A, const CIRCLE &B) {TYPE area = 0.0;const CIRCLE & M = (A.r > B.r) ? A : B;const CIRCLE & N = (A.r > B.r) ? B : A;TYPE D = Distance(Center(M), Center(N));if((D < M.r + N.r) && (D > M.r - N.r)) {TYPE cosM = (M.r * M.r + D * D - N.r * N.r) / (2.0 * M.r * D);TYPE cosN = (N.r * N.r + D * D - M.r * M.r) / (2.0 * N.r * D);TYPE alpha = 2.0 * ArcCos(cosM);TYPE beta = 2.0 * ArcCos(cosN);TYPE TM = 0.5 * M.r * M.r * Sin(alpha);TYPE TN = 0.5 * N.r * N.r * Sin(beta);TYPE FM = (alpha / 360.0) * Area(M);TYPE FN = (beta / 360.0) * Area(N);area = FM + FN - TM - TN;}else if(D <= M.r - N.r) {area = Area(N);}return area; } //判斷圓是否在矩形內(nèi)(不允許相切) bool IsInCircle(const CIRCLE &circle, const RECT &rect) {return (circle.x - circle.r > rect.a.x) &&(circle.x + circle.r < rect.b.x) &&(circle.y - circle.r > rect.a.y) &&(circle.y + circle.r < rect.b.y); } //判斷矩形是否在圓內(nèi)(不允許相切) bool IsInRect(const CIRCLE & circle, const RECT & rect) {POINT c,d;c.x = rect.a.x; c.y = rect.b.y;d.x = rect.b.x; d.y = rect.a.y;return (Distance( Center(circle) , rect.a ) < circle.r) &&(Distance( Center(circle) , rect.b ) < circle.r) &&(Distance( Center(circle) , c ) < circle.r) &&(Distance( Center(circle) , d ) < circle.r); } //判斷矩形是否與圓相離(不允許相切) bool Isoutside(const CIRCLE & circle, const RECT & rect) {POINT c,d;c.x = rect.a.x; c.y=rect.b.y;d.x = rect.b.x; d.y=rect.a.y;return ((Distance( Center(circle) , rect.a ) > circle.r) &&(Distance( Center(circle) , rect.b ) > circle.r) &&(Distance( Center(circle) , c ) > circle.r) &&(Distance( Center(circle) , d ) > circle.r) &&(rect.a.x > circle.x || circle.x > rect.b.x || rect.a.y > circle.y || circle.y > rect.b.y)) ||((circle.x - circle.r > rect.b.x) ||(circle.x + circle.r < rect.a.x) ||(circle.y - circle.r > rect.b.y) ||(circle.y + circle.r < rect.a.y)); } //三角形 struct TRIANGLE {POINT a, b, c;TRIANGLE() {};TRIANGLE(const POINT & aa, const POINT & bb, const POINT & cc):a(aa), b(bb), c(cc) {} }; //三角形標(biāo)準(zhǔn)化 TRIANGLE StdTRIANGLE(TYPE len1, TYPE len2, TYPE len3) //把3邊長(zhǎng)轉(zhuǎn)化成三角形 {POINT a, b, c; //已知邊長(zhǎng)后將其轉(zhuǎn)化成坐標(biāo)三角形a.x = a.y = 0.0;b.x = len1;b.y = 0.0;c.x = (len2 * len2 - len3 * len3 + len1 * len1)/len1/2.0;c.y = Sqrt(len2 * len2 - c.x * c.x);TRIANGLE temp(a,b,c);return temp; }; //三角形費(fèi)馬點(diǎn)(即三角形到三頂點(diǎn)之和最小的點(diǎn)) POINT fermentPONT(TRIANGLE &t) {POINT u, v;double step = fabs(t.a.x) + fabs(t.a.y) + fabs(t.b.x) + fabs(t.b.y) + fabs(t.c.x) + fabs(t.c.y);int i, j, k;u.x = (t.a.x + t.b.x + t.c.x)/3;u.y = (t.a.y + t.b.y + t.c.y)/3;while (step > EPS)for(k = 0; k < 10; step /= 2, ++k)for(i = -1; i <= 1; ++i)for(j = -1; j <= 1; ++j) {v.x = u.x + step * i;v.y = u.y + step * j;if(Distance(u,t.a) + Distance(u,t.b) + Distance(u,t.c) > Distance(v,t.a) + Distance(v,t.b) + Distance(v,t.c))u = v;}return u; }; //三角形重心(中線交點(diǎn)) POINT InCenter(const TRIANGLE &t) {return POINT((t.a.x + t.b.x + t.c.x)/3, (t.a.y + t.b.y + t.c.y)/3); } //三角形外心(中垂線交點(diǎn)) POINT CcCenter(const TRIANGLE &t) {POINT u, v;LINE A, B;A.a = t.a; A.b = t.b; B.a = t.b; B.b = t.c;TYPE A1, B1, C1;TYPE A2, B2, C2;Coefficient(A, B1, A1, C1);Coefficient(B, B2, A2, C2);B1 = -B1; B2 = -B2;C1 = -((A.a.x + A.b.x) * A1 + (A.a.y + B.a.y) * B1)/2;C2 = -((B.a.x + B.b.x) * A2 + (B.a.y + B.b.y) * B2)/2;POINT I(0, 0);I.x = - (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1);I.y = (A2 * C1 - A1 * C2) / (A1 * B2 - A2 * B1);return I; } //三角形內(nèi)心(角平分線交點(diǎn)) POINT incenter(const TRIANGLE &t) {LINE u, v;TYPE m, n;u.a = t.a;m = atan2(t.b.y - t.a.y, t.b.x - t.a.x);n = atan2(t.c.y - t.a.y, t.c.x - t.a.x);u.b.x = u.a.x + cos((m + n)/2);u.b.y = u.a.y + sin((m + n)/2);v.a = t.b;m = atan2(t.a.y - t.b.y, t.a.x - t.b.x);n = atan2(t.c.y - t.b.y, t.c.x - t.b.x);v.b.x = v.a.x + cos((m + n)/2);v.b.y = v.a.y + sin((m + n)/2);return Intersection(u, v); }; //三角形內(nèi)切圓面積 TYPE InArea(const TRIANGLE &t) {TYPE p, a, b, c, s;a = Distance(t.a, t.b);c = Distance(t.a, t.c);b = Distance(t.c, t.b);s = (a + b + c)/2;p = s * (s - a) * (s - b) * (s - c);s = p * PI/s/s;return s; } //三角形外接圓面積 TYPE OutArea(const TRIANGLE & t) {TYPE a,b,c;a = (t.a.x - t.b.x) * (t.a.x - t.b.x) + (t.a.y - t.b.y) * (t.a.y - t.b.y);b = (t.a.x - t.c.x) * (t.a.x - t.c.x) + (t.a.y - t.c.y) * (t.a.y - t.c.y);c = (t.c.x - t.b.x) * (t.c.x - t.b.x) + (t.c.y - t.b.y) * (t.c.y - t.b.y);a = PI * sqrt(c/(1-(a+b-c)*(a+b-c)/a/b/4));return a; } // 多邊形 ,逆時(shí)針或順時(shí)針給出x,y struct POLY {int n; //n個(gè)點(diǎn)TYPE *x, *y; //x,y為點(diǎn)的指針，首尾必須重合POLY(): n(0), x(NULL), y(NULL) {};POLY(int nn, const TYPE *xx, const TYPE *yy) {n = nn;x = new TYPE[n + 1];memcpy(x, xx, n*sizeof(TYPE));x[n] = xx[0];y = new TYPE[n + 1];memcpy(y, yy, n*sizeof(TYPE));y[n] = yy[0];} }; //多邊形頂點(diǎn) POINT Vertex(const POLY &poly, int idx) {// idx %= poly.n;return POINT(poly.x[idx], poly.y[idx]); } //多邊形的邊 SEG Edge(const POLY &poly, int idx) {// idx %= poly.n;return SEG(POINT(poly.x[idx], poly.y[idx]), POINT(poly.x[idx + 1], poly.y[idx + 1])); } //多邊形的周長(zhǎng) TYPE Perimeter(const POLY & poly) {TYPE p = 0.0;for (int i = 0; i < poly.n; ++i)p += Distance(Vertex(poly, i), Vertex(poly, i + 1));return p; } //求多邊形面積 TYPE Area(const POLY & poly) {if (poly.n < 3) return TYPE(0);double s = poly.y[0] * (poly.x[poly.n - 1] - poly.x[1]);for (int i = 1; i < poly.n; ++i) {s += poly.y[i] * (poly.x[i - 1] - poly.x[(i + 1) % poly.n]);}return s/2; } //多邊形的重心 POINT InCenter(const POLY & poly) {TYPE S,Si,Ax,Ay;POINT p;Si = (poly.x[poly.n - 1] * poly.y[0] - poly.x[0] * poly.y[poly.n - 1]);S = Si;Ax = Si * (poly.x[0] + poly.x[poly.n - 1]);Ay = Si * (poly.y[0] + poly.y[poly.n - 1]);for(int i = 1; i < poly.n; ++i){Si = (poly.x[i-1] * poly.y[i] - poly.x[i] * poly.y[i-1]);Ax += Si * (poly.x[i-1] + poly.x[i]);Ay += Si * (poly.y[i-1] + poly.y[i]);S += Si;}S = S * 3;return POINT(Ax/S,Ay/S); } //判斷點(diǎn)是否在多邊形上 bool IsOnPoly(const POLY &poly, const POINT &p) {for(int i = 0; i < poly.n; ++i){if( IsOnSeg(Edge(poly, i), p) ) return true;}return false; } //判斷點(diǎn)是否在多邊形內(nèi)部 bool IsInPoly(const POLY &poly, const POINT &p) {SEG L(p, POINT(Infinity, p.y));int count = 0;for (int i = 0; i < poly.n; i++) {SEG S = Edge(poly, i);if (IsOnSeg(S, p)) {return true; //如果想讓 在poly上則返回 true,}if(!IsEqual(S.a.y, S.b.y)) {POINT & q = (S.a.y > S.b.y)?(S.a):(S.b);if(IsOnSeg(L, q)) ++count;else if(!IsOnSeg(L, S.a) && !IsOnSeg(L, S.b) && IsIntersect(S, L)) {++count;}}}return (count % 2 != 0); } //判斷是否簡(jiǎn)單多邊形 bool IsSimple(const POLY & poly) {if(poly.n < 3) return false;SEG L1, L2;for(int i = 0; i < poly.n - 1; ++i) {L1 = Edge(poly, i);for(int j = i + 1; j < poly.n; ++j) {L2 = Edge(poly, j);if(j == i + 1) {if(IsOnSeg(L1, L2.b) || IsOnSeg(L2, L1.a)) return false;}else if(j == poly.n - i - 1) {if (IsOnSeg(L1, L2.a) || IsOnSeg(L2, L1.b)) return false;}else if(IsIntersect(L1, L2)) return false;} // for j} // for ireturn true; } // 點(diǎn)陣的凸包，返回一個(gè)多邊形(求法一，已測(cè)試)，不適用于點(diǎn)少于三個(gè)的情況 POLY ConvexHull(const POINT * set, int n) {POINT * points = new POINT[n];memcpy(points, set, n * sizeof(POINT));TYPE * X = new TYPE[n];TYPE * Y = new TYPE[n];int i, j, k = 0, top = 2;for(i = 1; i < n; ++i) {if ((points[i].y < points[k].y) || ((points[i].y == points[k].y) && (points[i].x < points[k].x)))k = i;}swap(points[0], points[k]);for(i = 1; i < n - 1; ++i) {k = i;for(j = i + 1; j < n; ++j) {if((Cross(points[j], points[k], points[0]) > 0) ||((Cross(points[j], points[k], points[0]) == 0) &&(Distance(points[0], points[j]) < Distance(points[0], points[k]))))k = j;}swap(points[i], points[k]);}X[0] = points[0].x; Y[0] = points[0].y;X[1] = points[1].x; Y[1] = points[1].y;X[2] = points[2].x; Y[2] = points[2].y;for(i = 3; i < n; i++) {while(Cross(points[i], POINT(X[top], Y[top]),POINT(X[top - 1], Y[top - 1])) >= 0)--top;++top;X[top] = points[i].x;Y[top] = points[i].y;}delete [] points;POLY poly(++top, X, Y);delete [] X;delete [] Y;return poly; } // 點(diǎn)陣的凸包，返回一個(gè)多邊形 (求法二，未測(cè)試) bool cmp(const POINT& aa, const POINT& bb) {if(aa.y != bb.y) return aa.y < bb.y;else return aa.x < bb.x; } POLY ConvexHull2(const POINT * set, int n) {// 不適用于點(diǎn)少于三個(gè)的情況POINT *a = new POINT[n];memcpy(a, set, n * sizeof(POINT));sort(a, a + n, cmp);TYPE * X = new TYPE[n];TYPE * Y = new TYPE[n];X[0] = a[0].x;Y[0] = a[0].y;X[1] = a[1].x;Y[1] = a[1].y;int i,j;for(i = 2, j = 2; i < n; ++i) {X[j] = a[i].x;Y[j] = a[i].y;while(((X[j]-X[j-1]) * (Y[j-1]-Y[j-2]) - (Y[j]-Y[j-1]) * (X[j-1]-X[j-2])) >= 0 && j > 1) {X[j-1] = X[j];Y[j-1] = Y[j];--j;}++j;}X[j] = a[n-2].x;Y[j] = a[n-2].y;++j;for(i = n - 3; i >= 0; --i) {X[j] = a[i].x;Y[j] = a[i].y;while(((X[j]-X[j-1]) * (Y[j-1]-Y[j-2]) - (Y[j]-Y[j-1]) * (X[j-1]-X[j-2])) >= 0) {X[j-1] = X[j];Y[j-1] = Y[j];--j;}++j;}delete [] a;POLY poly(j, X, Y);delete [] X;delete [] Y;return poly; }int main() {double a,b,c,l1,l2,l3,l4;POINT d;int n;scanf("%d",&n);while(n--){scanf("%lf%lf%lf",&a,&b,&c);TRIANGLE t=StdTRIANGLE(a,b,c);d=fermentPONT(t);//費(fèi)馬點(diǎn)l1=Distance(d,t.a)+Distance(d,t.b)+Distance(d,t.c);d=incenter(t);//內(nèi)心l2=Distance(d,t.a)+Distance(d,t.b)+Distance(d,t.c);d=InCenter(t);//重心l3=Distance(d,t.a)+Distance(d,t.b)+Distance(d,t.c);d=CcCenter(t);//外心l4=Distance(d,t.a)+Distance(d,t.b)+Distance(d,t.c);printf("%.3lf %.3lf %.3lf %.3lf\n",l1,l2,l3,l4);}return 0; }

其他2：

#include <iostream> #include <cmath> #include <cstdlib> #include <ctime> using namespace std;struct Tpoint{double x, y; };inline double length(Tpoint A, Tpoint B){double x = A.x - B.x, y = A.y - B.y;return sqrt(x * x + y * y); }inline double across(Tpoint A, Tpoint B, Tpoint C){double x1 = C.x - A.x, x2 = C.x - B.x;double y1 = C.y - A.y, y2 = C.y - B.y;return x1 * y2 - x2 * y1; }double distance(Tpoint A, Tpoint B, Tpoint C){return abs(across(A,B,C))/length(B,C); } //http://iask.sina.com.cn/b/14697945.html int main(){Tpoint A, B, C, cen;while(true){cin>>A.x>>A.y>>B.x>>B.y>>C.x>>C.y;/*srand(time(0));A.x = (double)rand()/rand(); A.y = (double)rand()/rand();B.x = (double)rand()/rand(); B.y = (double)rand()/rand();C.x = (double)rand()/rand(); C.y = (double)rand()/rand();cout<<A.x<<' '<<A.y<<' '<<B.x<<' '<<B.y<<' '<<C.x<<' '<<C.y<<endl;*/double a = length(B,C), b = length(C, A), c = length(A, B);if(abs(across(A, B, C)) < 0.01) {cout<<"Wrong Input!"<<endl; continue;}cen.x = (a * A.x + b * B.x + c * C.x)/(a + b + c);cen.y = (a * A.y + b * B.y + c * C.y)/(a + b + c);cout<<cen.x<<' '<<cen.y<<endl;cout<<distance(cen, A, B)<<' '<<distance(cen, B, C)<<' '<<distance(cen, C, A)<<endl;//cout<<abs(across(A, B, C))/(a + b + c)<<endl;system("pause");} }#include <iostream> #include <cmath> #include <cstdlib> #include <ctime> using namespace std;struct Tpoint{double x, y, z; };inline double length(Tpoint A, Tpoint B){double x = A.x - B.x, y = A.y - B.y, z = A.z - B.z;return sqrt(x * x + y * y + z * z); }inline double across(Tpoint A, Tpoint B, Tpoint C){Tpoint p1, p2;p1.x = C.x - A.x; p1.y = C.y - A.y; p1.z = C.z - A.z;p2.x = C.x - B.x; p2.y = C.y - B.y; p2.z = C.z - B.z;double x = p1.y * p2.z - p1.z * p2.y;double y = p1.z * p2.x - p1.x * p2.z;double z = p1.x * p2.y - p1.y * p2.x;//cout<<x<<' '<<y<<' '<<z<<endl;return sqrt(x * x + y * y + z * z); }double distance(Tpoint A, Tpoint B, Tpoint C){return abs(across(A,B,C))/length(B,C); }void neijieyuan(Tpoint A, Tpoint B, Tpoint C){Tpoint cen;srand(time(0));A.x = (double)rand()/rand(); A.y = -(double)rand()/rand();B.x = -(double)rand()/rand(); B.y = (double)rand()/rand();C.x = -(double)rand()/rand(); C.y = (double)rand()/rand();cout<<A.x<<' '<<A.y<<' '<<B.x<<' '<<B.y<<' '<<C.x<<' '<<C.y<<endl;double a = length(B,C), b = length(C, A), c = length(A, B);//cout<<a<<' '<<b<<' '<<c<<endl;if(abs(across(A, B, C)) < 0.01) {cout<<"Wrong Input!"<<endl;return;}cen.x = (a * A.x + b * B.x + c * C.x)/(a + b + c);cen.y = (a * A.y + b * B.y + c * C.y)/(a + b + c);cen.z = (a * A.z + b * B.z + c * C.z)/(a + b + c);cout<<cen.x<<' '<<cen.y<<' '<<cen.z<<endl;cout<<distance(cen, A, B)<<' '<<distance(cen, B, C)<<' '<<distance(cen, C, A)<<endl; }int main(){Tpoint A, B, C;while(true){//cin>>A.x>>A.y>>A.z>>B.x>>B.y>>B.z>>C.x>>C.y>>C.z;//cout<<across(A, B, C)<<endl;neijieyuan(A, B, C);//cout<<abs(across(A, B, C))/(a + b + c)<<endl;system("pause");} }＃include <stdio.h> ＃include <math.h>typedef struct TPoint{//平面點(diǎn)double x;double y;}TPoint;typedef struct TTriangle{TPoint t[2];}TTriangle;typedef struct TCircle{//圓double r;TPoint centre;}TCircle;double distance(const TPoint p1, const TPoint p2){//計(jì)算平面上兩個(gè)點(diǎn)之間的距離return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y)); }double triangleArea(const TTriangle t){//已知三角形三個(gè)頂點(diǎn)的坐標(biāo)，求三角形的面積return fabs(t.t[0].x * t.t[1].y + t.t[1].x * t.t[2].y + t.t[2].x * t.t[0].y- t.t[1].x * t.t[0].y - t.t[2].x * t.t[1].y - t.t[0].x * t.t[2].y) / 2; }TCircle circumcircleOfTriangle(const TTriangle t){//三角形的外接圓TCircle tmp;double a, b, c, c1, c2;double xA, yA, xB, yB, xC, yC;a = distance(t.t[0], t.t[1]);b = distance(t.t[1], t.t[2]);c = distance(t.t[2], t.t[0]);printf("a = %lf b = %lf, c = %lf\n", a, b, c);//根據(jù)S = a * b * c / R / 4;求半徑Rtmp.r = a * b * c / triangleArea(t) / 4;xA = t.t[0].x; yA = t.t[0].y;xB = t.t[1].x; yB = t.t[1].y;xC = t.t[2].x; yC = t.t[2].y;c1 = (xA * xA + yA * yA - xB * xB - yB * yB) / 2;c2 = (xA * xA + yA * yA - xC * xC - yC * yC) / 2;tmp.centre.x = (c1 * (yA - yC) - c2 * (yA - yB)) /(xA - xB) * (yA - yC) - (xA - xC) * (yA - yB);tmp.centre.y = (c1 * (xA - xC) - c2 * (xA - xB)) /(yA - yB) * (xA - xC) - (yA - yC) * (xA - xB);return tmp; }TCircle incircleOfTriangle(const TTriangle t){//三角形的內(nèi)接圓TCircle tmp;double a, b, c, angleA, angleB, angleC, p, p2, p3, f1, f2;double xA, yA, xB, yB, xC, yC;a = distance(t.t[0], t.t[1]);b = distance(t.t[1], t.t[2]);c = distance(t.t[2], t.t[0]);printf("a = %lf b = %lf, c = %lf\n", a, b, c);tmp.r = 2 * triangleArea(t) / (a + b +c);angleA = acos((b * b + c * c - a * a) / (2 * b * c));angleB = acos((a * a + c * c - b * b) / (2 * a * c));angleC = acos((a * a + b * b - c * c) / (2 * a * b));p = sin(angleA / 2);p2 = sin(angleB / 2);p3 = sin(angleC / 2);xA = t.t[0].x; yA = t.t[0].y;xB = t.t[1].x; yB = t.t[1].y;xC = t.t[2].x; yC = t.t[2].y;f1 = ((tmp.r / p2) * (tmp.r / p2) - (tmp.r / p) * (tmp.r / p) +xA * xA - xB * xB + yA * yA - yB * yB) / 2;f2 = ((tmp.r / p3) * (tmp.r / p3) - (tmp.r / p) * (tmp.r / p) +xA * xA - xC * xC + yA * yA - yC * yC) / 2;tmp.centre.x = (f1 * (yA - yC) - f2 * (yA - yB)) /((xA - xB) * (yA - yC) - (xA - xC) * (yA - yB));tmp.centre.y = (f1 * (xA - xC) - f2 * (xA - xB)) /((yA - yB) * (xA - xC) - (yA - yC) * (xA - xB));return tmp; }int main(){TTriangle t;TCircle circumcircle, incircle;while(scanf("%lf%lf%lf%lf%lf%lf",&t.t[0].x, &t.t[0].y, &t.t[1].x, &t.t[1].y,&t.t[2].x, &t.t[2].y) != EOF){incircle = incircleOfTriangle(t);circumcircle = circumcircleOfTriangle(t);printf("%lf %lf \n", incircle.r, circumcircle.r);printf("%.3lf\n", incircle.r / circumcircle.r);} return 0;}

　　

轉(zhuǎn)載于:https://www.cnblogs.com/oyking/p/3269194.html

總結(jié)

以上是生活随笔為你收集整理的算法模板の计算几何的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。