Universität Leipzig Fakultät für Mathematik und Informatik (Institut für Informatik)

Visualisierung von Nervenfaserflächen im menschlichen Gehirn aus DTI Daten.

Diplomarbeit

Leipzig, Januar 2009

vorgelegt von Schurade, Ralph Studiengang Informatik

Betreuer : Dr. Alfred Anwander (MPI für Kognitions- und Neurowissenschaften, Leipzig) Dr. Mario Hlawitschka (University of California, Davis) Prof. Dr. Gerik Scheuermann (Universität Leipzig) 1

2

Inhalt 1. Einleitung ........................................................................................................................5 2. Grundlagen ......................................................................................................................8 2.1 Allgemeines ......................................................................................................................... 8 2.2 Magnetresonanzbildgebung ............................................................................................... 9 2.3 Diffusions-Tensor-Bildgebung ........................................................................................... 10 2.3.1 Der Diffusionstensor und der Diffusionsellipsoid ...................................................... 10 2.3.2 Zweidimensionale Visualisierung der DTI-Resultate .................................................. 12 2.3.3 Dreidimensionale Rekonstruktion von Nervenfaserbündeln..................................... 15 2.3.4 Der Fact-Algorithmus ................................................................................................. 16 2.3.5 Probabilistische Traktographie................................................................................... 16 3. Methoden ..................................................................................................................... 18 3.1 kd-Baum ............................................................................................................................ 18 3.2 Marching Cubes................................................................................................................. 20 3.3 NURBS-Flächen .................................................................................................................. 23 3.3.1 NURBS-Kurven ............................................................................................................ 23 3.3.2 NURBS-Flächen ........................................................................................................... 24 3.4 LIC ...................................................................................................................................... 26 3.5 Loop Subdivision ............................................................................................................... 28 Regeln für innere ungerade Punkte: ................................................................................... 28 Regeln für innere gerade reguläre Punkte: ......................................................................... 29 Regeln für innere gerade irreguläre Punkte:....................................................................... 29 Regeln für ungerade Randpunkte ....................................................................................... 30 Regeln für gerade Randpunkte ........................................................................................... 30 Exakte Berechnung von Randpunkten ................................................................................ 30 4. Implementierung ........................................................................................................... 32 4.1 Allgemeines ....................................................................................................................... 32 4.2 Anzeige von multimodalen MR-Daten .............................................................................. 32 4.2.1 Shaderbasiertes Multitexturing ................................................................................. 33 4.2.2 Interpolation .............................................................................................................. 34 4.3 Fasern ................................................................................................................................ 36

3

4.3.1 Auswahl mittels Box ................................................................................................... 37 4.3. 2 kd-Baum .................................................................................................................... 40 4.3. 3 Auswahl über ROI ...................................................................................................... 41 4.3. 4 Farbgebung................................................................................................................ 43 4.3.5 Rendering ................................................................................................................... 44 4.3.6 Rendering der Fasern als Röhren ............................................................................... 45 4.3.7 Betrachtung des Speicherverbrauches ...................................................................... 46 4.4 Isoflächen .......................................................................................................................... 48 4.4.1 Verbesserung ............................................................................................................. 49 4.5 Schnittfläche...................................................................................................................... 52 4.5.1 Clipping einer Oberfläche ........................................................................................... 52 5. Ergebnisse ..................................................................................................................... 54 5.1 Beispiel einer virtuellen Klingler-Dissektion...................................................................... 54 5.2 Anwendungen ................................................................................................................... 58 5.2.1 Vergleichende Morphometrie eines Faserbündels für verschiedene Probanden ..... 58 5.2.2 Segmentierung von Unterkomponenten eines Faserbündels ................................... 59 5.2.3 Visualisierung von Eigenschaften der weißen Masse auf Faserbündeln ................... 59 5.2.4 Evaluierung von Ultra-Hochfeld 7-Tesla diffusionsgewichteten Daten ..................... 60 5.2.5 Analyse der anatomischen Verbindungen zwischen funktionell aktivierten Hirnregionen ....................................................................................................................... 61 5.2.7 Analyse anatomischer Eigenschaften der weißen Substanz mithilfe der LICTexturierung ........................................................................................................................ 63 5.2.8 Kombination mehrerer Faserbündel in der virtuellen Klingler Dissektion................. 64 6. Zusammenfassung und Ausblick ..................................................................................... 65 Anhang A Quelltextfragmente ........................................................................................... 67 A.1 Methode zum Erzeugen eines kd-Baums auf einem STL Vector ...................................... 67 A.2 Rendering der Linien mit Vertex Buffer Objects ............................................................... 68 A.3 Shader-Code für Pseudo-Röhren ...................................................................................... 69 Anhang B Beschreibung des Nutzerinterface ...................................................................... 71 Literaturverzeichnis ........................................................................................................... 76 Abbildungsverzeichnis ....................................................................................................... 79 Abkürzungsverzeichnis ...................................................................................................... 82 Danksagung....................................................................................................................... 83

4

Kapitel 1

Einleitung Unter den verschiedenen nichtinvasiven Bildgebungsverfahren zur Analyse des menschlichen Gehirns in vivo1 hat die diffusionsgewichtete Magnetresonanzbildgebung die einzigartige Eigenschaft Richtungsinformationen, wie die Faserrichtungen in der weißen Gehirnsubstanz, aufzunehmen. Die Visualisierung der Diffusionsrichtung in jedem Voxel2 ist ein effizientes Verfahren, um die lokale Faserrichtung zu analysieren. Leider ist es auf gegenwärtigen Desktopcomputern nicht möglich, Regionen größer als ca. 20 x 20 Voxel mit Eigenvektoren oder Tensorglyphen sinnvoll darzustellen. Außerdem kann die hauptsächliche Diffusionsrichtung von der gewählten Ebene abweichen, was einen falschen Eindruck der wirklichen Faserrichtung vermittelt, oder das willkürliche Platzieren der Glyphen erzeugt visuelle Artefakte, die zu Fehlinterpretationen führen [1]. Zur Visualisierung der Diffusionsrichtung auf einem kompletten Schnitt durch das Gehirn wird üblicherweise eine Farbkodierung verwendet. Das reduziert jedoch den Eindruck von gerichteten Fasern für das menschliche Auge. Ein besserer Eindruck der gerichteten Fasern wird mit einer LIC3-Textur auf der Schnittebene erzeugt. Eine lokale Abbildung der Oberflächentangente auf die Faserstruktur erlaubt visuellen Zugriff auf die dreidimensionale Faserstruktur. In der Vergangenheit wurde ein großer Teil der Forschung auf Algorithmen zur Rekonstruktion und Darstellung der Nervenfasern verwendet. Das Interesse galt dabei meist der Selektion [2, 3, 4], dem Clustering [5,6], dem Rendering als Linien [7], Röhren [8] oder Oberflächen [9], dem Rendering mit hoher Qualität [10] oder der animierten Visualisierung [11].

1 2

3

am lebenden Objekt Der Begriff Voxel wird meist in der 3D-Computergrafik verwendet und setzt sich aus den Wörtern volumetric und pixel zusammen line integral convolution

5

Für ein besseres Verständnis dieser Faserstrukturen ist es jedoch von Vorteil, sie in die umgebende Anatomie einzubetten, die zum Beispiel aus hochauflösenden T1gewichteten MRI gewonnen wird. Die Methode von Park u.a. [12] rendert Stromlinien und eine Rekonstruktionen der grauen Materie, ist aber auf Grund der Verwendung von rechenintensiver Segmentierung der Daten zu langsam, um die gewünschte Interaktivität zu erlauben. Die allgemein akzeptierte Vorgehensweise ist die Kombination der dreidimensionalen Anzeige der Fasern mit achsenparallelen Schichten der T1-Daten oder Volumengrafik. Zur Darstellung der exakten räumlichen Position der Faserbündel wäre es jedoch besser das Gehirn entlang des Verlaufs der Faserbündel zu schneiden. Dabei wird das Gewebe entfernt, das die Faserbündel verdeckt und benachbarte Strukturen sichtbar gemacht. Ein erster Versuch mit aus T1-Daten erzeugten dreidimensionalen Strukturen wurde von Schultz u.a. [13] unternommen. Deren Methode zielt aber auf die Verformung einer planaren Schnittfläche, in die die Fasern eingebettet werden. Die hier betrachtete Methode soll eine verformbare Schnittfläche erzeugen, die parallel zu selektierten Faserbündeln verläuft und mit der hauptsächlichen Faserrichtung texturiert wird. Dabei soll das gewünschte Ergebnis Schaubildern in medizinischen Fachbüchern wie [14] (siehe Abbildung 1.1) möglichst nahe kommen.

Abbildung 1.1: Abbildung aus E.Ludwig, J. Klingler: Atlas Cerebri Humani; Hemisphaerium sinistrum ab latere externo conspectum [14]

6

Die Vorgehensweise wurde dabei von Visualisierungstechniken für Vektorfelder inspiriert. Vor allen Dingen von Laramee u.a. [15] und Hotz u.a. [16], die Texturen auf Oberflächen zur Veranschaulichung von technischen Tensordaten benutzen. Die zu entwickelnde Software soll außerdem zur Analyse von multi-modalen MRIDaten im täglichen Gebrauch zum Einsatz kommen. Dazu wurden folgende Funktionen implementiert: Anzeige von verschiedenen MRI-Daten als 3D-Texturen auf achsenparallelen Navigationsschichten, verschiedene Modi des Überblendens mehrerer Texturen Anzeige von vorberechneten Fasern als Linien oder Röhren Interaktive Auswahl von Faserbündeln Erstellung von Iso-Oberflächen von skalaren Datensätzen, Texturierung der Oberflächen MRI-Daten wie bei den Navigationsschichten Erstellung einer verformbaren Schnittfläche entlang selektierter Faserbündel Texturierung der Oberflächen mit LIC-Textur

Da hauptsächlich Faserstrukturen FiberNavigator getauft.

analysiert

7

werden,

wurde

das

Programm

Kapitel 2

Grundlagen 2.1 Allgemeines Die Grundlagen der Diffusions-Tensor-Bildgebung (engl. diffusion tensor imaging, DTI) wurden schon 1965 gelegt, als Stejskal und Tanner [17] die Diffusionskonstante von Wasser mit einer Kernspin Magnet-Resonanz-Tomographie (MRT) und Magnetfeld Gradient Systemen ermittelten. Verglichen mit der Diffusionsmessung mittels chemischer Marker, hat diese Methode einige Vorteile. Zum Einen ist sie nichtinvasiv, zum Anderen misst sie molekulare Bewegung entlang einer willkürlich gewählten Richtung. Man kann Diffusion von links nach rechts, von oben nach unten, von vorn nach hinten oder entlang eines beliebigen Winkels messen. Wenn man frei bewegendes Wasser misst, bedeutet dies noch nicht viel, da Messungen aus jeder Richtung das gleiche Resultat liefern. Dies bezeichnet man als isotrope Diffusion. Die Situation ändert sich jedoch, wenn man biologisches Gewebe, wie Muskeln oder Gehirne, betrachtet. Dieses besteht aus Fasern mit zusammenhängender Ausrichtung. In solchen Systemen tendiert das Wasser entlang der Fasern zu diffundieren, die Diffusion wird anisotrop. Das bedeutet, dass die Resultate der Diffusionsmessungen sind nicht mehr gleich sind, wenn in verschiedenen Richtungen gemessen wird. Wird die Diffusion entlang der Faserrichtung gemessen ist die berechnete Diffusionskonstante am größten, während sie ihren kleinsten Wert bei Messungen senkrecht zur Faserrichtung erhält. Das Interessante an anisotroper Diffusion ist, dass man die zugrundeliegende Anatomie des Untersuchungsgegenstandes rekonstruieren kann. Für Muskelgewebe besteht üblicherweise a priori Wissen über die Faserrichtung. Was ist aber wenn diese Richtung unbekannt ist? Theoretisch sollte es möglich sein, die Faserrichtung aus mehreren Messungen aus unterschiedlichen Richtungen zu ermitteln. Die Faser sollte dabei entlang der Richtung mit der größten berechneten Diffusionskonstante ausgerichtet sein.

8

Diese ergibt sich aus dem diffusionsgewichteten Signal und einem Referenzsignal. Ein Beispiel hierfür ist die Messung der anisotropen Diffusion im Gehirn. Es ist bekannt, dass das Gehirn axonale4 Fasern enthält, deren Strukturen sehr kompliziert und noch nicht voll verstanden sind. Es wäre also sehr hilfreich, aus der anisotropen Diffusion die Faserstrukturen ableiten zu können. Dies ist leider nicht ohne weiteres möglich. Anders als bei Muskelgewebe, das aus Fasern mit einheitlicher Orientierung besteht, haben im Hirn verschiedene Regionen unterschiedliche Faserrichtungen. Die Struktur des Gehirns ist zu kompliziert, um mit simpler Messung der Diffusionskonstanten beschrieben zu werden, es sei denn man schneidet kleinere Regionen aus. Das ist natürlich am lebenden Menschen nicht möglich. Dieses Problem wurde gelöst, in dem Diffusionsmessung mit Magnetresonanzbildgebung verbunden wurde, mittels derer die Messung der Diffusionskonstante an jedem Voxel mit einer Auflösung von 1-3 Millimeter möglich ist. Diese Technik wird als diffusions MRT bezeichnet und seit den späten 1980er Jahren angewandt. Diffusions MRT lieferte eine Möglichkeit die Faserstruktur des Gehirns abzuschätzen. Eine wesentlich Hürde musste aber noch übersprungen werden: Wie lassen sich aus einer Serie von Diffusionsmessungen nutzbare quantifizierende Parameter ableiten um die Diffusionsanisotropie und Faserrichtungen zu beschreiben. Steht unendlich Zeit für Messungen zur Verfügung, könnte man aus tausenden Richtungen messen, aus denen dann die mit der größten Diffusionskonstante identifiziert werden kann. Natürlich sind die verfügbare Zeit und auch die Anzahl der Richtungen, aus denen man messen kann, endlich. Aus diesem Grund wurden mathematische Modelle nötig, mit denen sich aus diesen limitierten Messungen die gewünschten Werte berechnen lassen. In den frühen 1990er Jahren wurden einige Versuche unternommen. Ein Modell auf Basis von Tensorberechnungen ging als das meist akzeptierte hervor. Das auf dieser Methode basierende diffusions MRT wird Diffusions-Tensor-Bildgebung oder DTI genannt.

2.2 Magnetresonanzbildgebung MR-Bilder bestehen aus einzelnen Bildelementen (Voxel) mit unterschiedlicher Intensität (Helligkeit). Die wesentlichen Parameter sind dabei räumliche Auflösung und Kontrast. Die räumliche Auflösung beträgt in modernen Geräten 1-3 mm oder sogar noch weniger, womit sehr detailierte Aufnahmen der Gehirnstruktur möglich sind. MRKontrast basiert üblicherweise auf den unterschiedlichen Relaxationszeiten der Wassermoleküle im Gewebe, T1 oder T2, die zur Unterscheidung der unterschiedlichen Gehirnregionen zum Beispiel graue oder weiße Gehirnsubstanz benutzt werden. MR ist 4

