HDU 5230 (计算几何圆和多边形面积交)

最新推荐文章于 2021-05-17 08:48:37 发布

原创最新推荐文章于 2021-05-17 08:48:37 发布 · 504 阅读

1 ·

本内容遵循CC 4.0 BY-SA版权协议

平面几何专栏收录该内容

28 篇文章

订阅专栏

本文介绍了一种计算给定多边形区域内，满足特定条件的安全区域面积的方法。通过数学推导，将问题转换为求解一个圆与多边形的交集面积。提供了完整的C++代码实现，包括点、直线、圆等基本几何对象的操作。

题目链接:点击这里

题意:给出一个多边形区域, 给出两个点 $A$ , $B$ , 某一个点距离 $B$ 的距离为 $d_B$ , 距离 $A$ 的距离为 $d_A$ , 如果 $k\times d_A \leq d_B$ 一个点是安全的. 求安全的面积.

先假设两个点是 $(x_A, y_A)$ , $(x_B, y_B)$ , 假设点 $(x, y)$ 是安全的, 那么满足 $d \times (x - x A) 2 + (y - y a) 2 - - - - - - - - - - - - - - - - - \sqrt \leq (x - x B) 2 + (y - y B) 2 - - - - - - - - - - - - - - - - - - \sqrt$ $d\times \sqrt{(x-x_A)^2+(y-y_a)^2}\leq\sqrt{(x-x_B)^2+(y-y_B)^2}$

平方整理一下就是 $(d 2 - 1) x 2 - (2 d x A - 2 x B) x + d 2 x 2 A - x 2 B + (d 2 - 1) y - (2 d y A - 2 y B) y + d 2 y 2 A - y 2 B \leq 0$ $(d^2-1)x^2-(2dx_A-2x_B)x+d^2x_A^2-x_B^2+(d^2-1)y-(2dy_A-2y_B)y+d^2y_A^2-y_B^2\leq 0$

这个方程整理一下就是圆方程了. 所以所求的面积就是一个圆和一个多边形的面积交.

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <vector>
using namespace std;

const double eps = 1e-8;
const double INF = 1e20;
const double pi = acos (-1.0);

int dcmp (double x) {
    if (fabs (x) < eps) return 0;
    return (x < 0 ? -1 : 1);
}
inline double sqr (double x) {return x*x;}

//*************点
struct Point {
    double x, y;
    Point (double _x = 0, double _y = 0):x(_x), y(_y) {}
    void input () {scanf ("%lf%lf", &x, &y);}
    void output () {printf ("%.2f %.2f\n", x, y);}
    bool operator == (const Point &b) const {
        return (dcmp (x-b.x) == 0 && dcmp (y-b.y) == 0);
    }
    bool operator < (const Point &b) const {
        return (dcmp (x-b.x) == 0 ? dcmp (y-b.y) < 0 : x < b.x);
    }
    Point operator + (const Point &b) const {
        return Point (x+b.x, y+b.y);
    }
    Point operator - (const Point &b) const {
        return Point (x-b.x, y-b.y);
    }
    Point operator * (double a) {
        return Point (x*a, y*a);
    }
    Point operator / (double a) {
        return Point (x/a, y/a);
    }
    double len2 () {//返回长度的平方
        return sqr (x) + sqr (y);
    }
    double len () {//返回长度
        return sqrt (len2 ());
    }
    Point change_len (double r) {//转化为长度为r的向量
        double l = len ();
        if (dcmp (l) == 0) return *this;//零向量返回自身
        r /= l;
        return Point (x*r, y*r);
    }
    Point rotate_left () {//顺时针旋转90度
        return Point (-y, x);
    }
    Point rotate_right () {//逆时针旋转90度
        return Point (y, -x);
    }
    Point rotate (Point p, double ang) {//绕点p逆时针旋转ang
        Point v = (*this)-p;
        double c = cos (ang), s = sin (ang);
        return Point (p.x + v.x*c - v.y*s, p.y + v.x*s + v.y*c);
    }
    Point normal () {//单位法向量
        double l = len ();
        return Point (-y/l, x/l);
    }
};

