java - graficar - Área de intersección entre círculo y rectángulo
graficar circulos en java (7)
Estoy buscando una manera rápida de determinar el área de intersección entre un rectángulo y un círculo (necesito hacer millones de estos cálculos).
Una propiedad específica es que, en todos los casos, el círculo y el rectángulo siempre tienen 2 puntos de intersección.
A continuación se explica cómo calcular el área de superposición entre el círculo y el rectángulo donde el centro del círculo se encuentra fuera del rectángulo. Otros casos pueden reducirse a este problema.
El área se puede calcular integrando la ecuación de círculo y = sqrt [a ^ 2 - (xh) ^ 2] + k donde a es el radio, (h, k) es el centro del círculo, para encontrar el área bajo la curva. Puede usar la integración de computadora donde el área se divide en muchos rectángulos pequeños y calcular la suma de ellos, o simplemente usar la forma cerrada aquí.
Y aquí hay una fuente de C # que implementa el concepto anterior. Tenga en cuenta que hay un caso especial donde la x especificada se encuentra fuera de los límites del círculo. Solo uso una solución simple aquí (que no produce las respuestas correctas en todos los casos)
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;
}
Nota: Este problema es muy similar a uno en el problema de la ronda de clasificación de Google Code Jam 2008 : Swatter de mosca . También puede hacer clic en los enlaces de puntuación para descargar el código fuente de la solución.
Aquí hay otra solución para el problema:
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)));
}
Dados 2 puntos de intersección:
0 vértices están dentro del círculo: el área de un segmento circular
XXXXX -------------------
X X X X Circular segment
X X XX XX
+-X-------X--+ XXXXXXXX
| X X |
| XXXXX |
1 vértice está dentro del círculo: la suma de las áreas de un segmento circular y un triángulo.
XXXXX XXXXXXXXX
X X Triangle ->X _-X
X X X _- X
X +--X--+ X _- X <- Circular segment
X | X | X- XXX
XXXXX | XXXX
| |
2 vértices están dentro del círculo: la suma del área de dos triángulos y un segmento circular
XXXXX +------------X
X X | _--''/''X
X +--X--- Triangle->| _-- / X
X | X |_-- /XX <- Circular segment
X +-X---- +-------XX
XXXXX Triangle^
3 vértices están dentro del círculo: el área del rectángulo menos el área de un triángulo más el área de un segmento circular
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
Para calcular estas áreas:
La mayoría de los puntos que necesitarás usar se pueden encontrar al encontrar la intersección de una línea y un círculo
Las áreas que necesita calcular se pueden encontrar al calcular el área de un segmento circular y al área de un triángulo .
Puedes determinar si un vértice está dentro del círculo calculando si su distancia desde el centro es menor que el radio.
Espero que no sea una mala forma publicar una respuesta a una pregunta tan antigua. Revisé las soluciones anteriores y elaboré un algoritmo que es similar a la primera respuesta de Daniel, pero un poco más estricto.
En resumen, suponga que el área completa está en el rectángulo, reste los cuatro segmentos en los semiplanos externos, luego agregue cualquiera de las áreas en los cuatro cuadrantes externos, descartando los casos triviales en el camino.
pseudocde (mi código real es solo ~ 12 líneas ..)
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
Incidentalmente, la última fórmula para el área de un círculo contenido en un cuadrante plano se deriva fácilmente como la suma de un segmento circular, dos triángulos rectos y un rectángulo.
Disfrutar.
Gracias por las respuestas,
Olvidé mencionar que las estimaciones de área eran suficientes. Ese; Por eso, al final, después de ver todas las opciones, fui con la estimación de monte-carlo, donde genero puntos aleatorios en el círculo y pruebo si están en el cuadro.
En mi caso es probable que esto sea más eficaz. (Tengo una cuadrícula colocada sobre el círculo y tengo que medir la proporción del círculo que pertenece a cada una de las celdas de la cuadrícula).
Gracias
Me di cuenta de que esto fue respondido hace un tiempo, pero estoy resolviendo el mismo problema y no pude encontrar una solución viable que pudiera usar. Tenga en cuenta que mis cuadros están alineados con el eje , esto no está especificado por el OP. La solución a continuación es completamente general y funcionará para cualquier número de intersecciones (no solo dos). Tenga en cuenta que si sus cuadros no están alineados con el eje (pero los cuadros fijos con ángulos rectos, en lugar de los quads generales), puede aprovechar que los círculos sean redondos, rotar las coordenadas de todo para que el cuadro termine alineado con el eje y usa este código
Quiero usar la integración, me parece una buena idea. Empecemos por escribir una fórmula obvia para trazar un círculo:
x = center.x + cos(theta) * radius
y = center.y + sin(theta) * radius
^
|
|**### **
| #* # * * x
|# * # * # y
|# * # *
+-----------------------> theta
* # * #
* # * #
* #* #
***###
Esto es bueno, pero no puedo integrar el área de ese círculo sobre x o y ; Son cantidades diferentes. Solo puedo integrarme sobre el ángulo theta , produciendo áreas de porciones de pizza. No es lo que quiero. Intentemos cambiar los argumentos:
(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
Eso es lo que más me gusta. Ahora, dado el rango de x , puedo integrarme sobre y para obtener un área de la mitad superior de un círculo. Esto solo es válido para x en [center.x - radius, center.x + radius] (otros valores causarán salidas imaginarias) pero sabemos que el área fuera de ese rango es cero, por lo que se maneja fácilmente. Asumamos el círculo unitario para simplificar, siempre podemos volver a conectar el centro y el radio más adelante:
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
Esta función tiene una integral de pi/2 , ya que es la mitad superior de un círculo unitario (el área de la mitad de un círculo es pi * r^2 / 2 y ya hemos dicho unidad , lo que significa r = 1 ). Ahora podemos calcular el área de intersección de un semicírculo y un cuadro infinitamente alto, parados en el eje x (el centro del círculo también se encuentra en el eje x) al integrar sobre 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
Esto puede no ser muy útil, ya que las cajas infinitamente altas no son lo que queremos. Necesitamos agregar un parámetro más para poder liberar el borde inferior del cuadro infinitamente alto:
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
Donde h es la distancia (positiva) del borde inferior de nuestra caja infinita desde el eje x. La función de section calcula la posición (positiva) de la intersección del círculo unitario con la línea horizontal dada por y = h y podríamos definirla resolviendo:
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)
Ahora podemos poner las cosas en marcha. Entonces, cómo calcular el área de intersección de un cuadro finito que se interseca con un círculo unitario sobre el eje 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
Eso es bueno. Entonces, ¿qué tal un cuadro que no está por encima del eje x? Yo diría que no todas las cajas son. Surgen tres casos simples:
- el cuadro está sobre el eje x (use la ecuación anterior)
- el cuadro está debajo del eje x (voltea el signo de las coordenadas y y usa la ecuación anterior)
- el recuadro está intersecando el eje x (divida el recuadro en la mitad superior e inferior, calcule el área de ambos usando lo anterior y súmelos)
¿Suficientemente fácil? Vamos a escribir un código:
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);
}
Y algunas pruebas unitarias básicas:
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
La salida de esto es:
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
Lo que me parece correcto. Es posible que desee integrar las funciones si no confía en que su compilador lo haga por usted.
Esto usa 6 sqrt, 4 asin, 8 div, 16 mul y 17 agregados para cuadros que no se intersecan con el eje x y un doble de eso (y 1 agregado) para cuadros que sí lo hacen. Tenga en cuenta que las divisiones son por radio y podrían reducirse a dos divisiones y un montón de multiplicaciones. Si el recuadro en cuestión intersectaba el eje x pero no intersectaba el eje y, podría rotar todo en pi/2 y hacer el cálculo en el costo original.
Lo estoy usando como una referencia para depurar el rasterizador de círculos con antialias con precisión de subpíxel. Es lento como el infierno :), necesito calcular el área de intersección de cada píxel en el cuadro delimitador del círculo con el círculo para obtener alfa. Puedes ver que funciona y no parece que aparezcan artefactos numéricos. La siguiente figura es una gráfica de un grupo de círculos con un radio que aumenta de 0,3 px a aproximadamente 6 px, dispuesto en una espiral.
Quizás pueda usar la respuesta a esta pregunta , donde se pregunta el área de intersección entre un círculo y un triángulo. Divide tu rectángulo en dos triángulos y usa los algoritmos descritos allí.
Otra forma es dibujar una línea entre los dos puntos de intersección, esto divide su área en un polígono (3 o 4 bordes) y un segmento circular , para que ambos puedan encontrar bibliotecas más fácilmente o hacer los cálculos usted mismo.