2009-02-12 12 views
18

Comment calculer l'aire d'intersection entre un triangle (spécifié comme trois (X, Y) paires) et un cercle (X, Y, R)? J'ai fait quelques recherches en vain. C'est pour le travail, pas pour l'école. :)Calculer la zone d'intersection entre un cercle et un triangle?

Il ressemblerait à quelque chose comme ça en C#:

struct { PointF vert[3]; } Triangle; 
struct { PointF center; float radius; } Circle; 

// returns the area of intersection, e.g.: 
// if the circle contains the triangle, return area of triangle 
// if the triangle contains the circle, return area of circle 
// if partial intersection, figure that out 
// if no intersection, return 0 
double AreaOfIntersection(Triangle t, Circle c) 
{ 
... 
} 

Répondre

27

Si vous voulez une solution exacte (ou au moins aussi précis que vous pouvez obtenir en utilisant l'arithmétique à virgule flottante), alors cela va impliquer beaucoup de travail, car il y a tellement de cas à considérer.

je compte neuf cas différents (classés dans la figure ci-dessous par le nombre de sommets du triangle à l'intérieur du cercle, et le nombre d'arêtes du triangle qui se croisent ou sont contenus dans le cercle):

Nine cases for intersection: 1, 2. no vertices, no edges; 3. no vertices, one edge; 4. no vertices, two edges; 5. no vertices, three edges; 6. one vertex, two edges; 7. one vertex, three edges; 8. two vertices, three edges; 9. three vertices, three edges.

(Cependant, ce genre d'énumération des cas géométriques est bien connu pour être difficile, et il ne me surprendrait pas du tout si je manqué un ou deux!)

ainsi, l'approche est la suivante:

  1. Déterminez pour chaque sommet du triangle s'il se trouve à l'intérieur du cercle. Je vais supposer que vous savez comment faire cela.

  2. Déterminez pour chaque arête du triangle si elle coupe le cercle. (J'ai écrit une méthode here, ou voir n'importe quel livre de géométrie de calcul.) Vous devrez calculer le point ou les points d'intersection (le cas échéant) pour l'utilisation à l'étape 4.

  3. Déterminez lequel des neuf cas vous avoir.

  4. Calculez la zone de l'intersection. Les cas 1, 2 et 9 sont faciles. Dans les six cas restants, j'ai dessiné des lignes pointillées pour montrer comment diviser la zone d'intersection en triangles et circular segments sur la base des sommets d'origine du triangle, et sur les points d'intersection que vous avez calculés à l'étape 2.

Cet algorithme va être assez délicate et sujette à des erreurs qui affectent seulement l'un des cas, alors assurez-vous que vous avez des cas de test qui couvrent l'ensemble des neuf cas (et je suggère permutant les sommets des triangles de test aussi). Portez une attention particulière aux cas où l'un des sommets du triangle est sur le bord du cercle. Si vous n'avez pas besoin d'une solution exacte, alors pixelliser les chiffres et compter les pixels dans l'intersection (comme suggéré par quelques autres répondants) semble une approche beaucoup plus facile au code, et par conséquent moins sujette aux erreurs .

+0

+1 maths! On dirait que la solution exacte fonctionnerait beaucoup plus vite qu'une techinque rastérisée aussi. – Crashworks

+1

Je suis dûment impressionné par votre rigueur. –

+0

Notez que le moyen le plus simple de faire # 4 et # 5 est de prendre la zone du cercle et de soustraire les segments à l'extérieur du triangle (plutôt que d'additionner tous les soustriangles et segments à l'intérieur). Je suis vraiment impressionné, Gareth. – Crashworks

1

En supposant que vous parlez pixels entiers, pas réel, la mise en œuvre naïve serait de boucle à travers chaque pixel du triangle et vérifier la distance du centre du cercle contre son rayon.

Ce n'est pas une formule mignonne, ou particulièrement rapide, mais elle fait le travail.

1

essayer computational geometry

Note: ce n'est pas un problème trivial, j'espère que ce n'est pas devoir ;-)

0

Comment exact que vous devez être? Si vous pouvez approximer le cercle avec des formes plus simples, vous pouvez simplifier le problème. Il ne serait pas difficile de modéliser un cercle comme un ensemble de triangles très étroits se rencontrant au centre, par exemple.

1

Si vous avez un processeur graphique à votre disposition, vous pouvez utiliser this technique pour obtenir un nombre de pixels de l'intersection ..

0

Si seulement l'un des segments de ligne du triangle coupe le cercle, la solution pure mathématique ISN » T trop dur. Une fois que vous savez quand les deux points d'intersection sont, vous pouvez utiliser la formule de distance pour trouver la longueur de corde.

Selon these equations:

ϑ = 2 sin⁻¹(0.5 c/r) 
A = 0.5 r² (ϑ - sin(ϑ)) 

où c est la longueur de la corde, r est le rayon, θ devient l'angle au centre, et A est la surface. Notez que cette solution casse si plus de la moitié du cercle est coupée.

Cela ne vaut probablement pas la peine si vous avez juste besoin d'une approximation, car elle fait plusieurs suppositions sur ce à quoi ressemble l'intersection réelle.

1

Je pense que vous ne devriez pas approximer un cercle comme un ensemble de triangles, au lieu de cela vous pouvez approximer sa forme avec un polygone. L'algorithme naïf peut ressembler à:

  1. Convertissez votre cercle en polygone avec un nombre désiré de sommets.
  2. Calculer l'intersection de deux polygones (cercle converti et un triangle).
  3. Calculer le carré de cette intersection.

Vous pouvez optimiser cet algorithme en combinant l'étape 2 et l'étape 3 en fonction unique.

Lire ce lien:
Area of convex polygon
Intersection of convex polygons

0

Mon premier instinct serait de tout transformer de telle sorte que le cercle est centré sur l'origine, trans triangle en coordonnées polaires, et résoudre l'intersection (ou englobement) du triangle avec le cercle. Je n'ai pas encore travaillé sur papier, mais ce n'est qu'une intuition.

+0

Je suis en train d'examiner cette approche en ce moment ... dans le cas général, il y a une intégration assez laide impliquée. Je ne pense pas qu'il y aura une bonne formule simple qu'un ordinateur peut calculer. –

+2

Cela ressemble au genre de chose qui a dû être travaillé par un mathématicien du 19ème siècle, mais malheureusement, Google Scholar ne retourne pas si loin! =) – Crashworks

1

Étant donné que vos formes sont convexes, vous pouvez utiliser l'estimation de la surface de Monte Carlo. Dessinez une boîte autour du cercle et du triangle.

Choisissez des points aléatoires dans la boîte et tenez compte du nombre de tombes dans le cercle et du nombre de points dans le cercle et le triangle.

zone d'intersection ≅ Zone de cercle * # points en cercle et triangle/# points cercle

points d'arrêt lorsque le choix de la zone estimée ne change pas de plus d'un certain montant sur un certain nombre de tours , ou simplement choisir un nombre fixe de points en fonction de la surface de la boîte. L'estimation de la zone devrait converger assez rapidement, sauf si l'une de vos formes a très peu de surface.

Note: Voici comment déterminer si un point est dans un triangle: Barycentric coordinates

2

J'ai presque un an et demi de retard, mais je pensais que les gens seraient peut-être intéressés par code here que j'ai écrit et que je pense faire correctement. Regardez dans la fonction IntersectionArea près du fond. L'approche générale est de décoller le polygone convexe circonscrit par le cercle, puis de traiter les petits bouchons circulaires.

31

D'abord, je vais nous rappeler comment trouver l'aire d'un polygone. Une fois que nous avons fait cela, l'algorithme pour trouver l'intersection entre un polygone et un cercle devrait être facile à comprendre.

Comment trouver la surface d'un polygone

Examinons le cas d'un triangle, parce que toute la logique essentielle y apparaît. Supposons que nous avons un triangle dont les sommets (x1, y1), (x2, y2) et (x3, y3) que vous allez dans le sens antihoraire triangle, comme le montre la figure suivante: triangleFigure

Ensuite, vous pouvez calculer la surface par la formule

A = (x1 y2 + x2 y3 + x3 y1 - x2y1 - x3 y2 - x1y3)/2.

Pour comprendre pourquoi cette formule fonctionne, nous allons réarranger il est donc sous la forme

A = (x1 y2 - x2 y1)/2 + (x2 Y3 - x3 y2)/2 + (x3 y1 - x1y3)/2.

Maintenant, le premier terme est la zone suivante, ce qui est positif dans notre cas: enter image description here

Si on ne sait pas que la zone de la région verte est en effet (x1 y2 - x2 y1)/2 , puis lisez this.

Le deuxième terme est ce domaine, qui est à nouveau positif:

enter image description here

Et la troisième zone est représentée dans la figure suivante. Cette fois, la région est négative

enter image description here

L'ajout de ces trois jusqu'à nous obtenons l'image suivante

enter image description here

On voit que la zone verte qui était en dehors du triangle est annulé par la zone rouge , de sorte que la zone nette est juste la zone du triangle, et cela montre pourquoi notre formule était vraie dans ce cas.

Ce que j'ai dit ci-dessus était l'explication intuitive de la raison pour laquelle la formule de surface était correcte. Une explication plus rigoureuse consisterait à observer que lorsque nous calculons l'aire à partir d'un bord, la zone que nous obtenons est la même que celle que nous obtiendrions de l'intégration r^2dθ/2, donc nous intégrons efficacement r^2dθ/2 autour de la limite du polygone, et par le théorème de Stokes, cela donne le même résultat que l'intégration de rdrdθ sur la région délimitée par le polygone. Puisque l'intégration de rdrdθ sur la région délimitée par le polygone donne la surface, nous concluons que notre procédure doit correctement donner la zone.

zone de l'intersection d'un cercle avec un polygone

Maintenant, nous allons discuter de la façon de trouver la zone de l'intersection d'un cercle de rayon R avec un polygone comme indiqué dans la figure suivante:

enter image description here

Nous sommes intéressés à trouver la zone de la région verte. Nous pouvons, comme dans le cas du polygone unique, casser notre calcul en trouvant une zone pour chaque côté du polygone, puis ajouter ces zones.

Notre première zone ressemblera: enter image description here

La deuxième zone ressemblera enter image description here

Et la troisième zone sera enter image description here

Encore une fois, les deux premiers domaines sont positifs notre cas alors que le troisième sera négatif. Espérons que les annulations fonctionneront de sorte que la zone nette est en effet la zone qui nous intéresse. Voyons voir.

enter image description here

En effet, la somme des surfaces sera zone qui nous intéresse.

Encore une fois, nous pouvons donner une explication plus rigoureuse des raisons pour lesquelles cela fonctionne. Soit I la région définie par l'intersection et soit P le polygone. Ensuite, à partir de la discussion précédente, nous savons que nous voulons calculer l'intégrale de r^2dθ/2 autour de la limite de I. Cependant, cela est difficile à faire parce qu'il faut trouver l'intersection. Au lieu de cela, nous avons fait une intégrale sur le polygone. Nous avons intégré max (r, R)^2 dθ/2 sur la limite du polygone. Pour voir pourquoi cela donne la bonne réponse, définissons une fonction π qui prend un point en coordonnées polaires (r, θ) au point (max (r, R), θ). Il ne devrait pas être déroutant de se référer aux fonctions de coordonnées de π (r) = max (r, R) et π (θ) = θ. Nous avons ensuite intégré π (r)^2 dθ/2 sur la limite du polygone. Par contre, puisque π (θ) = θ, cela revient à intégrer π (r)^2 dπ (θ)/2 sur la limite du polygone. Maintenant, en faisant un changement de variable, nous obtenons la même réponse si nous avons intégré r^2 dθ/2 sur la limite de π (P), où π (P) est l'image de P sous π . En utilisant encore le théorème de Stokes, nous savons que l'intégration de r^2 dθ/2 sur la limite de π (P) nous donne l'aire de π (P). En d'autres termes, il donne la même réponse que l'intégration de dxdy sur π (P). En utilisant à nouveau un changement de variable, nous savons que l'intégration de dxdy sur π (P) revient à intégrer Jdxdy sur P, où J est le jacobien de π.

Maintenant, nous pouvons diviser l'intégrale de Jdxdy en deux régions: la partie dans le cercle et la partie à l'extérieur du cercle. Maintenant, π laisse des points dans le cercle seul, J = 1 là, donc la contribution de cette partie de P est la surface de la partie de P qui se trouve dans le cercle, c'est-à-dire la zone de l'intersection. La deuxième région est la région en dehors du cercle. Là J = 0 puisque π effondre cette partie jusqu'à la limite du cercle.

Ainsi, ce que nous calculons est en effet la zone de l'intersection.

Maintenant que nous sommes relativement sûrs que nous savons conceptuellement comment trouver la zone, parlons plus spécifiquement de la façon de calculer la contribution d'un seul segment. Commençons par regarder un segment dans ce que j'appellerai «géométrie standard». Il est montré ci-dessous.

enter image description here

En géométrie standard, le bord horizontalement en va de gauche à droite. Il est décrit par trois nombres: xi, la coordonnée x où le bord commence, xf, la coordonnée x où le bord se termine, et y, la coordonnée y du bord.

Maintenant nous voyons que si | y | R, comme sur la figure, le bord intersectera le cercle aux points (-xint, y) et (xint, y) où xint = (R^2-y^2)^(1/2). Ensuite, la zone que nous devons calculer est divisée en trois parties étiquetées sur la figure. Pour obtenir les aires des régions 1 et 3, on peut utiliser arctan pour obtenir les angles des différents points et ensuite assimiler la surface à R^2 Δθ/2. Ainsi, par exemple, nous définirons θi = atan2 (y, xi) et θl = atan2 (y, -xint). Ensuite, la zone de la région 1 est R^2 (θl-θi)/2. Nous pouvons obtenir la zone de la région 3 de manière similaire.

La zone de la région 2 est juste la zone d'un triangle. Cependant, nous devons faire attention au signe. Nous voulons que la zone montrée soit positive donc nous dirons que la zone est - (xint - (-xint)) y/2.

Une autre chose à garder à l'esprit est qu'en général, xi ne doit pas être inférieur à -xint et xf ne doit pas être supérieur à xint.

L'autre cas à prendre en compte est | y | > R. Ce cas est plus simple, car il n'y a qu'une seule pièce semblable à la région 1 de la figure.

Maintenant que nous savons comment calculer la surface à partir d'une arête en géométrie standard, il ne reste plus qu'à décrire comment transformer une arête quelconque en géométrie standard.

Mais ce n'est qu'un simple changement de coordonnées. Etant donné le sommet initial vi et le sommet final vf, le nouveau vecteur unitaire sera le vecteur unitaire de vi à vf. Alors xi est juste le déplacement de vi du centre du cercle pointé dans x, et xf est juste xi plus la distance entre vi et vf. Pendant ce temps, y est donné par le produit de coin de x avec le déplacement de vi du centre du cercle.

code

qui complète la description de l'algorithme, maintenant il est temps d'écrire un code. Je vais utiliser Java.

Tout d'abord, puisque nous travaillons avec des cercles, nous devrions avoir une classe de cercle

public class Circle { 

    final Point2D center; 
    final double radius; 

    public Circle(double x, double y, double radius) { 
     center = new Point2D.Double(x, y); 
     this.radius = radius; 
    } 

    public Circle(Point2D.Double center, double radius) { 
     this(center.getX(), center.getY(), radius); 
    } 

    public Point2D getCenter() { 
     return new Point2D.Double(getCenterX(), getCenterY()); 
    } 

    public double getCenterX() { 
     return center.getX(); 
    } 

    public double getCenterY() { 
     return center.getY(); 
    } 

    public double getRadius() { 
     return radius; 
    } 

} 

Pour les polygones, je vais utiliser la classe de java Shape. Shape s j'ai un PathIterator que je peux utiliser pour itérer à travers les bords du polygone.

Maintenant pour le travail actuel. Je vais séparer la logique d'itération à travers les bords, en mettant les arêtes en géométrie standard, etc., à partir de la logique de calcul de la surface une fois que cela est fait. La raison en est que vous pouvez dans l'avenir vouloir calculer quelque chose d'autre en plus ou en plus de la zone et que vous voulez pouvoir réutiliser le code ayant à traiter avec les itérations. J'ai donc une classe générique qui calcule une propriété de la classe T à propos de notre intersection de polygones.

public abstract class CircleShapeIntersectionFinder<T> {

Il dispose de trois méthodes statiques qui qu'aider Compute la géométrie:

private static double[] displacment2D(final double[] initialPoint, final double[] finalPoint) { 
    return new double[]{finalPoint[0] - initialPoint[0], finalPoint[1] - initialPoint[1]}; 
} 

private static double wedgeProduct2D(final double[] firstFactor, final double[] secondFactor) { 
    return firstFactor[0] * secondFactor[1] - firstFactor[1] * secondFactor[0]; 
} 

static private double dotProduct2D(final double[] firstFactor, final double[] secondFactor) { 
    return firstFactor[0] * secondFactor[0] + firstFactor[1] * secondFactor[1]; 
} 

Il y a deux champs d'instance, un Circle qui maintient juste une copie du cercle, et le currentSquareRadius, qui maintient une copie du rayon carré. Cela peut sembler étrange, mais la classe que j'utilise est en fait équipée pour trouver les zones d'une collection complète d'intersections cercle-polygone. C'est pourquoi je fais référence à l'un des cercles comme «courant».

private Circle currentCircle; 
private double currentSquareRadius; 

Vient ensuite la méthode de calcul de ce que nous voulons calculer:

public final T computeValue(Circle circle, Shape shape) { 
    initialize(); 
    processCircleShape(circle, shape); 
    return getValue(); 
} 

initialize() et getValue() sont abstraites. initialize() définirait la variable qui maintient un total de la zone à zéro, et getValue() retournerait simplement la zone. La définition de processCircleShape est

private void processCircleShape(Circle circle, final Shape cellBoundaryPolygon) { 
    initializeForNewCirclePrivate(circle); 
    if (cellBoundaryPolygon == null) { 
     return; 
    } 
    PathIterator boundaryPathIterator = cellBoundaryPolygon.getPathIterator(null); 
    double[] firstVertex = new double[2]; 
    double[] oldVertex = new double[2]; 
    double[] newVertex = new double[2]; 
    int segmentType = boundaryPathIterator.currentSegment(firstVertex); 
    if (segmentType != PathIterator.SEG_MOVETO) { 
     throw new AssertionError(); 
    } 
    System.arraycopy(firstVertex, 0, newVertex, 0, 2); 
    boundaryPathIterator.next(); 
    System.arraycopy(newVertex, 0, oldVertex, 0, 2); 
    segmentType = boundaryPathIterator.currentSegment(newVertex); 
    while (segmentType != PathIterator.SEG_CLOSE) { 
     processSegment(oldVertex, newVertex); 
     boundaryPathIterator.next(); 
     System.arraycopy(newVertex, 0, oldVertex, 0, 2); 
     segmentType = boundaryPathIterator.currentSegment(newVertex); 
    } 
    processSegment(newVertex, firstVertex); 
} 

Prenons une seconde pour regarder initializeForNewCirclePrivate rapidement. Cette méthode définit simplement les champs d'instance et permet à la classe dérivée de stocker n'importe quelle propriété du cercle. Sa définition est

private void initializeForNewCirclePrivate(Circle circle) { 
    currentCircle = circle; 
    currentSquareRadius = currentCircle.getRadius() * currentCircle.getRadius(); 
    initializeForNewCircle(circle); 
} 

initializeForNewCircle est abstraite et une mise en œuvre serait pour elle de stocker le rayon des cercles pour éviter d'avoir à faire des racines carrées. Quoi qu'il en soit, revenez à processCircleShape. Après avoir appelé initializeForNewCirclePrivate, nous vérifions si le polygone est null (que j'interprète comme un polygone vide), et nous retournons s'il est null. Dans ce cas, notre zone calculée serait zéro. Si le polygone n'est pas null alors nous obtenons le PathIterator du polygone. L'argument de la méthode getPathIterator que j'appelle est une transformation affine qui peut être appliquée au chemin. Je ne veux pas en appliquer un, alors je passe juste null.

Ensuite, je déclare le double[] s qui gardera la trace des sommets. Je dois me souvenir du premier sommet parce que le PathIterator ne me donne qu'une seule fois chaque sommet, donc je dois revenir en arrière après m'avoir donné le dernier sommet, et former un bord avec ce dernier sommet et le premier sommet.

La méthode currentSegment sur la ligne suivante place le sommet suivant dans son argument. Il renvoie un code qui vous indique quand il est hors des sommets. C'est pourquoi l'expression de contrôle de ma boucle while est ce qu'elle est.

La majeure partie du reste du code de cette méthode est une logique inintéressante liée à l'itération à travers les sommets.L'important est qu'une fois par itération de la boucle while j'appelle processSegment puis j'appelle processSegment à la fin de la méthode pour traiter le bord qui relie le dernier sommet au premier sommet.

Regardons le code pour processSegment:

private void processSegment(double[] initialVertex, double[] finalVertex) { 
    double[] segmentDisplacement = displacment2D(initialVertex, finalVertex); 
    if (segmentDisplacement[0] == 0 && segmentDisplacement[1] == 0) { 
     return; 
    } 
    double segmentLength = Math.sqrt(dotProduct2D(segmentDisplacement, segmentDisplacement)); 
    double[] centerToInitialDisplacement = new double[]{initialVertex[0] - getCurrentCircle().getCenterX(), initialVertex[1] - getCurrentCircle().getCenterY()}; 
    final double leftX = dotProduct2D(centerToInitialDisplacement, segmentDisplacement)/segmentLength; 
    final double rightX = leftX + segmentLength; 
    final double y = wedgeProduct2D(segmentDisplacement, centerToInitialDisplacement)/segmentLength; 
    processSegmentStandardGeometry(leftX, rightX, y); 
} 

Dans cette méthode, je mets en œuvre les étapes pour transformer un bord dans la géométrie standard comme décrit ci-dessus. D'abord je calcule segmentDisplacement, le déplacement du sommet initial au sommet final. Ceci définit l'axe des x de la géométrie standard. Je fais un retour anticipé si ce déplacement est nul. Ensuite, je calcule la longueur du déplacement, car cela est nécessaire pour obtenir le vecteur unité x. Une fois que j'ai cette information, je calcule le déplacement du centre du cercle au sommet initial. Le produit scalaire de ceci avec segmentDisplacement me donne leftX que j'avais appelé xi. Puis rightX, que j'avais appelé xf, est juste leftX + segmentLength. Enfin, je fais le produit coin pour obtenir y comme décrit ci-dessus.

Maintenant que j'ai transformé le problème en géométrie standard, ce sera facile à gérer. C'est ce que fait la méthode processSegmentStandardGeometry. Regardons le code

private void processSegmentStandardGeometry(double leftX, double rightX, double y) { 
    if (y * y > getCurrentSquareRadius()) { 
     processNonIntersectingRegion(leftX, rightX, y); 
    } else { 
     final double intersectionX = Math.sqrt(getCurrentSquareRadius() - y * y); 
     if (leftX < -intersectionX) { 
      final double leftRegionRightEndpoint = Math.min(-intersectionX, rightX); 
      processNonIntersectingRegion(leftX, leftRegionRightEndpoint, y); 
     } 
     if (intersectionX < rightX) { 
      final double rightRegionLeftEndpoint = Math.max(intersectionX, leftX); 
      processNonIntersectingRegion(rightRegionLeftEndpoint, rightX, y); 
     } 
     final double middleRegionLeftEndpoint = Math.max(-intersectionX, leftX); 
     final double middleRegionRightEndpoint = Math.min(intersectionX, rightX); 
     final double middleRegionLength = Math.max(middleRegionRightEndpoint - middleRegionLeftEndpoint, 0); 
     processIntersectingRegion(middleRegionLength, y); 
    } 
} 

Le premier if distingue les cas où y est assez petit pour que le bord peut recouper le cercle. Si y est grand et il n'y a aucune possibilité d'intersection, alors j'appelle la méthode pour gérer ce cas. Sinon, je gère le cas où l'intersection est possible.

Si l'intersection est possible, je calcule la coordonnée x de l'intersection, intersectionX, et je divise le bord en trois parties, qui correspondent aux régions 1, 2 et 3 de la figure de géométrie standard ci-dessus. Je gère d'abord la région 1.

Pour gérer la région 1, je vérifie si leftX est effectivement inférieur à -intersectionX sinon il n'y aurait pas de région 1. S'il y a une région 1, alors j'ai besoin de savoir quand elle se termine. Il se termine au minimum rightX et -intersectionX. Après avoir trouvé ces coordonnées x, je traite de cette région sans intersection.

je fais une chose semblable à gérer la région 3.

Pour la région 2, je dois faire une certaine logique pour vérifier que leftX et rightX ne fait suPPoRt une région entre -intersectionX et intersectionX. Après avoir trouvé la région, je ne ai besoin de la longueur de la région et y, donc je passe ces deux chiffres à une méthode abstraite qui gère la région 2.

Maintenant, regardons le code pour processNonIntersectingRegion

private void processNonIntersectingRegion(double leftX, double rightX, double y) { 
    final double initialTheta = Math.atan2(y, leftX); 
    final double finalTheta = Math.atan2(y, rightX); 
    double deltaTheta = finalTheta - initialTheta; 
    if (deltaTheta < -Math.PI) { 
     deltaTheta += 2 * Math.PI; 
    } else if (deltaTheta > Math.PI) { 
     deltaTheta -= 2 * Math.PI; 
    } 
    processNonIntersectingRegion(deltaTheta); 
} 
J'utilise simplement atan2 pour calculer la différence d'angle entre leftX et rightX. Ensuite, j'ajoute du code pour traiter la discontinuité dans atan2, mais cela est probablement inutile, car la discontinuité se produit soit à 180 degrés soit à 0 degré. Ensuite, je passe la différence d'angle sur une méthode abstraite.Enfin, nous avons juste des méthodes et getters abstraites:

protected abstract void initialize(); 

    protected abstract void initializeForNewCircle(Circle circle); 

    protected abstract void processNonIntersectingRegion(double deltaTheta); 

    protected abstract void processIntersectingRegion(double length, double y); 

    protected abstract T getValue(); 

    protected final Circle getCurrentCircle() { 
     return currentCircle; 
    } 

    protected final double getCurrentSquareRadius() { 
     return currentSquareRadius; 
    } 

} 

Maintenant regardons la classe extension, CircleAreaFinder

public class CircleAreaFinder extends CircleShapeIntersectionFinder<Double> { 

public static double findAreaOfCircle(Circle circle, Shape shape) { 
    CircleAreaFinder circleAreaFinder = new CircleAreaFinder(); 
    return circleAreaFinder.computeValue(circle, shape); 
} 

double area; 

@Override 
protected void initialize() { 
    area = 0; 
} 

@Override 
protected void processNonIntersectingRegion(double deltaTheta) { 
    area += getCurrentSquareRadius() * deltaTheta/2; 
} 

@Override 
protected void processIntersectingRegion(double length, double y) { 
    area -= length * y/2; 
} 

@Override 
protected Double getValue() { 
    return area; 
} 

@Override 
protected void initializeForNewCircle(Circle circle) { 

} 

}

Il a un champ area de garder une trace de la région. initialize définit la zone sur zéro, comme prévu. Lorsque nous traitons un arête qui ne se recoupe pas, nous augmentons la surface de R^2 Δθ/2 car nous avons conclu que nous devrions le faire plus haut. Pour un bord d'intersection, nous décrémenter la zone par y*length/2. C'est ainsi que les valeurs négatives pour y correspondent à des zones positives, comme nous l'avons décidé.

Maintenant, la meilleure chose est que si nous voulons garder une trace du périmètre, nous n'avons pas à faire beaucoup plus de travail. Je définissais une AreaPerimeter classe:

public class AreaPerimeter { 

    final double area; 
    final double perimeter; 

    public AreaPerimeter(double area, double perimeter) { 
     this.area = area; 
     this.perimeter = perimeter; 
    } 

    public double getArea() { 
     return area; 
    } 

    public double getPerimeter() { 
     return perimeter; 
    } 

} 

et maintenant nous avons juste besoin d'étendre notre classe abstraite en utilisant à nouveau AreaPerimeter comme le type.

public class CircleAreaPerimeterFinder extends CircleShapeIntersectionFinder<AreaPerimeter> { 

    public static AreaPerimeter findAreaPerimeterOfCircle(Circle circle, Shape shape) { 
     CircleAreaPerimeterFinder circleAreaPerimeterFinder = new CircleAreaPerimeterFinder(); 
     return circleAreaPerimeterFinder.computeValue(circle, shape); 
    } 

    double perimeter; 
    double radius; 
    CircleAreaFinder circleAreaFinder; 

    @Override 
    protected void initialize() { 
     perimeter = 0; 
     circleAreaFinder = new CircleAreaFinder(); 
    } 

    @Override 
    protected void initializeForNewCircle(Circle circle) { 
     radius = Math.sqrt(getCurrentSquareRadius()); 
    } 

    @Override 
    protected void processNonIntersectingRegion(double deltaTheta) { 
     perimeter += deltaTheta * radius; 
     circleAreaFinder.processNonIntersectingRegion(deltaTheta); 
    } 

    @Override 
    protected void processIntersectingRegion(double length, double y) { 
     perimeter += Math.abs(length); 
     circleAreaFinder.processIntersectingRegion(length, y); 
    } 

    @Override 
    protected AreaPerimeter getValue() { 
     return new AreaPerimeter(circleAreaFinder.getValue(), perimeter); 
    } 

} 

Nous avons une variable perimeter de garder une trace du périmètre, nous nous souvenons de la valeur du radius pour éviter devoir appeler Math.sqrt beaucoup, et nous déléguer le calcul de la zone à notre CircleAreaFinder. Nous pouvons voir que les formules pour le périmètre sont faciles.

Pour référence est ici le code complet de CircleShapeIntersectionFinder

private static double[] displacment2D(final double[] initialPoint, final double[] finalPoint) { 
     return new double[]{finalPoint[0] - initialPoint[0], finalPoint[1] - initialPoint[1]}; 
    } 

    private static double wedgeProduct2D(final double[] firstFactor, final double[] secondFactor) { 
     return firstFactor[0] * secondFactor[1] - firstFactor[1] * secondFactor[0]; 
    } 

    static private double dotProduct2D(final double[] firstFactor, final double[] secondFactor) { 
     return firstFactor[0] * secondFactor[0] + firstFactor[1] * secondFactor[1]; 
    } 

    private Circle currentCircle; 
    private double currentSquareRadius; 

    public final T computeValue(Circle circle, Shape shape) { 
     initialize(); 
     processCircleShape(circle, shape); 
     return getValue(); 
    } 

    private void processCircleShape(Circle circle, final Shape cellBoundaryPolygon) { 
     initializeForNewCirclePrivate(circle); 
     if (cellBoundaryPolygon == null) { 
      return; 
     } 
     PathIterator boundaryPathIterator = cellBoundaryPolygon.getPathIterator(null); 
     double[] firstVertex = new double[2]; 
     double[] oldVertex = new double[2]; 
     double[] newVertex = new double[2]; 
     int segmentType = boundaryPathIterator.currentSegment(firstVertex); 
     if (segmentType != PathIterator.SEG_MOVETO) { 
      throw new AssertionError(); 
     } 
     System.arraycopy(firstVertex, 0, newVertex, 0, 2); 
     boundaryPathIterator.next(); 
     System.arraycopy(newVertex, 0, oldVertex, 0, 2); 
     segmentType = boundaryPathIterator.currentSegment(newVertex); 
     while (segmentType != PathIterator.SEG_CLOSE) { 
      processSegment(oldVertex, newVertex); 
      boundaryPathIterator.next(); 
      System.arraycopy(newVertex, 0, oldVertex, 0, 2); 
      segmentType = boundaryPathIterator.currentSegment(newVertex); 
     } 
     processSegment(newVertex, firstVertex); 
    } 

    private void initializeForNewCirclePrivate(Circle circle) { 
     currentCircle = circle; 
     currentSquareRadius = currentCircle.getRadius() * currentCircle.getRadius(); 
     initializeForNewCircle(circle); 
    } 

    private void processSegment(double[] initialVertex, double[] finalVertex) { 
     double[] segmentDisplacement = displacment2D(initialVertex, finalVertex); 
     if (segmentDisplacement[0] == 0 && segmentDisplacement[1] == 0) { 
      return; 
     } 
     double segmentLength = Math.sqrt(dotProduct2D(segmentDisplacement, segmentDisplacement)); 
     double[] centerToInitialDisplacement = new double[]{initialVertex[0] - getCurrentCircle().getCenterX(), initialVertex[1] - getCurrentCircle().getCenterY()}; 
     final double leftX = dotProduct2D(centerToInitialDisplacement, segmentDisplacement)/segmentLength; 
     final double rightX = leftX + segmentLength; 
     final double y = wedgeProduct2D(segmentDisplacement, centerToInitialDisplacement)/segmentLength; 
     processSegmentStandardGeometry(leftX, rightX, y); 
    } 

    private void processSegmentStandardGeometry(double leftX, double rightX, double y) { 
     if (y * y > getCurrentSquareRadius()) { 
      processNonIntersectingRegion(leftX, rightX, y); 
     } else { 
      final double intersectionX = Math.sqrt(getCurrentSquareRadius() - y * y); 
      if (leftX < -intersectionX) { 
       final double leftRegionRightEndpoint = Math.min(-intersectionX, rightX); 
       processNonIntersectingRegion(leftX, leftRegionRightEndpoint, y); 
      } 
      if (intersectionX < rightX) { 
       final double rightRegionLeftEndpoint = Math.max(intersectionX, leftX); 
       processNonIntersectingRegion(rightRegionLeftEndpoint, rightX, y); 
      } 
      final double middleRegionLeftEndpoint = Math.max(-intersectionX, leftX); 
      final double middleRegionRightEndpoint = Math.min(intersectionX, rightX); 
      final double middleRegionLength = Math.max(middleRegionRightEndpoint - middleRegionLeftEndpoint, 0); 
      processIntersectingRegion(middleRegionLength, y); 
     } 
    } 

    private void processNonIntersectingRegion(double leftX, double rightX, double y) { 
     final double initialTheta = Math.atan2(y, leftX); 
     final double finalTheta = Math.atan2(y, rightX); 
     double deltaTheta = finalTheta - initialTheta; 
     if (deltaTheta < -Math.PI) { 
      deltaTheta += 2 * Math.PI; 
     } else if (deltaTheta > Math.PI) { 
      deltaTheta -= 2 * Math.PI; 
     } 
     processNonIntersectingRegion(deltaTheta); 
    } 

    protected abstract void initialize(); 

    protected abstract void initializeForNewCircle(Circle circle); 

    protected abstract void processNonIntersectingRegion(double deltaTheta); 

    protected abstract void processIntersectingRegion(double length, double y); 

    protected abstract T getValue(); 

    protected final Circle getCurrentCircle() { 
     return currentCircle; 
    } 

    protected final double getCurrentSquareRadius() { 
     return currentSquareRadius; 
    } 

Quoi qu'il en soit, c'est ma description de l'algorithme. Je pense que c'est bien parce que c'est exact et qu'il n'y a pas vraiment beaucoup de cas à vérifier.

+5

Réponse intense!Devrait l'avoir séparément sur un blog je pense – Lodewijk

+2

Je crois que la quantité de temps et d'efforts pour mettre cette réponse méritent une appréciation. Et voici le mien. Je vous remercie! – Alp

+0

Je souhaite que je pourrais upvote ceci deux fois. Merci! – Daerst

Questions connexes