圆与任意多边形交的长度

最新推荐文章于 2023-07-14 16:31:19 发布
PSU_VJD
最新推荐文章于 2023-07-14 16:31:19 发布
阅读量628
点赞数
分类专栏：计算几何
本文链接：https://blog.csdn.net/PSU_VJD/article/details/49150943
版权
计算几何专栏收录该内容
5 篇文章 0 订阅
订阅专栏
#include <bits/stdc++.h>
#define EXIT exit(0);
#define DEBUG puts("Here is a BUG");
#define what_is(x) cout << #x << " is " << x << endl;
using namespace std;

const int maxn = 1000 + 10;
const double eps = 1e-8;
const double PI = acos(-1.0);
const int offset=1e8;
int dcmp(double x){
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}

struct Point{
     double x, y;
     Point(double x = 0, double y = 0) : x(x), y(y) {}
     void read(){
        cin >> x >> y;
     }
     bool operator < (const Point &rhs) const{
          if(dcmp(x - rhs.x) != 0) return x < rhs.x;
          return y < rhs.y;
     }
};

typedef Point Vector;
Vector operator + (const Vector &A, const Vector &B){ return Vector(A.x + B.x, A.y + B.y); }
bool operator == (const Vector &A, const Vector &B){ return dcmp(A.x - B.x) == 0 && dcmp(A.y - B.y) == 0; }
Vector operator - (const Vector &A, const Vector &B){ return Vector(A.x - B.x, A.y - B.y); }
Vector operator * (const Vector &A, double b){ return Vector(A.x*b, A.y*b); }
Vector operator / (const Vector &A, double b){ return Vector(A.x/b, A.y/b); }
double Dot(Point a, Point b){ return a.x * b.x + a.y * b.y; }
double Length(const Vector &A){ return sqrt(Dot(A, A)); }
double Angle(Point A, Point B){ return acos(Dot(A, B) / Length(A) / Length(B)); }
double Dist2(Point A, Point B){ return (A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y); }
double Cross(Point A, Point B){ return A.x * B.y - A.y * B.x; }
bool cmp_(const Point& p1, const Point& p2){ return atan2(p1.y, p1.x) > atan2(p2.y, p2.x); }

struct Line{
     Point p;
     Vector v;
     double ang;
     Line(){}
     Line(Point p, Vector v) : p(p), v(v){
          ang = atan2(v.y, v.x);
     }
     bool operator < (const Line &L)const{
          return ang < L.ang;
     }
     Point point(double t){
          return p + v * t;
     }
};
struct Circle{
    Point c;
    double r;
    Circle() {}
    Circle(Point c, double r) : c(c), r(r) {}
    Point point(double a){
        return Point(c.x + cos(a) * r, c.y + sin(a) * r);
    }
    double ccfr(void) const {
        return 2 * PI * r;
    }
    void read(){
        c.read();
        cin >> r;
    }
};
#define Vector Point
typedef vector<Point> Polygon;
// int n;
Polygon pols[305]; //pWall
Circle center; //hit
Circle cir[305]; //cWall
vector<Point> intersection;//保存圆与多边形的交点
double segment, rad_len;

bool in(Point a, Point l, Point r){
     double max_x = max(l.x, r.x);
     double min_x = min(l.x, r.x);
     double max_y = max(l.y, r.y);
     double min_y = min(l.y, r.y);
     if(dcmp(max_x - a.x) >= 0 && dcmp(a.x - min_x) >= 0 && dcmp(max_y - a.y) >= 0 && dcmp(a.y - min_y) >= 0) return true;
     return false;
}

