All Classes Functions Variables Enumerations Pages
Berechnung der VESTA-Zyklen

Einführung FaceCenterPoints Zyklenberechnung VESTA-Zyklen Triangulierung Abstandsgesetz Ergebnisse Ausblick Anwendung Quellennachweise


Nachdem in der Methode fillup sämtliche Face Center Points ermittelt und ihre Lage zueinander durch die Vektoren xfIndex, yfIndex oder zfIndex bestimmt wurde, gelangen wir zu dem eigentlichen Kern des VESTA-Algorithmus: Benachbarte Face Center Points müssen jeweils so miteinander verbunden werden, dass sie ein Polygon bilden, das ein Oberflächenelement der gesuchten Isofläche darstellt. Dazu iteriert man über alle Face Center Points und bildet an jedem Face Center Point 4 Zyklen, indem man bei jedem Face Center Point in 4 Richtungen nach dem nächsten angrenzenden Face Center Point sucht. Das Ganze lässt sich am besten anhand eines Beispiels illustrieren:

Abb8.jpg
Abb. 8: Bildung eines Zyklus der Länge 5. Bildquelle von (b): Schlei, S. 3.

In Abbildung 8 wird gezeigt, wie ein solcher Zyklus sukzessive entsteht. Man sieht drei Würfel, die drei Zellen repräsentieren, in denen sich jeweils ein aktiver Voxel befindet (Abb. 8 (a)). Alle übrigen Voxel in der näheren Nachbarschaft sind inaktiv und wurden daher nicht dargestellt. Jede Würfelseite stellt somit eine Grenzfläche zwischen einem aktiven und einem inaktiven Voxel dar. Im Zentrum jeder dieser Grenzflächen liegt ein Face Center Point (Abb. 8 (b); Quelle: Schlei, S. 3 [1] ) - daher rührt auch der Name "Face Center Point" (Schlei nennt ihn in seinem Paper " Voxel Face Center ", vgl. Schlei, S. 3 [1] ). Nun beginnen wir mit dem FCP oben links und starten mit der Zyklenbildung, indem wir zunächst nach unten wandern (Abb. 8 (d), die Wahl dieser Richtung ist zunächst willkürlich). Dort treffen wir auf auf den nächsten FCP. Um die Richtung zu ermitteln, in die wir als nächste zu gehen haben, halten wir uns an das Schema, das unter (c) abgebildet ist: Kommt man von oben, so biegt man nach rechts ab, kommt man von rechts, so biegt man nach unten ab, kommt man von unten, biegt man nach links ab, und kommt man von links, dann biegt man nach oben ab (da die Grenzflächen auch horizontal verlaufen können, wurden die 4 Richtungen mit den Labeln n (north), e (east), s (south) und w (west) versehen). Das bedeutet nun, dass wir am zweiten FCP des Zyklus nach rechts abbiegen müssen. So gelangen wir dann zum dritten, vierten und fünften FCP, bevor wir im letzten Schritt wieder bei unserem Ausgangspunkt ankommen (Abb. 8 (e)). Damit haben wir einen Zyklus der Länge 5 konstruiert. Die 5 Face Center Points werden nun als Eckpunkte eines Polygons aufgefasst, und dieses Polygon stellt nun ein Oberflächenelement der Isofläche dar (Abb. 8 (f)). Um alle an einen Face Center Point angrenzenden Oberflächenelemente zu konstruieren und damit die vollständig geschlossene Isofläche zu generieren, muss man bei einem Face Center Point jeweils in 4 Richtungen starten, wie Abbildung 9 zeigt:

Abb9.png
Abb. 9 Vier Zyklen an einem Face Center Point

Allerdings sollte man bei einem Face Center Point nicht in eine Richtung zu starten, die bereits in einem anderen Zyklendurchlauf eingeschlagen worden ist, ansonsten würde man ja diesen Zyklus mehr als einmal durchlaufen und ein- und dasselbe Oberflächenelement mehrmals erstellen. Das wirkt sich nicht nur sehr negativ auf die Performance aus (ein Zyklus der Länge n würde genau n-mal berechnet), sondern kann auch zu fehlerhaften Oberflächen-Meshes führen: So wird in dieser Implementation zur graphischen Darstellung der Oberfläche die Bibliothek OpenMesh benutzt, und dort dürfen die gleichen Oberflächenelemente an denselben Vertices und Kanten nicht mehrfach eingefügt werden.

Disconnected Mode

