Domanda

Sto cercando un modo rapido per determinare l'area di intersezione tra un rettangolo e un cerchio (devo fare milioni di questi calcoli).

Una proprietà specifica è che in tutti i casi il cerchio e il rettangolo hanno sempre 2 punti di intersezione.

È stato utile?

Soluzione

Dati 2 punti di intersezione:

0 vertici è all'interno del cerchio: l'area di un segmento circolare

    XXXXX              -------------------
   X     X               X            X Circular segment
  X       X               XX        XX 
+-X-------X--+              XXXXXXXX 
|  X     X   |
|   XXXXX    |

1 vertice è all'interno del cerchio: la somma delle aree di un segmento circolare e un triangolo.

    XXXXX                   XXXXXXXXX
   X     X       Triangle ->X     _-X
  X       X                 X   _-  X 
  X    +--X--+              X _-   X <- Circular segment 
   X   | X   |              X-  XXX 
    XXXXX    |              XXXX
       |     | 

2 vertici sono all'interno del cerchio: la somma dell'area di due triangoli e un segmento circolare

    XXXXX                   +------------X
   X     X                  |      _--'/'X 
  X    +--X---    Triangle->|   _--   / X
  X    |  X                 |_--     /XX <- Circular segment
   X   +-X----              +-------XX
    XXXXX                 Triangle^

3 vertici sono all'interno del cerchio: l'area del rettangolo meno l'area di un triangolo più l'area di un segmento circolare

    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

Per calcolare queste aree:

Altri suggerimenti

Mi rendo conto che è stata data risposta qualche tempo fa, ma sto risolvendo lo stesso problema e non sono riuscito a trovare una soluzione praticabile pronta all'uso che potrei usare. Nota che le mie caselle sono asse allineato , questo non è del tutto specificato dall'OP. La soluzione di seguito è completamente generale e funzionerà per qualsiasi numero di intersezioni (non solo due). Nota che se le tue caselle non sono allineate agli assi (ma sono comunque caselle con angoli retti, piuttosto che quads ), puoi sfruttare i cerchi rotondi, ruotare le coordinate di tutto in modo che la casella finisca in asse e quindi usi questo codice.

Voglio usare l'integrazione - sembra una buona idea. Cominciamo con la scrittura di una formula ovvia per tracciare un cerchio:

x = center.x + cos(theta) * radius
y = center.y + sin(theta) * radius

^
|
|**###        **
| #*  #      *          * x
|#  *  #    *           # y
|#  *  #    *     
+-----------------------> theta
     *  #  *  # 
     *  #  *  #
      *  #*  #
       ***###

È carino, ma non riesco a integrare l'area di quella cerchia su x o y ; quelli sono quantità diverse. Posso solo integrarmi sopra l'angolo theta , producendo aree di fette di pizza. Non quello che voglio. Proviamo a cambiare gli argomenti:

(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

È più simile. Ora data la gamma di x , posso integrarmi su y , per ottenere un'area della metà superiore di un cerchio. Questo vale solo per x in [center.x - radius, center.x + radius] (altri valori causeranno output immaginari) ma sappiamo che l'area al di fuori di tale intervallo è zero, quindi può essere gestito facilmente. Supponiamo che il cerchio unitario sia semplice, possiamo sempre ricollegare il centro e il raggio in un secondo momento:

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

Questa funzione ha effettivamente un integrale di pi / 2 , poiché è una metà superiore di un cerchio unitario (l'area di un semicerchio è pi * r ^ 2/2 e abbiamo già detto unità , che significa r = 1 ). Ora possiamo calcolare l'area di intersezione di un semicerchio e una casella infinitamente alta, in piedi sull'asse x (il centro del cerchio si trova anche sull'asse x) integrando su 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

Questo potrebbe non essere molto utile, poiché le scatole infinitamente alte non sono ciò che vogliamo. Dobbiamo aggiungere un altro parametro per poter liberare il bordo inferiore della casella 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

Dove h è la distanza (positiva) del bordo inferiore della nostra scatola infinita dall'asse x. La funzione section calcola la posizione (positiva) dell'intersezione del cerchio unitario con la linea orizzontale data da y = h e potremmo definirla risolvendo:

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)

Ora possiamo iniziare. Quindi, come calcolare l'area di intersezione di una scatola finita che interseca un cerchio unitario sopra l'asse 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

