poj 2986 A Triangle and a Circle

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<cmath>
  4 #include<cstring>
  5 #include<vector>
  6 #include<algorithm>
  7 using namespace std;
  8 typedef long long LL;
  9 
 10 const double PI = acos(-1.0);
 11 const double EPS = 1e-8;
 12 
 13 inline int sgn(double x) {
 14     return (x > EPS) - (x < -EPS);
 15 }
 16 
 17 struct Point {
 18     double x, y;
 19     Point() {}
 20     Point(double x, double y): x(x), y(y) {}
 21     void read() {
 22         scanf("%lf%lf", &x, &y);
 23     }
 24     double angle() {
 25         return atan2(y, x);
 26     }
 27     Point operator + (const Point &rhs) const {
 28         return Point(x + rhs.x, y + rhs.y);
 29     }
 30     Point operator - (const Point &rhs) const {
 31         return Point(x - rhs.x, y - rhs.y);
 32     }
 33     Point operator * (double t) const {
 34         return Point(x * t, y * t);
 35     }
 36     Point operator / (double t) const {
 37         return Point(x / t, y / t);
 38     }
 39     double operator *(const Point &b)const
 40     {
 41         return x*b.x + y*b.y;
 42     }
 43     double length() const {
 44         return sqrt(x * x + y * y);
 45     }
 46     Point unit() const {            //单位向量
 47         double l = length();
 48         return Point(x / l, y / l);
 49     }
 50 };
 51 double cross(const Point &a, const Point &b) {
 52     return a.x * b.y - a.y * b.x;
 53 }
 54 double sqr(double x) {
 55     return x * x;
 56 }
 57 double dist(const Point &p1, const Point &p2) {
 58     return (p1 - p2).length();
 59 }
 60 double sdist(Point a,Point b){
 61     return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
 62 }
 63 //向量 op 逆时针旋转 angle
 64 Point rotate(const Point &p, double angle, const Point &o = Point(0, 0)) {
 65     Point t = p - o;
 66     double x = t.x * cos(angle) - t.y * sin(angle);
 67     double y = t.y * cos(angle) + t.x * sin(angle);
 68     return Point(x, y) + o;
 69 }
 70 Point line_inter(Point A,Point B,Point C,Point D){ //直线相交交点
 71         Point ans;
 72         double a1=A.y-B.y;
 73         double b1=B.x-A.x;
 74         double c1=A.x*B.y-B.x*A.y;
 75 
 76         double a2=C.y-D.y;
 77         double b2=D.x-C.x;
 78         double c2=C.x*D.y-D.x*C.y;
 79 
 80         ans.x=(b1*c2-b2*c1)/(a1*b2-a2*b1);
 81         ans.y=(a2*c1-a1*c2)/(a1*b2-a2*b1);
 82         return ans;
 83 }
 84 Point p_to_seg(Point p,Point a,Point b){        //点到线段的最近点
 85     Point tmp=p;
 86     tmp.x+=a.y-b.y;
 87     tmp.y+=b.x-a.x;
 88     if(cross(a-p,tmp-p)*cross(b-p,tmp-p)>0) return dist(p,a)<dist(p,b)?a:b;
 89     return line_inter(p,tmp,a,b);
 90 }
 91 void line_circle(Point c,double r,Point a,Point b,Point &p1,Point &p2){
 92     Point tmp=c;
 93     double t;
 94     tmp.x+=(a.y-b.y);//求垂直于ab的直线
 95     tmp.y+=(b.x-a.x);
 96     tmp=line_inter(tmp,c,a,b);
 97     t=sqrt(sqr(r)-sqr( dist(c,tmp)))/dist(a,b); //比例
 98     p1.x=tmp.x+(b.x-a.x)*t;
 99     p1.y=tmp.y+(b.y-a.y)*t;
100     p2.x=tmp.x-(b.x-a.x)*t;
101     p2.y=tmp.y-(b.y-a.y)*t;
102 }
103 struct Region {
104     double st, ed;
105     Region() {}
106     Region(double st, double ed): st(st), ed(ed) {}
107     bool operator < (const Region &rhs) const {
108         if(sgn(st - rhs.st)) return st < rhs.st;
109         return ed < rhs.ed;
110     }
111 };
112 struct Circle {
113     Point c;
114     double r;
115     vector<Region> reg;
116     Circle() {}
117     Circle(Point c, double r): c(c), r(r) {}
118     void read() {
119         c.read();
120         scanf("%lf", &r);
121     }
122     void add(const Region &r) {
123         reg.push_back(r);
124     }
125     bool contain(const Circle &cir) const {
126         return sgn(dist(cir.c, c) + cir.r - r) <= 0;
127     }
128     bool intersect(const Circle &cir) const {
129         return sgn(dist(cir.c, c) - cir.r - r) < 0;
130     }
131 };
132 void intersection(const Circle &cir1, const Circle &cir2, Point &p1, Point &p2) {   //两圆相交 交点
133     double l = dist(cir1.c, cir2.c);                            //两圆心的距离
134     double d = (sqr(l) - sqr(cir2.r) + sqr(cir1.r)) / (2 * l);  //cir1圆心到交点直线的距离
135     double d2 = sqrt(sqr(cir1.r) - sqr(d));                     //交点到 两圆心所在直线的距离
136     Point mid = cir1.c + (cir2.c - cir1.c).unit() * d;
137     Point v = rotate(cir2.c - cir1.c, PI / 2).unit() * d2;
138     p1 = mid + v, p2 = mid - v;
139 }
140 Point calc(const Circle &cir, double angle) {
141     Point p = Point(cir.c.x + cir.r, cir.c.y);
142     return rotate(p, angle, cir.c);
143 }
144 const int MAXN = 1010;
145 Circle cir[MAXN],cir2[MAXN];
146 bool del[MAXN];
147 int n;
148 double get_area(Circle* cir,int n) {            //多个圆的相交面积
149     double ans = 0;
150     memset(del,0,sizeof(del));
151     for(int i = 0; i < n; ++i) {
152         for(int j = 0; j < n; ++j) if(!del[j]) {                //删除被包含的圆
153             if(i == j) continue;
154             if(cir[j].contain(cir[i])) {
155                 del[i] = true;
156                 break;
157             }
158         }
159     }
160     for(int i = 0; i < n; ++i) if(!del[i]) {
161         Circle &mc = cir[i];
162         Point p1, p2;
163         bool flag = false;
164         for(int j = 0; j < n; ++j) if(!del[j]) {
165             if(i == j) continue;
166             if(!mc.intersect(cir[j])) continue;
167             flag = true;
168             intersection(mc, cir[j], p1, p2);                   //求出两圆的交点
169             double rs = (p2 - mc.c).angle(), rt = (p1 - mc.c).angle();
170             if(sgn(rs) < 0) rs += 2 * PI;
171             if(sgn(rt) < 0) rt += 2 * PI;
172             if(sgn(rs - rt) > 0) mc.add(Region(rs, PI * 2)), mc.add(Region(0, rt)); //添加相交区域
173             else mc.add(Region(rs, rt));
174         }
175         if(!flag) {
176             ans += PI * sqr(mc.r);
177             continue;
178         }
179         sort(mc.reg.begin(), mc.reg.end());                 //对相交区域进行排序
180         int cnt = 1;
181         for(int j = 1; j < int(mc.reg.size()); ++j) {
182             if(sgn(mc.reg[cnt - 1].ed - mc.reg[j].st) >= 0) {   //如果有区域可以合并，则合并
183                 mc.reg[cnt - 1].ed = max(mc.reg[cnt - 1].ed, mc.reg[j].ed);
184             } else mc.reg[cnt++] = mc.reg[j];
185         }
186         mc.add(Region());
187         mc.reg[cnt] = mc.reg[0];
188         for(int j = 0; j < cnt; ++j) {
189             p1 = calc(mc, mc.reg[j].ed);
190             p2 = calc(mc, mc.reg[j + 1].st);
191             ans += cross(p1, p2) / 2;                           //
192             double angle = mc.reg[j + 1].st - mc.reg[j].ed;
193             if(sgn(angle) < 0) angle += 2 * PI;
194             ans += 0.5 * sqr(mc.r) * (angle - sin(angle));      //弧所对应的的面积
195         }
196     }
197     return ans;
198 }
199 double two_cir(Circle t1,Circle t2){            //两个圆的相交面积
200     if(t1.contain(t2)||t2.contain(t1))    return PI * sqr(min(t2.r,t1.r));
201     if(!t1.intersect(t2)) return 0;
202     double ans=0,len=dist(t1.c,t2.c);
203     double x=(sqr(t1.r)+sqr(len)-sqr(t2.r))/(2*len);
204     double angle1=acos(x/t1.r),angle2=acos((len-x)/t2.r);
205     ans=sqr(t1.r)*angle1+sqr(t2.r)*angle2-len*t1.r*sin(angle1);    // 两个扇形 减去一个四边形面积
206     return ans;
207 }
208 double triangle_circle(Point a,Point b,Point o,double r){
209     double sign=1.0;
210     double ans=0;
211     Point p1,p2;
212     a=a-o;b=b-o;
213     o=Point(0,0);
214     if(sgn(cross(a,b))==0) return 0;
215     if(sdist(a,o)>sdist(b,o)){
216         swap(a,b);
217         sign=-1.0;
218     }
219     if(sdist(a,o)<r*r+EPS){  //a 在内， b 在外
220         if(sdist(b,o)<r*r+EPS) return cross(a,b)/2.0*sign;
221         line_circle(o,r,a,b,p1,p2);
222         if(dist(p1,b)>dist(p2,b)) swap(p1,p2);
223         double ans1=fabs(cross(a,p1));
224         double ans2=acos(p1*b/p1.length()/b.length())*r*r;
225         ans=(ans1+ans2)/2.0;
226         if(cross(a,b)<EPS&&sign>0.0||cross(a,b)>EPS&&sign<0.0) return -ans;
227         return ans;
228     }
229     Point tmp=p_to_seg(o,a,b);
230     if(sdist(o,tmp)>r*r-EPS){    //a,b所在直线与圆没有交点
231         ans=acos(a*b/a.length()/b.length())*r*r/2.0;
232         if(cross(a,b)<EPS&&sign>0.0||cross(a,b)>EPS&&sign<0.0) return -ans;
233         return ans;
234     }
235     line_circle(o,r,a,b,p1,p2);
236     double cm=r/(dist(a,o)-r);
237     Point m=Point((o.x+cm*a.x)/(1+cm),(o.y+cm*a.y)/(1+cm) );
238     double cn=r/(dist(o,b)-r);
239     Point n=Point( (o.x+cn*b.x)/(1+cn) , (o.y+cn*b.y)/(1+cn) );
240     double ans1 = acos( m*n/m.length()/n.length() )*r*r;
241     double ans2 = acos( p1*p2/p1.length()/p2.length() )*r*r-fabs(cross(p1,p2));
242     ans=(ans1-ans2)/2.0;
243     if(cross(a,b)<EPS&&sign>0.0||cross(a,b)>EPS&&sign<0.0) return -ans;
244     return ans;
245 }
246 int main() {
247     #ifndef ONLINE_JUDGE
248     freopen("input.txt","r",stdin);
249     #endif // ONLINE_JUDGE
250     Point a,b,c;
251     int cas=0;
252     Circle t1;
253     double d1,d2,d3,d4,d5,d6,d7,d8,d9;
254     while(scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf",&d1,&d2,&d3,&d4,&d5,&d6,&d7,&d8,&d9)!=EOF){
255         a.x=d1;a.y=d2;
256         b.x=d3;b.y=d4;
257         c.x=d5;c.y=d6;
258         t1.c.x=d7;t1.c.y=d8;t1.r=d9;
259         double ans=0;
260         ans+=triangle_circle(a,b,t1.c,t1.r);
261         ans+=triangle_circle(b,c,t1.c,t1.r);
262         ans+=triangle_circle(c,a,t1.c,t1.r);
263         printf("%.2f\n",fabs(ans)+EPS);
264     }
265 }