Axon: langer faserartiger Fortsatz der Nervenzelle

9

aber nicht in der Lage guten Kontrast innerhalb der weißen Substanz zu liefern. Diese erscheint als homogene Masse, egal wie hoch die räumliche Auflösung gewählt wird. Aus der MRT-Sicht erscheint die weiße Gehirnsubstanz wie homogene Flüssigkeit. Weiße Gehirnsubstanz besteht aus Axonen, die verschiedene Regionen des Gehirns verbinden. Diese Axone gruppieren sich mit anderen Axonen zu Faserbündeln. Der Durchmesser dieser Faserbündel kann je nach Zielregion mehrere Zentimeter betragen. Einige davon, wie das corpus callosum5, sind auch im konventionellen MRT sichtbar. Da die meisten dieser Bündel aber gleiche chemische Zusammensetzungen und damit gleiche T1- und T2- Relaxationszeiten haben, können sie nicht auf diese Art identifiziert werden. Es ist sehr schwer zu erkennen, wo sich interessante Abschnitte der weißen Gehirnsubstanz befinden und wie sie untereinander verbunden sind. Konventionelles MRT bietet auch keinen guten Kontrast für Strukturen innerhalb der grauen Gehirnsubstanz. Es ist unmöglich die Unterschiede entlang der kortikalen Schichten oder zwischen diesen zu erkennen. Andererseits ist es möglich, die Muster der Hirnwindungen als optische Hinweise für die Erkennung der Abgrenzungen der kortikalen Anatomie zu benutzen. Leider gibt es keine entsprechenden optischen Hinweise für Strukturen innerhalb der weißen Gehirnsubstanz.

2.3 Diffusions-Tensor-Bildgebung

2.3.1 Der Diffusionstensor und der Diffusionsellipsoid

MRT ist eine Technik, die Protonensignale von Wassermolekülen misst. Die Bilder daraus geben daher Wasserdichte und die Eigenschaften des beweglichen Wassers als Position im Raum an. Zusätzlich können mit MRT noch die lokalen chemischen und physikalischen Eigenschaften des Wassers gemessen werden. Zwei dieser Eigenschaften sind molekulare Diffusion und Fluss. MRT kann mit unterschiedlichen Methoden beide messen. Der Diffusionsprozess spiegelt die thermische Brownsche Bewegung wieder. Dieser Prozess kann mit dem einfachen Beispiel eines Tropfens Tinte, der auf ein Blatt Papier getropft wird erklärt werden. Normalerweise breitet sich dieser Tropfen kreisförmig aus 5

große, quer verlaufende Verbindung zwischen den beiden Hemisphären des Großhirns

10

und wird mit der Zeit größer. Das Ausmaß des Fleckes ist umso größer, je schneller die Diffusion ist. Da die Ausdehnung des Fleckes in jede Richtung gleich ist, spricht man hier von isotroper Diffusion. Wenn das Papier aber eine darunter liegende Faserstruktur hat, in der zum Beispiel dichte Fasern horizontal verlaufen und weniger dichte vertikal, wird der Tintenfleck eine eher ovale Form annehmen. Man spricht dann von anisotroper Diffusion. Ein ähnlicher Prozess findet im Gehirn statt, wo Wasser entlang der axonalen Fasern diffundiert. Würde man Tinte in die weiße Gehirnsubstanz injizieren, würde sich diese entlang der axonalen Fasern verteilen. Innerhalb der grauen Gehirnsubstanz würde die Verteilung eine zufälligere Form annehmen, da hier weniger oder keine gerichteten Faserstrukturen vorliegen.

Bei Diffusionsmessung mittels MRT kann der Grad der Wasserdiffusion entlang beliebig gewählter Achsen gemessen werden. Diffusionsmessungen können mit Bildgebung verbunden werden, aus der die effektive Diffusionskonstante an jedem Pixel gemessen wird, sogenannte Karten der Diffusionskonstanten.

Abbildung 2.1 Karten der Diffusionskonstante einer Schicht gemessen in 3 verschiedenen Richtungen. In diesen Bildern zeigt die Helligkeit das Ausmaß der Diffusion. Zum Beispiel hat das corpus callosum (gelber Pfeil) eine hohe Diffusionskonstante wenn die Messung horizontal erfolgt (erstes Bild), aber gering wenn die Messung entlang der vertikalen Achse (zweites Bild) oder rechtwinklig dazu erfolgt (drittes Bild).

Abbildung 2.1 zeigt in drei unterschiedlichen Richtungen gemessene Karten der Diffusionskonstanten derselben Schicht. Die Diffusion wurde entlang der Pfeilrichtungen gemessen, im dritten Bild senkrecht zur Ebene. Die Helligkeit der Pixel zeigt das Ausmaß der Diffusion an. Man sieht zum Beispiel im corpus callosum (gelber Pfeil) eine hohe Diffusionskonstante, wenn die Diffusion entlang der horizontalen Achse gemessen wird (linkes Bild), aber einen geringeren Wert entlang der anderen 11

Achsen. Wie man an diesen Bildern leicht sehen kann, ist die Diffusion im Gehirn anisotrop: das Ausmaß der Diffusion hängt besonders in der weißen Substanz von der Messungrichtung ab. Die Größe der Anisotropie kann aus mehreren Messungen entlang unterschiedlicher Achsen abgeleitet werden. Auf Grund der Komplexität der betrachteten Zytoarchitektur und nicht vermeidbarem Rauschen können die Messungen recht komplexe Muster der Anisotropie liefern. Ist jedoch eine parallele Faserorientierung innerhalb eines Voxels vorhanden, kann die Diffusion mit einem simplen 3D-Ellipsoiden beschrieben werden. In der Diffusions-Tensor-Bildgebung wird die Diffusion in einer Serie von Messungen entlang unterschiedlicher räumlicher Richtungen ermittelt. Aus diesen Messungen wird für jeden Voxel der Diffusionsellipsoid, der die durchschnittliche Form der Zytoarchitektur in diesem Voxel beschreibt, erzeugt. Während der Messungen und der Erzeugung des Ellipsoiden wird ein 3 x 3 Tensor verwendet, von welchem der Name Diffusions-Tensor-Bildgebung abgeleitet wurde. Der symmetrische 3 x 3 Diffusionstensor hat 6 unabhängige Elemente, die zur vollständigen Beschreibung Größe und Ausrichtung des Ellipsoiden nötig sind. Für die Erzeugung des Tensors sind daher mindestens sechs Messungen der Diffusionskonstanten entlang verschiedener Abbildung 2.2: Parameter zur räumlicher Richtungen nötig (Abb. 2.2). Beschreibung des Diffussionselipsoiden

Die drei Eigenwerte 𝜆1 , 𝜆2 , 𝜆3 definieren die Form und die Eigenvektoren 𝑣1 , 𝑣2 , 𝑣3 die Ausrichtung des Ellipsoiden.

2.3.2 Zweidimensionale Visualisierung der DTI-Resultate

Nachdem die sechs Tensorparameter ermittelt wurden, wird die Information üblicherweise wieder reduziert um visuell auswertbare Ergebnisse zu erhalten. Es ist fast unmöglich alle sechs Parameter in einem einzigen Bild intuitiv verständlich zu visualisieren. Die beiden gebräuchlichsten Verfahren sind dabei die Karte der Anisotropie und die Richtungskarte oder auch Farbkarte der Faserorientierung. Da der Begriff Farbkarte auch für die Abbildung eines skalaren Wertes mit einer Farbpalette benutzt wird, soll hier im Weiteren nur Richtungskarte verwendet werden und Farbkarte für die farbige Darstellung skalarer Daten wie Wahrscheinlichkeitskarten und fMRTDaten.

12

Die Karte der Anisotropie (Abbildung 2.3) zeigt die Grad der Gerichtetheit (Anisotropie) des Diffusionsellipsoiden. Bei dieser Darstellung erscheint die weiße Gehirnsubstanz heller als die graue und zeigt damit Regionen im Gehirn mit zusammenhängenden Faserstrukturen. Man nimmt an, dass Faktoren wie axonale Dichte, Myelination und Homogenität in der axonalen Ausrichtung den Grad der Diffusionsanisotropie beeinflussen. Für diese Karte wird die sogenannte Fraktionelle Anisotropie (FA), die auf [0,1] skaliert ist, berechnet.

𝐹𝐴 =

1

(𝜆 1 −𝜆 2 )2 +(𝜆 2 −𝜆 3 )2 +(𝜆 3 −𝜆 1 )2

2

𝜆 21 +𝜆 22 +𝜆 23

Abbildung 2.3: Karte der fraktionellen Anisotropie als Graustufenbild

13

(1)

Für die Richtungskarte (Abbildung 2.4) wird der Eigenvektor des Tensors mit dem größten Eigenwert benutzt, da angenommen wird, dass dieser die dominante Faserausrichtung repräsentiert. Dieser Vektor ist ein Einheitsvektor der die Bedingung 𝑣𝑥2 + 𝑣𝑦2 + 𝑣𝑧2 = 1 erfüllt.

Abbildung 2.4: Farbige Darstellung der Haupt- Eigenvektoren

Für die Einfärbung [18] wird eine 24-Bit Farbpalette benutzt. Die Werte 𝑣𝑥 x 255, 𝑣𝑦 x 255 und 𝑣𝑧 x 255 werden den Farbkanälen RGB zugewiesen. In dieser Darstellung erscheinen die Richtungen rechts-links Rot, anterior-posterior Grün und superiorinferior Blau. Auf diese Art kann jeder beliebige Winkel durch eine Mischung der Farbwerte ausgedrückt werden. Zum Beispiel werden Fasern, die in einem Winkel von 45 Grad zwischen rechts-links (rot) und aterior-posterior (grün) verlaufen, gelb (rot+grün) dargestellt.

Um die Darstellung der Richtungen in isotropen Gehirnregionen, in denen keine dominante Ausrichtung der Fasern existiert und die Messergebnisse auf Rauschen beruhen, zu unterdrücken, wird der 24-Bit Farbwert mit dem FA-Wert multipliziert.

14

2.3.3 Dreidimensionale Rekonstruktion von Nervenfaserbündeln

Als Traktografie [19] werden Verfahren bezeichnet, die den Verlauf größerer Nervenfaserbündel rekonstruieren. Üblich sind hierbei Darstellungen von HyperStromlinien, dreidimensionalen Linien, deren Verlauf der Richtung des größten Diffusions-Koeffizienten folgt. Die Tatsache, dass die Diffusions-Tensor-Bildgebung derzeit das einzige Verfahren ist, das eine nicht-invasive Darstellung der Nervenfaserbündel erlaubt, hat wesentlich zu ihrer Verbreitung beigetragen. Andererseits ist es aufgrund dessen schwierig zu überprüfen, inwiefern die Ergebnisse gängiger Traktografie-Verfahren mit dem tatsächlichen Verlauf der Nervenbahnen übereinstimmen. Erste Versuche der Validierung im Tierexperiment stützen jedoch die Vermutung, dass die Hauptrichtung der Diffusion die Ausrichtung kohärenter Nervenfasern anzeigt und weisen Übereinstimmungen zwischen nichtinvasiver Traktografie und nach dem Tod durchgeführten histologischen Untersuchungen nach. Bereiche, in denen Faserbündel sich auffächern oder kreuzen, werden von der DT-MRT jedoch nur sehr unzureichend erfasst. Zur Erzeugung der Linien sind zwei Ansätze verbreitet. Im ersten werden die Fasern aus einer vorbestimmten Region (region of interest – ROI) gestartet, im zweiten, dem sogenannten „brute force“ Ansatz wird aus jedem Voxel eine Linie gestartet. Dies ist zeitaufwendiger bietet aber eine wesentlich umfangreichere Darstellung der Fasern. Wir verwenden hier den zweiten Ansatz, da dieser die besseren Resultate liefert und mittels interaktiver Auswahlmöglichkeiten die Betrachtung der Fasern im gesamten Gehirn möglich ist. Zur Stromlinienberechnung wird der FACT6-Algorithmus [20,21] verwendet, der eine einfache, lineare Fortsetzung der Linie entsprechend des Winkels des größten Eigenvektors erzeugt. Um Fehler durch Rauschen und andere Störungen zu vermeiden, werden die Fasern nur in Voxeln mit einer FA größer einem Schwellwert gestartet und in Pixeln mit FA kleiner dieses Schwellwertes abgebrochen. Für den Schwellwert wird üblicherweise 0.2 gewählt. Der Winkel der Fortsetzung der Linie wird begrenzt, in dem das Skalarprodukt zweier aufeinanderfolgender Eigenvektoren gebildet wird und nur Werte größer 0.75 weiterbehandelt werden.

6

Fiber Assignment by Continuous Tracking

15

2.3.4 Der Fact-Algorithmus Ausgehend von einem Saatpunkt wird ein Polygon erzeugt, dessen Stützstellen genau auf den Voxelgrenzen liegen. Die Geradenstücke sind innerhalb der Voxel parallel zum entsprechenden Eigenvektor 𝑒1 des Diffusionstensors. Der Ablauf des Algorithmus ist im Folgenden zusammengefasst:

1. Wähle Saatpunkt 𝑝0 ∈ ℝ3 . Sei 𝑥0 ∈ ℕ3 der zu 𝑝0 korrespondierende Voxel und setze t := 0. 2. Bestimme Faserrichtung 𝑣𝑡+1 = 𝑒1 (𝑥𝑡 ) aus Diffusionstensor im Voxel𝑥𝑡 . 3. Berechne Skalar s, so dass 𝑠 × 𝑣𝑡+1 , 𝑣𝑡 > 0 und 𝑝𝑡+1 = 𝑝𝑡 + 𝑠𝑣𝑡+1 gerade auf der Voxelgrenze des Voxels 𝑥𝑡 liegt. 4. Bestimme 𝑥𝑡+1 das zu 𝑥𝑡 benachbarte Voxel, definiert durch die Voxelgrenze 𝑝𝑡+1 . 5. Teste Abbruchkriterium. Sind Abbruchkriterien erfüllt, gebe Polygon (𝑝0 , … , 𝑝𝑡 ) zurück. 6. Erhöhe Zähler um eins t := t+1 und gehe zu Schritt 2.

2.3.5 Probabilistische Traktographie

