圆形和矩形之间的交点区域

给定 2 个交叉点：

0 个顶点位于圆内：圆形线段的面积

    XXXXX              -------------------
   X     X               X            X Circular segment
  X       X               XX        XX 
+-X-------X--+              XXXXXXXX 
|  X     X   |
|   XXXXX    |

1 个顶点位于圆内：圆形线段和三角形的面积之和。

    XXXXX                   XXXXXXXXX
   X     X       Triangle ->X     _-X
  X       X                 X   _-  X 
  X    +--X--+              X _-   X <- Circular segment 
   X   | X   |              X-  XXX 
    XXXXX    |              XXXX
       |     | 

圆内有 2 个顶点：两个三角形和一条圆形线段的面积之和

    XXXXX                   +------------X
   X     X                  |      _--'/'X 
  X    +--X---    Triangle->|   _--   / X
  X    |  X                 |_--     /XX <- Circular segment
   X   +-X----              +-------XX
    XXXXX                 Triangle^

圆内有 3 个顶点：矩形的面积减去三角形的面积加上圆形线段的面积

    XXXXX
   X  +--X+             XXX
  X   |   X         -------XXX-----+ <- Triangle outside
 X    |   |X        Rect ''.  XXX  |
 X    +---+X                ''.  XX|  
 X         X                   ''. X <- Circular segment inside 
  X       X                       ^|X 
   X     X                         | X 
    XXXXX

要计算这些面积，请执行以下操作：

• 您需要使用的大多数点都可以通过查找直线和圆的交点来找到

• 通过计算圆形线段的面积并计算三角形的面积，可以找到需要计算的面积。

• 可以通过计算顶点与中心的距离是否小于半径来确定顶点是否在圆内。

我意识到这个问题在一段时间前得到了回答，但我正在解决同样的问题，我找不到一个可以使用的开箱即用的解决方案。请注意，我的框是轴对齐的，OP没有完全指定这一点。下面的解决方案是完全通用的，适用于任意数量的交叉路口（不仅仅是两个）。请注意，如果框不是轴对齐的（但仍然是具有直角的框，而不是常规四边形），则可以利用圆是圆形的，旋转所有内容的坐标，以使框最终成为轴对齐的，然后使用此代码。

我想使用集成 - 这似乎是一个好主意。让我们从编写一个明显的公式来绘制圆圈开始：

x = center.x + cos(theta) * radius
y = center.y + sin(theta) * radius

^
|
|**###        **
| #*  #      *          * x
|#  *  #    *           # y
|#  *  #    *     
+-----------------------> theta
     *  #  *  # 
     *  #  *  #
      *  #*  #
       ***###

这很好，但我无法将该圆的面积整合到 或 ;这些是不同的数量。我只能整合角度，产生比萨饼片的区域。不是我想要的。让我们尝试更改参数：xytheta

(x - center.x) / radius = cos(theta) // the 1st equation
theta = acos((x - center.x) / radius) // value of theta from the 1st equation
y = center.y + sin(acos((x - center.x) / radius)) * radius // substitute to the 2nd equation

那更像是这样。现在给定的范围，我可以积分，得到圆的上半部分的面积。这只适用于 in（其他值将导致虚数输出），但我们知道该范围之外的区域为零，因此可以轻松处理。为了简单起见，让我们假设单位圆，我们以后可以随时将中心和半径插回去：xyx[center.x - radius, center.x + radius]

y = sin(acos(x)) // x in [-1, 1]
y = sqrt(1 - x * x) // the same thing, arguably faster to compute http://www.wolframalpha.com/input/?i=sin%28acos%28x%29%29+

               ^ y
               |
            ***|***     <- 1
        ****   |   ****
      **       |       **
     *         |         *
    *          |          *
----|----------+----------|-----> x
   -1                     1

这个函数确实有一个积分，因为它是单位圆的上半部分（半圆的面积是，我们已经说过单位，这意味着）。现在，我们可以通过积分来计算站在x轴上的半圆和无限高的盒子的交点面积（圆的中心也位于x轴上）：pi/2pi * r^2 / 2r = 1y