Bevor wir uns damit befassen, wie man aus n vorgegebenen Punkten ein Oberflächenelement mittels Triangulierung erzeugt, wird im folgenden en detail erklärt, wie der Algorithmus von einem gegebenen Face Center Punkt zum nächsten Face Center Punkt eines Zyklus gelangt.
Wenn wir uns an einem Face Center Point befinden, kennen wir seine Position und seine Orientierung. Außerdem wissen wir, in welche Richtung wir uns bewegen müssen (North, East, South oder West). Nun gibt es stets drei Möglichkeiten, wo sich der nächste FCP befindet: entweder geradeaus oder in einer der beiden Richtungen, die vertikal zur Grenzflächen-Kante verlaufen, die wir passieren müssen, wenn wir uns in eine Richtung bewegen.

Abb10.png
Abb. 10: Die drei Abbiege-Möglichkeiten. Disconnected mode: (b),(c),(d),(e). Connected mode: (b),(c),(d),(f).

In Abbildung 10 ist ein Fall geschildert: Bewegen wir uns vom Face Center Point f nach rechts, dann kommt als nächster Face Center Point entweder ein Face Center Point an der Position x, an der Position y oder an der Position z in Frage. Ob an diesen Positionen ein FCP existiert - und wenn ja, welcher es ist -, darüber geben uns die Vektoren xfIndex oder zfIndex Auskunft (yfIndex spielt in diesem Beispiel keine Rolle). Nur wie entscheiden wir, welche dieser drei Positionen die richtige ist? Die Antwort hängt davon ab, ob Voxel A oder Voxel B (oder beide) aktiv oder inaktiv sind: Sind beide Voxel inaktiv, dann verbinden wir f mit x - der nächste FCP im Zyklus muss sich an Position x befinden (Abb. 10 (b)) -, sind beide Voxel aktiv (Abb. 10 (d)), dann befindet sich der nächste FCP an Position z. Ist nur Voxel A aktiv, Voxel B aber inaktiv, dann verbinden wir f mit y (Abb. 10 (c)). Ist umgekehrt Voxel B aktiv und Voxel A inaktiv, dann verbinden wir f mit x (Abb. 10 (e)). Eine weitere Möglichkeit (10 (f)) soll hier zunächst unbeachtet bleiben; sie wird weiter unten behandelt. Die verschiedenen Fälle lassen sich mit folgendem Pseudocode abdecken:

if (A ist inaktiv)
verbinde f mit x; neue Richtung <- north
else {
if (B ist inaktiv)
verbinde f mit y; neue Richtung <- north
else
verbinde f mit z; neue Richtung <- north
}

Wenn A inaktiv ist, dann bedeutet das aber nichts anderes, als dass der FCP an der Stelle x die Orientierung 2 hat, und wenn A aktiv und B inaktiv ist, dann bedeutet das nichts anderes, als dass der FCP an der Stelle y die Orientierung 4 hat (Abb. 10 (a)). Daher lässt sich der Code so umschreiben:

if (x.ori == 2)
verbinde f mit x; nextDirection = north
else {
if (y.ori == 4)
verbinde f mit y; nextDirection = north
else {
verbinde f mit z; nextDirection = north
}
}

Nun gilt dieser Pseudocode nur für den Fall, dass der Face Center Point f die Orientierung 4 aufweist (orientation = 4) und man nach rechts in Richtung East (direction = east) startet. Unter diesen Bedingungen wird der in Abbildung 10 demonstrierte Fall durch folgenden Code wiedergegeben, wobei die Variable position die aktuelle Position von f in dem zfIndex-Vektor beschreibt und die Face Center Points x,y und z über die Vektoren xfIndex und zfIndex mithilfe der Offsets in dem allFCP-Vektor gefunden werden:

if (orientation == 4) {
if (direction == east) {
if (allFCP.at(xfIndex[position]).ori == 2) {
nextFCP = xfIndex[position];
nextDirection = north;
}
else {
if (allFCP.at(zfIndex[position+1]).ori == 4) {
nextFCP = zfIndex[position+1];
nextDirection = north;
}
else {
nextFCP = xfIndex[position+edge2];
nextDirection = north;
}
}
}
}

Da f insgesamt 6 verschiedene Orientierungen aufweisen und man von einem Face Center Point in 4 verschiedene Richtungen gehen kann, müsste es von diesem Code-Abschnitt 6*4=24 Versionen geben - was aber kein guter Progranmmierstil wäre.
Führt man nun aber folgende Variablen ein: cube[7][4], indexStd[7][2], indexAlt[7], stdDirection[7][4], altDirection[7][4] und

