Algorithmes de recherche du volume et du centre de masse d'un polyèdre

Tout le monde connaît probablement cet algorithme, mais «les autorités me cachaient». J'ai trouvé sa description verbale sur la troisième page du moteur de recherche dans l'archive des traductions automatiques du forum anglophone. Il me semble que sa description détaillée (et avec le code) est digne de habrosta.

Ainsi, par exemple, vous devez générer des foules pour un jouet et quelque part dans le processus pour éliminer ceux qui ne sont pas debout. Pour ce faire, vous devez trouver le centre de masse de la foule (et c'est presque la même chose que trouver son volume) et vous assurer qu'il se trouve quelque part au-dessus des jambes de la foule.



Un mob est un polyèdre, pour simplifier, nous considérons qu'un polyèdre n'est composé que de triangles (dans l'algorithme à l'intérieur se trouve la formule de l'aire gaussienne , vous pouvez donc l'étendre pour n'importe quel polyèdre, mais pourquoi ...). De plus, le polyèdre ne devrait pas avoir d'auto-intersections et limiter le volume fermé, comme il sied à des polyèdres décents.


(enfin, comme ça)

Un petit UPD expliquant pourquoi sur le KDPV la foule de droite n'est pas OK, mais la gauche OK:
La bonne image n'est pas OK parce que la foule va tomber en avant, parce que son centre de masse s'étend au-delà de la zone d'appui. La zone de support d'un polygone se tenant sur la surface est définie comme le polygone minimum à l'intérieur duquel tous les points de la surface sont situés. Dans le cas de gauche, la zone du support est décalée vers le centre de masse et plus (parce que les pattes de dinosaure sont plus grandes), et dans l'image de droite, la zone elle-même est plus petite et plus proche de la queue.
Le rapport de la zone de référence au centre de masse sera quelque chose comme ceci:


Je vais commencer tout de suite avec le code de recherche de volume (Python, entrée - une liste de points et une matrice de transition):

du code
def RecSetDirsTriangles(para, Connects, TR): """ ,           """ #1.   , ,     for i in range(0,len(Connects)): if i != para[0] and i != para[1] and Connects[i][para[0]] and Connects[i][para[1]]: #  ! fl = 1 for T in TR: if i in T and para[1] in T and para[0] in T: fl = 0 #    break if fl: # ! TR += [(para[1],para[0],i)] Recc((para[0], i) , Connects, TR) Recc((i, para[1]) , Connects, TR) def FindV(dots, Connects): """ .   - dots     [x, y, z], Connects -  , Connects[i][j]=1      i, j,  =0 """ #1.      TR = [] for i in range(1,len(Connects)):#        - if Connects[i][0]: for j in range(i+1, len(Connects)): if Connects[0][j] and Connects[i][j]: TR += [(0,i,j)] break RecSetDirsTriangles((0,i),Connects, TR) break print(" : ", len(TR)) #2.        V = 0 for T in TR: ''' : x1y2 x2y3 x3y1 x2y1 x3y2 x1y3''' S = 0.5 * (dots[T[0]][0]*dots[T[1]][1] + dots[T[1]][0]*dots[T[2]][1] + dots[T[2]][0]*dots[T[0]][1] - dots[T[1]][0]*dots[T[0]][1] - dots[T[2]][0]*dots[T[1]][1] - dots[T[0]][0]*dots[T[2]][1]) #S   +  -    ,    V += S*(dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2])/3 #    ... return math.fabs(V) 


L'essence de l'algorithme est de considérer les volumes de figures qui forment les faces du polyèdre qui «tombent» sur le plan xy. Pour ce faire, vous devez connaître la zone de projection du triangle et le signe avec lequel ajouter le volume de la figure (prisme tronqué). En fait, si les triangles sont commandés à l'avance, le volume et le signe sont réduits à un seul calcul.

Par conséquent, la première chose qu'une fonction récursive collecte des triangles est son entrée. Assemble de telle manière que lorsque vous regardez à l'extérieur d'un polyèdre, les directions pour contourner les triangles sont les mêmes (idéalement dans le sens antihoraire; si vous prenez les directions dans le sens horaire, le résultat sera correct, mais négatif - par conséquent, le module de volume est donné au retour).

C'est très simple à réaliser - nous prenons un triangle (points a1, a2, a3), recherchons ses voisins et listons deux sommets correspondants dans l'ordre inverse (par exemple, comme ceci: a2, a1, b1).
Il s'avère quelque chose comme ceci:



Maintenant, si nous projetons un tel triangle sur le plan xy, alors l'ordre de traversée pour la projection du triangle "supérieur" coïncidera avec celui initialement sélectionné, et l'ordre de traversée pour la projection du triangle "inférieur" changera de direction. En conséquence, cela changera le signe et l'aire de ce triangle, calculés par la formule de Gauss. Ici, le triangle «inférieur» - un concept conditionnel - signifie que le volume immédiatement en dessous n'est pas inclus dans le volume du polyèdre. Le triangle «inférieur» d'un polyèdre non convexe peut être supérieur à celui «supérieur».

Après ces étapes préliminaires, pour calculer le volume total du polyèdre, il suffit d'ajouter (en tenant compte du signe, qui est obtenu "par lui-même") tous les volumes des prismes tronqués collectés sur les faces et les projections de ces faces sur le plan xy. Et les volumes de prismes sont considérés comme le produit de l'aire (gaussienne, avec un signe) et les coordonnées z moyennes arithmétiques des sommets du triangle.

Si le polyèdre coupe le plan xy, alors lors du calcul du volume, tous les signes s'annulent et le résultat reste correct (il suffit de prendre les hauteurs de prisme sans module).


(en quelque sorte le prisme tronqué "supérieur" ressemble)

Avec la recherche du centre de masse, tout est à peu près le même. De même, nous devons trouver les centres de masse pour chaque prisme tronqué et résumer en coordonnées, en multipliant par le volume du prisme (on suppose que la masse est répartie uniformément dans tout le volume et que l'un peut être remplacé par un autre). Pour trouver le centre de masse d'un prisme tronqué, il est nécessaire de calculer les centres de masse de deux tétraèdres (+1 fonction) et d'un prisme ordinaire. L'algorithme "ne se détériore pas non plus" si le polyèdre traverse le plan xy (et ici pourrait être une reproduction de Magritte).


(ces deux tétraèdres, marqués en rouge et rouge, ainsi qu'un prisme triangulaire (sous le tétraèdre rouge) forment le prisme tronqué souhaité. Nous devons trouver les centres de masse et les volumes des trois figures. Les désignations correspondent à peu près aux désignations du code)

Code qui compte ceci et cela:

un peu plus de code
 def RecSetDirsTriangles(para, Connects, TR): #1.   , ,     for i in range(0,len(Connects)): if i != para[0] and i != para[1] and Connects[i][para[0]] and Connects[i][para[1]]: #  ! fl = 1 for T in TR: if i in T and para[1] in T and para[0] in T: fl = 0 break if fl: # ! TR += [(para[1],para[0],i)] Recc((para[0], i) , Connects, TR) Recc((i, para[1]) , Connects, TR) def TetrV(mas):#dot1, dot2, dot3, dot4): """   """ M = np.zeros((3,3),float) for i in range(1,4): for j in range(0,3): M[i-1][j] = mas[i][j] - mas[0][j] #print(M) return math.fabs(np.linalg.det(M)/6) def FindVandCM(dots, Connects): """     """ #1.      TR = [] for i in range(1,len(Connects)): #        - if Connects[i][0]: for j in range(i+1, len(Connects)): if Connects[0][j] and Connects[i][j]: TR += [(0,i,j)] break RecSetDirsTriangles((0,i),Connects, TR) break print(" : ", len(TR)) #2.   ,          V = 0 CM = [0, 0, 0] for T in TR: ''' : x1y2 x2y3 x3y1 x2y1 x3y2 x1y3''' S = 0.5 * (dots[T[0]][0]*dots[T[1]][1] + dots[T[1]][0]*dots[T[2]][1] + dots[T[2]][0]*dots[T[0]][1] - dots[T[1]][0]*dots[T[0]][1] - dots[T[2]][0]*dots[T[1]][1] - dots[T[0]][0]*dots[T[2]][1]) #S   +  -    ,    V += S*(dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2])/3 #    ... #c        c1 = ((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0])/3, (dots[T[0]][1]+ dots[T[1]][1]+ dots[T[2]][1])/3) #    hm = min([dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]]) hM = max([dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]]) indM = [dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]].index(hM) indm = [dots[T[0]][2] , dots[T[1]][2] , dots[T[2]][2]].index(hm) V3 = S * hm if indM == indm: # ! CM[0] += V3*c1[0] CM[1] += V3*c1[1] CM[2] += V3*hm/2 continue L = [0,1,2] L.remove(indM) L.remove(indm) indmidle = L[0] dots1 = [dots[T[0]], dots[T[1]], dots[T[2]], (dots[T[indM]][0], dots[T[indM]][1] , hm)] #  V1 = TetrV(dots1) if S < 0: V1 = -V1 V2 = S * ( dots[T[indmidle]][2] - hm)/3 #V3 = S * hm CM[0] += V1*((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0] + dots[T[indM]][0])/4) + V2*((dots[T[0]][0] + dots[T[1]][0] + dots[T[2]][0] + dots[T[indmidle]][0])/4) + V3*c1[0] CM[1] += V1*((dots[T[0]][1] + dots[T[1]][1] + dots[T[2]][1] + dots[T[indM]][1])/4) + V2*((dots[T[0]][1] + dots[T[1]][1] + dots[T[2]][1] + dots[T[indmidle]][1])/4) + V3*c1[1] CM[2] += V1*((dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2] + hm)/4) + V2*((dots[T[0]][2] + dots[T[1]][2] + dots[T[2]][2] + hm)/4) + V3*hm/2 CM[0] = CM[0]/V CM[1] = CM[1]/V CM[2] = CM[2]/V return (math.fabs(V), CM) 


Un morceau de l'algorithme où les directions des triangles sont considérées et utilisées pour comprendre le volume externe et interne est un mouvement très fort, il peut être utilisé beaucoup lorsque vous travaillez avec des polyèdres. Par exemple, si vous avez besoin de calculer la direction des normales "out" - il suffit de connaître la direction "antihoraire" pour une face - et le tour est joué!


(devinez le film!)

Source: https://habr.com/ru/post/fr479290/


All Articles