void get_Line_CircleIntersection(Line L, Circle C, vector<Point>& sol, Point A, Point B){
     double t1, t2;
     double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
     double e = a * a + c * c , f = 2 * (a * b + c * d), g= b * b + d * d - C.r * C.r;
     double delta = f * f - 4 * e * g;
     Point cur;
     sol.push_back(A);
     sol.push_back(B);
     if(dcmp(delta) < 0) return;
     if(dcmp(delta) == 0){
          t1 = t2 = -f / (2 * e);
          cur = L.point(t1);
          if(in(cur, A, B)) sol.push_back(cur), intersection.push_back(cur);
          return;
     }
     t1 = (-f + sqrt(delta)) / (2 * e);
     cur = L.point(t1);
     if(in(cur, A, B)) sol.push_back(cur), intersection.push_back(cur);
     t2 = (-f - sqrt(delta)) / (2 * e);
     cur = L.point(t2);
     if(in(cur, A, B)) sol.push_back(cur), intersection.push_back(cur);
}
/*void init(){
     rad_len = 0;
     segment = 0;
     intersection.clear();//保存圆与多边形的交点
}*/
bool Point_in_Circle(Point A){
     if(dcmp(center.r * center.r - Dist2(A, center.c)) > 0) return true;
     return false;
}
void get_Circle_Segment_Intersection(Point a, Point b){
     vector<Point> save;
     get_Line_CircleIntersection(Line(a, a-b), center, save, a, b);
     int size_ = save.size();
     sort(save.begin(), save.end());
     for(int i = 0; i < size_ - 1; i++){
          Point mid = (save[i] + save[i+1]) / 2;
          if(Point_in_Circle(mid)) segment += Length(save[i] - save[i+1]);
     }
}
int Circle_in_Ploygon(const Point *arr,const int &len,const Point &A,int on_edge)
{
    Point q;
    int i=0,counter;
    while(i<len)
    {
        q.x=rand()+offset;//随机取一个足够远的点q
        q.y=rand()+offset;//以p为起点q为终点做射线L
        for(counter=i=0; i<len; i++) //依次对多边形的每条边进行考察
        {
            Point t1 = arr[i]-A;
            Point t2 = arr[(i+1)%len]-A;
            Point t3 = q-A;
            Point t4 = arr[i]-A;
            Point t5 = A-arr[i];
            Point t6 = arr[(i+1)%len]-arr[i];
            Point t7 = q-arr[i];
            if(fabs(Cross(t1,t2))<eps &&
                    (arr[i].x-A.x)*(arr[(i+1)%len].x-A.x)<eps && (arr[i].y-A.y)*(arr[(i+1)%len].y-A.y)<eps)
                return on_edge; //点p在边上,返回on_edge
            else if(fabs(Cross(t3,t4))<eps) break; //点arr[i]在射线pq上，停止本循环，另取q
            else if(Cross(t4,t3)*Cross(t2,t3)<-eps &&
                    Cross(t5,t6)*Cross(t7,t6)<-eps)
                counter++;
        }
    }
    return counter&1;
}
double Angle_(Point A){
    double ang = Angle(A, Point(1, 0));
    if(dcmp(Cross(Point(1, 0), A)) < 0) ang = - ang;
    return ang;
}
Point get_rad_midPoint(Point A, Point B){
    double rad_A_x = Angle_(A);
    double rad_B_x = Angle_(B);
    double rad_A_B = Angle(A, B);
    double rad_;
    if(dcmp(Cross(A, B) > 0)) rad_ = rad_A_x - (PI - rad_A_B / 2);
    else rad_ = rad_A_x - rad_A_B / 2;
    return Circle(Point(0, 0), center.r).point(rad_);
}
/*void input(){
    init();
    for(int i = 0; i < n; i++) {
        p[i].read();
    }
    center.read();
}
void solve_ploygon(){
    for(int i = 0; i < n; i++) {
        get_Circle_Segment_Intersection(p[i], p[(i+1)%n]);
    }
}
void solve_circle(){
    double sum_rad = 0;
    for(int i = 0; i < (int)intersection.size(); i++)
        intersection[i] = intersection[i] - center.c;
    sort(intersection.begin(), intersection.end(), cmp_);
    intersection.erase(unique(intersection.begin(), intersection.end()), intersection.end());
    int size_ = intersection.size();
    if(size_ <= 1 && Circle_in_Ploygon(p, n, center.point(0), 1) && Circle_in_Ploygon(p, n, center.point(1), 1)){
        rad_len = 2 * PI * center.r;
        return;
    }
    for(int i = 0; i < n; i++) p[i] = p[i] - center.c;
    for(int i = 0; i < size_; i++){
        int j = (i + 1) % size_;
        Point mid = get_rad_midPoint(intersection[i], intersection[j]);
        if(Circle_in_Ploygon(p, n, mid, 0)){
            double cur_rad = Angle(intersection[i], intersection[j]);
            if(dcmp(Cross(intersection[i], intersection[j])) > 0) cur_rad = 2 * PI - cur_rad;
            sum_rad += cur_rad;
        }
    }
    rad_len = sum_rad * center.r;
}*/
int cnt_cir = 0, cnt_pol = 0;

#define ChongHe 0
#define NeiHan 1
#define NeiQie 2
#define INTERSECTING 3
#define WaiQie 4
#define XiangLi 5