int offset[7][4][3] = {{{0,0,0},{0,0,0},{0,0,0},{0,0,0}}, // ori == 0 exisitert nicht
{{-edge2,-edge2,-edge2+edge1},{0,+1,+edge1},{0,+edge2,+edge1},{-1,-1,+edge1-1}}, // ori == 1
{{0,+edge1,+1},{-edge2,-edge2,+1-edge2},{-edge1,-edge1,+1-edge1},{0,+edge2,+1}}, // ori == 2
{{+edge2,+edge1,0},{-1+edge2,-1,-1},{-edge1+edge2,-edge1,-edge1},{+edge2,+1,0}}, // ori == 3
{{0,+edge1,+edge2},{0,+1,+edge2},{-edge1,-edge1,+edge2-edge1},{-1,-1,-1+edge2}}, // ori == 4
{{+1,+edge1,0},{+1,+edge2,0},{+1-edge1,-edge1,-edge1},{+1-edge2,-edge2,-edge2}}, // ori == 5
{{+edge1,+edge2,0},{+edge1,+1,0},{+edge1-edge2,-edge2,-edge2},{-1+edge1,-1,-1}}}; // ori == 6

dann lassen sich unter der Verwendung des Konstruktors von start alle Fälle mit folgendem Code abdecken:

if (allFCP.at(stdIndex[orientation][direction]->at(position+offset[orientation][direction][0])).ori == cube[orientation][direction] {
nextDirection = stdDir[orientation][direction];
nextFCP = stdIndex[orientation][direction]->at(position+offset[orientation][direction][0]);
return nextFCP;
}
if (allFCP.at(altIndex[orientation]->at(position+offset[orientation][direction][1])).ori == orientation) {
nextDirection = start((direction+3)%4);
nextFCP = altIndex[orientation]->at(position+offset[orientation][direction][1]);
return nextFCP;
}
nextDirection = altDir[orientation][direction];
nextFCP = stdIndex[orientation][dir]->at(position+offset[orientation][dir][2]);
return nextFCP;

Dieser Code ist der Kern des VESTA-Algorithmus in dieser Implementation und führt im Endergebnis zu einer geschlossenen Oberfläche.

Connected Mode

Wie bereits oben erwähnt, gibt es aber noch eine zweite Möglichkeit, die Berechnung des nächsten Face Center Points durchzuführen: Denn in dem Fall, wo Voxel B aktiv und Voxel A inaktiv, gibt es auch die Möglichkeit, den Face Center Point f mit z statt mit x zu verbinden (Abb. 10 (f)). Wenn wir uns für diese Variante entscheiden, dann erhalten wir folgenden Pseudocode:

if (B ist aktiv)
verbinde f mit z; neue Richtung <- north
else {
if (A ist aktiv)
verbinde f mit y; neue Richtung <- north
else
verbinde f mit x; neue Richtung <- north
}

Diese Abfrage lässt sich ganz analog zu oben umsetzen. Am Ende kommen wir wieder zu einer geschlossenen Oberfläche, deren Gestalt an bestimmten Stellen von der durch die erste Variante generierten abweicht, und zwar an den Stellen, wo sich zwei aktive Voxel nur diagonal an einer Kante berühren: Im zweiten Fall werden diese Voxel direkt mit Oberflächenelementen verbunden, im ersten Fall bleiben sie getrennt. Abbildung 11 zeigt diesen Unterschied:

Abb11.png
Abb. 11: Connected Mode (links) und Disconnected Mode (rechts)

Beide Abfrage-Varianten werden in dieser Implementation von der Methode nextFCP realisiert, der man mittel des Parameters connec mitteilt, mit welcher Abfrage-Variante die Oberfläche gezeichnet werden soll. Wichtig hierbei ist, dass diese beiden Varianten nicht willkürlich miteinander gemischt werden dürfen, da sonst Löcher in der Oberfläche entstehen können.
Wird eine Oberfläche mit der ersten Variante erzeugt, so sagt man, dass der Algorithmus im disconnected mode operiert, bei der zweiten Variante befindet man sich dagegen im connected mode. Man kann beide Varianten aber trotzdem kombinieren, allerdings darf dabei kein Fehler gemacht werden.

Mixed Mode

Der Unterschied zwischen den beiden Varianten der Zyklenbildung kommt nur dann zum Tragen, wenn zwei aktive Voxel diagonal benachbart und die an diese beiden Voxel angrenzenden Voxel inaktiv sind (Abbildung 12, Bildquelle: Schlei, S. 6 [1] ). Um zu verhindern, dass Fehler entstehen, wenn beide Varianten bei der Generierung der Isofläche verwendet werden, muss man sich an eine einfache Regel halten - nur dann können innerhalb eines Durchlaufs beide Algorithmen zur Bestimmung der jeweils nächsten Face Center Points benutzt werden. Der Schlüssel dazu liegt in dem Punkt, der genau in der Mitte der Kante liegt, an der sich die beiden aktiven Voxel berühren. Dieser Punkt wird von Schlei Point Of Ambiguity genannt und hat genau 4 Face Center Points als direkte Nachbarn.

Abb12.png
Abb. 12: Zwei aktive Voxel mit den dazugehörigen Zyklen, im Disconneced Mode (links), im Connected Mode (Mitte) und mit rotem Point Of Ambiguity (rechts). Bildquelle: Schlei, S. 6.

Die Regel, die es zu beachten gilt, lautet nun so:
Wenn man beim Berechnen eines Zyklus von einem dieser 4 Face Center Points in Richtung des Point Of Ambiguity läuft (Abb. 12 (c)), dann muss man sich für genau den Konnektivitäts-Modus entscheiden, den man auch bei den anderen drei Face Center Points wählt, wenn sich man bei der Zyklenbildung in Richtung Point of Ambiguity bewegt. Mit anderen Worten: Die 4 Zyklen, die durch diesen Point Of Ambiguity laufen, müssen an dieser Stelle alle mit derselben Abfrage-Variante berechnet werden: Entweder alle 4 Verbindungen zwischen diesen 4 Face Center Points werden im disconnected mode oder im connected mode berechnet - und das ist auch intuitiv klar, denn die beiden Voxel können nur entweder verbunden oder nicht verbunden sein.

Die einheitlich Behandlung bei der Berechnung dieser 4 Zyklen-Abschnitte lässt sich dadurch gewährleisten, dass man bei jedem Face Center Point noch die Information hinterlegt, ob in der unmittelbaren Nachbarschaft ein Point of Ambiguity liegt und - wenn ja - welcher Konnektivitätsmodus beim Starten in eine bestimmte Richtung (north, east, south, west) benutzt werden soll. Dass die jeweils vier einem Point of Ambiguity benachbarten Face Center Points mit demselben Konnektivitätmodus ausgestattet werden, wird in dieser Implementation folgendermaßen bewerkstelligt: Es wird erst einmal über alle Face Center Points iteriert. Dabei wird geprüft, ob zwei in x-,y- oder z-Richtung benachbarte Face Center Points entgegengesetzte Orientierung aufweisen, denn das bedeutet, dass zwei aktive Voxel diagonal benachbart sind und an dieselben zwei inaktiven Voxel grenzen. Falls dieser Fall vorliegt, werden die 4 dem Point Of Ambiguity benachbarten Face Center Points ermittelt und deren Konnektivitätsmodus einheitlich festgelegt. Der Konnektivitätsmodus dieser Face Center Points wird jeweils in einem struct fcpPOA abgespeichert, das in dem std::vector allFCPpoa abgelegt wird (die Information über den Konnektivitätsmodus des Face Center Points allFCP[n] wird hierbei naheliegenderweise in allFCPpoa[n] gespeichert).

Die Frage, die sich nun stellt, ist: Welchen Konnektivitätsmodus sollte man wählen? Schlei schlägt in seinem Paper folgendes sehr naheliegendes Entscheidungskriterium vor: Man nehme den Mittelwert der (Dichte-)Werte der 4 benachbarten Voxel und prüfe, ob dieser Wert unter dem Isowert oder über ihm liegt. Falls er unter ihm liegt, so entscheide man sich für den disconnected mode, falls das Gegenteil der Fall ist, entscheide man sich für den connected mode. Die Begründung dafür ist sehr einleuchtend: Offenbar ist dieser Mittelwert eine sehr vernünftige Schätzung für den (Dichte-)Wert an dem Point Of Ambiguity, denn dieser liegt ja genau in der Mitte der vier Voxel; und wenn der (Dichte-)Wert über dem Isowert liegt, dann befindet sich der Point Of Ambiguity offenbar innerhalb der Isofläche im Inneren des Objekts, was bedeutet, dass die beiden aktiven Voxel verbunden werden sollten. Aus diesem Grund wurde dem Vorschlag von Bernd Schlei in dieser Implementation gefolgt.

weiter


Einführung FaceCenterPoints Zyklenberechnung Vesta-Zyklen Triangulierung Abstandsgesetz Ergebnisse Ausblick Anwendung Quellennachweise