Im Gegensatz zu der oben beschriebenen Stromlinien-Methode werden bei den Wahrscheinlichkeitskarten [22, 23] nicht die tatsächlichen Faserverläufe rekonstruiert. Es wird für jedes Voxel eine Wahrscheinlichkeit berechnet, die beschreibt, ob eine Verbindung zu dem definierten Saatpunkt existiert. Prinzipiell ist die Frage, ob eine Verbindung zwischen einem Voxel und dem Saatpunkt existiert, mit „Ja“ oder „Nein“ zu beantworten. Aufgrund der im Verhältnis zum Durchmesser der Axone groben Auflösung der diffusionsgewichteten MRT und der indirekten und rauschbehafteten Messmethode sind die Informationen über den Verlauf der Fasern nicht eindeutig. Die Wahrscheinlichkeitskarten spezifizieren den Grad der Sicherheit, dass tatsächlich eine Verbindung zwischen Saatpunkt und Voxel existiert. Dies hat zwangsläufig zur Folge, dass die Sicherheit der Aussage mit zunehmender Distanz zum Saatpunkt abnimmt. Im 16

Allgemeinen ist diese distanzabhängige Verringerung der Sicherheit bei großen und stark ausgeprägten Nervenbahnen geringer als bei kleinen Bahnen. Die Interpretation von Wahrscheinlichkeitskarten gestaltet sich häufig aufgrund der Distanzabhängigkeit und der nicht sehr intuitiven Darstellung ohne Erfahrung als schwierig. Ein großer Vorteil ist die einfache Handhabung verglichen zu den Kurven oder Polygonen der Stromlinien-Methoden. Zum Beispiel können Gruppenstatistiken relativ einfach damit erstellt werden.

Abbildung 2.5: T1-Graustufenbild überlagert von einer Warscheinlichkeitskarte als Farbkarte für eine Startpunkt im corpus callosum

17

Kapitel 3

Methoden Zur Erzeugung einer virtuellen Klingler-Dissektion werden verschiedene, etablierte Methoden der Visualisierung kombiniert. Für die schnelle Faserselektion ist es notwendig eine Datenstruktur zu verwenden, die die Anzahl der nötigen Tests zur Entscheidung, ob eine Faser ausgewählt ist, auf ein Minimum reduziert. Der kd-Baum hat sich dafür als optimal erwiesen. In diesem Kapitel sollen die theoretischen Grundlagen dieser Methoden beschrieben werden. Diese Methoden sind: -

kd-Baum zur Beschleunigung de Faserauswahl Marching Cubes Algorithmus zur Erzeugung von Iso-Oberflächen Loop Subdivision zur Glättung und Verfeinerung von Dreiecksgittern NURBS-Darstellung für die Schnittfläche LIC-Texturierung

3.1 kd-Baum

Ein kd-Baum ist eine Datenstruktur zur binären Raumaufteilung von Punkten im kdimensionalen (hier 3-dimensionalen) Raum. Kd-Bäume eignen sich besonders gut für orthogonale Bereichsanfragen. Beim Aufbau des kd-Baums wird über die Ebenen rotiert. An der Wurzel sind die Punkte entlang der x-Achse sortiert, die Kinder entlang der y-Achse, deren Kinder entlang der z-Achse usw. Auf jeder Ebene ist die Wurzel, an der das Volumen geteilt wird, üblicherweise der Median. Dieses Verfahren erzeugt einen balancierten kd-Baum.

18

Da die Anzahl der Punkte in unserem Fall statisch ist, kann auf Funktionen zum Einfügen und Löschen von Punkten verzichtet und der Baum effizient in einem Array gespeichert werden. Die Wurzel des Baums ist das mittlere Element des Arrays, die Knoten jeweils die mittleren Elemente der Teil-Arrays. Dies führt zu einem homogenen kd-baum, bei dem jeder Knoten einen Datensatz beinhaltet. Inhomogene kd-Bäume hingegen speichern Daten nur an den Blättern.

Abbildung 3.1: Dreidimensionale Darstellung eines kd-Baumes

Abbildung 3.1 zeigt einen 3-dimensionalen kd-Baum. Die erste Teilung entlang der roten Ebene teilt das Volumen in zwei Unterbereiche. Diese werden entlang der Grünen und die daraus entstehenden entlang der blauen Ebenen geteilt. Da keine weiteren Teilungen vorgenommen werden, ist das Ergebnis ein Baum, der das Volumen in acht Bereiche teilt.

Die Erzeugung des Baumes erfolgt über einen rekursiven divide-and-conquer Algorithmus (siehe Anhang). Dieser Algorithmus ist sehr einfach zu parallelisieren, da nach anfänglichem Teilen die weitere Bearbeitung auf voneinander unabhängigen Bereichen erfolgt. Die Erzeugung des kd-Baums für ca. 29 Millionen Punkte benötigt auf einem System mit Intel Core 2 Quad 9300 Prozessor 6 Sekunden.

19

3.2 Marching Cubes

Der Algorithmus, wie von Lorensen und Cline [8] beschrieben, besteht aus zwei primären Bestandteilen. Im Ersten werden die Lage der Oberfläche ermittelt und die Dreiecke erzeugt und im Zweiten wird für jeden Dreieckspunkt die Normale berechnet. Marching Cubes benutzt einen divide-and-conquer Ansatz zur Ermittlung der Oberfläche in einem Würfel, der aus acht Punkten, von denen je 4 aus zwei benachbarten Schichten des Datensatzes kommen, erzeugt wird.

Abbildung 3.2: Veranschaulichung eines Würfels zwischen zwei Schichten

Der Algorithmus bestimmt wie der Würfel (cube) von der Oberfläche geschnitten wird und geht dann zum nächsten weiter (marching). Dafür wird jedem Punkt des Würfels eine 1 zugewiesen wenn der Wert größer oder gleich dem für die Oberfläche gewählten ist. Die Punkte befinden sich innerhalb bzw. auf der Oberfläche. Punkte mit einem Wert kleiner als dem gewählten erhalten eine 0 und befinden sich außerhalb. Die Oberfläche schneidet den Würfel, wenn mindestens ein Punkt außerhalb der Oberfläche ist und die anderen innerhalb. Mit dieser Feststellung wird die Topologie der Oberfläche innerhalb des Würfels bestimmt.

20

Abbildung 3.3 Nummerierung der Ecken und Kanten des Würfels

Da es 8 Eckpunkte und zwei Zustände, innerhalb und außerhalb, gibt, existieren nur 28 =256 verschiedene Möglichkeiten, wie die Oberfläche einen Würfel schneiden kann. Mittels Nummerierung wird eine Tabelle der geschnittenen Kanten für diese 256 Fälle erzeugt.

Diese 256 Fälle zu triangulieren ist aufwendig und fehleranfällig. Zwei verschiedene Arten der Symmetrie reduzieren das Problem von 256 auf 15 Muster. Erstens ist die Topologie unverändert, wenn die Zustände der Punkte, innerhalb oder außerhalb, invertiert werden. Dies reduziert die Anzahl der zu betrachteten Fälle auf 128. Die zweite Art der Symmetrie, Rotationssymmetrie, reduziert dies weiterhin auf 15 Basiskonfigurationen. Abbildung 3.4 zeigt die 15 möglichen Konfigurationen.

21

Abbildung 3.4: Darstellung der triangulierten Würfel

Das einfachste Muster, 0, entsteht wenn sich alle Punkte innerhalb oder außerhalb der Ebene befinden. Das nächste, 1, entsteht wenn die Oberfläche einen der Punkte von den anderen 7 abschneidet, was in einem Dreieck aus den drei Kantenschnittpunkten resultiert. Die anderen Muster erzeugen mehrere Dreiecke. Durch Permutation mittels Komplementbildung und Rotation entstehen die 256 möglichen Konfigurationen. Für jede Konfiguration wird ein Index erzeugt, der sich aus dem Status jedes einzelnen Punktes ergibt. Ausgehend von der Nummerierung der Punkte wie in Abbildung 3.3, enthält der 8-Bit Index ein Bit für jeden Eckpunkt des Würfels. Dieser Index fungiert als Pointer auf eine Kantentabelle, die für jede Konfiguration die Kantenschnittpunkte enthält. Mit diesem Index kann der Oberflächenverlauf entlang der Kanten linear interpoliert werden.

22

3.3 NURBS-Flächen Die NURBS-Darstellung für Kurven und Flächen ist auf Grund ihrer Flexibilität und Kompaktheit beliebt. Stückweise polynomiale Kurven beliebigen Grades mit beliebigen Anschlussbedingungen an den Segmenttrennstellen sind mit NURBS darstellbar. Auch nicht-polynomiale Primitive, wie zum Beispiel Kreise, Ellipsen, Kugeln, Kegel etc. können mit NURBS exakt beschrieben werden.

3.3.1 NURBS-Kurven Eine offene NURBS-Kurve der Ordnung k (vom Grad k - 1) ist gegeben durch m + 1 de BoorPunkte 𝑑𝑖 mit Gewichten 𝑤𝑖 , (i = 0, … , m) und durch einen monotonen Knotenvektor {𝑡0 ≤ 𝑡1 ≤ ⋯ ≤ 𝑡𝑚 +𝑘 }:

𝑓 𝑡 =

𝑚 𝑘 𝑖=0 𝑤 𝑖 𝑑 𝑖 𝑁𝑖 (𝑡) 𝑚 𝑘 𝑖=0 𝑤 𝑖 𝑁𝑖 (𝑡)

, 𝑡 ∈ [𝑡𝑘−1 , 𝑡𝑚 +1 ]

(2)

Dabei sind 𝑁𝑖𝑘 die Basis-Splines (B-Splines) mit lokalem Träger [𝑡𝑖 , 𝑡𝑖+𝑘 ], das heißt 𝑁𝑖𝑘 ist außerhalb dieses Intervalls Null. Eine B-Spline Kurve interpoliert die Randpunkte 𝑑0 bzw. 𝑑𝑚 , wenn die Randknoten die Vielfachheit k besitzen, das heißt 𝑡0 = 𝑡𝑘−1 𝑏𝑧𝑤. 𝑡𝑚 +1 = 𝑡𝑘+𝑚 . Zum Auswerten einer B-Spline-Kurve an der Stelle 𝑡 ∈ [𝑡𝑘−1 , 𝑡𝑚 +1 ] empfiehlt sich der de Boor-Algorithmus, der im Folgenden für rationale B-Splines beschrieben ist. Zuerst wird der Index r des Knotenintervalls in dem sich s befindet ermittelt, also 𝑡 ∈ [𝑡𝑟 , 𝑡𝑟+1 ]. An der Stelle t haben nur k de Boor-Punkte Einfluss auf die Kurve. Diese werden zunächst umbenannt:

𝑑𝑖0 ≔ 𝑑𝑖 , 𝑤𝑖0 ≔ 𝑤𝑖 , 𝑖 = 𝑟 − 𝑘 + 1, … , 𝑟 . Die Menge dieser Punkte wird nun in einem mehrstufigen Prozess jeweils um eins verkleinert, in dem auf jeder Kante des zugehörigen Polygonzuges ein neuer Punkt eingefügt wird und die alten Punkte eliminiert werden. Der letzte Punkt 𝑑𝑟𝑘−1 , der dabei übrig bleibt, ist der gesuchte Punkt auf der B-Spline Kurve. Dieser Unterteilungsprozess lautet wie folgt:

23

𝑗

𝛼𝑖 =

𝑡−𝑡 𝑖 𝑡 𝑖+𝑘−𝑗 −𝑡 𝑖

,

𝑗

𝑗

𝑗

𝑗

𝑗 −1

𝑗

𝑗 −1

𝑤𝑖 = 1 − 𝛼𝑖 𝑤𝑖−1 + 𝛼𝑖 𝑤𝑖 𝑑𝑖 = 1 − 𝛼𝑖

𝑗 −1

𝑤 𝑖−1 𝑗 𝑤𝑖

𝑗 −1

, 𝑗 −1

𝑗 𝑤𝑖 𝑗 𝑤𝑖

𝑑𝑖−1 + 𝛼𝑖

𝑗 −1

𝑑𝑖

,

(𝑖 = 𝑟 – 𝑘 + 𝑗 + 1, … , 𝑟; 𝑗 = 1, … 𝑘 − 1) 𝑓 𝑡 = 𝑑𝑟𝑘−1

(3)

Nicht nur die B-Spline Kurve, sonder auch die zugehörige Tangente kann mit dem de Boor-Algortihmus ermittelt werden: Der Vektor 𝑘−2 𝑑𝑟𝑘−2 − 𝑑𝑟−1

ist nämlich proportional zur ersten Ableitung der B-Spline-Kurve.

3.3.2 NURBS-Flächen Aus den univarianten Basisfunktionen 𝜑𝑖 einer Kurvendarstellung

𝑓 𝑡 =

𝑖 𝑐𝑖 𝜑𝑖 (𝑡)

(4)

( 𝑐𝑖 sind die Koeffizienten oder Kontrollpunkte) lässt sich unmittelbar eine Flächendarstellung gewinnen:

𝑓 𝑡, 𝑢 = =

𝑖 𝑐𝑖

𝑖

𝑗 𝑐𝑖𝑗 𝜑𝑖

𝑡 𝜑𝑗 𝑢

𝑢 𝜑𝑖 𝑡 𝑚𝑖𝑡 𝑐𝑖 𝑢 =

𝑗 𝑐𝑖𝑗 𝜑𝑗 (𝑢)

(5)

Man bildet also zunächst eine Kurvenschaar 𝑐𝑖 (𝑢) , welche für konstantes u die Koeffizienten 𝑐𝑖 einer Isoparameterlinie, das heißt einer Flächenkurve nach t, liefert. Analog kann man auch zuerst eine Kurvenschaar nach t aufstellen und dann eine einzelne Kurve nach u auswerten. Dieses Konzept der Tensorproduktflächen lässt sich auch mit NURBS konstruieren. Eine beidseitig offene NURBS-Fläche ist gegeben durch zwei Ordnungen k, l, eine (m + 1) x (n + 1)-Matrix gewichteter de Boor-Punkte 𝑑𝑖𝑗 , 𝑤𝑖𝑗 , (𝑖 = 0, … , 𝑚; 𝑗 = 0, … , 𝑛) und durch zwei Knotenvektoren {𝑡0 ≤ 𝑡1 ≤ ⋯ ≤ 𝑡𝑚 +𝑘 } und {𝑢0 ≤ 𝑢1 ≤ ⋯ ≤ 𝑢𝑛+𝑙 }. Daraus ergibt sich die NURBS-Flächendarstellung 24

𝑓 𝑡, 𝑢 =

𝑚 𝑛 𝑘 𝑖=0 𝑗 =0 𝑤 𝑖𝑗 𝑑 𝑖𝑗 𝑁𝑖 𝑚 𝑛 𝑘 𝑖=0 𝑗 =0 𝑤 𝑖𝑗 𝑁𝑖 𝑡

𝑡 𝑁𝑗𝑙 (𝑢 ) 𝑁𝑗𝑙 (𝑢)

, 𝑡, 𝑢 ∈ 𝑡𝑘−1 , 𝑡𝑚 +1 𝑥 𝑢𝑙−1 , 𝑢𝑛 +1 (6).