//两圆位置关系判定
int GetCircleLocationRelation(const Circle& A, const Circle& B) {
    double d = Length(A.c-B.c);
    double sum = A.r + B.r;
    double sub = fabs(A.r - B.r);
    if (dcmp(d) == 0) return dcmp(sub) != 0;
    if (dcmp(d - sum) > 0) return XiangLi;
    if (dcmp(d - sum) == 0) return WaiQie;
    if (dcmp(d - sub) > 0 && dcmp(d - sum) < 0) return INTERSECTING;
    if (dcmp(d - sub) == 0) return NeiQie;
    if (dcmp(d - sub) < 0) return NeiHan;
}

//两圆相交的面积
double GetCircleIntersectionArea(const Circle& A, const Circle& B) {
    int rel = GetCircleLocationRelation(A, B);
    if (rel < INTERSECTING) {
        double a = A.ccfr(), b = B.ccfr();
        // printf("%.10lf %.10lf %d\n", A.r, B.r, dcmp(a - b));
        if (dcmp(a - b) < 0) return 0;
        return b;
    }
    if (rel > INTERSECTING) return 0;
    double dis = Length(A.c - B.c);
    // double ang1 = acos((A.r*A.r + dis*dis - B.r*B.r) / (2.0*A.r*dis));
    double ang2 = acos((B.r*B.r + dis*dis - A.r*A.r) / (2.0*B.r*dis));
    // return ang1*A.r*A.r + ang2*B.r*B.r - A.r*dis*sin(ang1);
    return ang2 * B.r * 2;
}

bool isOk(Circle center, double sum) {
    intersection.clear();
    segment = 0;
    for(int i = 0; i < cnt_pol; i++) {
        int n = pols[i].size();
        for(int j = 0; j < n; j++) {
            get_Circle_Segment_Intersection(pols[i][j], pols[i][(j+1)%n]);
        }
    }
    double cir_cir = 0;
    for(int i = 0; i < cnt_cir; i++) {
        cir_cir += GetCircleIntersectionArea(center, cir[i]);
        // what_is(cir_cir);
    }
    // what_is(cir_cir);
    // what_is(sum);
    return dcmp(segment+cir_cir - sum) == 0;
}

double judge(Circle center, double sum) {
    intersection.clear();
    segment = 0;
    for(int i = 0; i < cnt_pol; i++) {
        int n = pols[i].size();
        for(int j = 0; j < n; j++) {
            get_Circle_Segment_Intersection(pols[i][j], pols[i][(j+1)%n]);
        }
    }
    double cir_cir = 0;
    for(int i = 0; i < cnt_cir; i++) {
        cir_cir += GetCircleIntersectionArea(center, cir[i]);
        // what_is(cir_cir);
    }
    // what_is(segment);
    // what_is(cir_cir);
    return segment + cir_cir;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("D:\\Documents\\Disk_Synchronous\\Programs\\Acm\\input.txt", "r", stdin);
#endif
    int T; cin >> T;
    for(int kase = 1; kase <= T; kase++) {
        printf("Case #%d: ", kase);
        int N;
        double sum;
        cin >> N >> sum;
        cnt_cir = cnt_pol = 0;
        double sum_ccrf = 0;
        for(int i = 0; i < N; i++) {
            char Type[10];
            scanf("%s", Type);
            if (Type[0] == 'C') {
                // DEBUG;
                cir[cnt_cir++].read();
                sum_ccrf += cir[cnt_cir-1].ccfr();
            }
            else {
                int n; scanf("%d", &n);
                pols[cnt_pol].clear();
                for(int i = 0; i < n; i++) {
                    Point tmp; tmp.read();
                    pols[cnt_pol].push_back(tmp);
                }
                for(int i = 0; i < n; i++) {
                    sum_ccrf += Length(pols[cnt_pol][i] - pols[cnt_pol][(i+1)%n]);
                }
                cnt_pol++;
            }
        }
        cin >> center.c.x >> center.c.y;
        // cout << sum_ccrf << ' ' << sum << endl;
        if (dcmp(sum_ccrf - sum) == 0) {
            puts("inestimable");
            continue;
        }
        if (dcmp(sum_ccrf - sum) < 0) {
            puts("impossible");
            continue;
        }
        double flag = -1;
        double l = 0, r = INT_MAX, mid;
        for(int i = 0; i < 150; i++) {
            mid = (l + r) * 0.5;
            // mid = 1.22;
            center.r = mid;
            // what_is(mid);
            double res = judge(center, sum);
            if (dcmp(res - sum) <= 0) {
                l = mid;
                if (dcmp(res - sum) == 0) flag = mid;
            }
            else r = mid;
        }
        center.r = l;
        // what_is(l);
        // puts("*******************************");
        if (dcmp(flag + 1) == 0) puts("impossible");
        else printf("%.2lf\n", flag);
    }
    return 0;
}