double cross (Point a, Point b) {//叉积

    return a.x*b.y-a.y*b.x;
}
double dot (Point a, Point b) {//点积

    return a.x*b.x + a.y*b.y;
}
double dis (Point a, Point b) {//两个点的距离

    Point p = b-a; return p.len ();
}
double rad_degree (double rad) {//弧度转化为角度

    return rad/pi*180;
}
double rad (Point a, Point b) {//两个向量的夹角

    return fabs (atan2 (fabs (cross (a, b)), dot (a, b)) );
}
bool parallel (Point a, Point b) {//向量平行
    double p = rad (a, b);
    return dcmp (p) == 0 || dcmp (p-pi) == 0;
}

//************直线 线段
struct Line {
    Point s, e;//直线的两个点
    Line () {}
    Line (Point _s, Point _e) : s(_s), e(_e) {}
    //一个点和倾斜角确定直线
    Line (Point p, double ang) {
        s = p;
        if (dcmp (ang-pi/2) == 0) {
            e = s + Point (0, 1);
        }
        else
            e = s + Point (1, tan (ang));
    }
    //ax+by+c=0确定直线
    Line (double a, double b, double c) {
        if (dcmp (a) == 0) {
            s = Point (0, -c/b);
            e = Point (1, -c/b);
        }
        else if (dcmp (b) == 0) {
            s = Point (-c/a, 0);
            e = Point (-c/a, 1);
        }
        else {
            s = Point (0, -c/b);
            e = Point (1, (-c-a)/b);
        }
    }
    void input () {
        s.input ();
        e.input ();
    }
    void adjust () {
        if (e < s) swap (e, s);
    }
    double length () {//求线段长度
        return dis (s, e);
    }
    double angle () {//直线的倾斜角
        double k = atan2 (e.y-s.y, e.x-s.x);
        if (dcmp (k) < 0) k += pi;
        if (dcmp (k-pi) == 0) k -= pi;
        return k;
    }
};

double point_to_line (Point p, Line a) {//点到直线的距离
    return fabs (cross (p-a.s, a.e-a.s) / a.length ());
}

double point_to_seg (Point p, Line a) {//点到线段的距离
    if (dcmp (dot (p-a.s, a.e-a.s)) < 0 || dcmp (dot (p-a.e, a.s-a.e)) < 0)
        return min (dis (p, a.e), dis (p, a.s));
    return point_to_line (p, a);
}

Point projection (Point p, Line a) {//点在直线上的投影
    return a.s + (((a.e-a.s) * dot (a.e-a.s, p-a.s)) / (a.e-a.s).len2() );
}
//***************圆
struct Circle {
    //圆心 半径
    Point p;
    double r;
    Circle () {}
    Circle (Point _p, double _r) : p(_p), r(_r) {}
    Circle (double a, double b, double _r) {
        p = Point (a, b);
        r = _r;
    }
    void input () {
        p.input ();
        scanf ("%lf", &r);
    }
    void output () {
        p.output ();
        printf (" %.2f\n", r);
    }
    bool operator == (const Circle &a) const {
        return p == a.p && (dcmp (r-a.r) == 0);
    }
    double area () {//面积
        return pi*r*r;
    }
    double circumference () {//周长
        return 2*pi*r;
    }
    bool operator < (const Circle &a) const {
        return p < a.p || (p == a.p && r < a.r);
    }
};

int relation (Point p, Circle a) {//点和圆的关系
    //0:圆外 1:圆上 2:圆内
    double d = dis (p, a.p);
    if (dcmp (d-a.r) == 0) return 1;
    return (dcmp (d-a.r) < 0 ? 2 : 0);
}

int relation (Line a, Circle b) {//直线和圆的关系
    //0:相离 1:相切 2:相交
    double p = point_to_line (b.p, a);
    if (dcmp (p-b.r) == 0) return 1;
    return (dcmp (p-b.r) < 0 ? 2 : 0);
}

int relation (Circle a, Circle v) {//两圆的位置关系
    //1:内含 2:内切 3:相交 4:外切 5:相离
    double d = dis (a.p, v.p);
    if (dcmp (d-a.r-v.r) > 0) return 5;
    if (dcmp (d-a.r-v.r) == 0) return 4;
    double l = fabs (a.r-v.r);
    if (dcmp (d-a.r-v.r) < 0 && dcmp (d-l) > 0) return 3;
    if (dcmp (d-l) == 0) return 2;
    if (dcmp (d-l) < 0) return 1;
}