Eine NURBS-Fläche kann an der Stelle (t,u) durch mehrmaliges Anwenden des de Boor-Algorithmus für Kurven ausgewertet werden. Zuerst werden die Indizes r,s der Knotenintervalle, in den sich t und u befinden ermittelt:

𝑡, 𝑢 ∈ [𝑡𝑟 , 𝑡𝑟+1 )𝑥[𝑢𝑠 , 𝑢𝑠+1 ) Daraus ergibt sich die 𝑘 𝑥 𝑙 Untermatrix

𝑑𝑖𝑗 , 𝑤𝑖𝑗 (𝑖 = 𝑟 − 𝑘 + 1, … , 𝑟; 𝑗 = 𝑠 − 𝑙 + 1, … , 𝑠) der de Boor-Punkte, welche einen Einfluss auf die Fläche an der Stelle (t,u) haben. Für jede Zeile dieser Untermatrix (mit konstantem Index i) liefert der de Boor-Algorithmus unter Berücksichtigung des Knotenvektors {u} den gewichteten Punkt 𝑑𝑖 (𝑢) = 𝑑𝑟𝑘−1 , 𝑤𝑖 𝑢 = 𝑤𝑟𝑘−1 . Diese Punkte gehören nun zu einer NURBS-Kurve f(t), welche auf der Fläche liegt und unter Berücksichtigung des Knotenvektors {t} mit dem de Boor-Algorithmus ausgewertet wird. Für das Rendering der NURBS-Fläche wird diese auf einem regulären, rechteckigen Gitter abgetastet, wobei jedes Rechteck in zwei Dreiecke zerlegt wird

25

3.4 LIC LIC [25], entwickelt 1993, steht für line integral convolution, was mit Linienintegralfaltung übersetzt werden kann. Die Grundidee von LIC entstand aus der Beobachtung eines Ölfilms auf Objekten im Windkanal. Das Öl bildet dabei Schlieren entlang der Strömungsrichtungen. LIC ist heute eine allgemein akzeptierte Methode zur Visualisierung von Richtungsinformationen auf Oberflächen. In der Originalmethode wird für jeden Pixel eine lokale Stromlinie berechnet. Dann wird für jedes Segment der Stromlinie, die innerhalb des Pixels liegt, das Integral eines Faltungskerns berechnet. Diese Integrale werden benutzt, um eine gewichtete Summe aller von dieser Stromlinie getroffenen Pixel zu berechnen. Die resultierende Summe ergibt den Ausgabewert für diesen Punkt. Üblicherweise wird weißes Rauschen als Eingabewert für alle Punkte verwendet.

Abbildung 3.5 Strömungsvisualisierung um bluntfin

Das ursprünglich texturbasierte Verfahren für statische zweidimensionale Felder wird für unseren Fall der Tensorfelder im dreidimensionalen Raum erweitert. Dabei werden nicht mehr parametrisierte Texturen verwendet, sondern ein geometriebasierter Ansatz. Hierfür wird die Oberfläche mittels Loop-Subdivision [26] soweit verfeinert, dass die Auflösung der Oberfläche gleich der Displayauflösung ist. Die Linien werden effizient durch konstante Interpolation in den Dreiecken integriert. Durch die kurze Integrationslänge und die hohe Auflösung sind Fehler vernachlässigbar. Die Richtung der Eigenvektoren ist dabei unwichtig, da nur die tangentiale Information genutzt wird.

26

Zusätzlich könnte die Länge der Stromlinien skaliert werden, in dem die fraktionelle Anisotropie als Parameter genutzt wird. Zuerst wird die Tensorinformation auf die Oberfläche abgebildet und der größte Eigenvektor gebildet. Das erzeugt einen Richtungsvektor pro Dreieck. Für die Dreiecksmenge wird außerdem eine Rauschtextur erzeugt. Als nächstes werden in zufällig gewählten Dreiecken Stromlinien gestartet und für jedes Dreieck die Anzahl der Stromlinien, die es treffen, notiert. Dreiecke, die nicht von einer bestimmten Anzahl Stromlienen getroffen werden, werden als zusätzliche Saatpunkte für Stromlinien genutzt. Die Stromlinien dienen als Filterkern auf der Rauschtextur. Durch die hohe Auflösung des Dreiecksgitters erzeugt dies ein optisch ausreichendes Ergebnis. Im Vergleich zu Verfahren mit parametrisierten Texturen erfordert dieses Vorgehen eine wesentlich höhere Anzahl an Dreiecken. Es ist jedoch keine globale Parametrisierung erforderlich, was es für jegliche glatte Oberflächen beliebiger Topologie anwendbar macht. Dies ermöglicht es, zusätzlich zur Schnittfläche auch kompliziertere Iso-Oberflächen zu texturieren, wie sie auf Wahrscheinlichkeitskarten erzeugt werden. (Abb. 3.6).

Abbildung 3.6 LIC auf Iso-Oberfläche einer Wahrscheinlichkeitskarte

27

3.5 Loop Subdivision Dieses Verfahren wurde von Charles Loop [26] entwickelt und ist eines der Einfachsten. Ein weiterer großer Vorteil besteht darin, dass dieses Verfahren mit Dreiecksnetzen arbeitet welche auch Grundlage der Rendering Pipeline heutiger Grafikkarten sind. Grundsätzlich wird jedes Dreieck so aufgeteilt, dass 4 weitere Dreiecke entstehen (siehe Abbildung 3.7).

Abbildung 3.7 Aufteilung der Dreiecke in vier neue Dreiecke

Wie diese Aufteilung vorzunehmen ist, wird in den folgenden Punkten erläutert. Wenn das Basisnetz in anderer Topologie vorliegt (z.B.: Vierecksnetz) muss es in einem ersten Schritt durch Teilen der Flächen in ein Dreiecksnetz umgewandelt werden.

Regeln für innere ungerade Punkte: Abbildung 3.8 zeigt wie ungerade innere Punkte (Kantenpunkte die in diesem Schritt neu erzeugt werden) generiert werden. Dabei wird pro Kante ein neuer Kantenpunkt aus angrenzenden Punkten entsprechend der Verteilung in Abbildung 3.8 generiert.

Abbildung 3.8 Regel für innere Punkte

28

Regeln für innere gerade reguläre Punkte: Abbildung 3.9 zeigt wie die neue Position von geraden Punkten (bereits bestehende Punkte) aus den angrenzenden Punkten errechnet wird, wenn es sich um einen inneren regulären Punkt handelt.

Abbildung 3.9 Regel für innere gerade reguläre Punkte

Regeln für innere gerade irreguläre Punkte: Abbildung 3.10 zeigt das grundsätzliche Schema für innere gerade irreguläre Punkte. Für β bei gegebener Valenz k schlägt Loop folgendes vor: 1 5 3 1 2   (  (  cos( )) 2 ). Diese Wahl führt nachweislich zu C 1 Stetigkeit. Es sind k 8 8 4 k aber auch andere Werte möglich. Warren schlug z.B.: folgende Belegung vor: Für k  3 : 3 3 ; Für k  3 :   . Dies zeigt auch, dass Subdivision keinesfalls eindeutig ist.  8k 16 Jedes Verfahren bietet in gewisser Hinsicht Möglichkeiten das Glätteverhalten zu beeinträchtigen.

Abbildung 3.10 Regel für innere gerade irreguläre Punkte

29

Regeln für ungerade Randpunkte Abbildung 3.11 zeigt die Regel für ungerade Randpunkte. Diese Regel wird auch zur Generierung von Falten (engl. creases) innerhalb des Netzes verwendet. Dabei werden innere Punkte als Randpunkte deklariert und entsprechend derer Vorschrift berechnet.

Abbildung 3.11 Regel für ungerade Randpunkte

Regeln für gerade Randpunkte Abbildung 3.12 zeigt die Regel für gerade Randpunkte. Wie im vorigen Abschnitt wird dies auch zur Erzeugung von Falten verwendet.

Abbildung 3.12 Regel für gerade Randpunkte

Exakte Berechnung von Randpunkten Die beiden vorherigen Schemata haben den Nachteil dass Randpunkte nicht von ihrer Valenz abhängig sind. Dies führt dazu, dass Punkte ab Valenz > 7 nicht mehr formal C 1 stetig sind. Theoretisch könnte das Problem dazu führen, dass verschiedene Netze entlang eines Randes nicht korrekt verbunden werden, da diese Schnittstelle nicht exakt gleich ist. wenn Randpunkte unterschiedliche Valenz haben. In der Praxis hat sich aber gezeigt, dass dies nicht von Bedeutung ist, da die Abweichung vom exakten Ergebnis in der Regel nicht wahrnehmbar ist. Wenn man aber formal korrekte C 1 stetige Kanten und Falten erhalten will, muss man das Schema noch etwas modifizieren und Regeln für innere ungerade Punkte, welche an einen Randpunkt angrenzen (siehe Abbildung 3.13) einführen. Bei Randpunkten mit einer Valenz von k  7 wird für angrenzende ungerade innere Punkte das Schema aus Abbildung 3.14 verwendet. Eine kompliziertere Variante 1 1 die zu besseren Ergebnissen führt ergibt sich wenn man die Werte und aus 2 4 1 1 2 1 1 2 Abbildung 3.14 durch  cos und  cos ersetzt. 4 4 k 1 2 4 k 1

30

Abbildung 3.13 Exakte Berechnung von Randpunkten

Abbildung 3.14 Exakte Berechnung von Randpunkten

31

Kapitel 4

Implementierung 4.1 Allgemeines Das Programm wurde unter Linux in C++ entwickelt. Zum Einsatz kommen das Open Source Windows-Toolkit wxWidgets, OpenGL [27], die OpenGL shading language GLSL und die OpenGL Extension Wrangler Library GLEW. wxWidgets ermöglicht plattformunabhängige GUI-Programmierung. Zum jetzigen Zeitpunkt läuft FiberNavigator auf verschiedenen Linux-Distributionen, Microsoft Windows und MacOS. GLEW ist eine plattformunabhängige Bibliothek, die ermittelt, welche OpenGLFunktionalität zur Verfügung steht und diese in einem einzelnen Header-File bereitstellt.

4.2 Anzeige von multimodalen MR-Daten Die Datensätze werden in drei Kategorien eingeteilt: -

T1, T2 gewichtete Daten Richtungs- bzw. Vektorinformationen funktionelle Daten, Wahrscheinlichkeitskarten

Die Anzeige der Daten erfolgt nach gängigen Konventionen, wie sie auch in vergleichbaren Programmen benutzt werden. T1, T2-Daten werden als Graustufenbilder angezeigt, Richtungsinformationen als farbige Bilder, und funktionelle Daten sowie Wahrscheinlichkeitskarten als Farbkarten. Als Standardansicht dienen achsenparallele Navigationsschichten. Diese dienen der räumlichen Orientierung im Datensatz und der exakten Positionierung von Objekten.

32

4.2.1 Shaderbasiertes Multitexturing Die Datensätze werden in OpenGL-3D-Texturen, die das gesamte Volumen umspannen, überführt. Alle Objekte können so texturiert werden, ohne dass Texturkoordinaten übergeben werden müssen. Dafür werden die Texturwerte für das jeweilige Voxel in einem Fragment-Shader ermittelt, welcher auch das Multi-Texturing, Clipping und Überblenden der Texturen übernimmt. Die Anzahl der möglichen Texturen wird dabei von der jeweils verwendeten Grafikkarte begrenzt, beträgt jedoch mindestens 10 auf gängigen Modellen. Über Schieberegler lassen sich ein Schwellwert für das Anzeigen und der Alphawert für das Überblenden für jede Textur einstellen.

Die Ermittlung des Farbwertes eines Fragments erfolgt iterativ. Hierbei wird die jeweilige Textur nur berücksichtigt wenn der Betrag des Farbwertes größer als der eingestellte Schwellwert ist. Ist dies der Fall, wird der Texturwert abhängig vom eingestellten Alphawert zum Farbwert des Fragments addiert.

Abbildung 4.1 Linkes Bild T1-Graustufenbild überlagert mit Wahrscheinlichkeitskarte, mittleres Bild FA als Graustufenbild, rechtes Bild FA als Farbkarte

Abbildung 4.1 zeigt die Darstellung verschiedener Datensätze auf einem axialen Schnitt. Links ist ein T1 Graustufenbild überlagert von einer Wahrscheinlichkeitskarte die durch eine Farbkarte repräsentiert wird. Diese Farbkarte entspricht der Aufnahme einer Wärmebildkamera. Hohe Werte werden gelb bis rot dargestellt, niedrige grün bis dunkelblau. Es kann aus mehreren vordefinierten Farbkarten gewählt oder eigene hinzugefügt werden.

33

Daneben ist eine Karte der fraktionellen Anisotropie als Graustufenbild und rechts als Farbkarte abgebildet.

Abbildung 4.2 Linkes Bild Eigenvektordatensatz farbkodiert, rechtes Bild Eigenvektordatensatz farbkodiert und mit FA multipliziert

In Abbildung 4.2 ist im linken Bild die farbige Abbildung der größten Eigenvektoren zu sehen. Da diese Vektoren die Bedingung 𝑣𝑥2 + 𝑣𝑦2 + 𝑣𝑧2 = 1 erfüllen erscheint diese Darstellung sehr hell. Für die Darstellung im rechten Bild wird die fraktionelle Anisotropie mit dem größten Eigenvektor des Voxels multipliziert. Voxel, in denen die Diffusion isotrop ist, erscheinen hier sehr dunkel bis schwarz, was Vektoren in Regionen in denen keine wesentliche Faserausrichtung vorliegt, verschwinden lässt.

4.2.2 Interpolation Interpolierte Texturen erzeugen einen glatten optischen Eindruck, der für die meisten Darstellungen wünschenswert ist. Jedoch ist es für die wissenschaftliche Auswertung oft erforderlich, genaue Informationen über die Werte in einzelnen Voxeln zu erhalten. Aus diesem Grund kann für jede Textur lineare Interpolation ein- oder ausgeschaltet werden (Abbildung 4.3).

34

Abbildung 4.3 Wahrscheinlichkeitskarte als Farbkarte interpoliert (links) und nicht interpoliert (rechts)

35

