圆与矩形的交面积

Question

我正在寻找一种快速方法来确定矩形和圆形之间的相交面积（我需要进行数百万次这样的计算）。

一个特定的属性是，在所有情况下，圆形和矩形始终有 2 个交点。

Answer 1

给定 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

要计算这些面积：

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

• 您需要计算的面积可以通过 计算圆弧段的面积 和 计算三角形的面积。

• 可以通过计算顶点到圆心的距离是否小于半径来判断顶点是否在圆内。

  可以通过计算顶点到圆心的距离是否小于半径来判断

Answer 2

我意识到这个问题不久前就得到了回答，但我正在解决同样的问题，但我找不到可以使用的开箱即用的可行解决方案。 请注意，我的框是轴对齐，OP并未完全指定这一点。 下面的解决方案是完全通用的，并且适用于任意数量的交叉点（不仅仅是两个）。 请注意，如果您的框不是轴对齐的（但仍然是直角框，而不是一般的四边形），您可以利用圆形的优势，旋转所有物体的坐标，使盒子最终轴对齐，然后然后使用此代码。

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

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

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

这很好，但我无法在 x 或 y 上积分该圆的面积； 这些是不同的数量。 我只能在 theta 角度上积分，产生披萨片区域。 不是我想要的。 让我们尝试改变一下论点：

(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

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

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

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

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

现在我们可以通过对y积分 可能不是很有用，因为无限高的盒子不是我们想要的。 我们需要再添加一个参数，以便能够释放无限高盒子的底部边缘：

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

其中 h 是无限盒子的底部边缘与 x 轴的（正）距离。 section 函数计算单位圆与由 y = 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 个相加，对于与 x 轴相交的框，使用该值的两倍（以及另外 1 个相加）。 请注意，除法是按半径进行的，可以减少为两次除法和一堆乘法。 如果相关框与 x 轴相交但不与 y 轴相交，则可以将所有内容旋转 pi/2 并按原始成本进行计算。

我使用它作为调试亚像素精度抗锯齿圆形光栅器的参考。 它太慢了:)，我需要计算圆的边界框中每个像素与圆的相交面积以获得alpha。 您可以看到它有效，并且似乎没有出现数字伪影。 下图是一堆半径从 0.3 px 增加到约 6 px 的圆，呈螺旋状排列。

Answer 3

我希望对这样一个老问题发表答案并不是坏事。 我查看了上述解决方案，并制定了一个与丹尼尔斯第一个答案类似的算法，但更严格一些。

简而言之，假设整个区域都在矩形中，减去外部半平面中的四个线段，然后添加四个外部象限中的任何区域，并在此过程中丢弃琐碎的情况。

伪cde（我的实际代码只有~12行..）

find the signed (negative out) normalized distance from the circle center
to each of the infinitely extended rectangle edge lines,
ie.
d_1=(xcenter-xleft)/r
d_2=(ycenter-ybottom)/r
etc

for convenience order 1,2,3,4 around the edge. If the rectangle is not
aligned with the cartesian coordinates this step is more complicated but
the remainder of the algorithm is the same

If ANY d_i <=- 1 return 0

if ALL d_i >=  1 return Pi r^2

this leave only one remaining fully outside case: circle center in
an external quadrant, and distance to corner greater than circle radius:

for each adjacent i,j (ie. i,j=1,2;2,3;3,4;4,1)
     if d_i<=0 and d_j <= 0 and d_i^2+d_j^2 > 1 return 0

now begin with full circle area  and subtract any areas in the
four external half planes

Area= Pi r^2
for each  d_i>-1
     a_i=arcsin( d_i )  #save a_i for next step
     Area -= r^2/2 (Pi - 2 a_i - sin(2 a_i)) 

At this point note we have double counted areas in the four external
quadrants, so add back in:

