圆与长方形相交的面积
Posted
技术标签:
【中文标题】圆与长方形相交的面积【英文标题】：Area of intersection between circle and rectangle 【发布时间】：2010-10-11 23:01:12 【问题描述】：我正在寻找一种快速的方法来确定矩形和圆形之间的相交面积（我需要进行数百万次这样的计算）。
一个特定的属性是，在所有情况下，圆和矩形总是有 2 个交点。
【问题讨论】：
他们只有2个交点吗？还是它们至少有 2 个交叉点？ 需要以平方为单位计算面积，还是返回一组定义面积的线段？ 如果一个在另一个内部，或者如果两者完全不重叠，则没有交点。如果圆与矩形的任意一条边相切，则只有一个交点。 您到底需要知道什么？如果是为了简单的比较，你可能会做的比你需要做的要少。 【参考方案1】：给定 2 个交点：
0 个顶点在圆圈内：circular segment 的面积
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^
三个顶点在圆内：矩形的面积减去三角形的面积加上圆弧的面积
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
要计算这些面积：
您需要使用的大部分积分都可以通过查找intersection of a line and a circle
您需要计算的区域可以通过computing the area of a circular segment和computing the area of a triangle找到。
您可以通过计算顶点到圆心的距离是否小于半径来确定顶点是否在圆内。
【讨论】：
【参考方案2】：我意识到这已在不久前得到解答，但我正在解决同样的问题，但我找不到可以使用的开箱即用的可行解决方案。请注意，我的盒子是axis aligned，OP 并没有完全指定。下面的解决方案是完全通用的，适用于任意数量的交叉点（不仅仅是两个）。请注意，如果您的框不是轴对齐的（但仍然是直角框，而不是一般的quads），您可以利用圆形的优势，旋转所有内容的坐标，使框最终轴对齐然后使用此代码。
我想使用集成 - 这似乎是个好主意。让我们从写一个明显的画圆公式开始：
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，我们已经说过unit，意思是r = 1)。现在我们可以通过对y进行积分来计算一个半圆和一个无限高的盒子在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
这可能不是很有用，因为无限高的盒子不是我们想要的。我们需要再添加一个参数，以便能够释放无限高盒子的底部边缘：
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 add，而对于与 x 轴相交的框则使用两倍（以及另外 1 个 add）。请注意，除法是按半径划分的，可以减少为两个除法和一堆乘法。如果有问题的框与 x 轴相交但不与 y 轴相交，您可以将所有内容旋转pi/2 并以原始成本进行计算。
我将它用作调试亚像素精度抗锯齿圆光栅化器的参考。它慢得要命:)，我需要计算圆的边界框中每个像素与圆的交点面积才能得到 alpha。您可以看到它有效，并且似乎没有出现数字伪影。下图是一堆半径从 0.3 px 增加到大约 6 px 的圆，呈螺旋状排列。
【讨论】：
【参考方案3】：我希望发布这样一个老问题的答案还不错。我查看了上述解决方案并制定了一个类似于丹尼尔斯第一个答案的算法，但更严格一些。
简而言之，假设整个区域都在矩形中，减去外部半平面中的四个段，然后添加四个外部象限中的任何区域，沿途丢弃琐碎的情况。
伪代码（我的实际代码只有~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
顺便说一句，平面象限中包含的圆的面积的最后一个公式很容易导出为圆段、两个直角三角形和一个矩形的总和。
享受吧。
【讨论】：
【参考方案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 Qualification round 问题中的一个非常相似：Fly Swatter。您也可以点击分数链接下载解决方案的源代码。
【讨论】：
【参考方案5】：谢谢你的回答，
我忘了提到面积估计就足够了。 那;这就是为什么最后，在查看了所有选项之后，我使用了蒙特卡罗估计，我在圆圈中生成随机点并测试它们是否在框中。
就我而言，这可能更高效。 （我在圆上放置了一个网格，我必须测量属于每个网格单元的圆的比率。）
谢谢
【讨论】：
啊，估计没问题有很大的不同：]【参考方案6】：也许您可以使用this question 的答案，其中询问圆形和三角形之间的相交面积。将您的矩形分成两个三角形并使用那里描述的算法。
另一种方法是在两个交点之间画一条线，这会将您的区域分成一个多边形（3 或 4 条边）和一个 circular segment，因为您应该能够更轻松地找到库或自己进行计算.
【讨论】：
【参考方案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)));
【讨论】：
以上是关于圆与长方形相交的面积的主要内容，如果未能解决你的问题，请参考以下文章