4.3 Fasern Kern des Programms ist die interaktive Auswahl interessanter Faserbündel. Wichtig für die Akzeptanz der Anwender ist die gefühlte Geschwindigkeit, mit der die Navigation im Datensatz und die Auswahl der Fasern erfolgt. Die zwei Hauptfaktoren sind dabei der Test, welche Linien von den gegenwärtig aktiven Auswahlwerkzeugen erfasst werden und das Rendern der Linien in OpenGL. Da ersteres allein auf der CPU berechnet wird, sind geeignete Datenstrukturen notwendig um diesen Test zu beschleunigen. Für das Rendern der Linien bietet OpenGL mit vertex buffer objects beziehungsweise vertex arrays Funktionalität, die die Möglichkeiten moderner Graphikhardware ausnutzt. Um die gewünschte Interaktivität zu erreichen, ist es nicht möglich eine Region of Interest (ROI) zu definieren und die davon ausgehenden Fasern in Echtzeit zu berechnen. Vielmehr wird mit einem geeigneten Tracking-Algorithmus eine globale Berechnung aller möglichen Fasern vorgenommen und als Menge von Punkten und die Fasern bildenden Punktindizes übergeben. Dazu wird in der weißen Gehirnsubstanz in jedem Voxel eine Stromlinie gestartet. Dieses Vorgehen erzeugt je nach Auflösung des Datensatzes eine große Menge Fasern. Im Fall des mir zur Verfügung stehenden Testdatensatzes sind dies ca. 400 000 Linien bestehend aus ca. 29 Millionen Punkten. Aus dieser Menge an Fasern werden mit geeigneten Selektionswerkzeugen die interessanten ausgewählt. Ein hybrider Ansatz aus rechteckigen Auswahlboxen [28] und beliebig geformten Regionen, gegeben durch funktionelle Datensätze und Wahrscheinlichkeitskarten, hat sich dabei als sinnvoll erwiesen. Sowohl die ROI als auch die top-level Auswahlboxen können mit weiteren Auswahlboxen logisch verknüpft werden. Zur Gewinnung der Daten wurden hochauflösende DTI Messungen mit einem Siemens Trio Scanner (1.7mm isotrop, 60 Richtungen, b=1000s/mm2, GRAPPA/2, NEX=3) durchgeführt. Nach Bewegungskorrektur und Registrierung auf die T1-Anatomie wurden die Diffusionstensoren berechnet. Aus diesen Daten werden die Fasern werden mit MedINRIA [29] berechnet. Die Anzahle der erzeugten Fasern ist abhängig von der gewählten Auflösung des DTI-Bildes. Bei einer Interpolation auf 1 mm entstehen die o.a. 400 000 Fasern. Mit einer gröberen Auflösung kann ein entsprechend ausgedünnter Datensatz erzeugt werden.

36

4.3.1 Auswahl mittels Box

Wichtig für die Auswahl mittels Boxen ist intuitive Positions- und Größenveränderung, klare visuelle Veranschaulichung der Funktion, sowie einfache Verknüpfung mittels logischer Funktionen. Die Boxen werden in Saat- oder Masterboxen und Verfeinerungsboxen unterschieden. Die hellblau gezeichneten Boxen, im Weiteren als Masterboxen bezeichnet, fungieren als Saatregion für die Linienselektion (Abbildung 4.4). Alle Linien, die über mindestens einen Punkt innerhalb des Volumens der Box verfügen, werden ausgewählt.

Abbildung 4.4 Faserselektion mit einer Masterbox

37

Für unsere Zwecke ist es jedoch nötig diese recht grobe Auswahl weiter zu verfeinern. Dazu können weitere Boxen mit einer Masterbox mittels der logischen Funktion UND und Nicht-UND verknüpft werden. Abbildung 4.5 zeigt die Verknüpfung der in Abbildung 4.4 gezeigten Masterbox mit einer UND-Box, um Fasern auszuwählen, die durch zwei bestimmte Region verlaufen, und Abbildung 4.6 eine weitere Verknüpfung mit einer Nicht-UND-Box um unerwünschte Fasern auszublenden.

Abbildung 4.5 Faserselektion mit Masterbox verfeinert durch eine UND-Box

Abbildung 4.6 Faserselektion mit Masterbox verfeinert durch eine UND-Box und Entfernung unerwünschter Fasern mit NICHT-Box

38

Beliebig viele Masterboxen und mit ihnen verknüpfte UND- und NICHT-Boxen können mittels ODER-Funktion verbunden werden. Die Boxen können mit Hilfe der Maus und der Tastatur frei im Bildraum bewegt und in ihrer Größe verändert werden.

Abbildung 4.7 Kombination mehrerer Auswahlboxen zur Auswahl unterschiedlicher Faserbündel

Die in Abbildung 4.7 gezeigte Selektion zweier Faserbündel mittels fünf miteinander verknüpfter Auswahlboxen entspricht der logischen Formel 𝑏1 ∧ 𝑏2 ∨ (𝑏3 ∧∼ 𝑏4 ∧ 𝑏5).

39

4.3. 2 kd-Baum

Nach dem Laden des Faserdatensatzes wird für die Punktmenge ein kd-Baum [2] aufgebaut. Dieser ermöglicht die schnelle Entscheidung ob sich ein Punkt in einem beliebigen rechtwinkligen, achsenparallelen Volumen befindet. Bei Erstellung, Lageoder Größenveränderung einer Auswahlbox wird für jeden Punkt des Faserdatensatzes getestet ob er sich innerhalb des Volumens der Box befindet. Ist dies der Fall, wird in einer invertierten Liste der Index der zu diesem Punkt gehörigen Linie nachgeschlagen und die Linie in einem Bitfeld als ausgewählt markiert.

Abbildung 4.8 Die besuchten Knoten des kd-Baum für die Auswahl des cingulum [2]

40

4.3. 3 Auswahl über ROI

Der Nachteil an der Auswahl mittels Boxen ist ihre starre geometrische Form, die es stellenweise doch recht umständlich gestaltet interessante Regionen zu modellieren. Desweiteren setzt sie a priori Wissen über den Faserverlauf voraus. Die Modellierung und Manipulation beliebig geformter Regionen in 3D gestaltet sich jedoch sehr aufwendig. Eine Möglichkeit die Auswahl auf beliebig geformte Regionen auszudehnen bieten funktionelle Datensätze, zum Beispiel aus Wahrscheinlichkeitskarten oder fMRT. In diesen Datensätzen ist jedem Voxel ein skalarer Wert zwischen 0 und 1 zugewiesen. Für jeden Punkt im Faserdatensatz wird ein Schwellwert vom Wert des zugehörigen Voxels abgezogen. Ist das Ergebnis größer als Null, wird die zugehörige Linie als aktiv markiert. Der Schwellwert kann interaktiv über einen Regler verändert werden, so dass nur Regionen mit einer gewissen Wahrscheinlichkeit oder einer besonders starken Aktivierung als Saatregionen der Linien gewählt werden.

Abbildung 4.9 Iso-Oberfläche auf Wahrscheinlichkeitskarte

41

Abbildung 4.9 zeigt eine Iso-Oberfläche auf einer Wahrscheinlichkeitskarte mit einem gewählten Schwellwert von 0,6. In Abbildung 4.10 werden die Regionen mit der stärksten Aktivierung in einem fMRT-Datensatz zur Auswahl herangezogen.

Abbildung 4.10 Faserselektion mit Wahrscheinlichkeitskarte

42

4.3. 4 Farbgebung Oftmals ist es wichtig zusätzlich zur räumlichen Lage noch andere Information darzustellen. Bevorzugtes Mittel ist dabei Farbgebung. Fasern können im FiberNavigator auf unterschiedliche Weise eingefärbt werden. Dies ist: -

Richtungsvektor zwischen Start- und Endpunkt (globale Richtung) (Abb.4.11) Tangentenvektor (lokale Richtungsfarbgebung) (Abb.4.12) manuelle Färbung selektierter Faserbündel (Abb.4.13) Farbkarte auf selektierten skalaren Daten (Abb.4.14)

Abbildung 4.11 Globale Färbung

Abbildung 4.12 Lokale Färbung

Abbildung 4.13 Manuelle Färbung

Abbildung 4.14 Färbung mit FA

43

4.3.5 Rendering Das Zeichnen der Linien erfolgt über den OpenGL-Befehl GL_LINE_STRIP. Um die Anzahl der Funktionsaufrufe zu vermindern und die Fähigkeiten moderner Grafikkarten auszunutzen werden Vertex Buffer Objects (VBO) [30] eingesetzt. Vertex Buffer Objects sind Teil des OpenGL Standards seit Version 1.5, welcher 2003 vorgestellt wurde. Dafür werden die verwendeten Felder für Positionen der Punkte, Normalen und Farben in den Grafikkartenspeicher kopiert. Der Befehl zum Zeichnen einer Linie benötigt dann nur den Startindex und die Anzahl der Punkte. Für ältere Grafikkarten, auf denen VBO’s nicht zur Verfügung stehen, wird auf Vertex Arrays zurückgegriffen. Hierbei werden die verwendeten Felder, für Lesezugriffe optimiert, im Hauptspeicher abgelegt. Der Funktionsaufruf zum Zeichnen einer Linie gleicht dem für VBO’s. Beleuchtung ist wichtig für die räumliche Wahrnehmung dreidimensionaler Strukturen. Die Abbildung 4.15 und 4.16 zeigen das gleiche Faserbündel unbeleuchtet und beleuchtet. Wie man im linken Bild sieht, erscheinen die unbeleuchteten Linien als homogene Masse, während im beleuchteten Fall Strukturen zu erkennen sind. Für Linien wird dabei die Standard-OpenGL-Beleuchtung verwendet. Zur Vereinfachung wird hier die Tangente der Linie an jedem Vertex als Normale übergeben. Die Lichtquelle befindet sich im Punkt des Betrachters, was man sich auch als Helmlampe vorstellen kann.

Abbildung 4.15 Unbeleuchtete Fasern

Abbildung 4.16 Beleuchtete Fasern

44

4.3.6 Rendering der Fasern als Röhren Während das Rendern der Fasern als Linien in vielen Fällen ausreichend ist, bietet die Darstellung als Röhren eine bessere räumliche Vorstellung und ein optisch ansprechenderes Ergebnis. Dies wird jedoch durch einen wesentlich höheren Aufwand an zu rendernder Geometrie und damit verbundener Anzahl von OpenGL-Befehlen erkauft. Eine Möglichkeit dies zu umgehen und trotzdem den Eindruck von Röhren zu vermitteln wird in [8] beschrieben.

Abbildung 4.17 Kompletter Faser-Datensatz als Röhren gerendert

Hierfür wird die Linie nun als GL_QUAD_STRIP gerendert. Jeder Punkt wird doppelt übergeben und erhält einen binären Parameter 1 oder -1. Zusätzlich ist für jeden Punkt die Linientangente verfügbar. Da diese für die Einfärbung mit der lokalen Richtungsinformation schon zur Verfügung steht, erzeugt dies keinen zusätzlichen Aufwand. In einem Vertex-Shader wird nun das Kreuzprodukt des Sichtvektors 𝑣𝑖 mit der Linientangente 𝑡𝑖 gebildet und mit dem binären Parameter und einem Offsetwert multipliziert. Diese Operation verschiebt den Punkt senkrecht zur Linientangente und dem Sichtvektor um einen festgelegten Wert. Dies erzeugt ein GL_QUAD_STRIP,

45

dessen Eckpunkte jeweils konstanten Abstand zur ursprünglichen Linie haben, also ein Band konstanter Breite. Um den Eindruck einer Röhre zu erzeugen, wird dieses Band mit einer 1D-Textur, deren Helligkeitsverlauf der Sinusfunktion entspricht, texturiert. Im Gegensatz zur ursprünglichen Beschreibung des Verfahrens in [8] wird in der hier verwendeten Implementierung die Textur nicht durch OpenGL übergeben sondern algorithmisch im Fragment-Shader erzeugt.

Abbildung 4.18 Konstruktion des Quadstrips, entsprechend der Linienpunkte werden Eckpunkte für den Quadstrip generiert

4.3.7 Betrachtung des Speicherverbrauches

Schnelles Rendern und Interaktion werden im Normalfall mit erhöhtem Speicherverbrauch erkauft. Anhand des Testdatensatzes mit ca. 400.000 Fasern und ca. 29 Millionen Punkten möchte ich kurz den Speicherverbrauch diskutieren. Hierbei gehe ich von 4 Byte pro int und float aus. Für das Rendern der Linien werden 3 Arrays aus jeweils 3 float-Werten, pro Punkt vorgehalten. Dies bedeutet 36 Byte pro Punkt und ca. 1.1 Gigabyte für den gesamten Datensatz. Der kd-Baum und die invertierte Liste der zu jeden Punkt gehörigen Linie werden jeweils in einem Integer-Array gespeichert, was zusätzlich ca. 230 MB verbraucht. Für den gesamten Datensatz, jede Auswahlbox und jede Gruppe verknüpfter 46

Boxen wird jeweils ein Bitfeld der Größe der Anzahl Linien angelegt. Dies verbraucht bei optimaler Speicherung ca. 50k Byte pro Bitfeld. Da im Gegensatz zu Bytedaten für den Zugriff auf die einzelnen Bit Maskierungs- und Verschiebeoperationen notwendig sind, ist hier sogar noch Optimierungspotenzial auf Kosten von zusätzlichem Speicherverbrauch vorhanden. Es ist zu sehen, dass allein der Faserdatensatz mit der nicht unüblichen Anzahl von 1015 Auswahlboxen circa 1.5 GB verbraucht. Sind noch weitere Datensätze geladen, werden leicht die in heutigen Arbeitsplatzrechnern verwendeten 2 GB erreicht. Durch Auslagerungsvorgänge ist dann meist kein angenehmes Arbeiten mehr möglich. Es empfiehlt sich daher die Positionierung der Auswahlboxen und Schnittebene mit einem ausgedünnten Faserdatensatz vorzunehmen.

47

4.4 Isoflächen Die Rekonstruktion einer Oberfläche aus den Skalardaten der T1- oder T2-Scans ist eine beliebte Technik. Hierbei ist zu beachten, dass natürlich nur drei Regionen sinnvolle Ergebnisse liefern. Diese sind die Übergänge von Umgebung zu Gesichtsgewebe, das Gehirn umgebende Flüssigkeit zu grauer Materie und graue Materie zu weißer Materie. Außerdem können Iso-Oberflächen aus Wahrscheinlichkeitskarten (Abbildung 4.19, links) erzeugt werden. Dies ergibt eine dreidimensionale Visualisierung der Lage der interessanten Regionen im Gehirn. So erzeugte Oberflächen aus Daten aus dem probabilistischen Tracking können mit der Faserrichtung texturiert werden und ergeben eine realistischere Ansicht der Faserbündel (Abb.4.19, rechts). Auch können so Regionen, die als Saatregionen für die Faserauswahl verwendet werden, sichtbar gemacht werden.

Abbildung 4.19 Iso-Oberfläche auf Wahrscheinlichkeitskarte eingefärbt mit Eigenvektorrichtung (links)und zusätzlich texturiert mit LIC (rechts)

Für die Erzeugung der Oberfläche kommt der 1987 von Lorensen und Cline vorgestellte Marching Cubes Algorithmus zum Einsatz [8]. Der zur Erstellung verwendete Iso-Wert kann über einen Schieberegler verändert werden. Die Lage und Größe der zu erzeugenden Iso-Fläche kann damit zum Beispiel den zur Anzeige verwendeten Schwellwerten angeglichen werden.