f(x): integral(sqrt(1 - x * x) * dx) = (sqrt(1 - x * x) * x + asin(x)) / 2 + C // http://www.wolframalpha.com/input/?i=sqrt%281+-+x*x%29
area(x0, x1) = f(max(-1, min(1, x1))) - f(max(-1, min(1, x0))) // the integral is not defined outside [-1, 1] but we want it to be zero out there

        ~         ~
        |      ^  |
        |      |  |
        |   ***|***     <- 1
        ****###|##|****
      **|######|##|    **
     *  |######|##|      *
    *   |######|##|       *
----|---|------+--|-------|-----> x
   -1   x0        x1      1

这可能不是很有用，因为无限高的盒子不是我们想要的。我们需要再添加一个参数，以便能够释放无限高的框的底部边缘：

g(x, h): integral((sqrt(1 - x * x) - h) * dx) = (sqrt(1 - x * x) * x + asin(x) - 2 * h * x) / 2 + C // http://www.wolframalpha.com/input/?i=sqrt%281+-+x*x%29+-+h
area(x0, x1, h) = g(min(section(h), max(-section(h), x1))) - g(min(section(h), max(-section(h), x0)))

        ~         ~
        |      ^  |
        |      |  |
        |   ***|***     <- 1
        ****###|##|****
      **|######|##|    **
     *  +------+--+      *   <- h
    *          |          *
----|---|------+--|-------|-----> x
   -1   x0        x1      1

其中，无限框的底部边缘与 x 轴的（正）距离。该函数计算单位圆与水平线相交的（正）位置，我们可以通过求解来定义它：hsectiony = h

sqrt(1 - x * x) = h // http://www.wolframalpha.com/input/?i=sqrt%281+-+x+*+x%29+%3D+h
section(h): (h < 1)? sqrt(1 - h * h) : 0 // if h is 1 or above, then the section is an empty interval and we want the area integral to be zero

               ^ y
               |  
            ***|***     <- 1
        ****   |   ****  
      **       |       **
-----*---------+---------*------- y = h
    *          |          *
----||---------+---------||-----> x
   -1|                   |1
-section(h)          section(h)

现在我们可以把事情做好了。那么如何计算一个有限盒子与 x 轴上方的单位圆相交的相交面积：

area(x0, x1, y0, y1) = area(x0, x1, y0) - area(x0, x1, y1) // where x0 <= x1 and y0 <= y1

        ~         ~                              ~         ~
        |      ^  |                              |      ^  |
        |      |  |                              |      |  |
        |   ***|***                              |   ***|*** 
        ****###|##|****                          ****---+--+****      <- y1
      **|######|##|    **                      **       |       **
     *  +------+--+      *   <- y0            *         |         *
    *          |          *                  *          |          *
----|---|------+--|-------|-----> x      ----|---|------+--|-------|-----> x
        x0        x1                             x0        x1


               ^
               |
            ***|***
        ****---+--+****      <- y1
      **|######|##|    **
     *  +------+--+      *   <- y0
    *          |          *
----|---|------+--|-------|-----> x
        x0        x1

那很好。那么，一个不在x轴上方的盒子呢？我会说不是所有的盒子都是。出现三种简单情况：

• 框位于 x 轴上方（使用上述等式）
• 框位于 x 轴下方（翻转 y 坐标的符号并使用上述等式）
• 盒子与x轴相交（将盒子分成上半部分和下半部分，使用上面的计算两者的面积并求和）

够简单吗？让我们编写一些代码：

float section(float h, float r = 1) // returns the positive root of intersection of line y = h with circle centered at the origin and radius r
{
    assert(r >= 0); // assume r is positive, leads to some simplifications in the formula below (can factor out r from the square root)
    return (h < r)? sqrt(r * r - h * h) : 0; // http://www.wolframalpha.com/input/?i=r+*+sin%28acos%28x+%2F+r%29%29+%3D+h
}

float g(float x, float h, float r = 1) // indefinite integral of circle segment
{
    return .5f * (sqrt(1 - x * x / (r * r)) * x * r + r * r * asin(x / r) - 2 * h * x); // http://www.wolframalpha.com/input/?i=r+*+sin%28acos%28x+%2F+r%29%29+-+h
}

