----- ====== Dritter Tag: Aktueller Code ====== **19.03.2020** [[ws1920:tag_02_code|<<]] | [[ws1920:scutoids|Home]] | [[ws1920:tag_03|Text]] Aktuelle Programmversion Stand 19:00 Uhr import bpy from random import random, randint from math import * import numpy as np import time from bpy_extras.mesh_utils import ngon_tessellate # # # # KLASSEN class Sphere(object): # Klasse, die spherische Objekte beschreibt, die für das Kugelkriterium # der Delaunay-Tatraedisierung notwendig sind pos = np.array([0,0,0]) # Mittelpunkt der Sphere radius = 0 # Radius der Sphere epsilon = 0.001 # epsilon-Bereich für das Kugelkriterium def __init__(self, pos, radius): # eine Sphere wird immer mit einer Position und einem Radius initialisiert self.pos = pos self.radius = radius def punktInKugel(self, p): # diese Funktion stellt fest, ob sich ein Punkt innerhalb der Sphere befindet # (Funktion des Kugelkriteriums) if (abstandPunkte(self.pos, p) < self.radius - self.epsilon): return True else: return False class Tetraeder(object): # Klasse, die Tetraedische Objekte beschreibt, für die Delaunay-Tetraedisierung verts = [] # vertecies (Eckpunkte) des Tetraeders faces = [] # faces (Flächen) des Tetraeders (enthält Flächen-Listen, # die aus Ketten von Indexen der vertecie-Liste bestehen) sp = np.array([0,0,0]) # Schwerpunkt (= Mittelpunkt der Umkugel) -> Diese Punkte bilden die vertecies für das Voronoi-Diagramm!! name = "tetraeder" # ... umkugel = Sphere(np.array([0,0,0]), 0) # Spherenobjekt, auf dessen Oberfläche alle vertecies des Tetraeders liegen def __init__(self, verts, name): # Ein Tetraeder wird immer mit 4 Eckpunkten, vier Flächen und einem Namen initialisiert. # Außerdem berechnet der Code hier automatisch die resultierende Umkugel und den daraus # resultierenden Schwerpunkt. self.name = name self.verts = verts self.faces = [(0,1,2),(1,2,3),(2,3,0),(3,0,1)] # jedes Tetraeder hat immer dieselben 4 faces... if (self.schwerpunkt()[0] != None): self.sp = self.schwerpunkt() self.umkugel = Sphere(self.sp,abstandPunkte(self.sp,self.verts[0])) else: self.sp = np.array([None,None,None]) def schwerpunkt(self): # Funktion errechnet den Schwerpunkt des Tetraeders, indem es zwischen je # drei Eckpunkten Dreiecke aufspannt, die Schwerpunkte dieser errechnet, dann die Flächennormalen # an diesen Punkten ansetzt und den Schnittpunkt der darasu resultierenden Geraden berechnet. # Dieser Punkt ist dann der Schwerpunkt des Tetraeders (bzw. der Mittelpunkt der Umkugel) v1 = np.subtract(self.verts[1],self.verts[0]) v2 = np.subtract(self.verts[2],self.verts[0]) v3 = np.subtract(self.verts[3],self.verts[0]) sp = tetraSchwerpunkt(self.verts[0],v1,v2,v3) return sp # gibt den Schwerpunkt des Tetraeders zurück... def selfRender(self, ker, col, tetra): # falls innerhalb der Umkugel des Tetraeders kein Zellfremder Kern enthalten ist, wird das # Tetraeder als Meshobjekt generiert, in dem 3D View angezeigt und der Liste "tetra" hinzugefügt. # diese verfügt somit zum Schluss über eine Ansammlung von Referenzen aller zur Delaunay-Tetraedisierung # gehörenden Tetraeder. frei = True for m in range(0, len(ker)): if (self.umkugel.punktInKugel(ker[m]) == True): frei = False if (frei == True): generateObj("Tetraeder", t.verts, t.faces, (0,0,0), col) tetra.append(self) class Gerade(object): # Klasse, die Geradenformeln mithilfe von Stützvektor/Ortsvektor OV und # Richtungsvektor RV beschreibt. Das Objektattribut Kerne beschreibt, durch welche # beiden Kerne die Gerade verläuft ov = (0,0,0) rv = (0,0,0) kerne = [0,0] def __init__(self,ov,rv,kerne = [0,0]): # Jede Gerade wird mit einem Ortsvektor und Richtungsvektor initialisiert. Da nicht alle Funktionen # Geraden durch Kerne ziehen, wird das Attribut standardmäßig auf 0,0 gessetzt self.ov = ov self.rv = rv self.kerne = kerne def punktAn(self,k): # Gibt die Koordinaten eines Punktes auf der Geraden an, der mit der Koordinatenform berechnet wird # (im Detail: Stützvektor + Richtungsvektor*k) punkt = np.add(self.ov,(self.rv[0]*k,self.rv[1]*k,self.rv[2]*k)) return punkt # # # # FUNKTIONEN def dreieckSchwerpunkt(pos,a,b): # Berechnet mithilfe von drei Punkten den Schwerpunkt eines Dreieckes # unter Berücksichtigung der dreidimensionalen Lage des Dreieckes g1 = Gerade(pos,a) g2 = Gerade(pos,b) c = kreuzProdukt(a,b) d = kreuzProdukt(a,c) e = kreuzProdukt(b,c) g3 = Gerade(g1.punktAn(.5), np.divide(d, 40)) g4 = Gerade(g2.punktAn(.5), np.divide(e, 40)) schnitt = schnittPunkt(g3,g4) return schnitt def tetraSchwerpunkt(pos, v1, v2, v3): # berechnet den Schwerpunkt des Tetraeders mit Hilfe von drei Richtungsvektoren, die gemeinsam # zwei Dreiecke aufspannen (ein Vektor liegt genau an der Grenze zwischen beiden Dreiecken) sp1 = dreieckSchwerpunkt(pos, v1, v2) sp2 = dreieckSchwerpunkt(pos, v1, v3) # Es wird nach None-Type untersucht, da "falsche Dreiecke" sich einschleichen könnten # Dies sind zum Beispiel jene Dreiecke, deren Innenfläche 0 beträgt - somit wäre eine Seite gleich 0 # und es existiert kein Richtungsvektor, der sie beschreiben könnte if (sp1[0] != None and sp2[0] != None): g1 = Gerade(sp1, kreuzProdukt(v1, v2)) g2 = Gerade(sp2, kreuzProdukt(v1, v3)) schnitt = schnittPunkt(g1, g2) return schnitt else: return np.array([None,None,None]) def schnittTest(g1,g2): # testet, ob zwei Geraden sich schneiden und gibt gegebenfalls den Faktor aus, der in die Koordinatenform # der einen Geraden zur Berechnung des Schnittpunktes eingesetzt werden kann a = [[g2.rv[0],-g1.rv[0]],[g2.rv[1],-g1.rv[1]]] b = [[g1.ov[0]-g2.ov[0]],[g1.ov[1]-g2.ov[1]]] c = [[g2.rv[2],-g1.rv[2]],[g2.rv[1],-g1.rv[1]]] d = [[g1.ov[2]-g2.ov[2]],[g1.ov[1]-g2.ov[1]]] e = [[g2.rv[2],-g1.rv[2]],[g2.rv[0],-g1.rv[0]]] f = [[g1.ov[2]-g2.ov[2]],[g1.ov[0]-g2.ov[0]]] #print("Det :",np.linalg.det(a)) if (np.linalg.det(a) != 0): return np.linalg.solve(a,b)[0][0] elif (np.linalg.det(c) != 0): #print(">> ERGEBNIS") return np.linalg.solve(c,d)[0][0] elif (np.linalg.det(e) != 0): #print(">> SINBEGRE") return np.linalg.solve(e,f)[0][0] else: return None def schnittPunkt(g1,g2): # setzt das Ergebnis des Schnitttests in die Koordinatengleichung der einen Gerade ein # (mithilfe der Funktion "Punkt an") und gibt den Schnittpunkt von g1 und g2 aus z = schnittTest(g1,g2) if (z != None): return g2.punktAn(z) else: return np.array([None,None,None]) def abstandPunkte(p1, p2): # gibt den Abstand zweier Punkte im dreidimensionalen Raum aus (Satz des Pythagoras) a = sqrt(((p1[0]-p2[0])**2) + ((p1[1]-p2[1])**2) + ((p1[2]-p2[2])**2)) return a def zeichnePunkt(pos, size = 0.2, name = "Punkt"): # zeichnet einen kleinen Würfel an der übergebenen Stelle, sodass dieser für Debugging-Zwecke # später leicht gefunden und nachvollzogen werden kann (wird in dieser Programmversion nicht mehr benötigt) # Größe kann angepasst werden size = size/2 verts = [(size,size,size),(-size,size,size),(-size,-size,size),(size,-size,size) ,(size,size,-size),(-size,size,-size),(-size,-size,-size),(size,-size,-size)] faces = [(0,1,2,3),(4,5,6,7),(0,1,5,4),(1,2,6,5),(2,3,7,6),(3,0,4,7)] generateObj(name,verts,faces,pos) def kreuzProdukt(a,b): # Bildet das Kreuzprodukt der beiden übergebenen Parameter c = (a[1]*b[2]-a[2]*b[1] , a[2]*b[0]-a[0]*b[2] , a[0]*b[1]-a[1]*b[0]) return c def genKerne(anzahl, x_max, y_max, z_max): # Generiert in der übergebenen Menge zufällige Punkte und gibt diese in einer Liste zurück. # Die drei anderen Parameter sind für die Begrenzung des Raumes, in dem sich die Punkte befinden sollen ker = [] for n in range(0, anzahl): x = randint(-x_max, x_max)*((random()/(1+random()))) y = randint(-y_max, y_max)*((random()/(1+random()))) z = randint(-z_max, z_max)*((random()/(1+random()))) if (len(ker) != 0): neu = True # verhindert, dass kein Punkt doppelt generiert wird... hätte man auch mit if element in liste # machen können, verbesserungswürdig for m in range(0, len(ker)): if (x == ker[m][0] and y == ker[m][1] and z == ker[m][2]): neu = False if (neu == True): ker.append([x,y,z]) zeichnePunkt(ker[n], 1) else: n = n-1 else: ker.append([x,y,z]) zeichnePunkt(ker[n], 1) return ker # AB HIER ERSTELLUNG DER ZELLWÄNDE def radianBetweenVectors(v1, v2, n): # berechnet den Zwischenwinkel zweier Vektoren unter Beachtung eines Normalenvektors # (später wird für n die Flächennormale der Zellwand eingesetzt, welche die Verbindung der # beiden Knotenpunkte ist, zwischen denen sich die Zellwand befindet). a = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] b = v1[0]*v2[1]*n[2] + v2[0]*n[1]*v1[2] + n[0]*v1[1]*v2[2] - v1[2]*v2[1]*n[0] - v2[2]*n[1]*v1[0] - n[2]*v1[1]*v2[0] rad = np.arctan2(b, a) if (rad < 0): rad = (2*np.pi) +rad deg = rad*180/np.pi #print(">> v02", rad, deg, "\n") return deg def polySortList(v_list, normal): # Funktion bekommt eine Liste von Punkten (mehr oder weniger in einer Ebene liegend) # übergeben und sortiert diese, sodass sie im Kreis herum nacheinander ind er Liste # angeordnet werden. Für die korrekte Drehrichtung braucht diese Funktion den Normalenvektor # der Ebene, um den sie dreht. # Könnte man wahrscheinlich deutlich kürzer schreiben, überarbeitungswürdig! # Fall die Vertexliste aus 3 Elementen besteht, handelt es sich um ein Dreieck und die Reihnefolge ist egal. if (len(v_list) == 3): return v_list geordnet = [] # Zielliste, die zum Schluss die geordneten vertecie-Koordinaten als numpy-arrays enthält start = v_list[0] # Funktion merkt sich, bei welchem Punkt sie anfängt -> weiß dann, wann sie einmal im Kreis herum iteriert ist aktu = start # aktu beinhaltet die Vertex-Daten des aktuellen Eckpunktes vv = np.subtract(v_list[0], v_list[1]) # Vergleichswinkel (zeigt immer aus dem Kreis raus), somit kann immer eine einzige Reihenfolge generiert werden geordnet.append(start) while True: # loop start values min_angle = 360 min_stelle = 0 # suche nach nächstem Punkt for v in range(0, len(v_list)): if (aktu[0] != v_list[v][0] or aktu[1] != v_list[v][1] or aktu[2] != v_list[v][2]): dir = np.subtract(v_list[v], aktu) deg = radianBetweenVectors(vv, dir, normal) if (deg < min_angle): min_angle = deg min_stelle = v geordnet.append(v_list[min_stelle]) vv = np.subtract(v_list[min_stelle], aktu) aktu = v_list[min_stelle] if (aktu[0] == start[0] and aktu[1] == start[1] and aktu[2] == start[2]): break return geordnet def genPoly(liste, normal): # generiert aus einer Liste von durcheinander angeordneten Polygon-Koordinaten ein Mesh-Flächenobjekt # VORSICHT: diese Funktion wird nur bei konvexen Polygonen zu sinnstiftenden Ergebnissen führen! geordnet = polySortList(liste, normal) f = [] # face-Datenliste, die angibt, welche Eckpunkte in welcher Reihenfolge verwendet werden sollen for i in range(0, len(geordnet)): f.append(i) # Eigenart der generateObj-Funktion, benötigt geschachtelte Liste [[faces]], um richtig zu funktionieren facelist = [] facelist.append(f) generateObj("Zellwand", geordnet, facelist) def generateObj(name, verts, faces, loc = np.array([0,0,0]), col = ( 1.0, 1.0, 1.0, 1.0 )): # generiert ein 3D-MeshObjekt mit den Übergebenen Parametern (Name, Eckpunkte, Seiten) # und gibt dieses im 3D Viewport aus mymesh = bpy.data.meshes.new(name) myobject = bpy.data.objects.new(name,mymesh) mymat = bpy.data.materials.new("Mesh_mat") mymat.diffuse_color = col mymesh.materials.append(mymat) myobject.location = loc bpy.context.collection.objects.link(myobject) mymesh.from_pydata(verts,[],faces) mymesh.update(calc_edges=True) # ZUORDNUNG DER WÄNDE ZU KOLLEKTIONEN (funkioniert lieder noch nicht komplett...) def createCollection(name): col = bpy.data.collections.new(name) bpy.context.scene.collection.children.link(col) return col def linkToCollection(obj, col): col.objects.link(obj) # LISTEN DURCHSUCHEN def inListeAnStelle(element, liste): # durchsucht eine Liste nach einem Vektor und gibt den entsprechenden Index aus, falls es enthalten ist if element in liste: for i in range(0, len(liste)): if (element == liste[i]): return i return None def inListeAnStelleNP(element, liste): # durchsucht eine Liste von np.arrays nach einem np.array und gibt die Stelle aus, # falls dieser enthalten ist for i in range(0, len(liste)): if (np.array_equal(element, liste[i])): return i return None # SCUTOID TEST!! def watertightTest(faces): # Testet, ob das Objekt mit der übergebenen Face-Liste watertight ist, also keinerlei Löcher aufweist. # wenn jede Edge mindestens zu genau zwei Flächen zugeordnet ist, ist eine Mesh watertight (in diesem Fall) for face in faces: for a in range(0, len(face)): if (a == len(face)-1): b = 0 else: b = a+1 enthalten = False for face02 in faces: if (face != face02): if ((face[a] in face02) and (face[b] in face02)): enthalten = True if (enthalten == False): print(">> ",a,b) print(face[a], face[b]) return False return True def grenzenFlächenAneinander(f1, f2): # Funktion bekommt zwei Flächen übergeben und ermittelt, ob diese direkt aneinander grenzen. for a in range(0, len(f1)): for b in range(a+1, len(f1)): print(">> ", a, b) print(" ", f1[a], f1[b]) if ((f1[a] in f2) and (f1[b] in f2)): return True def ScutoidTest(faces): # Diese Funktion testet, ob ein Meshobjekt ein Scutoid ist oder nicht. Hierzu bekommt es eine Liste aller # faces übermittelt und testet, ob es eine Kombination aus zwei Flächen gibt (diese beiden dürfen sich nicht berühren, # sind die Deckflächen des "Prismas"), an die jeweils beide jede andere Fläche anschließt # (bis auf eine Fläche, diese ist die Dreiecksfläche) # das kleinstmöglichste Scutoid hat eine Dreiecksgrundgläche und eine vierecksgrundfläche und hat somit 6 Seiten if (len(faces) < 6): return False if (watertightTest(faces) == True): for a in range(0, len(faces)): for b in range(a+1, len(faces)): scutellum = 0 #index = 0 if not (grenzenFlächenAneinander(faces[a], faces[b])): for g in faces: if not (grenzenFlächenAneinander(g, faces[a]) == True and grenzenFlächenAneinander(g, faces[b]) == True): #index = g scutellum = scutellum + 1 if (scutellum <= 1): #and len(faces[g]) == 3): return True return False # # # # CODE FÜR DIE ZELLGENERIERUNG col = createCollection('Delaunay') col = createCollection('Voronoi') col = createCollection('Kerne') x_max = 100 y_max = 100 z_max = 100 anzahlKerne = 25 #del_col = (1.0, 0.0, 0.0, 1) tetra = [] # Kernliste wird generiert ker = genKerne(anzahlKerne, x_max, y_max, z_max) # Generiert auf Basis der Kernliste die Delaunay-Tetraedisierung for i in range(0, len(ker)): for j in range(i+1, len(ker)): for k in range(j+1, len(ker)): for l in range(k+1, len(ker)): t = Tetraeder([ker[i], ker[j], ker[k], ker[l]],"Tetraeder") # random red color, weil es nice aussieht... green_tint = randint(0, 6)/100 blue_tint = randint(0,6)/100 del_col = (1.0, green_tint, blue_tint, 1.0) t.selfRender(ker, del_col, tetra) # erstellt die notwendigen Sammellisten, die die entsprechenden Information über ALLE finalen Zellen enthalten vertex = [[] for dummy in range(anzahlKerne)] # liste alle vertecies, index markiert die Zellnummer, zu der sie gehören faces = [[] for dummy in range(anzahlKerne)] # liste alle faces, index markiert die Zellnummer, zu der sie gehören # Generiert auf Basis der Delaunay-Tetraedisierung das Voronoi-Diagramm for u in range(0, len(ker)): for v in range(u+1, len(ker)): z_verts = [] # Alle Vertecies, die rund um die Verbindung beider Zellkerne liegen (Eckpunkte der Zelland zwischen ihnen) werden in z_verts abgespeichert for w in range(0, len(tetra)): if ((ker[u] in tetra[w].verts) and (ker[v] in tetra[w].verts)): z_verts.append(tetra[w].sp) if (len(z_verts) > 2): # checkt, ob eine "echte" Fläche aufgespannt wird normal = np.subtract(ker[u], ker[v]) geordnet = polySortList(z_verts, normal) temp_face1 = [] temp_face2 = [] for i in range(0, len(geordnet)): if (inListeAnStelleNP(geordnet[i], vertex[u]) != None): temp_face1.append(inListeAnStelleNP(geordnet[i], vertex[u])) else: vertex[u].append(geordnet[i]) temp_face1.append(inListeAnStelleNP(geordnet[i], vertex[u])) if (inListeAnStelleNP(geordnet[i], vertex[v]) != None): temp_face2.append(inListeAnStelleNP(geordnet[i], vertex[v])) else: vertex[v].append(geordnet[i]) temp_face2.append(inListeAnStelleNP(geordnet[i], vertex[v])) faces[u].append(temp_face1) faces[v].append(temp_face2) for i in range(0, anzahlKerne): if (ScutoidTest(faces[i]) == True): red_tint = randint(0, 20)/100 blue_tint = randint(0,25)/100 generateObj("Scutoid", vertex[i], faces[i], np.array([0,0,0]), ( red_tint, 1.0, blue_tint, 1.0 )) else: generateObj("Zelle", vertex[i], faces[i]) #def generateObj(name, verts, faces, loc = np.array([0,0,0]), col = ( 1.0, 1.0, 1.0, 1.0 )):