Powered by GitBook

圆

圆与线求交

算法思路

将线段 $AB$ 写成参数方程 $P=A+t(B-A)$ ，带入圆的方程，得到一个一元二次方程。解出 $t$ 就可以求得线段所在直线与圆的交点。如果 $0 \leq t \leq 1$ 则说明点在线段上。

模板代码

void circle_cross_line(Pt a, Pt b, Pt o, double r, Pt ret[], int &num) {
    double ox = o.x, oy = o.y, ax = a.x, ay = a.y, bx = b.x, by = b.y;
    double dx = bx-ax, dy = by-ay;
    double A = dx*dx + dy*dy;
    double B = 2*dx*(ax-ox) + 2*dy*(ay-oy);
    double C = sqr(ax-ox) + sqr(ay-oy) - sqr(r);
    double delta = B*B - 4*A*C;
    num = 0;
    if (sgn(delta) >= 0) {
        double t1 = (-B - Sqrt(delta)) / (2*A);
        double t2 = (-B + Sqrt(delta)) / (2*A);
        if (sgn(t1-1) <= 0 && sgn(t1) >= 0)
            ret[num++] = Pt(ax + t1*dx, ay + t1*dy);
        if (sgn(t2-1) <= 0 && sgn(t2) >= 0)
            ret[num++] = Pt(ax + t2*dx, ay + t2*dy);
    }
}

经典题目

POJ 1263 Reflections

圆与圆求交

// 计算圆a和圆b的交点，注意要先判断两圆相交
void circle_circle_cross(Pt ap, double ar, Pt bp, double br, Pt p[]) {
    double d = (ap - bp).norm();
    double cost = (ar*ar + d*d - br*br) / (2*ar*d);
    double sint = sqrt(1.0 - cost*cost);
    Pt v = (bp - ap) / (bp - ap).norm() * ar;
    p[0] = ap + rotate(v, cost, -sint);
    p[1] = ap + rotate(v, cost, sint);
}

经典题目

POJ 2546 Circular Area

圆与多边形交

算法思路

按照圆心为中心，将多边形三角剖分，并计算出每个三角形与圆交的面积后求和。求三角形与圆交面积的方法要按照情况讨论。

$AB$ 都在圆内，计算三角形 $OAB$ 的面积。
$A$ 在圆内 $B$ 不在圆内，这时计算一个三角形的面积和一个扇形的面积。
$B$ 在圆内 $A$ 不在圆内，这时计算一个三角形的面积和一个扇形的面积。
$AB$ 都不在圆内，如果 $AB$ 与圆无交点，则计算一个扇形的面积，否则计算两个扇形和一个三角形的面积。

模板代码

// 圆与直线相交
void circle_cross_line(Pt a, Pt b, Pt o, double r, Pt ret[], int &num) {
    double ox = o.x, oy = o.y, ax = a.x, ay = a.y, bx = b.x, by = b.y;
    double dx = bx-ax, dy = by-ay;
    double A = dx*dx + dy*dy;
    double B = 2*dx*(ax-ox) + 2*dy*(ay-oy);
    double C = sqr(ax-ox) + sqr(ay-oy) - sqr(r);
    double delta = B*B - 4*A*C;
    num = 0;
    if (sgn(delta) >= 0) {
        double t1 = (-B - Sqrt(delta)) / (2*A);
        double t2 = (-B + Sqrt(delta)) / (2*A);
        if (sgn(t1-1) <= 0 && sgn(t1) >= 0)
            ret[num++] = Pt(ax + t1*dx, ay + t1*dy);
        if (sgn(t2-1) <= 0 && sgn(t2) >= 0)
            ret[num++] = Pt(ax + t2*dx, ay + t2*dy);
    }
}

// 计算扇形面积
double sector_area(Pt a, Pt b, double r) {
    double theta = atan2(a.y, a.x) - atan2(b.y, b.x);
    while (theta <= 0) theta += 2*PI;
    while (theta > 2*PI) theta -= 2*PI;
    theta = min(theta, 2*PI - theta);
    return r*r*theta / 2;
}

double area(Pt res[], int n, double r) {
    double ans = 0.0;
    for (int i = 0; i < n; ++i) {
        Pt a = res[i], b = res[next(i)];
        Pt p[2];
        int num = 0;
        int ina = sgn(a.norm() - r) < 0;
        int inb = sgn(b.norm() - r) < 0;
        int s = sgn(det(a, b));

        double part = 0.0;

        if (ina) {
            if (inb) {
                part = Abs(det(a, b)) / 2.0;
            }else{
                circle_cross_line(a, b, Pt(0, 0), r, p, num);
                part = sector_area(b, p[0], r) + fabs(det(a, p[0])) / 2.0;
            }
        }else{
            if (inb) {
                circle_cross_line(a, b, Pt(0, 0), r, p, num);
                part = sector_area(p[0], a, r) + fabs(det(p[0], b)) / 2.0;
            }else{
                circle_cross_line(a, b, Pt(0, 0), r, p, num);
                if (num == 2) {
                    part = sector_area(a, p[0], r) + sector_area(p[1], b, r)
                        + fabs(det(p[0], p[1])) / 2.0;
                }else{
                    part = sector_area(a, b, r);
                }
            }
        }
        ans += s * part;
    }
    return ans;
}

const int MaxN = 55;
Pt res[MaxN];
double r;
int n;

int main() {
    while (cin >> r >> n) {
        for (int i = 0; i < n; ++i)
            cin >> res[i].x >> res[i].y;
        double ans = Abs(area(res, n, r));
        cout << setiosflags(ios::fixed) << setprecision(2) << ans << endl;
    }

    return 0;
}

results matching ""

No results matching ""