如何计算两个圆的交集面积？

Question

主题链接: https://codeforces.com/problemset/problem/600/D

对于这个问题，我在测试28上的答案是错误的，可能看起来像这样:

正确答案：119256.95877838134765625000

我的答案：120502.639190673828125

我猜这是由计算精度引起的，但我没有证据。也许算法本身有问题，请指出。

算法思路：

对于任何给定的两个圆，为了简化计算，我们可以将坐标原点平移到其中一个圆的中心，然后通过旋转将另一个圆旋转到x轴。为了计算圆的交集区域，在之前和之后是等价的，最后形成了一个紫色的圆和一个红色的圆。实际上，最终的交集面积等于两个扇形面积之和减去中间菱形的面积（下图，水平轴x，垂直轴y）。但在此之前，我们必须先计算两个圆的交点。

第一个圆的中心在开始时的坐标： 。

第二个圆的中心在开始时的坐标： 。

一系列变换后，第一个圆的中心坐标为： 。

一系列变换后，第二个圆的中心坐标为： 。

.

两个圆的方程被结合：

分别是第一个和第二个圆的半径，因此：

我们可以使用扇形面积公式：。 ，。

在这里，弧度值的正负会有问题（如下两个图），但可以证明它们可以被吸收在最终结果中。

最终结果是两个弧形面积之和减去中间菱形的面积。

我的代码：

# define MY_PI 3.14159265358979323846
using namespace std;

typedef long double ld;

ld pow2(ld x){return x * x;}
int main()
{
    cout.precision(25);
    ld x1, y1, r1, x2, y2, r2;   // (x1, y1) is the center point of the circle, r1 is the radius of the circle, the other is similar.
    cin >> x1 >> y1 >> r1 >> x2 >> y2 >> r2;

    ld c = sqrt(pow2(x2 - x1) + pow2(y2 - y1));
    if(r1 + r2 < c)         // when r1+r2 < c, there is no intersection
    {
        cout << 0 << endl;
        return 0;
    }
    if((max(r1, r2) - min(r1, r2)) >= c)  // one circle is inside the other circle.
    {
        cout << MY_PI * pow2(min(r1, r2)) << endl;
        return 0;
    }

    ld x = (pow2(r1) - pow2(r2)) / (2 * c) + c * 0.5;
    ld y = sqrt(pow2(r1) - pow2(x));
    ld angle1 = 2.0 * acos(x / r1);
    ld angle2 = 2.0 * acos((c - x) / r2);
    ld ar1 = angle1 * pow2(r1) * 0.5;
    ld ar2 = angle2 * pow2(r2) * 0.5;
    ld res = ar1 + ar2 - y * c;

    cout << res << endl;

    return 0;
}

- RecharBao

Answer 1

我使用Wolfram Alpha检查了您的公式，因此我认为您看到的是大消除法而不是小误差积累。要计算三角形面积，我建议使用William Kahan的公式，如下所示：

#include <algorithm>
#include <cmath>
#include <iostream>
#include <limits>
typedef long double ld;

struct Circle {
  ld x;
  ld y;
  ld r;
};

static std::istream& operator>>(std::istream& i, Circle& c) {
  return i >> c.x >> c.y >> c.r;
}

static ld Pi() { return std::acos(-1); }

static ld Square(ld x) { return x * x; }

static void SortDescending2(ld& a, ld& b) {
  if (a < b) {
    using std::swap;
    swap(a, b);
  }
}

static void SortDescending3(ld& a, ld& b, ld& c) {
  SortDescending2(a, b);
  SortDescending2(b, c);
  SortDescending2(a, b);
}

static ld KahanAreaOfTriangle(ld a, ld b, ld c) {
  SortDescending3(a, b, c);
  return 0.25 * std::sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) *
                          (a + (b - c)));
}

static ld AreaOfIntersection(const Circle& c1, const Circle& c2) {
  ld R = c1.r;
  ld r = c2.r;
  ld d = std::hypot(c2.x - c1.x, c2.y - c1.y);
  if (R + r <= d) {
    return 0;
  }
  SortDescending2(R, r);
  if (d <= R - r) {
    return Pi() * Square(r);
  }
  ld area = 2 * KahanAreaOfTriangle(R, r, d);
  ld y = area / d;
  ld A = std::asin(y / R) * Square(R);
  ld a = std::asin(y / r);
  if (std::hypot(r, d) <= R) {
    a = Pi() - a;
  }
  a *= Square(r);
  SortDescending2(A, a);
  return (A - area) + a;
}

int main() {
  Circle c1;
  Circle c2;
  if (std::cin >> c1 >> c2) {
    std::cout.precision(std::numeric_limits<ld>::max_digits10);
    std::cout << AreaOfIntersection(c1, c2) << "\n";
  }
}

- David Eisenstat

---

计算三角形面积的方法是在大约2000年前发现的，被称为Heron公式。 - user3386109

Answer 2

不要使用太多中间浮点变量来得到最终答案。当这些小的不准确性相加时，会导致巨大的不准确性，这在你的问题的预期输出和实际输出中可以清楚地看到。你可以做的是尽量减少这些中间浮点变量。例如

ld c = sqrt(pow2(x2 - x1) + pow2(y2 - y1));
if(r1 + r2 < c)         // when r1+r2 < c, there is no intersection
{
    cout << 0 << endl;
    return 0;
}

您可以通过以下方式消除 c 作为浮点变量的情况

int c = pow(x2-x1) + pow(y2-y1);
if (pow(r1+r2) > pow(c)) {
   // do your stuff
}

- dottomato

---

我注意到这个问题，但如何减少中间浮点变量呢？虽然你的if语句没有问题，但当你计算x和y时，必须有一些浮点中间变量，并且很难减少变量，无论它是左值还是右值。 - RecharBao