48

4.4.1 Verbesserung

Durch Messfehler und Rauschen im Datensatz entstehen bei der Erstellung der Oberfläche Artefakte, die sich als optisch störend erweisen (Abbildung 4.20). Um diese zu vermeiden wurden zwei einfache Verfahren implementiert.

Abbildung 4.20 Artefakte nach Erstellung einer Iso-Oberfläche auf T1-Daten

Im Ersten wird der Datensatz, auf dem die Oberfläche erzeugt wird, vor der Erstellung der Oberfläche mit einem 3x3x3-Medianfilter geglättet. Dies erzeugt ein optisch ansprechendes Ergebnis, in dem so gut wie keine Artefakte enthalten sind. Leider „verfälscht“ dieses Verfahren die Werte im Datensatz, so dass zwar optisch gute Bilder entstehen, die Oberfläche aber etwas von der realen Position abweicht (Abbildung 4.21).

49

Abbildung 4.21 Iso-Oberfläche auf T1-Daten, die vorher mit Medianfilter geglättet wurden

Im zweiten Verfahren werden nach der Oberflächenerstellung die Dreiecke entsprechend ihrer Zugehörigkeit zu zusammenhängenden Objekten geclustert und nur das größte Objekt beibehalten. Da dieses Verfahren die Nachbarschaftsinformation der Dreiecke benutzt, funktioniert es sehr gut für frei im Raum befindliche Artefakte, versagt aber bei Verunreinigungen in Oberflächennähe, wo Verbindungen zur Oberfläche bestehen.

50

Abbildung 4.22 Iso-Oberfläche auf T1-Daten mit Artefaktentfernung nach der Erstellung

51

4.5 Schnittfläche Für eine Klingler-Dissektion an einem Gehirn wird dieses an einer Ebene geschnitten und dann weitere Gehirnsubstanz entfernt, um die interessanten Faserbündel freizulegen. In unserem virtuellen Fall bedeutet dies, die initiale Schnittfläche wird durch den Nutzer vorgegeben und automatisch anhand der selektierten Faserbündel verformt. Da dies nicht in jedem Fall optimale Ergebnisse liefert ist es möglich die Schnittfläche dann noch manuell zu verändern. Für die Erzeugung der Schnittfläche wurde die NURBS-Repräsentation gewählt, da sie -

wenige Kontrollpunkte benutzt und damit leicht zu speichern ist, nachträgliches Editieren erlaubt, effiziente Methoden des Renderns existieren.

Die initiale Schnittebene wird durch die Position einer Navigationsschicht vorgegeben. Darauf werden in regelmäßigen Abständen de Boor Kontrollpunkte platziert. Nun wird in einem vorgegebenen Radius um jeden Kontrollpunkt die Punktmenge der selektierten Faserbündel analysiert und der Schwerpunkt dieser Punktmenge interpoliert. Die de Boor Kontrollpunkte werden dann auf den jeweiligen Schwerpunkt verschoben. Dies erzeugt eine glatte Oberfläche, die sich der selektierten Faserstruktur gut anpasst. Zur Optimierung des erreichten Ergebnisses ist es möglich, die Kontrollpunkte nachträglich zu verschieben und weitere einzufügen. Aus der Spline-Repräsentation der Schnittfläche wird durch regelmäßige Abtastung eine Vergitterung erzeugt. Um die gewünschte Interaktivität zu erreichen, wird in der Editierphase ein relativ grob aufgelöstes Dreiecksgitter verwendet. Die benötigte Verfeinerung für die LIC-Texturierung erfolgt durch mehrmaliges Anwenden des Loop Subdivision Algorithmus.

4.5.1 Clipping einer Oberfläche Um zusätzlich den Eindruck eines Schnittes durch ein Gehirn zu erzielen, werden IsoOberflächen an der Schnittfläche abgeschnitten (Abbildung 4.23). Während dies für ebene Flächen mittels OpenGL Befehlen einfach möglich ist, ist dies für beliebig verformte Oberflächen nicht der Fall. Von Vorteil ist jedoch, dass die Fläche eine generelle Ausrichtung entlang einer Achse hat und zwei Flächenpunkte auf dieser 52

Achse nicht übereinanderliegen können. Damit ist es möglich die Fläche in eine Tiefentextur zu rendern und diese dem Fragment-Shader für die Iso-Oberfläche zu übergeben. Der prüft, ob der Koordinatenteil dieser Achse über diesem Tiefenwert liegt und verwirft die entsprechenden Fragmente.

Abbildung 4.23 Clipping von Iso-Oberfläche an Schnittfläche

53

Kapitel 5

Ergebnisse 5.1 Beispiel einer virtuellen Klingler-Dissektion

Abbildung 5.1 Endergebnis der virtuellen Klingler-Dissektion

Abbildung 5.1 zeigt Ergebnis der Anwendung des Verfahrens zur Erzeugung einer virtuellen Klingler-Dissektion. In diesem Abschnitt sollen die Schritte, die zu diesem Endergebnis führen aufgezeigt werden. 54

Obwohl es durchaus möglich ist, eine LIC-texturierte Schnittfläche nur mit einem Eigenvektordatensatz zu erzeugen, empfiehlt es sich noch mindestens ein T1-Bild und einen Faserdatensatz zu laden. Erster Schritt ist die Auswahl des gewünschten Faserbündels. Abbildung 5.2 zeigt eine Auswahl der cortico-cerebellaren Fasern. Zwei Nicht-Boxen dienen der Entfernung ungewünschter Fasern.

Abbildung 5.2 Faserauswahl (cortico-cerebellaren Fasern)

Für den initialen Schnitt wird nun die sagittale Navigationsebene in die Nähe des gewählten Faserbündels bewegt. Abbildung 5.3 zeigt die Positionierung der Navigationsebene, so dass Fasern sowohl vor und auch hinter der Schnittebene liegen. Ausgehend von der Navigationsebene wird eine verformte Schnittfläche erzeugt, die bereits dicht an den ausgewählten Faserbündeln liegt (Abbildung 5.4).

55

Abbildung 5.3 Positionierung der Navigationsebene

Abbildung 5.4 Erzeugte Schnittfläche

Über einen Schieberegler kann die Lage der Kontrollpunkte weiter verändert werden, was die Verformung der Fläche global modifiziert.

Abbildung 5.5 Editiermodus für die Schnittfläche

56

Ist dies nicht ausreichend können die Punkte im Editiermodus separat verschoben oder zusätzliche hinzugefügt werden. Nachdem die Position der Schnittfläche fixiert wurde kann die LIC-Textur hinzugefügt werden (Abbildung 5.6). Dieser Vorgang dauert auf Grund von mehrmaligem Anwenden des Loop-Subdivision Algorithmus und der Stromlinienintegration ca. 20-30 Sekunden. Im rechten unteren Bereich wurde keine LIC-Textur gerendert, da keine Vektordaten für diese Region (Kleinhirn) vorliegen.

Abbildung 5.6 Verformte Schnittfläche texturiert mit T1 und LIC

57

5.2 Anwendungen 5.2.1 Vergleichende Morphometrie eines Faserbündels für verschiedene Probanden

Der exakte Verlauf der verschiedenen Faserbündel, deren Form und Stärke unterscheidet zwischen den verschiedenen Versuchspersonen. Eine Studie der normalen Variabilität eines Bündels erlaubt u.a. die Untersuchung von anormalen Situationen bei Patienten. Abbildung 5.7 zeigt den Vergleich des Verlaufs des Fasciculus Uncinatus (UNC) in sechs Versuchspersonen. Das Faserbündel wurde dabei jeweils mit zwei „UND“ verknüpften Auswahlboxen bestimmt.

Abb. 5.7 Vergleich des Verlaufs des Fasciculus Uncinatus (UNC) in sechs Versuchpersonen.

58

5.2.2 Segmentierung von Unterkomponenten eines Faserbündels

Durch die logische Verknüpfung mehrere Boxen lassen sich Unterkomponenten von Faserbündeln analysieren. Abbildung 5.8 zeigt die Zerlegung des Fasciculus Longitudinales Superior (SLF) in einen horizontalen (grün) und vertikalen (gelb) Anteil. Der horizontale Anteil wurde durch eine Auswahlbox im inferioren Parietallappen und im ventralen Prä-motorischen Kortex gewählt. Der vertikale Anteil entsteht durch Kombination der parietalen Box mit einer Auswahlbox im Schläfenlappen. Durch Verknüpfung beider Endboxen entsteht der durchgehende Anteil, der auch als Fasciculus Arcuatus (AF) bezeichnet wird.

Abb. 5.8 Visualisierung der Unterkomponenten des SLF vor einem sagittalen Schnittbild der farbkodierten und FA gewichteten Eigenvektorrichtung.

5.2.3 Visualisierung von Eigenschaften der weißen Masse auf Faserbündeln

Die Untersuchung der FA-Werte entlang eines Faserbündels kann z.b. zur Studie des Entwicklungsstadiums eines Faserbündels verwendet werden. Dazu wurden 7-jährige Kinder mit einer Gruppe von Erwachsenen verglichen und Bereiche entlang einzelner Faserbündel identifiziert, die reduzierte FA-Werte bei Kindern zeigten.

59

Abb. 5.9 Farbkodierte FA-Werte entlang der Fasern des internen Kapselsystems.

5.2.4 Evaluierung von Ultra-Hochfeld 7-Tesla diffusionsgewichteten Daten Die Entwicklung neuer Bildgebungssequenzen auf dem 7-Tesla MR-Tomographen erfordern eine schnelle Visualisierung einzelner Faserbündel zur Evaluierung der Bildqualität. Abbildung 5.10 zeigt ein Beispiel einer kombinierten Darstellung mehrerer Faserbündel in Kombination mit anatomischen Schnittbildern erzeugt aus 7-Tesla DTIDaten die mit hoher räumlicher Auflösung ( 1,4x1,4x1,4 mm3 ) aufgenommen wurden.

Abb. 5.10 Faserbündel aus 7 Tesla DTI Bildern.

60

5.2.5 Analyse der anatomischen Verbindungen zwischen funktionell aktivierten Hirnregionen

Die Darstellung der Faserverbindungen zwischen funktionell aktivierten Hirnregionen ermöglich die Identifikation von Netzwerken verknüpfter Regionen die potentiell in bei der Verarbeitung einer Aufgabe im Gehirn zusammenarbeiten. Die Abbildung 5.11. zeigt die Faserauswahl mit einer aus fMRT Daten erzeugten Auswahlregion (rote Isoflächen) die unter anderem den linken und rechten auditiven Kortex beinhaltet. Die Fasern die durch diese Regionen verlaufen werden zusammen mit einer transparenten Offsetfläche und einem axialen Schnitt gezeigt.

Abb. 5.11 Faserauswahl anhand von funktionell aktivierten Regionen ( rote Isoflächen)

5.2.6 Multimodale Visualisierung anhand von Iso-Oberflächen a) Funktionelle MRT Bilder müssen zusammen mit den anatomischen Strukturen visualisiert werden. Eine neue Möglichkeit ist die Darstellung der fMRT Daten als Farbtextur auf einer nach innen verschobenen Offsetfläche zur Hirnoberfläche, die mit der T1 texturiert ist. Eine zusätzliche Iso-Oberfläche der funktionell aktivierten Areale gibt zusätzliche Information über deren räumliche Anordnung und Ausdehnung. (Abbildung 5.12)

61

Abb 5.12 Akustisch stimulierte fMRT-Aktivierung sprachrelevanter Hirnareale

b) Mit probabilistischer Traktografie erzeugte Konnektivitätskarten (Traktogramme) können als Iso-Oberflächen dargestellt werden. Abb. zeigt die kombinierte Darstellung der Traktogramme von drei Subregionen im inferioren Parietallappen. Die Startregionen sind umrandet. Die Faserbündel durchdringen die Offsetfläche und erlauben eine genaue Lokalisation der kortikalen Projektionen (rechts). Eine transparente IsoOberfläche der T1-Anatomie zeigt den 3D-Verlauf der Faserbündel und deren räumliche Anordnung (Abbildung 5.13).

Abb. 5.13 Traktogramme von drei Startregionen im inferioren Parietallappen (umrandet).

62

5.2.7 Analyse anatomischer Eigenschaften der weißen Substanz mithilfe der LIC-Texturierung

Die virtuelle Klingler Dissektion erlaubt die Visualisierung von Strukturen der weißen Substanz auf nicht achsenparallelen gekrümmten Schnittflächen. Die LIC Texturierung der Oberfläche mit der Faserrichtung kombiniert mit deren Farbkodierung erlaubt die Untersuchung komplexer Fasersysteme. Abbildung 5.14 zeigt eine virtuelle Klingler Dissektion der sub-insularen weißen Substanz (rote Umrandung). Dabei zeigt sich eine Einteilung dieser Region in einen vorderen und hinteren Bereich mit unterschiedlichen Faserverläufen.

Abb. 5.14 Virtuelle und post-mortem Dissektion der sub-insularen weißen Substanz und deren Verbindungen.

63

5.2.8 Kombination mehrerer Faserbündel in der virtuellen Klingler Dissektion

Viele Faserbündel liegen in engem räumlichen Bezug und können bei der virtuellen Klingler Dissektion effizient kombiniert visualisiert werden. Abb 5.15. zeigt den Verlauf von drei Faserbündeln, die wichtige Aufgaben in der Sprachverarbeitung haben und den Schläfenlappen mit dem Frontalhirn verbinden. Die gute Korrespondenz zu der post-mortem Fotographie (Ludwig, 1956) zeigt die Robustheit der Methode.

Abb. 5.15 Klingler Dissektion (oben, Ludwig, 1956) und virtueller Schnitt des sprachrelevanten Fasersystems bestehend aus Fasciculus Longitudinales Superior (SLF, türkis, 1; 6); Fasciculus Uncinatus (UNC, gelb, 3); Fasciculus Occipitofrontalis Inferior (IFO, violett, 5)

64

Kapitel 6

