圆和矩形之间的交叉区域
use*_*424 23 java math intersection area
我正在寻找一种快速的方法来确定矩形和圆形之间的交叉区域(我需要进行数百万次这些计算).
一个特定的属性是,在所有情况下,圆和矩形总是有2个交点.
Dan*_*ant 65
给定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
要计算这些区域:
我意识到这已经回答了一段时间,但我正在解决同样的问题,我找不到一个可以使用的开箱即用的解决方案.请注意,我的框是轴对齐的,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,以获得圆的上半部分的区域.这仅适用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,我们已经说过单位,这意味着r = 1).现在我们可以计算一个半圆和一个无限高的盒子的交叉区域,站在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添加,对于那些盒子,它使用了两倍(还有1个添加).注意,除法是半径的,可以减少到两个除法和一系列乘法.如果相关的框与x轴相交但未与y轴相交,则可以旋转所有内容pi/2并以原始成本进行计算.
我正在使用它作为调试亚像素精确抗锯齿圆光栅化器的参考.地狱很慢:),我需要计算圆圈的边界框中每个像素与圆圈的交点区域得到alpha.您可以看到它的工作原理并且似乎没有出现数字工件.下图是一组圆的图,半径从0.3 px增加到约6 px,呈螺旋状排列.
我希望能找到这样一个老问题的答案并不是坏事.我查看了上面的解决方案并制定了一个类似于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
顺便提及,包含在平面象限中的圆的面积的最后公式容易地导出为圆形段,两个直角三角形和矩形的总和.
请享用.
| 归档时间： |
|
| 查看次数： |
24902 次 |
| 最近记录： |