float area(float x0, float x1, float h, float r) // area of intersection of an infinitely tall box with left edge at x0, right edge at x1, bottom edge at h and top edge at infinity, with circle centered at the origin with radius r
{
    if(x0 > x1)
        std::swap(x0, x1); // this must be sorted otherwise we get negative area
    float s = section(h, r);
    return g(max(-s, min(s, x1)), h, r) - g(max(-s, min(s, x0)), h, r); // integrate the area
}

float area(float x0, float x1, float y0, float y1, float r) // area of the intersection of a finite box with a circle centered at the origin with radius r
{
    if(y0 > y1)
        std::swap(y0, y1); // this will simplify the reasoning
    if(y0 < 0) {
        if(y1 < 0)
            return area(x0, x1, -y0, -y1, r); // the box is completely under, just flip it above and try again
        else
            return area(x0, x1, 0, -y0, r) + area(x0, x1, 0, y1, r); // the box is both above and below, divide it to two boxes and go again
    } else {
        assert(y1 >= 0); // y0 >= 0, which means that y1 >= 0 also (y1 >= y0) because of the swap at the beginning
        return area(x0, x1, y0, r) - area(x0, x1, y1, r); // area of the lower box minus area of the higher box
    }
}

float area(float x0, float x1, float y0, float y1, float cx, float cy, float r) // area of the intersection of a general box with a general circle
{
    x0 -= cx; x1 -= cx;
    y0 -= cy; y1 -= cy;
    // get rid of the circle center

    return area(x0, x1, y0, y1, r);
}

以及一些基本的单元测试：

printf("%f\n", area(-10, 10, -10, 10, 0, 0, 1)); // unit circle completely inside a huge box, area of intersection is pi
printf("%f\n", area(-10, 0, -10, 10, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2
printf("%f\n", area(0, 10, -10, 10, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2
printf("%f\n", area(-10, 10, -10, 0, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2
printf("%f\n", area(-10, 10, 0, 10, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2
printf("%f\n", area(0, 1, 0, 1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4
printf("%f\n", area(0, -1, 0, 1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4
printf("%f\n", area(0, -1, 0, -1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4
printf("%f\n", area(0, 1, 0, -1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4
printf("%f\n", area(-.5f, .5f, -.5f, .5f, 0, 0, 10)); // unit box completely inside a huge circle, area of intersection is 1
printf("%f\n", area(-20, -10, -10, 10, 0, 0, 1)); // huge box completely outside a circle (left), area of intersection is 0
printf("%f\n", area(10, 20, -10, 10, 0, 0, 1)); // huge box completely outside a circle (right), area of intersection is 0
printf("%f\n", area(-10, 10, -20, -10, 0, 0, 1)); // huge box completely outside a circle (below), area of intersection is 0
printf("%f\n", area(-10, 10, 10, 20, 0, 0, 1)); // huge box completely outside a circle (above), area of intersection is 0

其输出为：

3.141593
1.570796
1.570796
1.570796
1.570796
0.785398
0.785398
0.785398
0.785398
1.000000
-0.000000
0.000000
0.000000
0.000000

这在我看来是正确的。如果您不相信编译器会为您内联函数，则可能需要内联函数。

对于不与 x 轴相交的框，这将使用 6 sqrt、4 asin、8 div、16 mul 和 17 个 add，对于不与 x 轴相交的框，则使用该值的两倍（以及 1 个额外的添加）。请注意，除法是按半径划分的，可以减少到两个除法和一堆乘法。如果有问题的框与 x 轴相交，但未与 y 轴相交，则可以旋转所有内容，然后按原始成本进行计算。pi/2

我使用它作为调试亚像素精度抗锯齿圆光栅器的参考。它很慢，:)，我需要计算圆与圆的边界框中每个像素的交集面积才能得到alpha。你可以看到它的工作原理，似乎没有数字伪影出现。下图是一堆圆的图，半径从0.3 px增加到大约6 px，呈螺旋形布局。