算法模板の计算几何
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 }
2、旋转卡壳（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 }
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 }
4、半平面交+最远点对
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 }
其他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);}
//点
struct POINT {
TYPE x, y;
POINT(TYPE xx = 0, TYPE yy = 0): x(xx), y(yy) {}
};
// 两个点的距离
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) {}
};
//直线(两点式)
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) {}
};
//两点式化一般式
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 的系数
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);
}
//直线相交的交点
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;
}
//点到直线的距离
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);
}
//有公共端点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);
}
//判等(值，点，直线)
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);
}
// 判断点是否在线段上
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)));
}
//判断两条线断是否相交，端点重合算相交
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是否相交，端点重合算相交
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; // 左下点
POINT b; // 右上点
RECT() {}
RECT(const POINT &aa, const POINT &bb): a(aa), b(bb) {}
};
//矩形化标准
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;
}
//根据下标返回矩形的边
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);
}
//两个矩形的公共面积
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;
}
//两个圆的公共面积
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;
}
//判断圆是否在矩形内(不允许相切)
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);
}
//判断矩形是否在圆内(不允许相切)
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) {}
};
//三角形标准化
TRIANGLE StdTRIANGLE(TYPE len1, TYPE len2, TYPE len3) //把3边长转化成三角形
{
POINT a, b, c; //已知边长后将其转化成坐标三角形
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;
};
//三角形费马点(即三角形到三顶点之和最小的点)
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;
};
//三角形重心(中线交点)
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);
}
//三角形外心(中垂线交点)
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;
}
//三角形内心(角平分线交点)
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);
};
//三角形内切圆面积
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;
}
// 多边形 ,逆时针或顺时针给出x,y
struct POLY {
int n; //n个点
TYPE *x, *y; //x,y为点的指针，首尾必须重合
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];
}
};
//多边形顶点
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]));
}
//多边形的周长
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);
}
//判断点是否在多边形上
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;
}
//判断点是否在多边形内部
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);
}
//判断是否简单多边形
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 i
return true;
}
// 点阵的凸包，返回一个多边形(求法一，已测试)，不适用于点少于三个的情况
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;
}
// 点阵的凸包，返回一个多边形 (求法二，未测试)
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) {// 不适用于点少于三个的情况
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);//费马点
l1=Distance(d,t.a)+Distance(d,t.b)+Distance(d,t.c);
d=incenter(t);//内心
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
{
//平面点
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)
{
//计算平面上两个点之间的距离
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
double triangleArea(const TTriangle t)
{
//已知三角形三个顶点的坐标，求三角形的面积
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);
//根据S = a * b * c / R / 4;求半径R
tmp.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)
{
//三角形的内接圆
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;
}
浙公网安备 33010602011771号