Zone d'intersection entre cercle et rectangle
https://stackoverflow.com/questions/622287
-
05-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
Je cherche un moyen rapide de déterminer la zone d'intersection entre un rectangle et un cercle (je dois faire des millions de ces calculs).
Une propriété spécifique est que, dans tous les cas, le cercle et le rectangle ont toujours 2 points d'intersection.
La solution
Étant donné 2 points d'intersection:
0 sommet se trouve à l'intérieur du cercle: zone d'un segment circulaire
XXXXX -------------------
X X X X Circular segment
X X XX XX
+-X-------X--+ XXXXXXXX
| X X |
| XXXXX |
1 sommet se trouve à l'intérieur du cercle: somme des zones d'un segment circulaire et d'un triangle.
XXXXX XXXXXXXXX
X X Triangle ->X _-X
X X X _- X
X +--X--+ X _- X <- Circular segment
X | X | X- XXX
XXXXX | XXXX
| |
2 sommets se trouvent à l'intérieur du cercle: somme de l'aire de deux triangles et d'un segment circulaire
XXXXX +------------X
X X | _--'/'X
X +--X--- Triangle->| _-- / X
X | X |_-- /XX <- Circular segment
X +-X---- +-------XX
XXXXX Triangle^
3 sommets se trouvent à l'intérieur du cercle: surface du rectangle moins la surface d'un triangle plus la surface d'un segment circulaire
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
Pour calculer ces zones:
-
La plupart des points à utiliser peuvent être trouvés en trouvant l'intersection d'une ligne et d'un cercle
-
Les zones à calculer sont accessibles en en calculant la surface d'un segment circulaire et calculer l'aire d'un triangle .
-
Vous pouvez déterminer si un sommet est à l'intérieur du cercle en calculant si sa distance au centre est inférieure au rayon.
Autres conseils
Je me rends compte que la réponse a été donnée il y a un moment, mais je résous le même problème et je ne pouvais pas trouver de solution prête à l'emploi, utilisable. Notez que mes boîtes sont alignées sur l'axe , ce qui n'est pas tout à fait spécifié par l'OP. La solution ci-dessous est tout à fait générale et fonctionnera pour n’importe quel nombre d’intersections (pas seulement deux). Notez que si vos boîtes ne sont pas alignées les unes sur les autres (mais restent sur des cases à angles droits plutôt que sur des quads ), vous pouvez tirer parti des cercles arrondis, faire pivoter les coordonnées de tout pour que la boîte se termine alignée et puis utiliser ce code.
Je souhaite utiliser l'intégration - cela semble être une bonne idée. Commençons par écrire une formule évidente pour tracer un cercle:
x = center.x + cos(theta) * radius
y = center.y + sin(theta) * radius
^
|
|**### **
| #* # * * x
|# * # * # y
|# * # *
+-----------------------> theta
* # * #
* # * #
* #* #
***###
C'est bien, mais je ne parviens pas à intégrer la zone de ce cercle sur x ou y ; ce sont des quantités différentes. Je ne peux intégrer que sur l'angle thêta , donnant des zones de tranches de pizza. Pas ce que je veux. Essayons de changer les arguments:
(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
C'est plus comme ça. Maintenant, étant donné la plage de x , je peux intégrer plus de y , pour obtenir une zone de la moitié supérieure d'un cercle. Ceci ne vaut que pour x dans [center.x - radius, center.x + radius] (les autres valeurs provoqueront des sorties imaginaires) mais nous savons que la zone en dehors de cette plage est zéro, donc cela est facilement manipulé. Supposons que l’unité soit un cercle pour plus de simplicité, nous pouvons toujours brancher le centre et le rayon plus tard:
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
Cette fonction a en effet l'intégrale de pi / 2 , puisqu'il s'agit de la moitié supérieure d'un cercle unitaire (la surface d'un demi-cercle est pi * r ^ 2/2 et nous avons déjà dit unité , ce qui signifie r = 1 ). Maintenant, nous pouvons calculer l'aire d'intersection d'un demi-cercle et d'une boîte infiniment haute, debout sur l'axe des x (le centre du cercle se trouve également sur l'axe des x) en intégrant sur 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
Cela peut ne pas être très utile, car des boîtes infiniment hautes ne sont pas ce que nous voulons. Nous devons ajouter un paramètre supplémentaire afin de pouvoir libérer le bord inférieur de la boîte infiniment haute:
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
Où h est la distance (positive) du bord inférieur de notre boîte infinie par rapport à l'axe des x. La fonction section calcule la position (positive) de l'intersection du cercle unité avec la ligne horizontale donnée par y = h et nous pourrions le définir en résolvant:
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)
Maintenant, nous pouvons faire avancer les choses. Alors, comment calculer l’aire d’intersection d’une boîte finie intersectant un cercle unitaire au-dessus de l’axe des 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
C'est bien. Alors que diriez-vous d'une boîte qui n'est pas au-dessus de l'axe des x? Je dirais que toutes les boîtes ne sont pas. Trois cas simples se présentent:
- la boîte est au-dessus de l'axe x (utilisez l'équation ci-dessus)
- la boîte est en dessous de l'axe x (retournez le signe des coordonnées y et utilisez l'équation ci-dessus)
- la case coupe l’axe des x (divisez la case en deux parties, calculez l’aire des deux en utilisant ce qui précède et additionnez-les)
Assez facile? Écrivons du code:
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);
}
Et quelques tests unitaires de base:
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
Le résultat est:
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
Ce qui me semble correct. Vous pouvez vouloir intégrer les fonctions si vous ne faites pas confiance à votre compilateur pour le faire à votre place.
Ceci utilise 6 sqrt, 4 asin, 8 div, 16 mul et 17 add pour les cases qui ne coupent pas l'axe des x et un double de cela (et une autre addition) pour les cases qui le font. Notez que les divisions sont par rayon et pourraient être réduites à deux divisions et à un tas de multiplications. Si la zone en question coupait l'axe des x mais ne coupait pas l'axe des y, vous pouvez tout faire pivoter de pi / 2 et effectuer le calcul dans le coût d'origine.
Je l'utilise comme référence pour déboguer le rasterisateur de cercle antialiasé précis au sous-pixel. C'est lent comme l'enfer :), je dois calculer l'aire d'intersection de chaque pixel dans la boîte englobante du cercle avec le cercle pour obtenir un alpha. Vous pouvez en quelque sorte voir que cela fonctionne et qu'aucun artefact numérique ne semble apparaître. La figure ci-dessous est un graphique représentant un groupe de cercles dont le rayon augmente de 0,3 à environ 6 px, disposés en spirale.
J'espère que sa forme n'est pas mauvaise pour poster une réponse à une question aussi ancienne. J'ai examiné les solutions ci-dessus et mis au point un algorithme similaire à celui de la première réponse de Daniels, mais un peu plus serré.
En résumé, supposons que toute la surface soit dans le rectangle, soustrayez les quatre segments des demi-plans externes, puis ajoutez les zones des quatre quadrants externes, en éliminant les cas triviaux.
pseudocde (mon code actuel n’est que ~ 12 lignes ..)
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
Incidemment, cette dernière formule pour l'aire d'un cercle contenu dans un quadrant plan est facilement calculée en tant que somme d'un segment circulaire, de deux triangles rectangles et d'un rectangle.
Profitez.
Voici comment calculer la zone de chevauchement entre le cercle et le rectangle où le centre du cercle se trouve en dehors du rectangle. D'autres cas peuvent être réduits à ce problème.
La surface peut être calculée en intégrant l’équation du cercle y = sqrt [a ^ 2 - (xh) ^ 2] + k où a est le rayon, (h, k) est le centre du cercle, pour trouver l'aire sous la courbe. Vous pouvez utiliser une intégration informatique où la zone est divisée en plusieurs petits rectangles et en calculer la somme, ou simplement utiliser un formulaire fermé ici.
Et voici une source C # implémentant le concept ci-dessus. Notez qu'il existe un cas spécial où le x spécifié est en dehors des limites du cercle. Je viens d'utiliser une solution de contournement simple ici (qui ne produit pas les bonnes réponses dans tous les cas)
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;
}
Remarque: ce problème est très similaire à celui de Confiture de code Google. Qualification 2008 problème: Fly Swatter . Vous pouvez également cliquer sur les liens de la partition pour télécharger le code source de la solution.
Merci pour les réponses,
J'ai oublié de mentionner que les estimations de surface suffisaient. Cette; C’est pourquoi, après avoir examiné toutes les options, j’ai opté pour l’estimation monte-carlo où je génère des points aléatoires dans le cercle et vérifie s’ils sont dans la boîte.
Dans mon cas, cela est probablement plus performant. (J'ai une grille placée sur le cercle et je dois mesurer le rapport du cercle appartenant à chacune des cellules de la grille.)
Merci
Peut-être pouvez-vous utiliser la réponse à cette question , où la zone d'intersection entre un cercle et un triangle est posée. Divisez votre rectangle en deux triangles et utilisez les algorithmes qui y sont décrits.
Vous pouvez également tracer une ligne entre les deux points d’intersection, ce qui divisera votre zone en un polygone (3 ou 4 arêtes) et un segment circulaire , vous devriez pouvoir trouver les bibliothèques plus facilement ou faire les calculs vous-même.
Voici une autre solution au problème:
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)));
}
