我正在寻找一种快速的方法来确定矩形和圆形之间的相交区域（我需要进行数百万次这样的计算）。
一个特定的属性是，在所有情况下，圆和矩形总是有 2 个交点。
问问题
30302 次
7 回答
64
给定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
要计算这些区域：
您需要使用的大部分点都可以通过找到一条线和一个圆的交点来找到
您可以通过计算顶点到中心的距离是否小于半径来确定顶点是否在圆内。
于 2009-03-07T19:11:31.800 回答
14
我意识到这已在不久前得到解答，但我正在解决同样的问题，但我找不到可以使用的开箱即用的可行解决方案。请注意，我的框是轴对齐的，OP 并没有完全指定。下面的解决方案是完全通用的，适用于任意数量的交叉点（不仅仅是两个）。请注意，如果您的框不是轴对齐的（但仍然是直角框，而不是一般的quads），您可以利用圆形的优势，旋转所有内容的坐标，使框最终轴对齐，然后使用此代码。
我想使用集成——这似乎是个好主意。让我们从写一个明显的画圆公式开始：
x = center.x + cos(theta) * radius
y = center.y + sin(theta) * radius
^
|
|**### **
| #* # * * x
|# * # * # y
|# * # *
+-----------------------> theta
* # * #
* # * #
* #* #
***###
这很好，但我无法将那个圆圈的面积整合到xor 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，得到一个圆的上半部分的面积。这仅适用于xin [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，我们已经说过单位，意思是r = 1）。现在我们可以通过积分 over 计算半圆和无限高的盒子的交点面积，该盒子位于 x 轴上（圆心也位于 x 轴上）y：
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 的圆，呈螺旋状排列。
于 2015-09-21T15:25:56.527 回答
5
我希望发布这样一个老问题的答案还不错。我查看了上述解决方案并制定了一个类似于 Daniels 第一个答案的算法，但更严格一些。
简而言之，假设整个区域都在矩形中，减去外部半平面中的四个段，然后添加四个外部象限中的任何区域，沿途丢弃琐碎的情况。
pseudocde（我的实际代码只有~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
顺便说一句，平面象限中包含的圆的面积的最后一个公式很容易导出为圆段、两个直角三角形和一个矩形的总和。
享受。
于 2012-08-18T13:25:25.213 回答
3
以下是如何计算圆心位于矩形外的圆和矩形之间的重叠区域。其他情况可以归结为这个问题。
面积可以通过积分圆方程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 资格赛中的一个问题非常相似：Fly Swatter。您也可以点击分数链接下载解决方案的源代码。
于 2009-03-07T20:18:05.120 回答
2
感谢您的回答，
我忘了提到面积估计就足够了。那; 这就是为什么最后，在查看了所有选项之后，我使用了蒙特卡罗估计，我在圆圈中生成随机点并测试它们是否在框中。
就我而言，这可能更高效。（我在圆上放置了一个网格，我必须测量属于每个网格单元的圆的比率。）
谢谢
于 2009-03-08T16:53:03.747 回答
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)));
}
于 2010-08-16T06:53:27.303 回答