6. Zusammenfassung und Ausblick Die effiziente Implementierung und die intuitiven Werkzeuge zur Navigation in Datensätzen und Auswahl von Faserbündeln ermöglicht es schon auf einem heute üblichen Standard-PC detailierte Gehirnstrukturen zu erforschen. Die Kombination von global implementierter Faserkonstruktion und effizienten Auswahlwerkzeugen und leicht verständlichen booleschen Operationen führt zu so noch nicht vorhandenen Möglichkeiten der Darstellung von Faserbündeln. Wissenschaftler am MPI-CBS in Leipzig benutzen das Tool bereits zur Selektion von Faserbündeln in unterschiedlichen Datensätzen aus Gehirnscans. Das Einfärben der Faserbündel mit gemessenen Parametern der weißen Gehirnsubstanz wie der fraktionellen Anisotropie ermöglicht quantitative Darstellung der Gewebeeigenschaften und statistische Auswertungen der Gehirnstrukturen. Die multimodale Integration von DTI, MRT und fMRT Bildern komplettiert das System. Die Berechnung einer Fläche, welche den ausgewählten Nervenbahnen folgt, erzeugt eine virtuelle Klingler-Dissektion. Die Methode deformiert eine Schnittfläche so, dass diese lokal parallel zu den gewählten Faserbündeln verläuft. Mittels dieser Schnittfläche, ist es möglich, die die Fasern umgebenden anatomische Strukturen und ihre räumlichen Verhältnisse zu zeigen. Typische Faserbündel im Gehirn sind in ihrem zentralen Teil kompakt und fächern sich an den Enden in den kortikalen Regionen auf. Die Schnittfläche wird daher lokal an den berechneten Schwerpunkt der selektierten Fasern angeschmiegt. Die Gehirnsubstanz vor der Fläche wird abgeschnitten wie dies auch in einem realen Schnitt durch das Gehirn zur Erzeugung eines medizinischen Präparates erfolgen würde. Das Angleichen der Schnittfläche an die selektierten Fasern erzeugt immer ein interpretierbares Ergebnis. Die Qualität des Resultates hängt natürlich von der Qualität der Selektion der Faserbündel ab. Durch die interaktive Möglichkeit der Optimierung der Schnittfläche ist dies jedoch für den Anwender leicht verständlich und mit wenig Übung möglich. Die Texturierung der Schnittfläche mit der Faserrichtung bietet dem Anwender einen guten Eindruck des Faserverlaufs auf der Fläche. Sie ermöglicht einen intuitiven Eindruck der Faserstrukturen und bietet ein besseres Verständnis der sich auffächernden Faserstrukturen, die mit normaler Faserdarstellung schwieriger zu analysieren sind. Die breite Auswahl an Anwendungsszenarien zeigt, dass durch die Kombination verschiedener Visualisierungsmethoden Ergebnisse erreicht werden, die weit über die virtuelle Klingler-Dissektion hinausgehen. Für die Zukunft ist eine Weiterentwicklung 65

der Software geplant um weitere Datenmodalitäten zu unterstützen. Fokus wird dabei wie bisher auf Interaktivität auf handelsüblichen Computern aus dem Consumerbereich liegen. Das Verfahren der virtuellen Klingler-Dissektion wurde zu den Konferenzen EuroVis 2009 [30] und Human Brain Mapping 2009 eingereicht [31].

66

Anhang A

Quelltextfragmente A.1 Methode zum Erzeugen eines kd-Baums auf einem STL Vector // left, right boundaries of tree // axis: 0-2 alignment of points for sort void KdTree::buildTree(int left, int right, int axis) { // abort condition, we're at a leaf if (left >= right) return; // find the root node for this branch int div = ( left+right )/2; // find median std::nth_element( m_tree.begin()+left, m_tree.begin()+div, m_tree.begin()+right, lessy( m_pointArray, axis ) ); // build sub trees for booth sides buildTree(left, div-1, (axis+1)%3); buildTree(div+1, right, (axis+1)%3); } struct lessy { float const * const data; const int pos; lessy( float const * const data( data ),pos( pos ) { }

data, const int pos ):

bool operator()( const wxUint32& a, const wxUint32&b ) const { return data[ 3*a+pos ] < data[ 3*b+pos ]; } };

Diese rekursive Funktion erzeugt eine als kd-Baum verwendbare Sortierung eines Indexarrays. Die Funktion wird mit den linken und rechten Grenzen und einer Startachse aufgerufen. Die eigentliche Arbeit übernimmt der STL-Befehl nth_element. Dieser erwartet drei Iteratoren sowie bei Bedarf eine benutzerdefinierte Vergleichsfunktion. nth_element sucht das gewählte n-te Element und sortiert das 67

gesamte Array partiell, so daß alle Elemente vor dem gewählten n-ten Element kleiner sind und die Elemente dahinter größer. Da in unserem Fall die Sortierung der Punkte jeweils abwechselnd über die drei Achsen im Raum erfolgt, wird eine benutzerdefinierte Vergleichsfunktion benötigt und für jeden rekursiven Aufruf die aktuelle Achse übergeben.

A.2 Rendering der Linien mit Vertex Buffer Objects Das Zeichnen der Stromlinien geschieht mittels Vertex Buffer Objects. Das folgende Quelltextfragment zeigt die dafür nötigen Funktionsaufrufe. Die drei Arrays m_pointArray, m_colorArray und m_normalArray enthalten die Punktinformationen, Farbwerte und Normalen. glGenBuffers(3, m_bufferObjects); glBindBuffer(GL_ARRAY_BUFFER, m_bufferObjects[0]); glBufferData(GL_ARRAY_BUFFER, sizeof(GLfloat)*m_countPoints*3, m_pointArray, GL_STATIC_DRAW ); glBindBuffer(GL_ARRAY_BUFFER, m_bufferObjects[1]); glBufferData(GL_ARRAY_BUFFER, sizeof(GLfloat)*m_countPoints*3, m_colorArray, GL_STATIC_DRAW ); glBindBuffer(GL_ARRAY_BUFFER, m_bufferObjects[2]); glBufferData(GL_ARRAY_BUFFER, sizeof(GLfloat)*m_countPoints*3, m_normalArray, GL_STATIC_DRAW );

glEnableClientState(GL_VERTEX_ARRAY); glEnableClientState(GL_COLOR_ARRAY); glEnableClientState(GL_NORMAL_ARRAY); glBindBuffer(GL_ARRAY_BUFFER, m_bufferObjects[0]); glVertexPointer(3, GL_FLOAT, 0, 0); glBindBuffer(GL_ARRAY_BUFFER, m_bufferObjects[1]); glColorPointer (3, GL_FLOAT, 0, 0); glBindBuffer(GL_ARRAY_BUFFER, m_bufferObjects[2]); glNormalPointer (3, GL_FLOAT, 0, 0); for ( int i = 0 ; i < m_countLines ; ++i ) { if (m_inBox[i] == 1) glDrawArrays(GL_LINE_STRIP, getStartIndexForLine(i), getPointsPerLine(i)); } glDisableClientState(GL_VERTEX_ARRAY); glDisableClientState(GL_COLOR_ARRAY); glDisableClientState(GL_NORMAL_ARRAY);

68

A.3 Shader-Code für Pseudo-Röhren vertex shader varying varying varying varying

float tangent_dot_view; vec3 tangentR3; float s_param; vec4 myColor;

