圆和矩形之间的交叉区域

我正在寻找一种快速的方法来确定矩形和圆形之间的交叉区域(我需要进行数百万次这些计算)。

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

给定2个交点:

0个顶点位于圆内: 圆弧段的区域

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

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

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

2个顶点位于圆内:两个三角形和圆形区域的总和

  XXXXX +------------X XX | _--'/'XX +--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| XX ''. X <- Circular segment inside XX ^|XXX | X XXXXX 

要计算这些区域:

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

  • 您需要计算的区域可以通过计算圆形 区域的面积和计算三角形的面积来找到 。

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

我希望能找到这样一个老问题的答案并不坏。 我查看了上面的解决方案并制定了一个类似于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 

顺便提及,包含在平面象限中的圆的面积的最后公式容易地导出为圆形段,两个直角三角形和矩形的总和。

请享用。

我意识到这已经回答了一段时间,但我正在解决同样的问题,我找不到一个可以使用的开箱即用的解决方案。 请注意,我的框是轴对齐的 ,OP没有明确指出。 下面的解决方案是完全通用的,适用于任意数量的交叉点(不仅仅是两个)。 请注意,如果您的方框不是轴对齐的(但仍然是直角的方框,而不是一般的四边形 ),您可以利用圆形的圆形,旋转所有内容的坐标,使框最终轴对齐, 然后使用此代码。

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

 x = center.x + cos(theta) * radius y = center.y + sin(theta) * radius ^ | |**### ** | #* # * * x |# * # * # y |# * # * +-----------------------> theta * # * # * # * # * #* # ***### 

这很好,但是我无法将该圆的面积整合到xy ; 那些是不同的数量。 我只能整合角度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 )。 现在我们可以通过积分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); } 

还有一些基本的unit testing:

 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添加,对于那些盒子,它使用了两倍(和另外1个添加)。 请注意,除法是半径的,可以减少为两个除法和一系列乘法。 如果相关的框与x轴相交但未与y轴相交,则可以通过pi/2旋转所有内容并以原始成本进行计算。

我正在使用它作为调试亚像素精确抗锯齿圆光栅化器的参考。 地狱很慢:),我需要计算圆圈的边界框中每个像素与圆圈的交点区域得到alpha。 您可以看到它有效并且似乎没有出现数字工件。 下图是一系列圆的图,半径从0.3 px增加到约6 px,呈螺旋状排列。

界

以下是如何计算圆和矩形之间的重叠区域,其中圆的中心位于矩形之外。 其他情况可以减少到这个问题。

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

替代文字

这是一个实现上述概念的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 - (hx)*(hx); 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 - (hx)*(hx)) * (hx); double area = 0.5 * Math.Atan(a1 / ((hx)*(hx) - a*a))*a*a - 0.5 * a1 + k * x; return area; } 

注意:此问题与Google Code Jam 2008资格赛问题非常相似: Fly Swatter 。 您也可以单击分数链接以下载解决方案的源代码。

谢谢你的回答,

我忘了提到区域估计就足够了。 那; 为什么最后,在查看了所有选项之后,我选择了monte-carlo,我在圆圈中生成随机点并测试它们是否在框中。

在我的情况下,这可能更高效。 (我在圆圈上放置了一个网格,我必须测量属于每个网格单元的圆的比例。)

谢谢

也许你可以使用这个问题的答案,其中要求圆和三角形之间的交叉区域。 将矩形拆分为两个三角形,并使用其中描述的算法。

另一种方法是在两个交叉点之间绘制一条线,这会将您的区域分成多边形(3或4条边)和圆形线段 ,因为您应该能够更轻松地找到库或自己进行计算。

以下是该问题的另一种解决方案:

 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))); }