È carino. Che ne dici di una scatola che non si trova sopra l'asse x? Direi che non tutte le scatole lo sono. Si presentano tre casi semplici:

  • la casella si trova sopra l'asse x (utilizzare l'equazione sopra)
  • la casella si trova sotto l'asse x (capovolgi il segno delle coordinate y e usa l'equazione sopra)
  • la casella interseca l'asse x (dividere la casella nella metà superiore e inferiore, calcolare l'area di entrambi usando quanto sopra e sommarli)

Abbastanza facile? Scriviamo un po 'di codice:

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 alcuni test unitari di 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

L'output di questo è:

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

Il che mi sembra corretto. Potresti voler incorporare le funzioni se non ti fidi del compilatore che lo fa per te.

Questo utilizza 6 sqrt, 4 asin, 8 div, 16 mul e 17 aggiunte per le scatole che non intersecano l'asse x e un doppio di quello (e 1 aggiunta in più) per le scatole che lo fanno. Si noti che le divisioni sono per raggio e potrebbero essere ridotte a due divisioni e un gruppo di moltiplicazioni. Se la casella in questione interseca l'asse x ma non interseca l'asse y, è possibile ruotare tutto di pi / 2 ed eseguire il calcolo nel costo originale.

Lo sto usando come riferimento al rasterizzatore di cerchi antialiasing con precisione subpixel. È lento come l'inferno :), ho bisogno di calcolare l'area di intersezione di ciascun pixel nel riquadro di delimitazione del cerchio con il cerchio per ottenere l'alfa. Puoi vedere che funziona e che non sembrano apparire artefatti numerici. La figura sotto è una trama di un gruppo di cerchi con il raggio che aumenta da 0,3 px a circa 6 px, disposti a spirale.

 cerchi

Spero che non sia una cattiva forma pubblicare una risposta a una domanda così vecchia. Ho esaminato le soluzioni di cui sopra e ho elaborato un algoritmo simile alla prima risposta di Daniels, ma un po 'più stretto.

In breve, supponiamo che l'intera area sia nel rettangolo, sottrarre i quattro segmenti nei semipiani esterni, quindi aggiungere qualsiasi area nei quattro quadranti esterni, scartando casi banali lungo la strada.

pseudocde (il mio codice attuale è solo ~ 12 righe ..)

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

Per inciso, l'ultima formula per l'area di un cerchio contenuta in un quadrante piano viene prontamente derivata come la somma di un segmento circolare, due triangoli rettangoli e un rettangolo.

Enjoy.

Di seguito viene spiegato come calcolare l'area di sovrapposizione tra cerchio e rettangolo in cui il centro del cerchio si trova all'esterno del rettangolo. Altri casi possono essere ridotti a questo problema.

L'area può essere calcolata integrando l'equazione del cerchio y = sqrt [a ^ 2 - (xh) ^ 2] + k dove a è raggio, (h, k) è centro del cerchio, per trovare l'area sotto la curva. È possibile utilizzare l'integrazione del computer in cui l'area è divisa in molti piccoli rettangoli e calcolarne la somma, oppure utilizzare semplicemente il modulo chiuso qui.

alt text

Ed ecco una fonte C # che implementa il concetto sopra. Si noti che esiste un caso speciale in cui la x specificata si trova al di fuori dei limiti del cerchio. Uso semplicemente una soluzione alternativa qui (che non produce le risposte corrette in tutti i casi)

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: questo problema è molto simile a quello in Google Code Jam Turno di qualificazione 2008 problema: Fly Swatter . Puoi fare clic sui collegamenti dei punteggi per scaricare anche il codice sorgente della soluzione.

Grazie per le risposte,

Ho dimenticato di menzionare che le stime dell'area erano sufficienti. Quello; Ecco perché alla fine, dopo aver esaminato tutte le opzioni, sono andato con la stima monte-carlo in cui ho generato punti casuali nel cerchio e testato se sono nella casella.

Nel mio caso questo è probabilmente più performante. (Ho una griglia posizionata sopra il cerchio e devo misurare il rapporto del cerchio appartenente a ciascuna delle celle della griglia.)

Grazie

Forse puoi usare la risposta a questa domanda , in cui viene posta l'area di intersezione tra un cerchio e un triangolo. Dividi il rettangolo in due triangoli e usa gli algoritmi qui descritti.

Un altro modo è quello di tracciare una linea tra i due punti di intersezione, questo divide la tua area in un poligono (3 o 4 bordi) e un segmento circolare , per entrambi dovresti essere in grado di trovare le biblioteche più facilmente o fare i calcoli da solo.

Ecco un'altra soluzione al 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)));
}
Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top