计算多边形与圆的面积交

基本思路：将多边形的面积拆成三角形有向面积之和，然后分别求每个小三角形与圆的面积交

求三角形与圆的面积交

积分法

• 梯形积分
• 辛普森积分

不管是什么积分大概三角形的个数一多就不行了…

强搞法

一开始的想法：

三角形和圆的面积交肯定可以分解为一个凸多边形+一些弓形的面积和

凸多边形的顶点可以通过算三角形有哪些顶点在圆内及圆上 + 三角形的边与圆的交点得出

如果顶点个数为1：
1.圆与三角形不相交
2.相交部分为整个圆

如果顶点个数为2：
1.相交部分是优弧所对弓形
2.相交部分是劣弧所对弓形
3.相交部分为整个圆

如果顶点个数为3:
首先将这些顶点进行极角排序，凸多边形的面积可以轻易的算出

分别枚举这些顶点，如果data[i], data[i+1]所组成的线段在三角形内部，则这里肯定有一个弓形需要计算进来

如上图的话有凸多边形有两条边在三角形内，计算一下这两条边所对弓形就好了。

然后我天真地以为如果凸多边形有边在三角形内部的话，其所对的弓形一定是劣弧所在的弓形！！！！之后我还傻乎乎的把整个算法都实现出来了…

直到….我发现了如下的数据…

这尼玛的!!!想要判断的话也不是不可以…只是好麻烦啊啊啊

进化了之后的想法：

基本参考http://www.cnblogs.com/lxglbk/archive/2012/08/12/2634192.html

注意其中的第二种情况：

另外两个顶点有1个在圆内（上）,另外1个在圆外，求得直线与圆一个交点后求一个三角形面积+上一个扇形面积。>

这段话不完全对，也不能够说是错的。比如

这里的话相交面积就是一个扇形，这是因为直线与圆的一个交点与三角形顶点重合，使得三角形面积为0，实现的时候注意一下这里。

他的方法是针对于三角形中一个顶点有在圆心，对于一般的三角形，可以把这个三角形在划分成3个其中一个顶点在圆心的三角形求解。

情况无外乎3种：

3个三角形全为正

1正2负

2正1负

我是分这三种情况讨论的，可能还有别的方法？

100204F - Little Mammoth