void main() { // make clipping planes work gl_ClipVertex = gl_ModelViewMatrix * gl_Vertex; vec3 tangent; float thickness = 0.005; tangentR3 = gl_Normal; tangent vector // transform our tangent vector tangent = (gl_ModelViewProjectionMatrix * vec4(gl_Normal,0.)).xyz; // store texture coordinate for shader s_param = gl_MultiTexCoord0.x; vec3 v = (gl_ModelViewMatrix * vec4(0., 0., -1., 0.)).xyz; vec3 offsetNN = cross(v, normalize(tangent.xyz)); vec3 offset = normalize(offsetNN); tangent_dot_view = length(offsetNN); // transform position to eye space vec4 pos = ftransform(); // add offset in y-direction (eye-space) pos.xyz += offset * (s_param * thickness); myColor = gl_Color; // store final position gl_Position = pos; }

69

fragment shader varying varying varying varying uniform

vec3 tangentR3; // Tangent vector in world space float s_param; // s parameter of texture [-1..1] float tangent_dot_view; vec4 myColor; bool globalColor;

/* * simple fragment shader that does rendering of tubes with diffuse illumination */ void main() { vec3 color; if (globalColor) color = abs(normalize(myColor.rgb)); else color = abs(normalize(tangentR3)); vec3 view = vec3(0., 0., -1.); float view_dot_normal = sqrt(1. - s_param * s_param) + .1; float s = (s_param + 1.) * 0.5; //< normalize in [0..1] gl_FragColor.rgb = clamp(view_dot_normal * (color + 0.15 * pow( view_dot_normal, 10.) * pow(tangent_dot_view, 10.) * vec3(1., 1., 1.)), 0., 1.); //< set the color of this fragment (i.e. pixel) gl_FragColor.a = 1.; }

70

Anhang B

Beschreibung des Nutzerinterface

Abbildung B.1 Ansicht der Benutzeroberfläche

Die Oberfläche (Abbildung B.1) ist in verschiedene Bereiche aufgeteilt. Den größten Anteil nimmt das Hauptanzeigefenster (1) ein. Standardmäßig sind hier die drei Navigationsschichten, texturiert mit den geladenen Datensätzen, zu sehen. Die drei Ansichten (2-4) bieten schnelle und genaue Navigation im Datensatz. Das Baum-Widget (5) zeigt allgemeine Informationen zu den geladenen Datensätzen, bietet Shortcuts zum Laden der Datensätze und listet die Auswahlboxen für Fasern und Kontrollpunkte der Schnittfläche. Über kontextsensitive Menüs können Boxen gelöscht, neue hinzugefügt und die Anzeigemodi der Boxen geändert werden. Das List-Widget (6) zeigt alle geladenen Datensätzen und die Einstellungen für diese. 71

Rotation der Ansicht Bewegen der Maus bei gedrückter linker Maustaste Zoom Der Zoom lässt sich mittels Mausrad oder der „+“- und „-„-Taste verändern. Verschieben der Navigationsschichten Die Navigationsschichten können auf verschiedene Art bewegt werden. -

durch Anfassen mit der rechten Maustaste und Verschieben der Maus in die gewünschte Richtung durch Klick mit der linken Maustaste in eine der drei Navigationsansichten (2-4) durch Veränderung der Schieberegler unter den Navigationsansichten

Auswahlboxen Die Auswahlboxen werden im Baum-Widget (5) aufgelistet.(Abbildung B.2)

Abbildung B.2 Listenansicht der Auswahlboxen

Masterboxen werden hellblau, UND-Boxen grün und NICHT-Boxen rot hinterlegt und auch in dieser Farbe im Hauptanzeigefenster angezeigt. Bei Rechtsklick auf einen Eintrag öffnet sich ein Menu, über welches die Boxen umbenannt, sichtbar/unsichtbar geschaltet, aktiv/inaktiv geschaltet und gelöscht werden können. Durch Linksklick auf einen Eintrag im Baum-Widget oder Rechtsklick auf eine Box im Hauptanzeigefenster wird diese Box ausgewählt. Zur Verdeutlichung wird die Transparenz der Box verringert.

72

Folgende Tastatureingaben gelten für die ausgewählte Box: Cursor links/rechts Cursor oben/unten Bild auf/ab Strg+ Cursor links/rechts Strg+ Cursor oben/unten Bild auf/ab Shift + Cursor links/rechts Shift + Cursor oben/unten Shift + Bild auf/ab Shift+ Strg + Cursor links/rechts Shift+Strg + Cursor oben/unten Shift+ Strg + Bild auf/ab Entf Pos1 -

Verschiebung auf X-Achse um 1 Voxel Verschiebung auf Y-Achse um 1 Voxel Verschiebung auf Z-Achse um 1 Voxel Größenänderung auf X-Achse um 1 Voxel Größenänderung auf Y-Achse um 1 Voxel Größenänderung auf Z-Achse um 1 Voxel Verschiebung auf X-Achse um 5 Voxel Verschiebung auf Y-Achse um 5 Voxel Verschiebung auf Z-Achse um 5 Voxel Größenänderung auf X-Achse um 5 Voxel Größenänderung auf Y-Achse um 5 Voxel Größenänderung auf Z-Achse um 5 Voxel Löschen der Box klebt die Box an das Fadenkreuz, dies ermöglicht die genaue Positionierung über die Navigationssichten (2-4), ein nochmaliges Betätigen der Pos1-Taste löst diese Fixierung wieder

Die Boxen können auch im Hauptanzeigefenster mit der rechten Maustaste angefasst und bei gedrückter rechter Maustaste verschoben werden. Wird dabei die Strg-Taste gehalten, bewirkt ein ziehen der Maus eine Größenänderung der Box.

Liste der geladenen Datensätze

Abbildung B.3 List-Widget für Datensätze mit Buttons zur Veränderung der Position der Datensätze in der List und Schiebreglern für Schwellwert und Alphawert

73

In der Liste (Abbildung B.3) werden alle geladenen Texturen, triangulierte Oberflächen und Faserdatensätze angezeigt. Die Liste bietet außerdem Shortcuts zu über Menüs wählbare Funktionen. Über das Augen-Icon in der ersten Spalte kann das Rendering für alle Objekte ein- oder ausgeschaltet werden. Ein Doppelklick auf den Namen (2. Spalte) schaltet für Texturen die lineare Interpolation ein bzw. aus. Der nichtinterpolierte Modus wird durch einen * hinter dem Namen gekennzeichnet. Mit den Buttons „up“ und „down“ wird die Position des ausgewählten Objektes in der Liste verändert. Dies ist unter anderem wichtig für das Überblenden der Texturen, wo es auf die Reihenfolge ankommt.

Beschreibung der Toolbar-Buttons (Abbildung B.4)

Abbildung B.4 Toolbar für die am häufigsten benutzten Funktion

1

Einzelnen Datensatz oder abgespeicherte Szene laden

2

Szene speichern; speichert die Pfade und Dateinamen der gegenwärtig geladenen Datensätze, Position und Status aller Auswahlboxen, sowie die Stützpunkte für die Schnittfläche

3

Anzeige des axialen Schnittes ein- oder ausschalten

4

Anzeige des koronaren Schnittes ein- oder ausschalten

5

Anzeige des sagittalen Schnittes ein- oder ausschalten

6

Ein- oder Ausschalten des Clippings am Datensatzrand

7

Neue Auswahlbox am Schnittpunkt der Navigationsebenen erzeugen

8

Anzeige aller Auswahlboxen ein- oder ausschalten, die Boxen sind für die Auswahl der Fasern immer noch aktiv

9

Ausgewählte Box aktiv oder inaktiv schalten

10

Neue Schnittfläche an Position der sagittalen Navigationsebene erzeugen.

11

Modus zum manuellen Hinzufügen von Stützpunkten für die Schnittfläche einoder ausschalten

74

12

Randpunkte der Schnittfläche nach links verschieben

13

Randpunkte der Schnittfläche nach rechts verschieben

14

Ausgewähltes Objekt einfärben, dies kann ein selektiertes Faserbündel oder eine triangulierte Oberfläche sein

15

Beleuchtung der Fasern ein- oder ausschalten

16

Iso-Oberfläche für ausgewählten Skalardatensatz erzeugen

75

Literaturverzeichnis [1] Kindlmann G., Westin C.-F.: Diffusion tensor visualization with glyph packing. In Proceedings of IEEE Visualization’ 06 (Los Alamitos, CA, USA, October 2006), Gröller E., Pang A., Silva C. T., Stasko J., van Wijk J., (Eds.), IEEE Computer Society Press, pp. 1329–1335. [2] Blaas J., Botha C. P., Vos F. M., Post F. H.: Fast and reproducible fiber bundle selection in DTI visualization. In Proceedings of IEEE Visualization 2005 (Los Alamitos, CA, USA, 2005), Silva C. T., Gröller E., Rushmeier H., (Eds.), IEEE Computer Society, IEEE Computer Society Press, pp. 59–64. [3] Akers D., Sherbondy A., Mackenzie R., Dougherty R., Wandell B.: Exploration of the brain’s white matter pathways with dynamic queries. In Proceedings of IEEE Visualization’04 (Los Alamitos, CA, USA, 2004), Rushmeier H., Turk G., van Wijk J. J., (Eds.), IEEE Computer Society Press,pp. 377–384. [4] Wedeen V., Wang R., Schmahmann J., Benner T., Tseng W., Dai G., Pandya D., Hagmann P., d’Arcuil H., de Crespigny A.: Diffusion spectrum magnetic resonance imaging (dsi) tractography of crossing fibers. NeuroImage, 4 (Jul 2008), 1267–1277. [5] Moberts B., Vilanova A., van Wijk J. J.: Evaluation of fiber clustering methods for diffusion tensor imaging. In Proceedings of IEEE Visualization 2005 (Los Alamitos, CA, USA, 2005), Silva C. T., Gröller E., Rushmeier H., (Eds.), IEEE Computer Society, IEEE Computer Society Press, pp. 65–72. [6] Enders F., Sauber N., Merhof D., Hastreiter P., Nimsky C., Stamminger M.: Visualization of white matter tracts with wrapped streamlines. In Proceedings of IEEE Visualization 2005 (Los Alamitos, CA, USA, 2005), Silva C. T., Gröller E., Rushmeier H., (Eds.), IEEE Computer Society, IEEE Computer Society Press, pp. 51–58. [7] Zhukov L., Barr A. H.: Oriented tensor reconstruction: Tracing neural pathways from diffusion tensor MRI. In Proceedings of IEEE Visualization ’02 (Los Alamitos, CA, 2002), IEEE Computer Society, pp. 387–394. [8] Merhof D., Enders M. S. F., Nimsky C., Hastreiter P., Greiner G.: Hybrid visualization for white matter tracts using triangle strips and point sprites. IEEE Transactions on Visualization and Computer Graphics 12 (Oct. 2006), 1181–1188.

76

[9] Zhang S., Demiralp ´C Agatay., Laidlaw D. H.: Visualizing diffusion tensor MR images using streamtubes and streamsurfaces. IEEE Transactions on Visualization and Computer Graphics TVCG 9 (October-December 2003), 454–462. [10] Mallo O., Peickert R., Sigg C., Sadlo F.: Illuminated lines revisited. In Proceedings of IEEE Visualization 2005 (Los Alamitos, CA, USA, 2005), Silva C. T., Gröller E., Rushmeier H., (Eds.), IEEE Computer Society, IEEE Computer Society Press, pp. 19– 26. [11] Kondratieva P., Krüger J., Westermann R.: The application of GPU particle tracing to diffusion tensor field visualization. In Proceedings of IEEE Visualization 2005 (Los Alamitos, CA, USA, 2005), Silva C. T., Gröller E., Rushmeier H., (Eds.), IEEE Computer Society, IEEE Computer Society Press, pp. 73–78. [12] Park H.-J., Kubicki M., Westin C.-F., Talos I.-F., Brun A., Peiper S., Kikinis R., Jolesz F., Mccarley R., Shenton M.: Method for combining information from white matter tracking and gray matter parcellation. American Journal of Neuroradiology (2004), 1318–1324. [13] Schultz T., Sauber N., Anwander A., Theisel H., Seidel H.-P.: Virtual klingler dissection: Putting fibers into context. Computer Graphics Forum 27, 3 (2008), 1063– 1070. [14] Ludwig E., Klingler J.: Atlas cerebri humani: the inner structure of the brain demonstrated on the basis of macroscopical preparations. Boston : Little, Brown, c1956., 1956. [15] Laramee R. S., Jobard B., Hauser H.: Image space based visualization of unsteady flow on surfaces. Proceedings of IEEE Visualization 2003 (2003). [16] Hotz I., Feng L., Hagen H., Hamann B., Joy K. I.: Tensor field visualization using a metric interpretation. In Visualization and Processing of Tensor Fields (2006), Weickert J., Hagen H., (Eds.), Springer–Verlag Berlin Heidelberg, pp. 269–281. [17] Stejskal E.O., Tanner, J.E.: Spin Diffusion Measurements - Spin echoes in the presence of a time-dependent field gradient. In: Journal of Chemical Physics. Melville 42.1965, 288−292. [18] C. Westin, S. Maier, H. Mamata, A. Nabavi, F. Jolesz, and R. Kikinis. Processing and visualization for diffusion tensor mri. Medical Image Analysis, 6(2):93.108, 2002. [19] Mori S. and van Zijl P.C.: Fiber tracking: principles and strategies - a technical review. NMR Biomed, 15:468.480, 2002.

77

[20] Mori S., Crain B.J., Chacko V.P., and van Zijl P.C.: Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann Neurol, 45(2):265.269, 1999. [21] Mori S., Wakana S., Nagae-Poetscher L.M., van Zijl P.C.: MRI Atlas of the Human white Matter, Baltimore, MD, USA, 2005. [22] Anwander, A.; Tittgemeyer, M.; von Cramon, D.Y. ; Friederici, A.D.; Knösche, T.R.: Connectivity-Based Parcellation of Broca's Area,Cerebral Cortex, Volume 17, Number 4, (2007), pp. 816-825. [23] Kreher, B.W.: Detektion von Hirnnervenfasern auf der Basis diffusionsgewichteten Magnetresonanzdaten, Dissertationsschrift, Freiburg, 2007.

von

[24] Lorensen W. E., Cline H. E., „Marching Cubes: A High Resolution 3D Surface Construction Algorithm“, ACM Computer Graphics Vol. 21 No. 4 (SIGGRAPH 1987 Proceedings). [25] Cabral B., Leedom L. C.: Imaging Vector Fields Using Line Integral Convolution. In: Proceedings of ACM SIGGRAPH '93, Aug 2-6, Anaheim, California, pp. 263-270, 1993. [26] Loop C.: Smooth Subdivision Surfaces Based on Triangles. Master’s thesis, University of Utah, Salt Lake City, UT, USA, 1987. [27] D. Schreiner. OpenGL Reference Manual: The Official Reference Document to OpenGL, Version 1.2. Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1999. [28] Sherbondy A., Akers D., Mackenzie R., Dougherty R., Wandell B.: Exploring Connectivity of the Brain’s White Matter with Dynamic Queries, IEEE Transactions on Visualization and Computer Graphics (2005). [29] Toussaint N., Souplet J.-C., Fillard P.: MedINRIA: DT-MRI processing and visualization software. Proceedings of MICCAI’07, Workshop on Interaction in medical image analysis and visualization (2007). [30] nVidia Corporation: nVidia Whitepaper Using Vertex Buffer Objects (VBOs), Santa Clara, California, 2003. [31] Hlawitschka M., Anwander A., Schurade R.: FiberNavigator: Dissect the embedding Anatomy, in Review for EuroVis 2009, Berlin. [32] Anwander A, Schurade R., Hlawitschka M.: White Matter Imaging with Virtual Klingler Dissection, submitted to Human Brain Mapping 2009, San Francisco (CA), 2009.

78

Abbildungsverzeichnis 1.1 Abbildung aus Atlas Cerebri Humani ............................................................................... 6

2.1 Karten der Diffusionskonstanten .................................................................................. 11 2.2 Parameter zur Beschreibung des Diffusionsellipsoiden ................................................. 12 2.3 Karte der fraktionellen Anisotropie .............................................................................. 13 2.4 Farbige Darstellung der größten Eigenvektoren ............................................................ 14 2.5 T1-Graustufenbild überlagert von Wahrscheinlichkeitskarte ........................................ 17

3.1 Dreidimensionale Darstellung eines kd-Baumes ........................................................... 19 3.2 Veranschaulichung eines Würfels zwischen zwei Schichten ........................................... 20 3.3 Nummerierung der Ecken und Kanten des Würfels ....................................................... 21 3.4 Darstellung der triangulierten Würfel ........................................................................... 22 3.5 Strömungsvisualisierung um bluntfin ............................................................................ 26 3.6 LIC auf Iso-Oberfläche einer Wahrscheinlichkeitskarte .................................................. 27 3.7 Aufteilung der Dreiecke in vier neue Dreiecke .............................................................. 28 3.8 Regel für innere Punkte ................................................................................................ 28 3.9 Regel für innere gerade reguläre Punkte ...................................................................... 29 3.10 Regel für innere gerade irreguläre Punkte ................................................................... 29 3.11 Regel für ungerade Randpunkte ................................................................................. 30 3.12 Regel für gerade Randpunkte ..................................................................................... 30 3.13 Exakte Berechnung von Randpunkten ......................................................................... 31 3.14 Exakte Berechnung von Randpunkten ......................................................................... 31

4.1 Darstellung verschiedene Datensätze .......................................................................... 33 4.2 Eigenvektordatensatz farbkodiert und FA-gewichtet ..................................................... 34 4.3 Wahrscheinlichkeitskarte interpoliert und nicht interpoliert ......................................... 35

79

4.4 Faserselektion mit einer Auswahlbox ........................................................................... 37 4.5 Faserauswahl mit Master- und UND-Box....................................................................... 38 4.6 Faserselektion mit mehreren Boxen ............................................................................. 38 4.7 Kombination mehrerer Auswahlboxen zur Auswahl unterschiedlicher Faserbündel ....... 39 4.8 Die besuchten Knoten des kd-Baum für die Auswahl des cingulum ............................... 40 4.9 Iso-Oberfläche auf Wahrscheinlichkeitskarte ................................................................ 41 4.10 Faserselektion durch fMRI Daten ............................................................................... 42 4.11 globale Farbgebung .................................................................................................... 43 4.12 lokale Farbgebung ...................................................................................................... 43 4.13 manuelle Färbung ...................................................................................................... 43 4.14 Färbung mit FA........................................................................................................... 43 4.15 unbeleuchtete Fasern ................................................................................................. 44 4.16 beleuchtete Fasern..................................................................................................... 44 4.17 kompletter Faser-Datensatz als Röhren gerendert....................................................... 45 4.18 Konstruktion des Quadstrips ...................................................................................... 46 4.19 Iso-Oberfläche auf Wahrscheinlichkeitskarte mit Eigenvektorrichtung gefärbt ............ 48 4.20 Artefakte nach Erstellung der Iso-Oberfläche ............................................................. 49 4.21 Iso-Oberfläche auf mit Medianfilter geglättetem Datensatz ........................................ 50 4.22 Artefaktentfernung nach Erstellung der Oberfläche .................................................... 51 4.23 Clipping von Iso-Oberfläche an Schnittfläche .............................................................. 53

5.1 Endergebnis der virtuellen Klingler-Disektion................................................................ 54 5.2 Faserauswahl ............................................................................................................... 55 5.3 Positionierung der Navigationsebene ........................................................................... 56 5.4 Erzeugte Schnittfläche .................................................................................................. 56 5.5 Editiermodus für Schnittfläche ..................................................................................... 56 5.6 Verformte Schnittfläche, texturiert mit T1 und LIC ....................................................... 57 5.7 Vergleich des Verlaufs des Fasciculus Uncinatus in 6 Versuchspersonen ........................ 58 5.8 Visualisierung der Unterkomponenten des SLF ............................................................. 59 5.9 Farbkodierte FA-Werte entlang der Fasern des internen Kapselsystems ........................ 60 5.10 Faserbündel aus 7 Tesla DTI Bildern ............................................................................ 60 5.11 Faserauswahl anhand von funktionell aktivierten Regionen ........................................ 61

80

5.12 Akustisch stimulierte FMRI-Aktivierung sprachrelevanter Hirnareale ........................... 62 5.13 Traktogramme von drei Startregionen in inferioren Parietallappen ............................. 62 5.14 Virtuelle und post-mortem Dissektion der sub-insularen weißen Substanz .................. 63 5.15 Klingler Dissektion und virtueller Schnitt des sprachrelevanten Fasersystems .............. 64

B.1 Ansicht der Benutzeroberfläche ................................................................................... 71 B.2 Listenansicht der Auswahlboxen ................................................................................. 72 B.3 List-Widget für Datensätze ........................................................................................... 73 B.4 Toolbar ........................................................................................................................ 74

81

Abkürzungsverzeichnis AF

-

Fasciculus Arcuatus

DTI

-

Diffusion Tensor Imaging

FA

-

Fraktionelle Anisotropie

FACT

-

Fiber Assignment by Continuous Tracking

GUI

-

Graphical User Interface

fMRT

-

funktionelle Magnet- Resonanz-Tomographie

IFO

-

Fasciculus occipitofrontalis inferior

LIC

-

Line Integral Convolution

MRT

-

Magnet-Resonanz-Tomographie

ROI

-

Region of Interest

SLF

-

Fasciculus Longitudinales Superior

UNC

-

Fasciculus Uncinatus

VBO

-

Vertex Buffer Objects

82

Danksagung An dieser Stelle möchte ich mich bei all jenen bedanken, die mich während und bei der Arbeit an dieser Diplomarbeit unterstützt haben. Insbesondere gilt mein Dank:

Dr. Alfred Anwander für die sehr unkomplizierte Art und Weise der Betreuung meiner Arbeit sowie die endlosen Stunden der Diskussion und Erörterung. Dr. Mario Hlawitschka für Design- und Programmiertipps und Aufmunterungen der Art: „das ist doch in 2-3 Stunden implementiert“. Bei Benutzung der nächstgrößeren Zeiteinheit stimmte dies sogar meist. Tobias Göbel für das unermüdliche Betatesting im Rahmen seiner Doktorarbeit. Viele seiner Anregungen zur Funktionalität und Benutzerinterface sind in das Programm eingeflossen. Prof. Gerik Scheuerman und Dr. Alexander Wiebel stellvertretend für das gesamte FAnToM-Team für die Erlaubnis Code-Fragmente aus FAnToM benutzen zu dürfen. Christian Heine für den Tipp zu nth_element. Dieser erst machte den kd-Baum zu dem eleganten und schnellen Konstrukt, der er jetzt ist.

83

Erklärung

Ich versichere, dass ich die vorliegende Arbeit selbständig und nur unter Verwendung der angegebenen Quellen und Hilfsmittel angefertigt habe, insbesondere sind wörtliche oder sinngemäße Zitate als solche gekennzeichnet. Mir ist bekannt, dass Zuwiderhandlung auch nachträglich zur Aberkennung des Abschlusses führen kann.

Leipzig, den 13. Januar 2008 (Ralph Schurade)

84