円と長方形の交差領域
https://stackoverflow.com/questions/622287
-
05-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian質問
長方形と円の交点の領域をすばやく決定する方法を探しています（これらの計算を何百万も行う必要があります）。
特定のプロパティは、すべての場合において、円と長方形には常に2つの交点があるということです。
解決
交差点を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つの頂点は円の中にあります：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によって完全には指定されていません。以下のソリューションは完全に一般的であり、2つだけでなく、任意の数の交差点で機能します。ボックスが軸に沿っていない場合（ただし、一般的な 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軸上にある）の上に立って、半円と無限に高いボックスの交差面積を計算できるようになりました。 / p>
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
これは非常に便利ではないかもしれません。無限に長いボックスは私たちが望んでいるものではないからです。無限の高さのボックスの下端を解放できるようにするには、もう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軸の上にないボックスはどうでしょうか？すべてのボックスが揃っているわけではありません。 3つの単純なケースが発生します。
- ボックスは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の追加を使用し、ボックスの場合はその2倍（さらに1追加）を使用します。除算は半径によるものであり、2つの除算と乗算の束に減らすことができることに注意してください。問題のボックスがx軸と交差したがy軸と交差しなかった場合、すべてを pi / 2 だけ回転させ、元のコストで計算できます。
私は
このような古い質問への回答を投稿するのが悪くないことを願っています。上記のソリューションを検討し、ダニエルズの最初の答えに似たアルゴリズムを作成しましたが、かなり厳密です。
要するに、領域全体が長方形の中にあると仮定し、外部半平面の4つのセグメントを差し引いてから、4つの外部象限の領域を追加し、途中で些細なケースを破棄します。
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
ちなみに、平面象限に含まれる円の面積の最後の式は、円形セグメント、2つの直角三角形、1つの長方形の合計として容易に導き出されます。
お楽しみください。
以下は、円の中心が長方形の外側にある、円と長方形の間の重なり合う面積を計算する方法です。他のケースはこの問題に減らすことができます。
面積は、円方程式 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予選ラウンドの問題：フライスワッター。スコアリンクをクリックして、ソリューションのソースコードもダウンロードできます。
回答をありがとう、
地域の推定値で十分だったことに言及するのを忘れました。 それ;結局、すべてのオプションを確認した後、モンテカルロ推定を使用して、円の中にランダムなポイントを生成し、それらがボックス内にあるかどうかをテストしました。
私の場合、これはおそらくより高性能です。 （円の上にグリッドを配置し、各グリッドセルに属する円の比率を測定する必要があります。）
ありがとう
おそらくへの回答を使用できますこの質問では、円と三角形の交差領域が尋ねられます。長方形を2つの三角形に分割し、そこに記載されているアルゴリズムを使用します。
別の方法は、2つの交差点の間に線を引くことです。これにより、エリアがポリゴン（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)));
}