#include <cmath>#include <vector>#include <cstdio>#include <cstring>#include <iostream>#define x first#define y second#define mk make_pair#define pb push_backusing namespace std;typedef pair<double, double>II;const int SIZE = 1010;const double pi = 3.14159265357;const double Eps = 1e-6;struct Line {    II st, ed;};II data[SIZE];struct Circle {    II p0;    double R;}cir;int N;double cross(II p0, II p1, II p2) {    return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);}Line mk_Line(II st, II ed) {    Line now;    now.st = st;    now.ed = ed;    return now;}int pd(double x) {    if (x > Eps)        return 1;    if (x < -Eps)       return -1;    return 0;     }II operator +(II A, II B){return mk(A.x + B.x, A.y + B.y);}II operator -(II A, II B){return mk(A.x - B.x, A.y - B.y);}II operator *(II A, double k){return mk(A.x * k, A.y * k);}II operator -(II A){return mk(-A.x, -A.y);}double mymin(double x, double y){return x < y ? x : y;}double mymax(double x, double y){return x > y ? x : y;}double calc_length(II p0) {    return sqrt(p0.x * p0.x + p0.y * p0.y);}double FABS(double x){return x > 0 ? x : -x;}bool in_segment(Line it, II it1) {    if (pd(cross(it.st, it.ed, it1)) != 0)      return false;    if (pd(mymin(it.st.x, it.ed.x) - it1.x) <= 0 && pd(it1.x - mymax(it.st.x, it.ed.x)) <= 0 &&        pd(mymin(it.st.y, it.ed.y) - it1.y) <= 0 && pd(it1.y - mymax(it.st.y, it.ed.y)) <= 0)   return true;    return false;}vector<II> calc_segment_circle_point(Line it, Circle it2) {    double x1 = it.st.x, y1 = it.st.y;    double x2 = it.ed.x, y2 = it.ed.y;    double x0 = it2.p0.x, y0 = it2.p0.y;    double b = 2 * ((x2 - x1) * (x1 - x0) + (y1 - y0) * (y2 - y1));    double a = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);    double c = (x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0) - it2.R * it2.R;    double delta = b * b - 4 * a * c;    vector<II>res;          res.clear();    if (pd(delta) < 0)      return res;    double k1 = (-b - (sqrt(delta + Eps))) / (2 * a);    II jiji;    jiji = it.st + (it.ed - it.st) * k1;    if (in_segment(it, jiji))       res.pb(jiji);    if (pd(delta) == 0)     return res;    double k2 = (-b + (sqrt(delta + Eps))) / (2 * a);    jiji = it.st + (it.ed - it.st) * k2;    if (in_segment(it, jiji))       res.pb(jiji);    return res;}double calc_pai_circle(II p1, II p2, Circle it) {    double L = calc_length(p2 - p1);    double ang = acos((2 * it.R * it.R - L * L) / (2 * it.R * it.R));    double area1 = pi * it.R * it.R * (ang / pi / 2);    //double area2 = L * sqrt((calc_length(p1 - it.p0) * calc_length(p1 - it.p0)) - L * L / 4) / 2;    //printf("%lf\n", area1);    return area1;}bool operator!=(II A, II B) {    if (pd(A.x - B.x) == 0 && pd(A.y - B.y) == 0)   return false;    return true;}double calc(Line it1, Circle it2) {    if (pd(cross(it1.st, it1.ed, it2.p0)) == 0)      return 0.0;    if (pd(calc_length(it1.ed - it2.p0) - it2.R) <= 0)        swap(it1.st, it1.ed);    if (pd(calc_length(it1.ed - it2.p0) - it2.R) <= 0)        return FABS(cross(it2.p0, it1.st, it1.ed)) / 2;    if (pd(calc_length(it1.st - it2.p0) - it2.R) <= 0) {        vector<II>res = calc_segment_circle_point(it1, it2);        Line it3;   it3.st = it2.p0;    it3.ed = it1.ed;        II p1 = it1.st;        for (int i = 0; i < (int)res.size(); i++) {            if (res[i] != it1.st && res[i] != it1.ed)                p1 = res[i];        }        res = calc_segment_circle_point(it3, it2);        II p2;        for (int i = 0; i < (int)res.size(); i++) {            if (res[i] != it1.ed && res[i] != it2.p0)                p2 = res[i];        }        return FABS(cross(it1.st, p1, it2.p0)) / 2 +                calc_pai_circle(p1, p2, it2);    }    vector<II>res = calc_segment_circle_point(it1, it2);    Line it3;   it3.st = it1.st;    it3.ed = it2.p0;    Line it4;   it4.st = it1.ed;    it4.ed = it2.p0;    vector<II>res1 = calc_segment_circle_point(it3, it2);    vector<II>res2 = calc_segment_circle_point(it4, it2);    if ((int)res.size() < 2)        return calc_pai_circle(res1[0], res2[0], it2);    return calc_pai_circle(res1[0], res[0], it2) +           FABS(cross(res[0], res[1], it2.p0)) / 2 +           calc_pai_circle(res[1], res2[0], it2);}II a[6];int sig[6];double calc_area(II p0, II p1, II p2, Circle it) {    a[1] = p0;  a[2] = p1;  a[3] = p2;  a[4] = p0;  a[5] = p1;    int ope = -1;    for (int i = 1; i <= 3; i++) {        if (pd(cross(a[i], it.p0, a[i + 1]) * cross(a[i], it.p0, a[i + 2])) <= 0 &&            pd(cross(a[i + 1], a[i + 2], a[i]) * cross(a[i + 1], a[i + 2], it.p0)) <= 0)            ope = i;    }    if (ope == 3)        ope = 0;    for (int i = 1; i <= 3; i++)        sig[i] = 1;    if (ope == -1) {        double maxarea = 0.0;        for (int i = 1; i <= 3; i++) {            if (FABS(cross(a[i], a[i + 1], it.p0) * 0.5) > maxarea) {                maxarea = FABS(cross(a[i], a[i + 1], it.p0) * 0.5);                ope = i;            }        }        if (pd(maxarea - FABS(cross(a[1], a[2], a[3]) * 0.5)) >= 0) {            for (int i = 1; i <= 3; i++)                sig[i] = -1;            sig[ope] = 1;        }    }    else    sig[ope + 1] = -1;    double res = 0.0;    for (int i = 1; i <= 3; i++) {        res += sig[i] * calc(mk_Line(a[i], a[i + 1]), it);        //printf("%lf\n", sig[i] * calc(mk_Line(a[i], a[i + 1]), it));    }    return res;}void Read() {    scanf("%d", &N);    for (int i = 1; i <= N; i++) {        scanf("%lf%lf", &data[i].x, &data[i].y);    }}double Solve() {    N = 4;    double ans = 0.0;    for (int i = 2; i <= N - 1; i++) {        ans += calc_area(data[1], data[i], data[i + 1], cir);    }    return FABS(ans);}int main() {    freopen("mammoth.in", "r", stdin);    freopen("mammoth.out", "w", stdout);    scanf("%lf%lf", &cir.p0.x, &cir.p0.y);    scanf("%lf", &cir.R);    scanf("%lf%lf", &data[1].x, &data[1].y);        scanf("%lf%lf", &data[3].x, &data[3].y);    data[2].x = data[1].x;  data[2].y = data[3].y;    data[4].x = data[3].x;  data[4].y = data[1].y;    //Read();    printf("%.6f\n", Solve() + Eps);    return 0;}