?? convex_hull.cpp
字號:
#include <iostream>
#include <cmath>
using namespace std;
#define EPS 1e-10
#define MAXN 1000
struct point
{
double x, y;
int id;
point() {}
point(double nx, double ny) : x(nx), y(ny) {}
point initi() const
{
double l = hypot(x, y);
return point(x / l, y / l);
}
point rotate(double ang) const
{
return point(x * cos(ang) - y * sin(ang), x * sin(ang) + y * cos(ang));
}
};
point p[MAXN]; int pnum;
point h[MAXN]; int chn;
int bt[MAXN], tp[MAXN];
int bt_top, tp_top;
struct circle
{
double x, y;
double r;
int id;
};
circle c[MAXN]; int cnum;
circle t[MAXN];
void process()
{
int i, j, k = 0;
for(i = 0; i < cnum; i++)
{
for(j = 0; j < cnum; j++)
{
if(i == j) continue;
double dist = hypot(c[i].x - c[j].x, c[i].y - c[j].y);
if(c[j].r >= c[i].r && c[j].r - c[i].r >= dist) break;
}
if(j == cnum && c[i].r > EPS) t[k++] = c[i];
}
//printf("%d\n", k);
cnum = k;
for(i = 0; i < cnum; i++)
{
c[i] = t[i];
c[i].id = i;
//printf("%lf %lf %lf\n", c[i].x, c[i].y, c[i].r);
}
//putchar('\n');
}
void add_point()
{
pnum = 0;
int i, j;
for(i = 0; i < cnum; i++)
{
for(j = i + 1; j < cnum; j++)
{
circle c1 = c[i], c2 = c[j];
if(c1.r < c2.r) swap(c1, c2);
double d = hypot(c1.x - c2.x, c1.y - c2.y);
double h = c1.r - c2.r;
double theta = fabs(acos(h / d));
point p1, p2;
point vt, vl, vr;
vt.x = c2.x - c1.x, vt.y = c2.y - c1.y;
vt = vt.initi();
vl = vt.rotate(theta), vr = vt.rotate(-theta);
p1.x = c1.x + c1.r * vl.x;
p1.y = c1.y + c1.r * vl.y;
p1.id = c1.id;
p2.x = c2.x + c2.r * vl.x;
p2.y = c2.y + c2.r * vl.y;
p2.id = c2.id;
p[pnum++] = p1;
p[pnum++] = p2;
p1.x = c1.x + c1.r * vr.x;
p1.y = c1.y + c1.r * vr.y;
p1.id = c1.id;
p2.x = c2.x + c2.r * vr.x;
p2.y = c2.y + c2.r * vr.y;
p2.id = c2.id;
p[pnum++] = p1;
p[pnum++] = p2;
}
}
/*
for(i = 0; i < pnum; i++)
{
printf("%lf %lf %d\n", p[i].x, p[i].y, p[i].id);
}
putchar('\n');
*/
}
bool cmp(const point &p1, const point &p2)
{
return (p1.x < p2.x) || (fabs(p1.x - p2.x) < EPS && p1.y < p2.y);
}
int dcmp(double d)
{
if(fabs(d) < EPS) return 0;
return d > 0 ? 1 : -1;
}
int det(double x1, double y1, double x2, double y2)
{
return dcmp(x2 * y1 - x1 * y2);
}
int cross(point a, point b, point c)
{
return det(b.x - a.x, b.y - a.y, c.x - b.x, c.y - b.y);
}
void convex_hull()
{
int i;
sort(p, p + pnum, cmp);
bt[0] = 0, tp[0] = 0, bt[1] = 1, tp[1] = 1;
bt_top = 1, tp_top = 1;
for(i = 0; i < pnum; i++)
{
if(cross(p[bt[1]], p[0], p[i]) <= 0) bt[1] = i;
if(cross(p[tp[1]], p[0], p[i]) >= 0) tp[1] = i;
}
for(i = bt[bt_top] + 1; i < pnum; i++)
{
while(1)
{
if(cross(p[bt[bt_top - 1]], p[bt[bt_top]], p[i]) >= 0) bt_top--;
else
{
bt[++bt_top] = i;
break;
}
}
}
for(i = tp[tp_top] + 1; i < pnum; i++)
{
while(1)
{
if(cross(p[tp[tp_top - 1]], p[tp[tp_top]], p[i]) <= 0) tp_top--;
else
{
tp[++tp_top] = i;
break;
}
}
}
chn = 0;
for(i = 0; i <= bt_top; i++) h[chn++] = p[bt[i]];
for(i = tp_top - 1; i >= 1; i--) h[chn++] = p[tp[i]];
}
void calc_total_area()
{
h[chn] = h[0];
double area = 0;
int i;
/*
for(i = 0; i < chn; i++)
{
printf("%lf %lf %d\n", h[i].x, h[i].y, h[i].id);
}
putchar('\n');
*/
for(i = 0; i < chn; i++)
{
area += h[i].x * h[i + 1].y - h[i].y * h[i + 1].x;
}
area = fabs(area / 2.0);
//printf("area = %lf\n", area);
for(i = 0; i < chn; i++)
{
if(h[i].id == h[i + 1].id)
{
int d = h[i].id;
point ct = point(c[d].x, c[d].y);
int f = cross(h[i], h[i + 1], ct);
double dist = hypot(h[i].x - h[i + 1].x, h[i].y - h[i + 1].y);
double theta = fabs(2 * asin(dist / (2 * c[d].r)));
//printf("%lf\n", theta);
double arc_area = c[d].r * c[d].r * (theta - sin(theta)) * .5;
if(f > 0)
{
//printf("YES\n");
arc_area = M_PI * c[d].r * c[d].r - arc_area;
}
area += arc_area;
}
}
printf("%lf\n", area);
}
int main(void)
{
int i;
while(scanf("%d", &cnum) != EOF)
{
for(i = 0; i < cnum; i++)
{
scanf("%lf %lf %lf", &c[i].x, &c[i].y, &c[i].r);
}
process();//去掉包含和半徑為0的圓
if(cnum == 1)
{
printf("%lf\n", M_PI * c[0].r * c[0].r);
continue;
}
add_point();
convex_hull();
calc_total_area();
}
return 0;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -