我正在寻找一种快速的方法来确定矩形和圆之间相交区域的面积（我需要执行数百万次这样的计算）。
一个特定的属性是，在所有情况下，圆和矩形总是有2个交点。
圆和矩形的交集面积
25
- user67424
4
1它们只有两个交点吗？还是至少有两个交点？ - bigmonachus
你需要计算面积（以平方单位）还是返回定义该区域的线段集合？ - Leonard
如果一个图形在另一个图形内部，或者两个图形根本不重叠，则它们没有交点。如果圆和矩形的任何一条边相切，那么它们只有一个交点。 - duffymo
你需要了解什么？如果只是简单的比较，你可能不需要进行完整的操作来得到精确的答案。 - Thorbjørn Ravn Andersen
7个回答
65
给定两个交点：
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
| |
两个顶点在圆内: 两个三角形的面积和一个圆形切片的面积之和。
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
要计算这些区域:
- Daniel L
15
我知道这个问题早就有答案了，但我正在解决同样的问题，而且我找不到可以直接使用的解决方案。请注意，我的框是与坐标轴对齐的轴对齐盒子，这并没有被提问者明确说明。下面的解决方案是完全通用的，并且适用于任意数量的交叉点（不仅限于两个）。请注意，如果您的框不是轴对齐的（但仍然是带有直角而不是一般四边形的框），则可以利用圆形的圆形特征，将所有内容的坐标旋转使框变成轴对齐的，然后再使用这段代码。
我想使用积分-那似乎是个好主意。让我们先写一个明显的公式来绘制一个圆：
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进行积分，从而得到上半圆的面积。这仅适用于x在[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
这个函数的积分确实为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轴上方（使用上面的等式）
- 盒子在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轴相交但不与y轴相交的框，这使用了其双倍（和1个更多的add）。请注意，除法是通过半径进行的，可以将其减少到两个除法和一堆乘法。如果所讨论的框与x轴相交但未与y轴相交，则可以通过旋转pi/2来进行计算并在原始成本中完成计算。
我将其用作参考来调试亚像素精确抗锯齿圆形光栅化程序。它非常慢 :), 我需要计算圆的边界框中每个像素与圆的相交面积以获取alpha值。您可以大致看到它起作用，并且似乎没有出现任何数值问题。下面的图形是一堆半径从0.3 px到约6 px逐渐增加的圆，排列成螺旋形状。
- the swine
5
我希望回答这样一个老问题不会给人留下坏的形象。我查看了上述解决方案并制定了一种算法，它类似于丹尼尔的第一个答案，但更为紧凑。
简而言之，假设整个区域在矩形中，减去在外部半平面中的四个分段，然后添加任何四个外部象限中的面积，在此过程中舍弃微不足道的情况。
伪代码（我的实际代码仅有约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
顺便提一下，在一个平面象限内的圆面积公式可以很容易地推导出来，它是由一个圆弓形、两个直角三角形和一个矩形组成的。
享受吧。
- agentp
3
以下是如何计算圆和矩形重叠区域的方法，其中圆心位于矩形外部。其他情况可以归结为这个问题。
可以通过积分圆方程y = sqrt[a^2 - (x-h)^2] + k（其中a是半径，(h,k)是圆心）来计算曲线下面的面积。您可以使用计算机积分，将该区域分成许多小矩形并计算它们的总和，或者在此处使用闭式公式。
以下是实现上述概念的C#源代码。请注意，有一种特殊情况，指定的x位于圆的边界之外。我在这里只使用了一个简单的解决方法（并不在所有情况下产生正确的答案）。
可以通过积分圆方程y = sqrt[a^2 - (x-h)^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资格赛中的一个问题非常相似：Fly Swatter。您可以单击得分链接下载解决方案的源代码。
- m3rLinEz
2
非常感谢您的解答，
我忘了提到区域估计就足够了。这就是为什么最后在查看所有选项后，我选择了 Monte Carlo 估计，在圆形内生成随机点并测试它们是否在框内。
在我的情况下，这可能更有效率。（我在圆上放置了一个网格，并且必须测量属于每个网格单元的圆的比例。）
谢谢！
- user67424
1
啊，估计结果正确确实有很大的影响:] - Daniel LeCheminant
1
这里有另一个解决问题的方案：
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)));
}
- Bassam Alugili
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
- 相关问题
- 16 两个圆的交集面积
- 3 快速计算三角形和单位正方形的交集面积
- 4 Ruby中的矩形交集
- 22 计算圆和三角形之间的交集面积？
- 116 重叠圆的组合面积
- 4 计算多边形的最小面积矩形
- 28 矩形交集面积
- 61 两个矩形重叠的面积是多少？
- 3 计算平面与矩形棱柱的交集面积
- 4 3D中的圆与圆的交集