for each adjacent i,j
   if  d_i < 1 and   d_j < 1  and d_i^2+d_j^2 < 1
       Area += r^2/4 (Pi- 2 a_i - 2 a_j -sin(2 a_i) -sin(2 a_j) + 4 sin(a_i) sin(a_j))

return Area

顺便说一句，平面象限中包含的圆面积的最后一个公式很容易导出为圆弧段、两个直角三角形和一个矩形的总和。

享受。

Answer 4

下面介绍如何计算圆与矩形的重叠面积，其中圆心位于矩形之外。 其他情况可以简化为这个问题。

面积可以通过积分圆方程 y = sqrt[a^2 - (xh)^2] + k 来计算，其中 a 是半径，(h,k) 是圆心，找到曲线下面积。 您可以使用计算机集成，将区域划分为许多小矩形并计算它们的总和，或者在这里仅使用封闭形式。

这里是实现上述概念的 C# 源代码。 请注意，有一种特殊情况，即指定的 x 位于圆的边界之外。 我在这里只是使用一种简单的解决方法（在所有情况下都不会产生正确的答案）

public static void RunSnippet()
{
    // test code
    double a,h,k,x1,x2;
    a = 10;
    h = 4;
    k = 0;
    x1 = -100;
    x2 = 100;

    double r1 = Integrate(x1, a, h, k);
    double r2 = Integrate(x2, a, h, k);

    Console.WriteLine(r2 - r1);

}

private static double Integrate(double x, double a,double h, double k)
{
    double a0 = a*a - (h-x)*(h-x);

    if(a0 <= 0.0){
        if(k == 0.0)
            return Math.PI * a * a / 4.0 * Math.Sign(x);
        else
            throw new Exception("outside boundaries");
    }

    double a1 = Math.Sqrt(a*a - (h-x)*(h-x)) * (h-x);
    double area = 0.5 * Math.Atan(a1 / ((h-x)*(h-x) - a*a))*a*a - 0.5 * a1 + k * x;
    return area;
}

注意：此问题与 Google Code Jam 2008 资格赛问题：苍蝇拍。 您也可以单击分数链接下载解决方案的源代码。

Answer 5

感谢您的回答，

我忘了提及面积估计就足够了。
那; 这就是为什么最后，在查看了所有选项之后，我采用了蒙特卡罗估计，在圆圈中生成随机点并测试它们是否在框中。

就我而言，这可能性能更高。 （我在圆上放置了一个网格，我必须测量属于每个网格单元的圆的比率。）

谢谢

Answer 6

也许您可以使用这个问题的答案< /a>，询问圆和三角形之间的交点面积。 将矩形分成两个三角形并使用此处描述的算法。

另一种方法是在两个交点之间画一条线，这会将您的区域分割为一个多边形（3 或 4 条边）和一个 圆弧段，您应该能够更轻松地找到库或自己进行计算。

Answer 7

这是该问题的另一种解决方案：

public static bool IsIntersected(PointF circle, float radius, RectangleF rectangle)
{

        var rectangleCenter = new PointF((rectangle.X +  rectangle.Width / 2),
                                         (rectangle.Y + rectangle.Height / 2));

        var w = rectangle.Width  / 2;
        var h = rectangle.Height / 2;

        var dx = Math.Abs(circle.X - rectangleCenter.X);
        var dy = Math.Abs(circle.Y - rectangleCenter.Y);

        if (dx > (radius + w) || dy > (radius + h)) return false;


        var circleDistance = new PointF
                                 {
                                     X = Math.Abs(circle.X - rectangle.X - w),
                                     Y = Math.Abs(circle.Y - rectangle.Y - h)
                                 };


        if (circleDistance.X <= (w))
        {
            return true;
        }

        if (circleDistance.Y <= (h))
        {
            return true;
        }

        var cornerDistanceSq = Math.Pow(circleDistance.X - w, 2) + 
                    Math.Pow(circleDistance.Y - h, 2);

        return (cornerDistanceSq <= (Math.Pow(radius, 2)));
}