LA 7072 Signal Interference 计算几何 圆与多边形的交

题意：

给出一个\(n\)个点的简单多边形，和两个点\(A, B\)还有一个常数\(k(0.2 \leq k < 0.8)\)。

点\(P\)满足\(\left | PB \right | \leq k \left | PA \right |\)，求点\(P\)构成的图形与多边形相交的面积。

分析：

推导圆的公式：

高中应该做过这样的题目，我们很容易知道这是一个圆。

下面推导一下圆的方程：

设\(A(x_1, y_1),B(x_2,y_2),P(x,y)\)，根据条件列出等式：

\(\sqrt{(x-x_2)^2+(y-y_2)^2}=k\sqrt{(x-x_1)^2+(y-y_1)^2}\)

\(\Rightarrow (x-x_2)^2+(y-y_2)^2=k^2 \left [ (x-x_1)^2+(y-y_1)^2 \right ]\)

展开整理得：

\((1-k^2)x^2-2(x_2-k^2x_1)x+(x_2^2-k^2x_1^2)+(1-k^2)y^2-2(y_2-k^2y_1)y+(y_2^2-k^2y_1^2)=0\)

\(\Rightarrow x^2+y^2-2\frac{x_2-k^2x_1}{1-k^2}x-2\frac{y_2-k^2y_1}{1-k^2}y+\frac{x_2^2+y_2^2-k^2(x_1^2+y_1^2)}{1-k^2}=0\)

对于圆的一般式\(x^2+y^2+Dx+Ey+F=0\)，可以用配方法知道圆心为\((-\frac{D}{2},-\frac{E}{2})\)，半径为\(\frac{\sqrt{D^2+E^2-4F}}{2}\)

计算多边形与圆的面积交

根据以往的经验，多边形的问题都是要分解成若干三角形的。

对于多边形的每条边，我们可以先求该线段与圆心构成三角形与圆的面积交。

会有下面几种情况：

两个点都在圆的内部，此时相交的面积为三角形的面积。



一个点在圆内，一个点在圆外，相交的面积可以分解为一个扇形和一个三角形的面积



两个点都在圆外，而且线段和圆没有交点或者相切，相交的面积为一个扇形



两个点都在圆外，而且线段和圆有两个交点，相交的面积可以分解为两个扇形和一个三角形



我们求出每个三角形与圆的面积交以后，根据叉积的正负来计算面积的面积的正负，最后再取个绝对值就是整个多边形与圆的面积交。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <map>
#include <cmath>
//#define DEBUG
using namespace std;

const double eps = 1e-10;
const double PI = acos(-1.0);
const double TOW_PI = PI * 2.0;

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

struct Point
{
double x, y;

void read() { scanf("%lf%lf", &x, &y); }

void print() { printf("(%.2f %.2f)", x, y); }

Point(double x = 0, double y = 0):x(x), y(y) {}
};

typedef Point Vector;

Point operator + (const Point& A, const Point& B) {
return Point(A.x + B.x, A.y + B.y);
}

Point operator - (const Point& A, const Point& B) {
return Point(A.x - B.x, A.y - B.y);
}

Point operator * (const Point& A, double p) {
return Point(A.x * p, A.y * p);
}

Point operator / (const Point& A, double p) {
return Point(A.x / p, A.y / p);
}

bool operator < (const Point& A, const Point& B) {
return A.x < B.x || (A.x == B.x && A.y < B.y);
}

bool operator == (const Point& A, const Point& B) {
return A.x == B.x && A.y == B.y;
}

double Dot(const Vector& A, const Vector& B) {
return A.x * B.x + A.y * B.y;
}

double Cross(const Vector& A, const Vector& B) {
return A.x * B.y - A.y * B.x;
}

double TriangleArea(Point A, Point B, Point C) {
return fabs(Cross(B-A, C-A)) / 2.0;
}

double Length(const Vector& A) {
return sqrt(Dot(A, A));
}

double Angle(Vector A, Vector B) {
return acos(Dot(A, B) / Length(A) / Length(B));
}

struct Line
{
Point p;
Vector v;

Line() {}
Line(Point p, Vector v):p(p), v(v) {}

Point point(double t) { return p + v * t; }
};

struct Circle
{
Point c;
double r;
};

int getSegCircleIntersection(Line L, Circle C, vector<Point>& sol) {
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;
double t1, t2;
int ans = 0;

if(dcmp(delta) < 0) return 0;

if(dcmp(delta) == 0) {
t1 = t2 = -f / (2 * e);
if(dcmp(t1) >= 0 && dcmp(t1 - 1) <= 0) {
ans++;
sol.push_back(L.point(t1));
}
return ans;
}

t1 = (-f - sqrt(delta)) / (2 * e);
t2 = (-f + sqrt(delta)) / (2 * e);
if(t1 > t2) swap(t1, t2);
if(dcmp(t1) >= 0 && dcmp(t1 - 1) <= 0) {
ans++;
sol.push_back(L.point(t1));
}
if(dcmp(t2) >= 0 && dcmp(t2 - 1) <= 0) {
ans++;
sol.push_back(L.point(t2));
}
return ans;
}

bool inCircle(Circle C, Point p) {
return dcmp(C.r - Length(C.c - p)) >= 0;
}

double IntersectionArea(Circle C, Point A, Point B) {
Line L(A, B - A);
int cnt = 0;
bool inA, inB;
if(inA = inCircle(C, A)) cnt++;
if(inB = inCircle(C, B)) cnt++;

if(cnt == 2) return TriangleArea(C.c, A, B);

if(cnt == 1) {
vector<Point> q;
getSegCircleIntersection(L, C, q);
if(inB) swap(A, B);
double theta = Angle(q[0]-C.c, B-C.c);
return C.r*C.r*theta/2 + TriangleArea(C.c, A, q[0]);
}

vector<Point> q;
int sz = getSegCircleIntersection(L, C, q);
if(sz <= 1) {
double theta = Angle(A-C.c, B-C.c);
return C.r*C.r*theta/2;
}

double theta = Angle(A-C.c, q[0]-C.c) + Angle(B-C.c, q[1]-C.c);
return C.r*C.r*theta/2 + TriangleArea(q[0], q[1], C.c);
}

int n;
double k;
vector<Point> poly;
Point A, B;
Circle C;

int main()
{

int kase = 1;
while(scanf("%d%lf", &n, &k) == 2) {
poly.resize(n);
for(int i = 0; i < n; i++) poly[i].read();
A.read(); B.read();
poly.push_back(poly[0]);

double D = 2.0 * (k*k*A.x-B.x) / (1.0-k*k);
double E = 2.0 * (k*k*A.y-B.y) / (1.0-k*k);
double F = (B.x*B.x+B.y*B.y-k*k*(A.x*A.x+A.y*A.y)) / (1.0-k*k);

C.c = Point(-D/2.0, -E/2.0);
C.r = sqrt(D*D+E*E-4.0*F) / 2.0;

#ifdef DEBUG
printf("Circle:");
C.c.print();
printf(" Radius : %.2f\n", C.r);
#endif

double ans = 0;
for(int i = 0; i < n; i++) {
int sign;
if(dcmp(Cross(poly[i]-C.c, poly[i+1]-C.c)) > 0) sign = 1;
else sign = -1;
ans += sign * IntersectionArea(C, poly[i], poly[i+1]);
}

printf("Case %d: %.10f\n", kase++, fabs(ans));
}

return 0;
}