int circle_intersection (Circle a, Circle v, Point &p1, Point &p2) {//两个圆的交点
    //返回交点个数 交点保存在引用中
    int rel = relation (a, v);
    if (rel == 1 || rel == 5) return 0;
    double d = dis (a.p, v.p);
    double l = (d*d + a.r*a.r - v.r*v.r) / (2*d);
    double h = sqrt (a.r*a.
                     r - l*l);
    Point tmp = a.p + (v.p-a.p).change_len (l);
    p1 = tmp + ((v.p-a.p).rotate_left ().change_len (h));
    p2 = tmp + ((v.p-a.p).rotate_right ().change_len (h));
    if (rel == 2 || rel == 4) return 1;
    return 2;
}

int line_cirlce_intersection (Line v, Circle u, Point &p1, Point &p2) {//直线和圆的交点
    //返回交点个数 交点保存在引用中
    if (!relation (v, u)) return 0;
    Point a = projection (u.p, v);
    double d = point_to_line (u.p, v);
    d = sqrt (u.r*u.r - d*d);
    if (dcmp (d) == 0) {
        p1 = a, p2 = a;
        return 1;
    }
    p1 = a + (v.e-v.s).change_len (d);
    p2 = a - (v.e-v.s).change_len (d);
    return 2;
}

double circle_traingle_area (Point a, Point b, Circle c) {//圆心三角形的面积
    //a.output (), b.output (), c.output ();
    Point p = c.p; double r = c.r; //cout << cross (p-a, p-b) << endl;
    if (dcmp (cross (p-a, p-b)) == 0) return 0;
    Point q[5];
    int len = 0;
    q[len++] = a;
    Line l(a, b);
    Point p1, p2;
    if (line_cirlce_intersection (l, c, q[1], q[2]) == 2) {
        if (dcmp (dot (a-q[1], b-q[1])) < 0) q[len++] = q[1];
        if (dcmp (dot (a-q[2], b-q[2])) < 0) q[len++] = q[2];
    }
    q[len++] = b;
    if (len == 4 && dcmp (dot (q[0]-q[1], q[2]-q[1])) > 0)
        swap (q[1], q[2]);
    double res = 0;
    for (int i = 0; i < len-1; i++) {
        if (relation (q[i], c) == 0 || relation (q[i+1], c) == 0) {
            double arg = rad (q[i]-p, q[i+1]-p);
            res += r*r*arg/2.0;
        }
        else {
            res += fabs (cross (q[i]-p, q[i+1]-p))/2;
        }
    }
    return res;
}
double area_polygon_circle (Circle c, Point *p, int n) {//多边形和圆交面积
    double ans = 0;
    for (int i = 0; i < n; i++) {
        int j = (i+1)%n;
        if (dcmp (cross (p[j]-c.p, p[i]-c.p)) >= 0)
            ans += circle_traingle_area (p[i], p[j], c);
        else
            ans -= circle_traingle_area (p[i], p[j], c);
    }
    return fabs (ans);
}

#define maxn 511
Point p[maxn];
int n;
double R, k;
Point A, B;

void solve () {
    double aa, bb, pp, qq;
    double gg = 1-k*k;
    aa = (2*B.x-2*A.x*k*k)/gg;
    bb = (B.x*B.x-k*k*A.x*A.x)/gg;
    pp = (2*B.y-2*A.y*k*k)/gg;
    qq = (B.y*B.y-k*k*A.y*A.y)/gg;
    Point O (aa/2, pp/2);
    double R = sqrt (-bb+aa*aa/4 - qq+pp*pp/4);
    Circle C (O, R);
    double ans = area_polygon_circle (C, p, n);
    printf("%.10f\n",ans);
}

int main() {
    int kase= 0 ;
    while(scanf ("%d%lf", &n, &k) == 2) {
        printf ("Case %d: ", ++kase);
        for (int i = 0; i < n; i++) {
            p[i].input ();
        }
        A.input (), B.input ();
        solve ();
    }
    return 0;
}