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

Question

主题link：https://codeforces.com/problemset/problem/600/D

对于这个问题，我在 test28 上的答案是错误的，它可能看起来像这样：

正确 answer:119256.95877838134765625000

我的回答：120502.639190673828125

我猜是计算精度造成的，但我没有证据。可能算法本身有问题，请指出。

算法思路：

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

开头第一个圆心的坐标： $(x_1, y_1)$ .

开头第二个圆心的坐标： $(x_2, y_2)$ .

一系列变换后第一个圆的圆心坐标： $(0, 0)$ .

一系列变换后的第二个圆的圆心坐标： $(c, 0)$ ,

$c = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}$ .

两个圆的方程合并：

$\left\{\begin{matrix} x^2 + y^2 = r_1^2 \ (x-c)^2 + y^2 = r_2^2 \end{matrix}\right.$

$r_1, r_2$ 分别是第一个和第二个圆的半径，所以：

$x = \frac{r_1^2 - r_2^2 + c ^ 2}{2c}$

$y = \pm \sqrt{r_1^2 - x^2}$

我们可以使用扇区面积公式： $area = \frac{1}{2} angle * r^2$ ,

$angle1 = 2.0 * arcos(x / r1)$ , $angle2 = 2.0 * arcos((c - x) / r2)$ .

这个地方会出现弧度正负值的问题（下图两张），但可以证明在最后的结果中是可以吸收的

最后的结果是两条圆弧的面积之和减去中间菱形的面积

我的代码：

# 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;
}

Answer 1

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

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
}

Answer 2

我使用 Wolfram Alpha 检查了你的公式，所以我认为你看到 catastrophic cancellation rather than an accumulation of small errors. I'd use William Kahan's formula 来计算三角形面积，如下所示：

#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";
  }
}

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

How to calculate the area of two circles' intersection?

c++

algorithm

geometry