Área de interseção entre círculo e retângulo
https://stackoverflow.com/questions/622287
-
05-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPergunta
Estou procurando uma maneira rápida de determinar a área de interseção entre um retângulo e um círculo (preciso fazer milhões desses cálculos).
Uma propriedade específica é que, em todos os casos, o círculo e o retângulo sempre têm 2 pontos de interseção.
Solução
Dado 2 pontos de interseção:
0 vértices está dentro do círculo: a área de um segmento circular
XXXXX -------------------
X X X X Circular segment
X X XX XX
+-X-------X--+ XXXXXXXX
| X X |
| XXXXX |
1 vértice está dentro do círculo: a soma das áreas de um segmento circular e um 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ão dentro do círculo: a soma da área de dois triângulos e um 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ão dentro do círculo: a área do retângulo menos a área de um triângulo mais a área de um 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:
A maioria dos pontos que você precisa usar pode ser encontrada ao encontrar o interseção de uma linha e um círculo
As áreas que você precisa calcular podem ser encontradas por computando a área de um segmento circular e computando a área de um triângulo.
Você pode determinar se um vértice está dentro do círculo calculando se sua distância do centro for menor que o raio.
Outras dicas
Sei que isso foi respondido há um tempo atrás, mas estou resolvendo o mesmo problema e não consegui encontrar uma solução viável fora da caixa que eu poderia usar. Observe que minhas caixas são eixo alinhado, isso não é especificado pelo OP. A solução abaixo é completamente geral e funcionará para qualquer número de cruzamentos (não apenas dois). Observe que se suas caixas não são alinhadas ao eixo (mas ainda as caixas com ângulos retos, em vez de geral Quads), você pode aproveitar os círculos redondos, girar as coordenadas de tudo para que a caixa acabe alinhada e alinhada e então Use este código.
Quero usar a integração - isso parece uma boa ideia. Vamos começar a escrever uma fórmula óbvia para planejar um círculo:
x = center.x + cos(theta) * radius
y = center.y + sin(theta) * radius
^
|
|**### **
| #* # * * x
|# * # * # y
|# * # *
+-----------------------> theta
* # * #
* # * #
* #* #
***###
Isso é bom, mas não consigo integrar a área daquele círculo sobre x ou y; Essas são quantidades diferentes. Eu só posso integrar o ângulo theta, produzindo áreas de fatias de pizza. Não o que eu quero. Vamos tentar mudar os 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
É mais parecido. Agora dado o alcance de x, Eu posso integrar y, para obter uma área da metade superior de um círculo. Isso só vale para x dentro [center.x - radius, center.x + radius] (Outros valores causarão saídas imaginárias), mas sabemos que a área fora desse intervalo é zero, de modo que é tratado facilmente. Vamos assumir o círculo unitário para simplificar, sempre podemos conectar o centro e o raio de volta mais tarde:
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 função realmente tem uma parte integrante de pi/2, uma vez que é uma metade superior de um círculo unitário (a área de meio círculo é pi * r^2 / 2 E já dissemos unidade, que significa r = 1). Agora, podemos calcular a área de interseção de um meio círculo e uma caixa infinitamente alta, de pé no eixo X (o centro do círculo também está no eixo X) integrando-se 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
Isso pode não ser muito útil, pois caixas infinitamente altas não são o que queremos. Precisamos adicionar mais um parâmetro para poder liberar a borda inferior da caixa infinitamente alta:
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
Onde h é a distância (positiva) da borda inferior da nossa caixa infinita do eixo x. o section A função calcula a posição (positiva) da interseção do círculo unitário com a linha horizontal dada por y = h E poderíamos defini -lo resolvendo:
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)
Agora podemos fazer as coisas que vão. Então, como calcular a área de interseção de uma caixa finita que cruzava um círculo unitário acima do eixo 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
Muito legal. Então, que tal uma caixa que não está acima do eixo x? Eu diria que nem todas as caixas são. Três casos simples surgem:
- A caixa está acima do eixo x (use a equação acima)
- A caixa está abaixo do eixo x (gire o sinal de coordenadas y e use a equação acima)
- A caixa está cruzando o eixo x (divida a caixa na metade superior e inferior, calcule a área de ambos usando o acima e soma -os)
Bastante fácil? Vamos escrever algum 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);
}
E alguns testes básicos de unidade:
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
A saída disso é:
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
O que me parece correto. Você pode embrulhar as funções se não confiar no seu compilador para fazer isso por você.
Isso usa 6 sqrt, 4 asin, 8 div, 16 mul e 17 adicionam para caixas que não cruzam o eixo x e um dobro disso (e 1 mais add) para caixas que o fazem. Observe que as divisões são por raio e podem ser reduzidas para duas divisões e um monte de multiplicações. Se a caixa em questão cruzou o eixo x, mas não cruzou o eixo y, você poderia girar tudo por pi/2 e faça o cálculo no custo original.
Estou usando-o como uma referência para depixar o Subpixel-Precise AntialieSed Circle Rasterizer. É lento como o inferno :), preciso calcular a área de interseção de cada pixel na caixa delimitadora do círculo com o círculo para obter alfa. Você pode ver que ele funciona e nenhum artefato numérico parece aparecer. A figura abaixo está um gráfico de um monte de círculos com o raio aumentando de 0,3 px para cerca de 6 px, dispostos em espiral.
Espero que não seja uma forma ruim postar uma resposta para uma pergunta tão antiga. Eu olhei as soluções acima e trabalhei um algoritmo que é semelhante à primeira resposta de Daniels, mas um pouco mais apertado.
Em resumo, suponha que a área completa esteja no retângulo, subtraia os quatro segmentos nos meios planos externos e adicione todas as áreas nos quatro quadrantes externos, descartando casos triviais ao longo do caminho.
pseudocde (meu código real é apenas ~ 12 linhas ..)
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
Aliás, a última fórmula para a área de um círculo contida em um quadrante de avião é prontamente derivada como a soma de um segmento circular, dois triângulos direito e um retângulo.
Apreciar.
A seguir, é como calcular a área de sobreposição entre o círculo e o retângulo, onde o centro do círculo fica fora do retângulo. Outros casos podem ser reduzidos a esse problema.
A área pode ser calculada integrando a equação do círculo y = sqrt [a^2 - (xh)^2] + k onde a é raio, (h, k) é o Circle Center, para encontrar a área sob curva. Você pode usar a integração do computador onde a área é dividida em muitos pequenos retângulo e calculando a soma deles ou apenas usar o formulário fechado aqui.
E aqui está uma fonte C# implementando o conceito acima. Observe que há um caso especial em que o X especificado está fora dos limites do círculo. Eu apenas uso uma solução alternativa simples aqui (que não está produzindo as respostas corretas em todos os 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;
}
Observação: Este problema é muito semelhante a um em Google Code Jam 2008 Rodada de qualificação problema: Mata-moscas. Você também pode clicar nos links de pontuação para baixar o código -fonte da solução.
Obrigado pelas respostas,
Esqueci de mencionar que as estimativas de área foram suficientes. Este; Por que, no final, depois de olhar para todas as opções, fui com a estimativa de Monte-Carlo, onde gero pontos aleatórios no círculo e testei se eles estão na caixa.
No meu caso, isso provavelmente é mais desempenho. (Eu tenho uma grade colocada sobre o círculo e tenho que medir a proporção do círculo pertencente a cada uma das células da grade.)
Obrigado
Talvez você possa usar a resposta para essa questão, onde é perguntada a área de interseção entre um círculo e um triângulo. Divida o retângulo em dois triângulos e use os algoritmos descritos lá.
Outra maneira é desenhar uma linha entre os dois pontos de cruzamento, isso divide sua área em um polígono (3 ou 4 arestas) e um segmento circular, para ambos, você deve encontrar bibliotecas com mais facilidade ou fazer os cálculos.
Aqui está outra solução para o 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)));
}
