Algorithmen zum Suchen des Volumens und des Massenschwerpunkts eines Polyeders

Wahrscheinlich kennt jeder diesen Algorithmus, aber "die Behörden haben sich vor mir versteckt". Ich fand seine verbale Beschreibung auf der dritten Seite der Suchmaschine im Archiv der automatischen Übersetzungen des englischsprachigen Forums. Es scheint mir, dass seine detaillierte Beschreibung (und mit dem Code) Habrosta verdient.

So müssen Sie zum Beispiel Mobs für ein Spielzeug generieren und irgendwo im Prozess diejenigen ausmerzen, die nicht auf den Beinen stehen. Um dies zu tun, müssen Sie den Massenschwerpunkt des Pöbels finden (und dies ist fast dasselbe wie das Volumen) und sicherstellen, dass er sich irgendwo über den Beinen des Pöbels befindet.



Ein Mob ist ein Polyeder. Der Einfachheit halber glauben wir, dass ein Polyeder nur aus Dreiecken besteht (der Algorithmus enthält die Gaußsche Flächenformel im Inneren, sodass Sie sie für jedes Polyeder erweitern können, aber warum ...). Darüber hinaus sollte das Polyeder keine Selbstüberschneidungen aufweisen und das geschlossene Volumen begrenzen, wie dies für anständige Polyeder angemessen ist.


(na so)

Ein kleines UPD, das erklärt, warum auf dem KDPV der rechte Mob nicht OK ist, aber der linke OK:
Das richtige Bild ist nicht in Ordnung, weil der Mob nach vorne fallen wird, weil Sein Massenschwerpunkt erstreckt sich über den Stützbereich hinaus. Die Auflagefläche eines auf der Oberfläche stehenden Polygons ist definiert als das minimale Polygon, innerhalb dessen sich alle Punkte auf der Oberfläche befinden. Im linken Fall ist der Bereich der Stütze zum Massenmittelpunkt und mehr verschoben (weil die Dinosaurierpfoten größer sind), und im rechten Bild ist der Bereich selbst kleiner und näher am Schwanz.
Das Verhältnis der Referenzfläche zum Massenmittelpunkt ist ungefähr so:


Ich beginne sofort mit dem Volume-Suchcode (Python, Eingabe - eine Liste von Punkten und eine Übergangsmatrix):

etwas 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) 


Die Essenz des Algorithmus besteht darin, die Volumen der Figuren zu betrachten, die die Flächen des Polyeders bilden, die auf die xy-Ebene "fallen". Dazu müssen Sie die Projektionsfläche des Dreiecks und das Zeichen kennen, mit dem Sie das Volumen der Figur addieren (abgeschnittenes Prisma). Wenn Dreiecke im Voraus bestellt werden, werden Volumen und Vorzeichen auf eine einzige Berechnung reduziert.

Das erste, woraus eine rekursive Funktion Dreiecke sammelt, ist die Eingabe. Wird so zusammengebaut, dass die Dreiecksrichtungen beim Blick nach außen auf ein Polyeder gleich sind (idealerweise gegen den Uhrzeigersinn; wenn Sie die Richtungen im Uhrzeigersinn nehmen, ist das Ergebnis korrekt, aber negativ - daher wird der Volumenmodul für den Rücklauf angegeben).

Dies ist sehr einfach zu erreichen - nehmen Sie ein Dreieck (Punkte a1, a2, a3), suchen Sie nach seinen Nachbarn und listen Sie zwei übereinstimmende Eckpunkte in umgekehrter Reihenfolge auf (zum Beispiel so: a2, a1, b1).
Es stellt sich so etwas heraus:



Wenn wir nun ein solches Dreieck auf die xy-Ebene projizieren, stimmt die Traversalreihenfolge für die Projektion des „oberen“ Dreiecks mit der ursprünglich ausgewählten überein, und die Traversalreihenfolge für die Projektion des „unteren“ Dreiecks ändert seine Richtung. Infolgedessen ändert sich das Vorzeichen und die Fläche dieses Dreiecks, berechnet nach der Gauß-Formel. Das „untere“ Dreieck - ein bedingtes Konzept - bedeutet hier, dass das unmittelbar darunter liegende Volumen nicht im Volumen des Polyeders enthalten ist. Das "untere" Dreieck eines nicht konvexen Polyeders kann höher sein als das "obere".

Nach diesen vorbereitenden Schritten müssen Sie zur Berechnung des Gesamtvolumens des Polyeders (unter Berücksichtigung des Vorzeichens, das "von selbst" erhalten wird) nur alle Volumina der auf den Flächen und Projektionen dieser Flächen in der xy-Ebene gesammelten abgestumpften Prismen addieren. Und die Volumina der Prismen werden als Produkt der Fläche (Gaußsch, mit Vorzeichen) und des arithmetischen Mittels der Z-Koordinaten der Eckpunkte des Dreiecks betrachtet.

Wenn das Polyeder die xy-Ebene schneidet, heben sich bei der Berechnung des Volumens alle Vorzeichen auf und das Ergebnis bleibt korrekt (Sie müssen nur die Prismenhöhen ohne Modul angeben).


(irgendwie sieht das "obere" abgeschnittene Prisma so aus)

Bei der Suche nach dem Schwerpunkt ist alles in etwa gleich. In ähnlicher Weise müssen wir die Massenschwerpunkte für jedes abgestumpfte Prisma finden und koordinativ zusammenfassen und mit dem Volumen des Prismas multiplizieren (es wird angenommen, dass die Masse gleichmäßig über das Volumen verteilt ist und eines durch ein anderes ersetzt werden kann). Um den Massenschwerpunkt eines abgestumpften Prismas zu ermitteln, müssen die Massenschwerpunkte von zwei Tetraedern (+1 Funktion) und einem gewöhnlichen Prisma berechnet werden. Der Algorithmus "verschlechtert sich auch nicht", wenn das Polyeder die xy-Ebene kreuzt (und hier könnte es sich um eine Reproduktion von Magritte handeln).


(Diese beiden rot und rot markierten Tetraeder bilden zusammen mit einem dreieckigen Prisma (unterhalb des roten Tetraeders) das gewünschte abgeschnittene Prisma. Wir müssen die Massen- und Volumenschwerpunkte aller drei Figuren finden. Die Bezeichnungen entsprechen in etwa den Bezeichnungen im Code.)

Code, der dies und das zählt:

ein bisschen mehr 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) 


Ein Teil des Algorithmus, bei dem die Richtungen der Dreiecke berücksichtigt und zum Verstehen des externen und internen Volumens verwendet werden, ist ein sehr wichtiger Schritt. Er kann häufig verwendet werden, wenn Sie mit Polyedern arbeiten. Wenn Sie zum Beispiel die Richtung der Normalen "aus" berechnen müssen, ist es ausreichend, die Richtung "gegen den Uhrzeigersinn" für ein Gesicht zu kennen - und voila!


(Errate den Film!)

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


All Articles