19.03.2020

<< | Home | 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 )):