GPU-unterstütztes Radiosity. Entwicklung, Probleme, Perspektiven

Inhalt

1 Einleitung
1.1 Zielsetzung
1.2 Gliederung der Arbeit

2 Grundlagen des Radiosity-Verfahrens
2.1 Zugrunde gelegte Theorie
2.2 Basisalgorithmus
2.3 Bestimmung der Formfaktoren
2.3.1 Halbkugelverfahren
2.3.2 Halbwürfelverfahren
2.3.3 Approximation mittels Kreisscheiben
2.4 Verbesserung der Leistungsfähigkeit
2.4.1 Hierarchisches Radiosity
2.4.2 Schrittweise Verfeinerung

3 GPGPU
3.1 Grundlagen der GPU-Programmierung
3.1.1 Grafikpipeline
3.1.2 Programmierbare Hardware
3.1.3 GPU-Programmiermodell
3.1.4 Ablaufsteuerung von GPU-Programmen
3.2 Programmiersysteme
3.2.1 Shader-Hochsprachen
3.2.2 GPGPU-spezifische Sprachen und Bibliotheken
3.2.3 Debugging-Tools
3.3 Einführung in die RapidMind-Entwicklungsplattform .
3.3.1 Value-Typ
3.3.2 Array-Typ
3.3.3 Program-Typ
3.3.4 Teilweise Berechnung, Programmalgebra und kollektive Operationen
3.3.5 Parallelisierung eines Beispiels
3.4 Zusammenfassung

4 Mathematische Hilfsmittel
4.1 Iteratives Lösen linearer Gleichungssysteme
4.1.1 Beschreibung klassischer Iterationsverfahren
4.1.2 Jacobi-Verfahren
4.1.3 Gauß-Seidel-Verfahren
4.1.4 Konvergenzkriterium
4.1.5 Parallele Realisierung des Jacobi-Verfahrens
4.2 Ein Schnitttest zwischen Strecke und Dreieck
4.2.1 Adaption des Strahl-Dreieck-Schnitttests von Möller
4.2.2 Vorschläge zur Verbesserung der Leistungsfähigkeit
4.3 Octrees

5 Entwurf und Implementierung
5.1 Entwurf
5.2 Besonderheiten der Implementierung auf der GPU
5.2.1 Verwendete Entwicklungsumgebung
5.2.2 Beschränkte Anzahl an Schleifeniterationen
5.2.3 Beschränkte Größe von Texturen, Arrays und Eingabedaten
5.2.4 Partitionierung des Jacobi-Verfahrens
5.3 Szenenbeschreibung
5.3.1 MDL
5.3.2 RML
5.4 Grafische Benutzeroberfläche

6 Vergleich von CPU- und GPU-Implementierung
6.1 Ergebnisse der Zeitmessungen
6.1.1 Laufzeitvergleich des Radiosity-Verfahrens
6.1.2 Bewertung der Vorschläge zur Beschleunigung
6.1.3 Vergleich von Kreisscheibenapproximation und Halbwürfelverfahren
6.1.4 Zusammenfassung
6.2 Probleme der GPU-Programmierung

7 Resümee und Ausblick
A Gerenderte Szenen
B Spezifikationsschemata und Beispiele zu den Szenenbeschreibungsformaten
B.1 MDL-Schema
B.2 RML-Schema
B.3 Beispiele zu den Szenenbeschreibungsformaten
C Beispiel eines generierten GLSL-Quelltextes

Quellen

Abbildungen

1.1 Das Radiosity-Verfahren im Einsatz

2.1 Reflexion nach Lambertschem Kosinusgesetz

2.2 Der Radiosity-Basisalgorithmus

2.3 Verschmelzen benachbarter Flächenelemente

2.4 Energieaustausch zwischen differentiellen Flächenelementen

2.5 Halbkugelverfahren

2.6 Halbwürfelverfahren

2.7 Herleitung des Delta-Formfaktors

2.8 Formfaktorapproximation mittels Kreisscheiben

2.9 Adaptive Unterteilung beim hierarchischen Radiosity

3.1 Vergleich des Wachstums der Leistungsfähigkeit von CPUs und GPUs

3.2 Vereinfachter Aufbau einer programmierbaren Grafikpipeline

3.3 Berechnungsmuster im Vergleich

4.1 Vorberechnung zum Strecke-Dreieck-Schnitttest mithilfe einer Kugel

4.2 Vorberechnung zum Strecke-Dreieck-Schnitttest mithilfe von drei Kugeln . .

4.3 Octree-Darstellung

5.1 Ausschnitt aus dem UML-Klassendiagramm von Radioso

5.2 Die grafische Benutzeroberfläche von Radioso

6.1 Vergleich der Sichtbarkeitenbestimmung mit Kugelvorberechnungen (CPU)

6.2 Vergleich der Sichtbarkeitenbestimmung mit Octree-Unterstützung (CPU) .

6.3 Laufzeitvergleich der Verfahren zur Formfaktorbestimmung

6.4 Qualitativer Vergleich der Verfahren zur Formfaktorbestimmung

A.1 Zu den Laufzeitvergleichen herangezogene Testszenen

Tabellen

6.1 Laufzeitvergleich des Radiosity-Verfahrens

6.2 Laufzeitvergleich der Sichtbarkeitenbestimmung auf der GPU

Quelltexte

4.1 Pseudocode des angepassten Algorithmus von Möller

5.1 Paralleles Jacobi-Verfahren

5.2 Partitioniertes paralleles Jacobi-Verfahren

B.1 Beipiel einer Szenenbeschreibung in MDL

B.2 Beipiel einer Szenenbeschreibung in RML

Abkürzungen

Abbildung in dieser Leseprobe nicht enthalten

1 Einleitung

Die Verwendung computergenerierter Bilder hat in den letzten Jahren in immer mehr Be- reichen des täglichen Lebens Einzug gehalten. Eine sehr bedeutende Rolle spielen diese vor allem in der Film- und Werbeindustrie, wo durch synthetisch erzeugte Objekte im Zusam- menspiel mit realen Objekten für den Betrachter die Grenzen zwischen Realität und Illusi- on zu verschmelzen scheinen. Aber auch im Fahrzeugbau und in Architekturbüros sind mit dem Computer erzeugte Modelle heute nicht mehr wegzudenken. Sie erreichen eine Form der Veranschaulichung, die der von Modellen aus Gips, Holz, Pappe oder Kunststoff weit überlegen ist. Weitere Anwendungsgebiete finden sich in der Medizin (z. B. bei der Erstel- lung dreidimensionaler Modelle von Organen mithilfe der Computertomographie), bei Fahr- und Flugsimulatoren und natürlich in Computerspielen. Daneben gewinnen die Verfahren der Computergrafik auch in wissenschaftlichen Zweigen zunehmend an Bedeutung. So lassen sich mit ihrer Hilfe beispielsweise chemische und biologische Strukturen visualisieren, physikalische Modelle simulieren oder geographische und geologische 3D-Karten erzeugen.

Ein wesentliches Ziel ist dabei das möglichst realistische Aussehen der dargestellten Ob- jekte. Um dies zu erreichen, kommen häufig komplexe mathematische Verfahren zum Ein- satz. Eine besondere Bedeutung kommt dabei dem Licht und dessen Ausbreitung in drei- dimensionalen Szenen zu. Hierbei ist zwischen lokalen und globalen Beleuchtungsverfahren zu unterscheiden. Ein globales Beleuchtungsverfahren, das sich durch seine sehr gute diffuse Ausleuchtung dreidimensionaler Szenen auszeichnet, ist das Radiosity-Verfahren. Dieses wird häufig in Kombination mit anderen Verfahren wie dem Raytracing oder Photon Mapping verwendet. Abbildung 1.1 zeigt einige mithilfe des Radiosity-Verfahrens generierte Bilder.

Ein großer Nachteil des Radiosity-Verfahrens ist sein extrem hoher Aufwand, der dazu führt, dass ein Einsatz des Verfahrens bei sehr komplexen Szenen nur selten in Frage kommt. Trotz stetig steigender Rechenleistung ist die Anwendung des Verfahrens auf Standard-PCs bislang unzumutbar. Bereits lange Zeit vor der Entwicklung programmierbarer Grafikpipeli-

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 1.1: Das Radiosity-Verfahren im Einsatz. Quellen: (a) und (b) vom Cornell University Program of Computer Graphics [Car], (c) von Jacobson [Jac].

1 Einleitung

neelemente wurde deshalb versucht, das Radiosity-Verfahren mithilfe der Grafikhardware zu beschleunigen, vgl. [CHH03]. Jedoch verblieb ein Großteil der Berechnung dabei stets auf der CPU. Dies änderte sich mit der Einführung programmierbarer Grafikprozessoren (GPUs). Im Jahr 2003 veröffentlichten Carr et al. in ihrer Arbeit [CHH03] erste Algorithmen für eine vollständige Umsetzung des Radiosity-Verfahrens auf der Grafikkarte. In den darauffolgenden Jahren präsentierten sowohl Coombe et al. [CHL04, CH05] als auch Borsdorf et al. [BDS05] ihre GPU-Entwicklungen des Radiosity-Verfahrens mit beachtlichen Leistungssteigerungen gegenüber herkömmlichen CPU-Implementierungen. Darüber hinaus gelang Berenguier [Ber] mithilfe programmierbarer Grafikhardware die Realisierung einer, zumindest für kleine Sze- nen echtzeitfähigen Radiosity-Anwendung.

Einen aktuellen Ansatz der GPU-Programmierung verfolgt der Bereich des sogenannten GPGPU (General-Purpose computation on the Graphics Processing Unit). Dabei soll die programmierbare Grafikhardware auf möglichst einfachem Wege auch allgemeinen Berechnungen zugänglich gemacht werden. Speziell hierfür entwickelte Werkzeuge erleichtern die Programmierung und verbergen grafikspezifische Details. Es bietet sich an, den GPGPUAnsatz auch für das Radiosity-Verfahren zu nutzen, da dieses keine reine Grafikanwendung darstellt, sondern größtenteils aus allgemeinen Algorithmen besteht.

1.1 Zielsetzung

Im Rahmen dieser Diplomarbeit soll untersucht werden, inwieweit sich die Leistungsfähigkeit des Radiosity-Verfahrens mithilfe programmierbarer Grafikhardware gegenüber herkömmlichen Implementierungsmethoden verbessern lässt. Dabei soll die Umsetzung des Verfahrens mit den Möglichkeiten des GPGPU-Ansatzes eruiert werden.

1.2 Gliederung der Arbeit

Zur Einführung in die Thematik werden zunächst die Grundlagen des Radiosity-Verfahrens näher erläutert. Anschließend gewährt das dritte Kapitel einen Einblick in die GPU-Program- mierung aus der Sicht des GPGPU-Ansatzes. Es werden verschiedene Werkzeuge der GPU- Programmierung vorgestellt und eine Einführung in die für die Implementierung verwende- te RapidMind-Entwicklungsplattform gegeben. Beschreibungen zu den für die Umsetzung benötigten mathematischen Techniken finden sich im vierten Kapitel. Das darauffolgende Kapitel widmet sich den Besonderheiten des im Rahmen dieser Diplomarbeit entwickelten Programms Radioso. Danach werden in einem weiteren Kapitel die CPU- und die GPU- Implementierung miteinander verglichen. Diese Gegenüberstellung erstreckt sich über die Auswertung der in diesem Zusammenhang durchgeführten Zeitmessungen und aufgetretenen Probleme. Schließlich zieht das letzte Kapitel ein Resümee über die vorliegende Arbeit und gibt einen Ausblick auf zukünftige Entwicklungen in diesem Gebiet.

2 Grundlagen des Radiosity-Verfahrens

Die Ursprünge des Radiosity-Verfahrens gehen zurück auf eine Forschungsarbeit von Goral et al. [GTGB84]. Diese entwickelten 1984 an der Cornell-Universität eine auf dem Prinzip der Wärmeübertragung von Strahlen basierende Methode, die heute als klassisches Radiosity- Verfahren bekannt ist. Dabei handelt es sich um ein globales Beleuchtungsverfahren für drei- dimensionale Szenen, welches sich vor allem durch seine sehr guten Eigenschaften bezüglich diffuser Reflexion auszeichnet und sich deshalb besonders gut zur Darstellung von Innen- raumszenen eignet. Im Gegensatz zu lokalen Beleuchtungsverfahren, die nur einzelne Punkte der Szene betrachten, berücksichtigen globale Verfahren auch die Wechselwirkungen der in der Szene befindlichen Objekte und erreichen damit ein weitaus realistischeres Erscheinungsbild. Als bekanntestes globales Beleuchtungsverfahren ist neben dem Radiosity das Raytracing zu nennen, bei dem mittels Strahlverfolgung der Weg von Lichtstrahlen durch die Szene simuliert wird.

Beide Verfahren unterscheiden sich durch verschiedene Stärken und Schwächen. Während das Raytracing auch die Darstellung von Spiegelungen und Brechung ermöglicht, eignet sich das Radiosity-Verfahren besser zur Visualisierung von diffuser Reflexion und weichen Schat- ten. Ein weiterer Vorteil des Radiosity-Verfahrens ist seine Unabhängigkeit vom Betrachter- standpunkt. Dadurch ist es möglich, sich durch eine Szene zu bewegen, ohne die gesamte Berechnung wiederholen zu müssen. Als größter Nachteil beider Verfahren gilt ihr hoher Berechnungsaufwand, der beim Radiosity sogar dazu führt, dass dieses ohne spezielle Be- schleunigungstechniken nicht nutzbar ist. Darüber hinaus handelt es sich beim Radiosity um ein sehr speicherintensives Verfahren.

Den Ausgangspunkt dieses Kapitels bilden einige theoretische Aspekte zum Radiosity-Ver- fahren. Darauf aufbauend folgt ein Überblick über den Radiosity-Basisalgorithmus in seiner Gesamtheit, bevor anschließend verschiedene Techniken zur Berechnung von Formfaktoren, dem aufwendigsten Teil des Verfahrens, erläutert werden. In einem letzten Abschnitt wer- den zudem die beiden bedeutendsten Ansätze zur Verbesserung der Leistungsfähigkeit des Radiosity vorgestellt.

2.1 Zugrunde gelegte Theorie

Die Grundlage des Radiosity-Verfahrens, vgl. [Rau93, Wat02, CW93, FDFH90], bildet der Energieerhaltungssatz , welcher besagt, dass die Gesamtenergie eines abgeschlossenen Systems konstant ist. Das abgeschlossene System wird dabei durch den Raum der darzustellenden Sze- ne repräsentiert, als Energie wird ausschließlich die Energie des Lichtes betrachtet. Diese wird von den Lichtquellen in Form von Photonen abgestrahlt. Beim Radiosity sind, im Gegensatz zu anderen Beleuchtungsverfahren wie beispielsweise dem Raytracing, Lichtquellen norma- le Objekte mit Ausdehnung, welche sich gegen die anderen Objekte der Szene nur dadurch abgrenzen, dass sie Licht aussenden.

Da nach obigem Erhaltungssatz die Gesamtenergie einer Szene konstant bleibt, muss die abgestrahlte Energie, um nicht verloren zu gehen, beim Auftreffen auf ein anderes Objekt reflektiert werden. Hierzu gehorchen die Oberflächen der Objekte dem Lambertschen Kosinusgesetz, sodass man perfekte diffuse Reflexion erhält. Nach Lambert ist die Strahlungsdichte der von einem Flächenelement ausgehenden Strahlung proportional zum Kosinus des Winkels φ zwischen der Flächennormale und der Richtung dieser Strahlung (vgl. [PL05, Wag93, Els93]). Die diffus abgegebene Strahlung ist damit senkrecht zur Fläche maximal und verringert sich mit größer werdendem Winkel φ (siehe Abbildung 2.1).

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.1: Reflexion nach Lambertschem Kosinusgesetz. Unabhängig vom Winkel des einfallenden Lichtstrahls ist die Reflexion an einer Lambertschen Fläche senkrecht zu dieser maximal und verringert sich mit größer werdendem Winkel zwischen Flächennormale n und ausgehendem Strahl. Quelle: Nach [Hau96].

Mithilfe dieser Annahmen kann die allgemeine Rendergleichung der Computergrafik in die Radiosity-Gleichung überführt werden. Bei der allgemeinen Rendergleichung handelt es sich um eine Integralgleichung zur Beschreibung der Energieerhaltung bei der Ausbreitung von Lichtstrahlen in einer Szene. Laut [FDFH90] können alle Algorithmen zur globalen Beleuch- tung auf diese Gleichung zurückgeführt werden. Da die resultierende Radiosity-Gleichung in ihrer stetigen Form jedoch nicht direkt zur Berechnung der gesuchten Radiosity-Werte genutzt werden kann, führt man eine Diskretisierung der Oberfläche der Szene durch. Eine Herleitung der diskreten Form der Radiosity-Gleichung erfolgt im nächsten Abschnitt.

2.2 Basisalgorithmus

Grundsätzlich lässt sich das Radiosity-Verfahren in vier aufeinander folgende Schritte glie- dern (siehe Abbildung 2.2). Zunächst werden die Oberflächen der in der Szene befindlichen Objekte in kleine Flächenelemente unterteilt. Für jedes dieser als ebenes Polygon behandelten Flächenelemente soll im weiteren Verlauf ein konstanter Strahlungswert berechnet werden, der alle anderen Flächenelemente der Szene berücksichtigt. Um dies zu gewährleisten, wer- den in einem zweiten Schritt die sogenannten Formfaktoren ermittelt, welche zur Lösung des Radiosity-Gleichungssystems im sich anschließenden dritten Schritt benötigt werden. Als Ergebnis des dritten Schrittes erhält man die gesuchten Strahlungswerte, mithilfe derer im abschließenden vierten Schritt die Szene gerendert wird. Bei einem Wechsel der Position des Betrachters muss lediglich der letzte Schritt wiederholt werden, da sich dabei die Strahlungs- werte nicht ändern.

Der Strahlungswert B eines Flächenelements, in der Physik auch als Strahlungsdichte be- zeichnet, ist definiert als die von diesem pro Zeit- und Flächeneinheit abgestrahlte Energie und wird daher auch in W/m[2] angegeben. Dem Energieerhaltungssatz zufolge setzt sich die von einem Flächenelement i abgestrahlte Energie aus emittierter und aus reflektierter Energie zusammen, wobei nur Lichtquellen einen von Null verschiedenen emittierten Anteil haben.

2.2 Basisalgorithmus

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.1: Reflexion nach Lambertschem Kosinusgesetz. Unabhängig vom Winkel des einfallenden Lichtstrahls ist die Reflexion an einer Lambertschen Fläche senkrecht zu dieser maximal und verringert sich mit größer werdendem Winkel zwischen Flächennormale n und ausgehendem Strahl. Quelle: Nach [Hau96].

Abbildung 2.2: Der Radiosity-Basisalgorithmus lässt sich in vier aufeinander folgende Schritte (grün) gliedern. Da die berechneten Strahlungswerte von der Be- trachterposition unabhängig sind, muss bei einem Wechsel nur das Rendern wiederholt werden.

Der reflektierte Teil der Energie kann von allen anderen Flächenelementen der Szene abstammen. Er steigt proportional zur Größe des diffusen Reflexionskoeffizienten ρi, welcher die reflektierende Eigenschaft der Oberfläche von i beschreibt. [Rau93]

Bezeichnen nun Bi und Bj die von den Flächenelementen i und j abgegebenen Strahlungen in Energie pro Zeit- und Flächeneinheit, Ai und Aj ihre Flächeninhalte sowie Ei die von i emittierte Strahlung ebenfalls in Energie pro Zeit- und Flächeneinheit, so erhält man für die gesamte von Flächenelement i abgegebene Strahlung

(2.1)

Abbildung in dieser Leseprobe nicht enthalten

Fji benennt dabei den Formfaktor, der angibt, welcher Anteil der vom gesamten Flächenelement j abgestrahlten Energie direkt auf das gesamte Flächenelement i trifft. Er stellt genau wie der diffuse Reflexionskoeffizient ρ eine dimensionslose Größe dar. Eine genaue Definition sowie die Bestimmung der Formfaktoren folgen in Abschnitt 2.3. Damit entspricht FjiBj Aj der Strahlungsleistung, die vom gesamten Element j auf das gesamte Element i ausgeübt wird. Die Summe dieser Strahlungsleistungen über alle n Flächenelemente der Szene bildet dann die gesamte auf i einfallende Strahlung. [Rau93, Wat02]

Um die Gleichung (2.1) weiter zu vereinfachen, formt man diese zunächst mittels Division durch Ai zu die Strahlungsleistung, die das gesamte Flächenelement j verlässt und eine Flächeneinheit von i erreicht. Des Weiteren lässt sich auf Grundlage der Berechnung der Formfaktoren die allgemeine Reziprozitätsbeziehung nachweisen (vgl. Abschnitt 2.3), mithilfe derer sich (2.2) leicht zu ∑

(2.3)

Abbildung in dieser Leseprobe nicht enthalten

(2.4)

Abbildung in dieser Leseprobe nicht enthalten

2 Grundlagen des Radiosity-Verfahrens

(a)

Abbildung in dieser Leseprobe nicht enthalten

(b)

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.3: Verschmelzen benachbarter Flächenelemente: In 1273 Flächenelemente unterteilte Cornell-Box-Szene (a) ohne Verschmelzen und (b) mit Verschmel- zen der Flächenelemente gerendert.

umstellen lässt. Letztere wird als Radiosity-Gleichung bezeichnet und ergibt in Matrixschreib- weise

(2.5)

Abbildung in dieser Leseprobe nicht enthalten

En

Sind die Formfaktoren Fij für alle i, j = 1, . . . , n bekannt, so lassen sich bei gegebenen Reflexionskoeffizienten ρi und gegebenen emittierten Energien Ei durch Lösen dieses linearen Gleichungssystems die gesuchten Strahlungswerte B1, . . . , Bn ermitteln. Üblicherweise muss dies für verschiedene Wellenlängen durchgeführt werden, da Reflexionskoeffizient und emittierte Energie von dieser abhängen. Bei Verwendung von RGB-Farbwerten sind beispielsweise drei Berechnungen notwendig. [Rau93]

Anschließend findet in einem letzten Schritt das Rendern statt. Um sichtbare Übergänge zwischen benachbarten Flächenelementen zu vermeiden, sollte man bei der Ermittlung der Intensitätswerte darauf achten, dass zu einer Oberfläche gehörige Elemente miteinander ver- schmolzen werden. In [CG85] wird hierzu ein Interpolationsverfahren ähnlich der Gouraud- Schattierung[1] vorgeschlagen. Eine Beispielszene, welche einmal mit und einmal ohne Ver- schmelzen der Flächenelemente gerendert wurde, ist in Abbildung 2.3 zu sehen.

2.3 Bestimmung der Formfaktoren

Formfaktoren repräsentieren die geometrische Beziehung zwischen je zwei Flächenelementen. Sie sind die zentrale Komponente des Radiosity-Verfahrens und zugleich elementar, will man

die Ausbreitung des Lichtes in einer Szene verstehen. Der Formfaktor FAiAj (kurz: Fij ) be- schreibt den Energieaustausch von Flächenelement Ai zu Flächenelement Aj . Er gibt an,

[1]Gouraud-Schattierung oder auch Gouraud Shading ist ein Verfahren der Computergrafik zum Schattieren von Polygonflächen, benannt nach seinem Entwickler Henri Gouraud.

2.3 Bestimmung der Formfaktoren

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.4: Energieaustausch zwischen den differentiellen Flächenelementen dAi und dAj . Quelle: Nach [CG85].

welcher Anteil der von Ai abgegebenen Energie Aj direkt erreicht. Seiner Berechnung liegt ausschließlich die Geometrie der Szene (Abstand und Lage der Flächen zueinander) zugrunde. Bei statischen Szenen braucht jeder Formfaktor deshalb nur einmal berechnet werden. Cohen und Wallace [CW93], Rauber [Rau93], Watt [Wat02] sowie Foley et al. [FDFH90] erklären den Begriff des Formfaktors auf ähnliche Weise.

Darüber hinaus können laut Watt [Wat02] und Rauber [Rau93] zwei weitere Eigenschaften der Formfaktoren direkt aus ihrer Definition abgeleitet werden. Zum einen gilt Fii = 0, da sich ein Flächenpatch bzw. eine konvexe Oberfläche nie selbst direkt bestrahlt, zum anderen folgt aus dem Energieerhaltungssatz, dass die Summe einer beliebigen Zeile von Formfaktoren gleich Eins ist.

Die Beschreibung der geometrischen Beziehung zweier Flächenelemente durch den Formfaktor wurde aus der technischen Physik übernommen, wo man mit ihm dieÜbertragung von Wärmestrahlung erklärt. Um für den Formfaktor der Flächenelemente Ai und Aj eine Berechnungsgleichung zu finden, betrachtet man zunächst die differentiellen Teile dAi und dAj dieser Elemente. Siegel et al. leiten für den Formfaktor zwischen zwei differentiellen Flächenelementen in [SHL91] folgenden Zusammenhang her[2]:

(2.6)

Abbildung in dieser Leseprobe nicht enthalten

Dabei sind φi und φj die jeweiligen Winkel zwischen Verbindungsgerade und Flächennormale und r der Abstand von dAi nach dAj (siehe Abbildung 2.4). Durch Integration über das Flächenelement Aj erhält man den Teil der von dAi abgestrahlten Energie, der das gesamte

(2.7)

Abbildung in dieser Leseprobe nicht enthalten

[2]Dieser steht nicht im Einklang mit der Formulierung in [Rau93], wo im Nenner zusätzlich der Faktor 2 hinzugefügt ist. Als Begründung dient dort die Oberfläche einer Halbkugel auf Grundlage der Feststellung, dass die Intensität der ausgesandten Strahlung nach jeder Richtung des Raumes gleich groß ist (vgl. hierzu auch [Wag93]). Meiner Meinung nach ist dieser Ansatz hier nicht zutreffend, da die Formfaktoren zur Ermittlung von Strahlungsdichten und nicht von Intensitäten verwendet werden. Die Strahlungsdichte jedoch gehorcht dem Lambertschen Kosinusgesetz, was keine Deutung mittels Halbkugeloberfläche zulässt. Damit einher gehen u. a. auch die Interpretationen von [FDFH90] und [Wat02], bei denen der Faktor 2 ebenso nicht vorkommt.

Zur Berechnung des Formfaktors zwischen den endlichen Flächenelementen Ai und Aj integriert man nun zusätzlich über Ai:

(2.8)

Abbildung in dieser Leseprobe nicht enthalten

Schließlich fügt man noch eine binäre Funktion { Vis(dAi, dAj ) = 0, zwischen dAi und dAj liegt ein verdeckendes Objekt 1, sonst ein, um eine mögliche Verdeckung der Flächenelemente zu berücksichtigen, vgl. [Rau93, FDFH90]:

(2.9)

Abbildung in dieser Leseprobe nicht enthalten

[2]

Abbildung in dieser Leseprobe nicht enthalten

In Abschnitt 2.2 wurde bereits die Reziprozität der Formfaktoren erwähnt. Zum Nachweis dieser Wechselseitigkeit leitet [SHL91] den Formfaktor zwischen Aj und Ai auf ähnliche Weise wie (2.9) her. Da Vis(dAi, dAj ) = Vis(dAj , dAi) gilt und die Integrale nach dem Satz von Fubini in ihrer Reihenfolge vertauscht werden dürfen, ergibt sich hierfür folgende Darstellung:

[2]

Abbildung in dieser Leseprobe nicht enthalten

(2.10)

Abbildung in dieser Leseprobe nicht enthalten

Aus den identischen Doppelintegralen und den gleichen Sichtbarkeiten in (2.9) und (2.10) resultiert die Reziprozitätsbeziehung für Formfaktoren zwischen endlichen Flächenelementen:

(2.11)

Abbildung in dieser Leseprobe nicht enthalten

Cohen und Wallace stellen in [CW93] eine Reihe von analytischen und numerischen Techni- ken zur Lösung der Gleichung (2.9) vor. Neben verbreiteten Vorgehensweisen beschreiben sie dabei u. a. auch Monte-Carlo-Methoden[3]. Die folgenden Unterabschnitte widmen sich aus- gehend vom Halbkugelverfahren dem in der Praxis häufig verwendeten Halbwürfelverfahren[4] sowie einem Algorithmus, welcher Formfaktoren mithilfe von Kreisscheiben approximiert. Bei allen drei Ansätzen handelt es sich um numerische Vereinfachungen. Analytische Lösungen des obigen Doppelintegrals beschränken sich meist nur auf bestimmte Formen von Flächen- elementen und sind extrem aufwendig.

2.3.1 Halbkugelverfahren

Ziel des Halbkugelverfahrens ist das Finden einer näherungsweisen numerischen Lösung der Gleichung (2.9), um die analytische Berechnung des enthaltenen Doppelintegrals zu vermei- den. Die folgende Beobachtung führt zu einer ersten Vereinfachung: Wenn die beiden Flächen- elemente Ai und Aj nicht teilweise voreinander verdeckt sind und ihr Abstand im Vergleich

[3]Eine Monte-Carlo-Methode ist ein stochastisches Verfahren, bei dem durch Zufallsexperimente versucht wird, nicht oder nur aufwendig analytisch lösbare Probleme numerisch zu lösen. Als Rechtfertigung dient dabei vor allem das Gesetz der großen Zahlen.

[4]Das Halbwürfelverfahren ist unter anderem in den Programmen Blender [ble], tucan radiosity [awa] und RadiosGL [har] implementiert.

2.3 Bestimmung der Formfaktoren

Abbildung 2.5: Halbkugelverfahren: nach Nusselts Analogon ist der Formfaktor des diffe- rentiellen Flächenelements dAi zum endlichen Element Aj proportional zur doppelten Projektion von Aj auf die Grundfläche einer über dAi errichteten Einheitshalbkugel. Quelle: Nach [CG85].

zu ihren Flächeninhalten groß ist, kann das innere Integral als nahezu konstant betrach- tet werden, da der Winkel φi für verschiedene differentielle Flächenelemente von Aj nahezu gleich ist. Das äußere Integral vereinfacht sich damit zu einer Multiplikation mit Eins, sodass sich der Formfaktor zweier Flächenelemente gut mit der Gleichung (2.7) approximieren lässt, vergleiche [CG85, Rau93]:

Abbildung in dieser Leseprobe nicht enthalten

(2.12)

Abbildung in dieser Leseprobe nicht enthalten

dAi sei dabei in der Mitte von Ai positioniert. Ist eine der obigen Bedingungen nicht erfüllt, unterteilt man die Flächenelemente solange in kleinere, bis die gewünschte Genauigkeit der Näherung (2.12) wiederhergestellt ist.

Nusselt entwickelte ein geometrisches Analogon zur Lösung des verbliebenen Integrals, vgl. [SHL91]. Demnach errichtet man mittig über dem differentiellen Flächenelement dAi eine Halbkugel mit dem Einheitsradius. Anschließend wird Aj auf deren Oberfläche projiziert. Der Formfaktor von dAi zur Fläche Aj ergibt sich dann als das Verhältnis zwischen der Projektion der entstandenen Fläche auf die Grundfläche der Halbkugel (siehe Abbildung 2.5) und dem Flächeninhalt dieser Grundfläche. Sei Aij die Fläche dieser Projektion auf die Grundfläche der Halbkugel. Mit der Gleichung (2.12) gilt dann für den Formfaktor FAiAj :

Abbildung in dieser Leseprobe nicht enthalten

(2.13)

wobei πr[2] bzw. π die Fläche des Einheitskreises darstellt.

Zur einfacheren Handhabung führt man eine Diskretisierung der Halbkugeloberfläche in hinreichend kleine Oberflächenstücke durch, vgl. [CG85]. Für jedes dieser Oberflächenstücke wird dann der sogenannte Delta-Formfaktor ermittelt. Er repräsentiert einen Formfaktor- wert für den Raumwinkel, den das zugehörige Stück Halbkugeloberfläche abdeckt. Um einen Formfaktor zwischen zwei Flächenelementen Ai und Aj zu berechnen, projiziert man Aj auf die über Ai errichtete diskretisierte Halbkugel und addiert die Delta-Formfaktoren der

Abbildung in dieser Leseprobe nicht enthalten

Oberflächenstücke der Halbkugel, welche von der Projektion überdeckt werden. Jedes Oberflächenstück liefert dabei nur einen Beitrag zum Formfaktor des nächstgelegenen Flächenelements, dessen Projektion dieses Stück überlappt. Mögliche Verdeckungen der Flächenelemente werden dadurch implizit berücksichtigt. Ein weiterer Vorteil ist, dass man die DeltaFormfaktoren nur einmal berechnen muss, da sie lediglich von der Unterteilung der Halbkugeloberfläche abhängen und eine einmal diskretisierte Halbkugel wieder verwendet werden kann. Dennoch ist das Halbkugelverfahren nicht praktikabel. Die Schwierigkeit besteht in der Erzeugung gleich großer Oberflächenstücke sowie dem Berechnen linearer Koordinaten zur Beschreibung von Punkten auf der Halbkugeloberfläche.

2.3.2 Halbwürfelverfahren

Das von Cohen und Greenberg entwickelte Halbwürfelverfahren, vergleiche [Rau93, FDFH90, CW93], verfolgt den gleichen Ansatz wie das Halbkugelverfahren. Jedoch verwendet es zur Vereinfachung anstelle einer Halbkugel einen Halbwürfel. Dieser entsteht aus einem Würfel der Kantenlänge 2 mit dem Mittelpunkt in dAi, sodass eine Hälfte des Würfels unter dem Horizont von dAi verschwindet, vgl. Abbildung 2.6. Zur Berechnung des Formfaktors FAj Aj projiziert man Aj auf die Flächen des diskretisierten Halbwürfels und bildet die Summe über die DeltaFormfaktoren der überdeckten Pixel. Ein Pixel bezeichnet dabei einen der quadratischen Bereiche, in die man die Würfelflächen bei der Diskretisierung unterteilt. Beschreibt ΔFq den Delta-Formfaktor zum überdeckten Pixel q, so gilt demnach:

(2.14)

Abbildung in dieser Leseprobe nicht enthalten

mit R als Anzahl der Pixel, die von der Projektion von Aj überdeckt werden. Keiner dieser Pixel wird dabei von der Projektion einer näher zum Halbwürfel liegenden anderen Fläche überdeckt. Zudem muss man beachten, dass die Pixel der Oberfläche des Halbwürfels je nach Lage unterschiedlich große Raumwinkel abdecken und somit verschiedene Delta-Formfaktoren haben. Bleibt noch zu klären, wie man die Delta-Formfaktoren bestimmt.

Dazu verwendet man ein rechtshändiges Koordinatensystem mit Ursprung in dAi, dessen z-Achse senkrecht zur Deckfläche des Halbwürfels verläuft. Betrachtet man einen Pixel q auf der Deckfläche des Halbwürfels (vergleiche Abbildung 2.7a), so gilt nach Pythagoras für den Abstand r des Mittelpunkts von q zum Koordinatenursprung:

Abbildung in dieser Leseprobe nicht enthalten

Dabei habe der Mittelpunkt von q die Koordinaten (xq , yq , 1). Des Weiteren lässt sich für die

Winkel φi und φj in Abbildung 2.7a Folgendes feststellen:

Abbildung in dieser Leseprobe nicht enthalten

(2.16)

Den Delta-Formfaktor für den Pixel q erhält man mithilfe des Flächeninhalts Δq von q und der Gleichung ([2].[6]) für Formfaktoren zwischen differentiellen Flächenelementen:

(2.17)

Abbildung in dieser Leseprobe nicht enthalten

2.3 Bestimmung der Formfaktoren

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.6: Halbwürfelverfahren: der Formfaktor des differentiellen Flächenelements dAi zum endlichen Element Aj wird durch Projektion von Aj auf einen über dAi errichteten Halbwürfel bestimmt. Quelle: Nach [CG85].

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.7: Herleitung des Delta-Formfaktors für (a) einen Pixel auf der Deckfläche des Halbwürfels sowie für (b) einen Pixel auf einer der Seitenflächen. Quelle: Nach [CG85].

Dieser Zusammenhang lässt sich mittels (2.15) und (2.16) zu

Abbildung in dieser Leseprobe nicht enthalten

umformen. Damit können die Delta-Formfaktoren für die Pixel der Deckfläche des Halbwürfels berechnet werden. Für die Pixel der Seitenflächen, vgl. zum Beispiel [FDFH[90]], gilt gemäß Abbildung [2].[7]b:

(2.19)

Abbildung in dieser Leseprobe nicht enthalten

r.

Durch eine analoge Herleitung vereinfacht sich (2.17) damit zu

Abbildung in dieser Leseprobe nicht enthalten

für einen Pixel der Seitenflächen senkrecht zur x-Achse und zu

Abbildung in dieser Leseprobe nicht enthalten

für einen Pixel der Seitenflächen senkrecht zur y-Achse. Aus Symmetriegründen braucht nicht zu jedem Pixel der Delta-Formfaktor bestimmt werden. Laut [FDFH[90]] genügt es diesen für ein Achtel der Pixel der Deckfläche sowie ein Viertel der Pixel einer Seitenfläche zu berechnen. Wie beim Halbkugelverfahren müssen auch hier die Delta-Formfaktoren nur einmal ermittelt werden. Sie können bei der Berechnung der Formfaktoren nach Gleichung ([2].[14]) stets wieder verwendet werden.

2.3.3 Approximation mittels Kreisscheiben

Eine praxisrelevante Alternative zum Lösen des verbliebenen Integrals in der Gleichung (2.12) zur näherungsweisen Berechnung der Formfaktoren bietet die Approximation mittels Kreis- scheiben. Sie wurde von Wallace et al. [WEH89] ursprünglich für ein Radiosity-Verfahren mit schrittweiser Verfeinerung (weitere Informationen dazu folgen in Abschnitt 2.4) entwickelt, ist aber keinesfalls auf ein solches Vorgehen beschränkt. Zwei wichtige Vorteile dieser Me- thode sind ihre gute Parallelisierbarkeit, was dem Einsatz auf dem Grafikprozessor entgegen kommt, sowie die Verringerung von Artefakten bei nah nebeneinander liegenden größeren Flächenelementen, vgl. [CHL04].

Ziel ist also erneut die Bestimmung des Formfaktors FdAiAj zwischen dem differentiellen Flächenelement dAi und der endlichen Fläche Aj als Näherung für FAiAj . Den Ausgangspunkt bildet die Beobachtung, dass die Integration mittels uniformer Unterteilung der Fläche Aj durchgeführt werden kann. Dazu unterteilen Wallace et al. Aj gleichmäßig in eine Menge von m kleineren Flächen ΔAj,k mit k = 1, . . . , m und ermitteln FdAiAj als Summe über die Teilfaktoren FdAiΔAj,k . Um eine kostenintensive Berechnung dieser Teilfaktoren zu vermeiden, approximiert man ΔAj,k mit einer Kreisscheibe gleichen Flächeninhalts. Der Formfaktor von einem differentiellen Flächenelement zu einer direkt gegenüber liegenden Kreisscheibe mit gleicher Orientierung und der Ausdehnung ΔAj,k ergibt sich als

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.8: Formfaktorapproximation zwischen differentiellem Flächenelement dAi und endlichem Flächenelement Aj mittels Kreisscheiben. Quelle: Nach [CW93].

Eine Verallgemeinerung wird durch Einbeziehung der Kosinus der jeweiligen Winkel zwischen Flächennormale und Verbindungsgerade erreicht. Näherungsweise erhält man damit

Abbildung in dieser Leseprobe nicht enthalten

für den genannten Formfaktor in Abhängigkeit der Lage der beiden Flächen zueinander. Bei Verwendung von m Kreisscheiben der GrößeAjm giltdementsprechendfürdengesuchten Formfaktor zwischen differentieller Fläche und dem Element Aj :

(2.24)

Abbildung in dieser Leseprobe nicht enthalten

unter Berücksichtigung möglicher Verdeckungen durch andere Flächenelemente der Szene. Abbildung 2.8 veranschaulicht den geometrischen Aufbau. Bei fein unterteilten Szenen erzielt der Algorithmus bereits mit m = 1 gute bis sehr gute Ergebnisse, vgl. Abschnitt 6.1.3. Die hier beschriebene zusätzliche Unterteilung der Flächenelemente zur Lösung des Integrals kann so eingespart werden. [CW93, WEH89]

2.4 Verbesserung der Leistungsfähigkeit

Als größter Nachteil des Radiosity-Verfahrens gilt sein enormer Verbrauch an Rechenleistung und Speicherplatz. Es existiert eine Reihe von Ansätzen, die versuchen diesen zu reduzieren. Im Folgenden werden die zwei wichtigsten Methoden kurz erläutert. Dabei wird neben der hierarchischen Technik auf ein Verfahren mit schrittweiser Verfeinerung eingegangen. Für de- taillierte Informationen sowie weitere Ansätze sei hiermit auf [CW93] und [Rau93] verwiesen.

2.4.1 Hierarchisches Radiosity

Hintergrund des hierarchischen Radiosity-Verfahrens ist die Beobachtung, dass die Anzahl der zu berechnenden Formfaktoren quadratisch mit der Anzahl der Flächenelemente in der Szene wächst. Eine grobe Segmentierung der Objekte bringt daher eine erhebliche Einsparung an Laufzeit. Andererseits nimmt mit feinerer Unterteilung die Genauigkeit der errechneten

Abbildung 2.9: Adaptive Unterteilung beim hierarchischen Radiosity. In Regionen mit großen Strahlungsunterschieden sind die Flächen feiner untergliedert als in Regionen mit weniger Strahlungsschwankungen. Dies wird hier an der sehr feinen Unterteilung der Bereiche nahe der Schattengrenzen deutlich. Quelle: [Pow].

Strahlungswerte zu, was eine bessere Qualität des erzeugten Bildes zur Folge hat. Eine fei- ne Unterteilung ist dabei vor allem in Regionen mit großen Strahlungsunterschieden wie beispielsweise an Schattengrenzen von Bedeutung. In Bereichen mit geringen Strahlungsun- terschieden ist dies dagegen weniger relevant. Es bietet sich also ein Verfahren mit adaptiver Unterteilung an, das in Regionen mit großen Strahlungsunterschieden sehr fein und ansonsten gröber arbeitet, vgl. Abbildung 2.9. Hierbei besteht jedoch das Problem, dass vor der Berech- nung der Strahlungswerte nicht bekannt ist, welche Bereiche große Strahlungsunterschiede aufweisen. Deshalb berechnet man die Strahlungswerte zunächst für eine grobe Unterteilung und unterteilt die Flächenelemente in Regionen mit großen Strahlungsunterschieden dann weiter. Die Formfaktoren der nicht weiter unterteilten Flächenelemente werden dabei nicht erneut berechnet. Damit müssen auch die Strahlungswerte nur für die neu entstandenen Flächenelemente neu bestimmt werden. Diese können mithilfe der Formfaktoren der Unte- relemente aus den Strahlungswerten der Anfangsunterteilung berechnet werden, sodass die Gleichung (2.5) nicht neu gelöst werden braucht. Die Unterteilung in Unterelemente wird bis zur Erreichung der gewünschten Genauigkeit wiederholt. Als nachteilig kann sich eine zu grobe Anfangsunterteilung auswirken, da mit ihr die Gefahr wächst, dass beispielsweise die Schatten von kleineren Objekten nicht bemerkt werden. [Rau93, FDFH90]

2.4.2 Schrittweise Verfeinerung

Laut [Rau93] und [FDFH90] hat das Radiosity-Verfahren in der bisher beschriebenen Form für die praktische Verwendung zwei wesentliche Nachteile. Zum einen kann das Ergebnis erst nach Bestimmung aller Formfaktoren dargestellt werden, was das Konstruieren komplexer Szenen erschwert, da nach jeder vorgenommenen Modifikation gewartet werden muss, bis der gesamte Berechnungsprozess erneut durchlaufen wurde. Hier wäre eine Vorgehensweise wünschenswert, die schnell eine erste Näherungsdarstellung liefert und diese mit fortschreiten- der Berechnungsdauer schrittweise verbessert. Zum anderen resultiert aus der Notwendigkeit der Bestimmung aller Formfaktoren ein Speicherplatzbedarf, der quadratisch mit der An-

2.4 Verbesserung der Leistungsfähigkeit

zahl der Flächenelemente wächst. Dies führt bei detailreichen Szenen schnell dazu, dass das Verfahren auf handelsüblichen Rechnern nicht mehr anwendbar ist.

Eine Abhilfe schafft die in [CCWG88] vorgestellte Methode der schrittweisen Verfeinerung, die auch unter Progressive Refinement bzw. Progressive Radiosity bekannt ist. Durch eine on-the-fly-Berechnung der Formfaktoren kommt diese mit linearem Speicherplatzbedarf aus. Zudem erhält man bereits nach relativ kurzer Zeit eine Näherungsdarstellung, welche mit je- dem weiteren Berechnungsschritt verbessert wird. Dazu kehrt man den bisherigen Ansatz der Berechnung der Strahlungswerte um. Anstatt zu betrachten, wie viel Strahlung ein Flächen- element aufnimmt, verfolgt man hier die von diesem Flächenelement abgegebene Strahlung und verteilt diese in der Umgebung. Es kann gezeigt werden, dass sich die Radiosity-Gleichung (2.4) zu diesem Zweck nach

(2.25)

Abbildung in dieser Leseprobe nicht enthalten

umformen lässt. Bei Verwendung eines der iterativen Verfahren aus Abschnitt 4.1 enthält die Darstellung der Szene bereits nach einem kompletten Iterationsschritt die vollständige direk- te Beleuchtung. Mit jedem weiteren Iterationsschritt wird die Szene durch die zunehmende indirekte Beleuchtung immer besser ausgeleuchtet. Durch geeignete Wahl des als nächstes be- trachteten Flächenelements (shooter ) genügen sogar wenige Lösungsschritte zur Darstellung einer ersten Näherung. Als Lösungsschritt wird dabei die Verteilung der noch nicht abge- gebenen Strahlung eines Flächenelements bezeichnet. Jeder Iterationsschritt besteht aus der Durchführung eines solchen Lösungsschritts für alle Flächenelemente der Szene. Ersetzt man in Gleichung (2.25) Fji mithilfe der Reziprozitätsbeziehung (2.3) durch FijAi, so genügt zurAj Bestimmung der für einen Lösungsschritt benötigten Formfaktoren die einmalige Anwendung des Halbwürfelverfahrens auf das aktuell betrachtete Flächenelement i, vergleiche [Rau93].

3 GPGPU

Angetrieben von der Computerspieleindustrie entwickelte sich die Grafikhardware innerhalb der letzten Jahre zu einer sehr leistungsfähigen Komponente heutiger Rechnersysteme. Auf den Grafikkarten verbaute Prozessoren (GPUs) übersteigen das Leistungsvermögen der Pro- zessoren von Mainboards (CPUs) inzwischen um ein Vielfaches (siehe Abbildung 3.1). Zu- dem sind Grafikkarten meist günstiger zu haben als vergleichbare CPUs. Beispielsweise lag der Preis eines 3.0 GHz Intel Core2 Quad (QX6850) im August 2007 bei $ 1100. Dieser besitzt eine Spitzenleistung von 96 GFLOPS sowie eine Speicherbandbreite von 21 GB/s. Die zu diesem Zeitpunkt aktuelle Grafikkarte Geforce 8800 GTX von NVIDIA kostete dage- gen nur $ 550 und übertrifft den Intel-Prozessor beim erreichten Leistungsniveau deutlich. Laut GPUBench-Projekt der Stanford-Universität beträgt die theoretisch mögliche Verarbei- tungsgeschwindigkeit 330 GFLOPS und die Speicherbandbreite 55.2 GB/s. Dabei ist jedoch zu berücksichtigen, dass diese empirisch bestimmten Werte aufgrund der Architektur nur bei datenparallelen Prozessen erreicht werden können. [Houa]

Mit der Entwicklung programmierbarer Grafikhardware kam zunehmend der Wunsch auf, die hohe Leistungsfähigkeit der GPU auch für andere Anwendungsgebiete als ausschließlich der Darstellung von 3D-Szenen zu nutzen. Unter der Bezeichnung GPGPU entstand eine recht lebhafte und stetig wachsende Community der Anhänger dieser Bewegung.[5] GPGPU steht hierbei für

”General-PurposecomputationontheGraphicsProcessingUnit“(engl.)und meint damit die Zweckentfremdung der GPU als Koprozessor für allgemeine Aufgaben. Die Vielfalt der GPGPU-Programme reicht dabei von der Umsetzung globaler Beleuchtungsverfahren wie dem Raytracing oder dem Radiosity-Verfahren über numerische und physikalische Berechnungen bis hin zum Datamining, vgl. [OLG+07, gpgb].

In den folgenden Abschnitten werden einige Grundlagen der GPU-Programmierung erörtert sowie diverse GPU-spezifische Programmiersysteme vorgestellt. Darüber hinaus wird eine Einführung in die für diese Arbeit verwendete Entwicklungsplattform gegeben.

[5]Eine Homepage der GPGPU-Community ist unter [gpga] zu finden.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 3.1: Vergleich des Wachstums der Leistungsfähigkeit von CPUs und GPUs. Das Leistungsvermögen von GPUs erhöhte sich im Vergleich zu dem von CPUs in den letzten Jahren erheblich. Quelle: Nach [Houa].

3.1 Grundlagen der GPU-Programmierung

Die folgenden Unterabschnitte skizzieren angelehnt an [OLG+07] die Entwicklung der GPU. In diesem Zusammenhang werden neben der Architektur der Grafikhardware auch relevan- te Modelle zur Implementierung entsprechender Software behandelt. Aufgrund der für diese Arbeit verfügbaren Hardware beschränken sich die folgenden Aussagen auf Grafikkartenge- nerationen bis zur GeForce 7 Serie von NVIDIA und vergleichbarer Grafikhardware.

3.1.1 Grafikpipeline

Der Bereich interaktiver 3D-Grafik unterscheidet sich von allgemeinen Anwendungsfeldern durch einige wesentliche Merkmale. Insbesondere verlangen interaktive 3D-Anwendungen nach sehr hohen Verarbeitungsgeschwindigkeiten, besitzen aber zugleich ein hohes Maß an Parallelisierbarkeit. Die Entwicklung einer speziellen Hardware, welche sich die natürliche Parallelität der Applikation zunutze macht, erlaubt eine höhere Leistungsfähigkeit von Grafikanwendungen als dies traditionelle Mikroprozessoren ermöglichen.

Heutige Grafikchiphersteller strukturieren ihre Grafikberechnungen alle auf ähnliche Weise in der sogenannten Grafikpipeline (auch Rendering-Pipeline). Diese ist so designt, dass sie durch parallele Ausführung hohe Verarbeitungsgeschwindigkeiten erreicht. Sie ist in verschie- dene Stufen gegliedert, die von allen geometrischen Primitiven vor der Darstellung auf dem

Bildschirm durchlaufen werden. Eine vereinfachte Pipeline besteht aus den Stufen Operations“ (Operationen auf Eckpunkten),

”Vertex

mitive und Clipping[6]),

”PrimitiveAssembly“(ZusammensetzenderPri- ”Rasterization“(KonvertierenderGeometriedateninBildpunkte)und

”FragmentOperations“(OperationenaufBildpunkten),wobeidiegebräuchlichenenglischen Termini beibehalten wurden. Abbildung [3].[2] zeigt eine solche Pipeline im Überblick. Jede der

Pipelinestufen ist als separates Stück Hardware auf der GPU implementiert und kann so parallel zu den anderen Stufen ausgeführt werden. Eine solche Organisation wird auch als taskparallel bezeichnet. Für tiefgreifende Informationen über die Grafikpipeline wird in [KF[05]] und [MM[05]] die GPU der NVIDIA GeForce [6] Serie näher beschrieben.

3.1.2 Programmierbare Hardware

Neben der zunehmenden Leistungsfähigkeit der Grafikhardware setzte man sich zum Ziel, den Realismus gerenderter Bilder mit jeder neuen GPU-Generation zu steigern. Bei der oben beschrieben Grafikpipeline handelte es sich ursprünglich um eine sogenannte Fixed-Function- Pipeline mit fest verdrahteten Operationen in jeder Stufe. Der Erfolg von Rendersystemen wie RenderMan von Pixar zeigte allerdings die Vorteile einer flexibleren Gestaltung. Anstelle die Beleuchtungs- und Schattierungsoperationen auf einige wenige zu beschränken, führte RenderMan auf jedem Primitiv ein benutzerdefiniertes Shader-Programm aus und erzielte damit beeindruckende visuelle Ergebnisse.

Innerhalb der letzten Jahre formten die Grafikkartenhersteller die Fixed-Function-Pipeline deshalb zu einer flexibleren, programmierbaren Pipeline um. Diese Bemühungen konzentrier- ten sich zunächst auf zwei der Pipelinestufen: die Vertex und die Fragment Operations. Die Stufe der Vertex Operations führte in der Fixed-Function-Pipeline Operationen auf Eckpunk- ten durch, z. B. Transformationen und Lichtberechnungen. In der programmierbaren Pipeline

[6]Unter Clipping versteht man das Abschneiden nicht sichtbarer Teile grafischer Primitive am Rand eines gewünschten Bildschirmausschnitts.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 3.2: Vereinfachter Aufbau einer programmierbaren Grafikpipeline. Die gebräuchlichen englischen Bezeichnungen wurden übernommen. Bei moder- nen Grafikkarten sind Vertex- und Fragment-Prozessor durch den Benutzer programmierbar. Quelle: [Stu07].

wurden diese Operationen gegen ein benutzerdefiniertes Vertex-Programm ausgetauscht. Auf ähnliche Weise ersetzt ein vom Benutzer definiertes Fragment-Programm die Funktionen der Stufe der Fragment Operations. Hier wurden während der Fixed-Function-Ära die Farben der einzelnen Pixel bestimmt, weshalb anstelle von Fragment auch häufig der Begriff Pixel ver- wendet wird. Abbildung 3.2 zeigt die beschriebene Grafikpipeline mit den programmierbaren Shader-Einheiten.

Jede neue Hardwaregeneration erweiterte die Funktionalität und die Zugänglichkeit der beiden programmierbaren Stufen. So war mit der Einführung der ersten programmierbaren Stufe im Jahre 1999 zunächst nur eine begrenzte Kombination zwischen Textur- und inter- polierten Farbwerten zur Berechnung einer Pixelfarbe möglich. Im Jahr 2002 leitete die ATI Radeon 9700 einen weiteren Wechsel ein. Von nun an konnten auch Fließkommaoperationen innerhalb der Fragment-Pipeline durchgeführt werden. Den entscheidenden Schritt zur Ver- wendung der GPU für allgemeine Berechnungen ermöglichte aber erst die Einführung einer vollständig programmierbaren Hardware zusammen mit einer Assemblersprache zur Spezifi- zierung von Programmen, welche auf Eck- bzw. Bildpunkten arbeiteten. Diese programmier- bare Shader-Hardware wurde explizit so designt, dass sie mehrere datenparallele Primitive zur selben Zeit verarbeitet. Sowohl für den Vertex- als auch für den Fragment-Shader exis- tiert ein Standard, welcher im Jahr 2006 bereits in seiner dritten Überarbeitung erschienen ist. OpenGL bietet in seiner API Erweiterungen für beide Standards an. Der Befehlssatz jeder Stufe ist wie bei CPUs beschränkt. Dabei existieren hauptsächlich grafikspezifische ma- thematische Befehle. Zuletzt wurde die Befehlsmenge um einige Operationen zur Steuerung von Programmabläufen erweitert.

Im Allgemeinen übergibt man den programmierbaren Pipelinestufen eine begrenzte An- zahl an Vektoren von je vier 32-Bit Fließkommazahlen. Als Ergebnis des Vertex-Prozessors wird eine begrenzte Anzahl solcher Vektoren an den Rasterizer weitergegeben, welcher diese Werte plus interpolierte Werte für alle dazwischen liegenden Bildpunkte an den Fragment- Prozessor weiterreicht. Letzterer liefert schließlich bis zu vier Vektoren pro Pixel, bei de- nen es sich typischerweise um Farbwerte handelt. Jede programmierbare Stufe kann auf die Register aller Primitive zugreifen, jedoch ohne ihre Inhalte zu ändern. Zudem können pro Pri- mitiv weitere Register ausgelesen und beschrieben werden. Alle Shader-Einheiten unterliegen Beschränkungen bezüglich der Anzahl an Input- und Output-Werten, Konstanten, Registern sowie Instruktionen. Mit jeder Überarbeitung des Vertex- und Pixel-Shader-Standards erhöht sich auch die Anzahl an möglichen Werten.

Für gewöhnlich bestehen GPUs aus mehreren Vertex- und Fragment-Prozessoren. Beispiels- weise besitzt eine ATI Radeon X1900 XTX 8 Vertex- sowie 48 Fragment-Prozessoren. Die Möglichkeit Daten aus dem Texturspeicher zu lesen, erlaubt den Fragment-Prozessoren die Durchführung eines Gather[7]. Jedoch muss die Ausgabeadresse (Texturkoordinate) eines Frag- ments stets vor der Ausführung des Fragment-Programms feststehen und kann zudem nicht von diesem geändert werden. Ein Scatter[8] ist damit im Fragment-Shader nicht möglich. Ver- tex-Prozessoren erwarben dagegen erst kürzlich die Fähigkeit auf Texturen zuzugreifen. Hinzu kommt, dass sie die Positionen der übergebenen Eckpunkte verändern können und so direkt beeinflussen, wo im Bild die entsprechenden Pixel gezeichnet werden. Vertex-Prozessoren sind damit in der Lage, sowohl Gather- als auch Scatter-Operationen durchzuführen. Bedauerli- cherweise kann ein Scatter der Vertex-Prozessoren im Zusammenhang mit Speicherzugriffen und dem Rasterizer zu Problemen in darauffolgenden Pipelinestufen führen. In Kombination mit der niedrigeren Anzahl dieser Prozessoren schränkt dies die Brauchbarkeit eines Scatters auf heutigen GPUs stark ein.

3.1.3 GPU-Programmiermodell

Wie bereits erwähnt, stellt Grafikhardware eine Alternative für Anwendungen dar, die hohe Verarbeitungsgeschwindigkeiten und Bandbreiten benötigen. GPUs erreichen ihre enorme Leistungsfähigkeit vor allem durch Datenparallelität. Diese Parallelität verlangt nach einem Programmiermodell, das sich vom traditionellen sequentiellen Programmiermodell unterscheidet. Aus diesem Grund soll das Programmiermodell der GPU nun näher erläutert werden. Dabei finden sowohl die Terminologie der Grafik-API als auch die des Stream-Processing Anwendung, da beide in der Literatur recht häufig vorkommen.

Das Modell des Stream-Processing bindet Parallelität und Kommunikationsmuster direkt an die Anwendung. Dies geschieht durch die Strukturierung der Daten in sogenannte Streams (Datenströme) und dem Ausdrücken der Berechnung als Kernel. Streams werden in diesem Kontext z. B. durch Arrays, gefüllt mit Elementen vom gleichen Datentyp, repräsentiert. Kernels bezeichnen abgeschlossene Folgen von Berechnungsschritten, welche die Elemente der Eingabeströme unabhängig voneinander bearbeiten und zumeist Ausgabeströme als Er- gebnisse zurückgeben. Dabei werden für jedes Stream-Element exakt die selben Operationen durchgeführt.

Da Szenen im Allgemeinen aus mehr Fragmenten als Eckpunkten bestehen, ist in heutigen GPUs der Fragment-Prozessor die Pipelinestufe mit der höchsten Leistungsfähigkeit. Ein typisches GPGPU-Programm verwendet deshalb auch hauptsächlich den Fragment-Prozessor als Rechenmaschine. Derartige Programme weisen folgende Struktur auf:

1. Zunächst ermittelt der Programmierer die datenparallelen Operationen seiner Applika- tion, um die Anwendung anschließend in unabhängige parallele Teile zu zerlegen. Jedes dieser Teile stellt einen Kernel dar und wird nun als Fragment-Programm implemen- tiert. Ein- und Ausgabewerte der Fragment-Programme sind ein oder mehrere Daten-

[7] Kommunikationsoperation auf Parallelrechnern, bei der ein Prozessor von jedem anderen Prozessor eine Nachricht erhält.

[8] Kommunikationsoperation auf Parallelrechnern, bei der ein Prozessor an jeden anderen Prozessor eine evtl. unterschiedliche Nachricht schickt.

Arrays, welche als Texturen im Grafikspeicher abgelegt werden. In der Terminologie des Stream-Processing sind die Daten in den Texturen demnach Streams und ein FragmentProgramm wird jeweils parallel auf den Stream-Elementen aufgerufen.

2. Um ein Fragment-Programm aufrufen zu können, muss die Größe des Ausgabe-Stre- ams bereits spezifiziert sein. Der Programmierer erledigt dies durch das Übergeben von Eckpunkten an die Grafikpipeline. Ein typischer GPGPU-Aufruf gleicht dabei einem parallel zur Bildebene orientierten Rechteck. Die Größe des Output-Arrays stimmt dann mit der davon überdeckten rechteckigen Pixelregion überein. Man beachte, dass sich GPUs durch die Verarbeitung zweidimensionaler Arrays auszeichnen, bei eindimensionalen Arrays hingegen Beschränkungen unterliegen.

3. Der Rasterizer generiert für jeden Bildpunkt des Rechtecks ein Fragment und produziert so Tausende bis Millionen Fragmente.

4. Jedes der generierten Fragmente wird nun vom aktiven Programm des Fragment-Pro-

zessors verarbeitet. Man beachte, dass dabei für jedes Fragment der gleiche Kernel ausgeführt wird. Das Fragment-Programm kann zwar via Texturzugriff jede beliebige Stelle des globalen Texturspeichers auslesen, im Gegensatz dazu aber nur an die vom Rasterizer vorgegebene, zum Fragment korrespondiere Stelle im Framebuffer schreiben. Der Bereich jeder Eingabetextur (Stream), auf dem die Berechnung stattfindet, wird durch die Texturkoordinaten der zuvor übergebenen Eckpunkte festgelegt und für die einzelnen Fragmente interpoliert. Texturkoordinaten können zudem unabhängig für je- de Eingabetextur spezifiziert und auch on-the-fly vom Fragment-Programm berechnet werden. Dies ermöglicht den oben genannten beliebigen Speicherzugriff.

5. Als Resultat des Fragment-Programms erhält man zu jedem Fragment einen Wert bzw.

Vektor. Diese Werte können entweder Endergebnis der Anwendung sein oder auch als Textur gespeichert und für weitere Berechnungen verwendet werden. Letzteres erlaubt auch komplexere Berechnungen, welche möglicherweise einige oder gar dutzende solcher Durchläufe durch die Pipeline benötigen. Man spricht in diesem Zusammenhang auch von Multipass.

Während die Komplexität eines einzelnen Durchlaufs durch die Grafikpipeline beschränkt sein kann, zum Beispiel durch die Anzahl der Befehle, die Anzahl der erlaubten Ausgabe- werte oder durch eine begrenzte Komplexität der Steuerstrukturen, erlaubt die Verwendung mehrerer Durchläufe die Implementierung von Programmen beliebiger Komplexität. Peercy et al. demonstrieren dies in [POAU00], wo sie anhand einer um einige Fließkommaopera- tionen erweiterten OpenGL-Simulation zeigen, dass sich mit genug Durchläufen beliebige RenderMan-Shader sogar auf der Fixed-Function-Pipeline umsetzen lassen.

3.1.4 Ablaufsteuerung von GPU-Programmen

Strukturen zur Steuerung von Programmabläufen sind im Gebiet computergestützter Berech- nungen fundamentale Prinzipien. Sowohl Schleifen als auch bedingte Sprünge (Verzweigun- gen/Branching) zählen heute zu solch grundlegenden Konzepten, dass es sehr entmutigend sein kann, Software für Plattformen zu entwickeln, die diese nur in beschränktem Umfang unterstützen. Die neuesten Grafikprozessoren erlauben Branching innerhalb von Vertex- und Fragment-Programmen in verschiedenen Formen, doch ihre stark parallele Bauweise erfordert

Umsicht bei der Art der Umsetzung. Dieser Abschnitt betrachtet einige der Einschränkungen bedingter Sprünge auf heutigen GPUs und beschreibt außerdem eine Reihe von Techniken zur Realisierung von Schleifen und Verzweigungen in GPGPU-Programmen. Detailliertere Informationen zur Steuerung von Programmabläufen auf GPUs sind in [HB05] zu finden.

Hardwaremechanismen

Auf heutigen GPUs existieren drei grundlegende Implementierungen des datenparallelen Branching: Predication, MIMD-Branching und SIMD-Branching.

Architekturen, die ausschließlich Predication unterstützen, besitzen keine wirklich daten- parallelen Branch-Befehle. Stattdessen berechnet die GPU beide Zweige einer if-Anweisung und verwirft anschließend eines der Ergebnisse in Abhängigkeit der if-Bedingung. Der Nach- teil einer solchen Architektur ist, dass eine Berechnung beider Zweige sehr kostenintensiv sein kann. Dennoch unterstützen nicht alle modernen GPUs echtes datenparalleles Branching. Er- laubt die Ziel-GPU nur Predication zur Steuerung des Programmablaufs, so generieren die Compiler der Shader-Hochsprachen wie Cg oder OpenGL Shading Language automatisch die entsprechenden Assemblerbefehle.

Nach der Klassifizierung von Flynn können die Prozessoren einer Multiple-Instruction-Mul- tiple-Data-Architektur (MIMD) zur selben Zeit unterschiedliche Aufgaben bewältigen. Die Prozessoren einer so designten Grafikhardware mit Branching-Unterstützung sind demnach in der Lage zeitgleich verschiedene Verzweigungspfade eines GPU-Programms unabhängig voneinander zu durchlaufen. In einer Single-Instruction-Multiple-Data-Architektur (SIMD) müssen dagegen alle aktiven Prozessoren zur selben Zeit den selben Rechenschritt ausführen. Laut [OLG+07] ist die MIMD-Bauweise in neueren Grafikkarten bisher nur bei Vertex-Prozes- soren vorzufinden, z. B. in der GeForce 6 und 7 Serie sowie in einigen Quadro-Modellen von NVIDIA. Eine Klassifizierung der Fragment-Prozessoren gestaltet sich etwas schwie- riger. Das Verhalten entspricht gewissermaßen einem Single-Program-Multiple-Data-Modell (SPMD), was bedeutet, dass die parallel arbeitenden Threads verschiedene Verzweigungs- pfade verfolgen können. Betrachtet man Architektur und Leistungsfähigkeit von Fragment- Prozessoren heutiger GPUs, so bearbeiten diese die Fragmente in SIMD-Gruppen. Ist in einer SIMD-Gruppe die Branch-Bedingung für alle Fragmente identisch, so braucht nur der verwen- dete Zweig berechnet werden. Werten allerdings ein oder mehrere Prozessoren die Branch- Bedingung verschieden aus, so müssen für alle Fragmente beide Zweige berechnet und die nicht benötigten Ergebnisse entsprechend verworfen werden. Dies hat zur Folge, dass Diver- genzen in den Branch-Bedingungen gleichzeitig bearbeiteter Fragmente zu Leistungseinbußen führen können.

Verlagerung des Branching in frühere Pipelinestufen

Da explizites Branching auf der GPU zu oben genannten Leistungseinbußen führen kann, ist es sinnvoll mehrere Techniken zur Reduzierung der Branching-Kosten zu kennen. Eine sehr nützliche Strategie ist es, Entscheidungen über den weiteren Programmablauf in eine frühere Pipelinestufe zu verlagern, wo sie effizienter berechnet werden können.

Statische Branch-Auflösung Wie auf der CPU, ist es auch auf der GPU günstig, Verzwei- gungen innerhalb von Schleifen zu vermeiden. Bestimmt man beispielsweise eine partielle Ableitung einer diskretisierten Fläche (Gitter), so zerlegt eine effiziente Implementierung die

3.2 Programmiersysteme

Berechnung in mehrere Schleifen: eine Schleife über das Innere des Gitters ohne Einbeziehung der Punkte am Rand sowie eine oder mehrere Schleifen über die Randpunkte. Die resultie- renden Schleifen dieser statischen Branch-Auflösung beinhalten effizienten Quelltext ohne Sprünge. Bezogen auf das Modell des Stream-Processing gleicht diese Technik dem Unter- teilen eines Streams in mehrere Substreams. Auf der GPU wird die Berechnung auf mehrere Fragment-Programme verteilt, sodass innere Punkte und Randpunkte von unterschiedlichen Programmen bearbeitet werden. Das innere Programm wird dabei auf die inneren Pixel eines im Grafikspeicher abgelegten Rechtecks angewendet, das Randprogramm auf die Pixel der Strecken zwischen den Eckpunkten dieses Rechtecks. Die statische Branch-Auflösung wird darüber hinaus von Goodnight et al. [GWL+05], Harris und James [HJ03] sowie Lefohn et al. [LKHW03] diskutiert.

Vorberechnung Beim vorherigen Verfahren war das Branch-Ergebnis über einen großen Bereich an Eingabewerten konstant. In ähnlicher Weise bleibt das Ergebnis eines Branch manchmal für eine gewisse Zeitdauer bzw. eine gewisse Anzahl an Iterationen einer Berech- nung unverändert. Im Fall der Vorberechnung kann ein Branch nur aufgelöst werden, wenn die nächsteÄnderung des Branch-Ergebnisses bekannt ist. Diese Information speichert man für die Verwendung über viele darauffolgende Iterationen. Ein Auswerten der Branch-Bedingung während der Berechnung ist nun nicht mehr notwendig. Dies kann zu hohen Leistungssteige- rungen führen.

Z-Culling Vorberechnete Ergebnisse von Branch-Bedingungen können bereits eine Pipeli- nestufe früher verwendet werden, um mithilfe eines weiteren GPU-Features unnötige Arbeit ganz auszulassen. Moderne Grafikkarten besitzen verschiedene integrierte Funktionen zur Ver- meidung des Schattierens nicht sichtbarer Pixel. Eine dieser Funktionen ist das sogenannte Z-Culling. Dabei werden die Tiefenwerte (Z-Werte) eines eingehenden Pixelblocks mit den Tiefenwerten des zugehörigen Pixelblocks im Z-Puffer verglichen. Bestehen eingehende Pixel den Tiefentest nicht, werden sie verworfen, bevor ihre Farbwerte im Fragment-Prozessor be- rechnet werden. Auf diese Weise werden nur Fragmente weiter verarbeitet, die den Tiefentest bestehen. Unnötige Arbeit wird so eingespart und die Anwendung läuft schneller. In einer Strömungssimulation können zum Beispiel Punkte, die ein festes Hindernis darstellen, mit einem Z-Wert gleich Null maskiert werden, sodass alle strömungsspezifischen Berechnungen für diese Punkte ausgelassen werden. Sind solche Hindernisse ziemlich groß, kann so durch das Nicht-Verarbeiten dieser Punkte eine Menge Arbeit eingespart werden. Sander et al. beschreiben diese Technik in [STM04] zusammen mit einer weiteren Z-Culling-Beschleuni- gungstechnik für Strömungssimulationen, Harris und Buck stellen in [HB05] entsprechenden Pseudocode bereit und Purcell et al. verwenden in [PBMH02] Z-Culling zur Beschleunigung von GPU-Raytracing.

Erfolgreiches Programmieren benötigt unabhängig von der Zielplattform mindestens folgende drei Basiskomponenten: eine Hochsprache zur Quellcodeentwicklung, eine Debugging-Umge- bung sowie Profiling-Tools. CPU-Programmierer können auf eine große Anzahl an bewährten Programmiersprachen, Debuggern und Profilern zurückgreifen. Im Gegensatz dazu haben

GPU-Programmierer nur eine kleine Menge an Programmiersprachen, Debuggern und Profilern zur Auswahl. [OLG+07]

Dieser Abschnitt stellt die bekanntesten Shader-Hochsprachen, einige GPGPU-spezifische Sprachen und Bibliotheken sowie für die GPU-Programmierung verfügbare Debugging-Tools vor. Da Profiler sehr architekturabhängig sind, unterliegt ihre Entwicklung eher den Grafikkartenherstellern. Sie sollen an dieser Stelle deshalb nicht diskutiert werden.

3.2.1 Shader-Hochsprachen

Viele Hochsprachen zur GPU-Programmierung haben eines gemeinsam: sie folgen der Idee, dass die GPU Bilder generiert. Deshalb werden diese auch häufig als Shading-Languages bezeichnet. In diesem Sinne handelt es sich um Hochsprachen, die ein Shader-Programm in einen Vertex- und einen Fragment-Shader kompilieren, um so vom Programm beschriebene Bilder zu erzeugen. Dabei müssen stets beide Shader generiert werden. Als bedeutendste Vertreter sind HLSL, GLSL und Cg zu nennen. Üblicherweise werden Shading-Languages während derÜbersetzung in Zwischensprachen überführt, vgl. [Stu07, OLG+07].

HLSL

Die High Level Shading Language (HLSL) wurde von Microsoft entwickelt und ist nur unter Microsoft Windows mit DirectX nutzbar. Da in der Spieleindustrie bevorzugt DirectX verwendet wird, orientieren sich viele Grafikkartenhersteller an HLSL.

Wie bereits oben angedeutet, kommen bei der Übersetzung Zwischensprachen zum Einsatz. Dabei handelt es sich um spezielle Assemblersprachen für Vertex- und Fragment-Shader, deren Fähigkeiten im sogenannten Shader-Model-Standard beschrieben werden. Mit jeder neuen Version des Shader-Model wird der Funktionsumfang der Assemblersprachen erweitert (siehe Abschnitt 3.1.2). Da die Zwischensprachen der OpenGL Shading Language (GLSL) keinen eigenen Standard besitzen, bezieht man sich auch dort auf das Shader-Model. [Stu07]

GLSL

Die OpenGL Shading Language ist seit der Version 2.0 im OpenGL-Standard integriert und kann ausschließlich mit der OpenGL-API verwendet werden. Sie basiert auf der Programmiersprache C und wurde vom OpenGL-ARB (Architectural Review Board) entwickelt, um Programmierern mehr Kontrolle über die Grafikpipeline zu geben, ohne dabei Assembleroder maschinenspezifische Sprachen verwenden zu müssen. Die Verwendung mit OpenGL 1.5 wird durch einige OpenGL-Erweiterungen ermöglicht. [Stu07, ope]

Cg

Cg steht für

”CforGraphics“undwurdevonNVIDIAentwickelt,umGPU-Anwendungen

weniger plattform- und API-abhängig zu machen. Die Sprache ist an C angelehnt und wird mit dem Compiler des Cg-Toolkits übersetzt. DieÜbersetzung erfolgt dabei in eine andere Shader-Sprache, welche von der verwendeten API unterstützt werden muss. Dementsprechend unterscheidet sich auch der zur Verfügung stehende Funktionsumfang.

Cg bietet den Vorteil, dass Effekte einfacher für verschiedene Hardware implementiert werden können. Die zusätzliche Zwischenschicht verlangsamt jedoch den Kompiliervorgang und stellt zudem eine weitere potentielle Fehlerquelle dar. [Stu07, Kil03]

3.2.2 GPGPU-spezifische Sprachen und Bibliotheken

Shader-Hochsprachen, wie die oben genannten HLSL, GLSL und Cg, erschweren das Entwi- ckeln von GPGPU-Anwendungen, da sie die Programmierer zwingen in Begriffen der Compu- tergrafik zu denken, obwohl das zu entwickelnde Programm womöglich konzeptionell nichts damit zu tun hat. Beispielsweise muss, um Daten aus dem Grafikspeicher auszulesen, auf eine Textur zugegriffen werden. Auch die Implementierung des Radiosity-Verfahrens würde sich sehr kompliziert gestalten, da die zugrunde gelegten Algorithmen größtenteils allgemeiner Natur sind. Im Folgenden wird deshalb eine Reihe von Systemen vorgestellt, welche GPGPU- Funktionalität bereitstellen und zugleich GPU-spezifische Details vor dem Programmierer verbergen.

Accelator

Bei Accelator handelt es sich um ein von Microsoft Research entwickeltes System für das .NET-Framework. Es stellt ein datenparalleles Programmiermodell in einer Bibliothek zur Verfügung, welche über traditionelle imperative Programmiersprachen verwendbar ist. Für die Nutzung der GPU werden datenparallele Arrays als eigener Datentyp eingeführt. Hinzu kommen spezielle Funktionen, die u. a. die explizite Kommunikation zwischen daten- parallelen und normalen Arrays erlauben. Accelator übersetzt datenparallele Operationen zur Ausführungszeit in GPU-Fragment-Shader oder in CPU-Code. [OLG+07, Houb]

BrookGPU

Brook erweitert den ANSI-C-Standard um einige Konzepte des Stream-Processing. Einen Compiler sowie ein Laufzeitsystem mit der GPU als Zielarchitektur stellt die BrookGPU- Software bereit. Diese verfolgt die Idee, die GPU als Datenflusskoprozessor zu nutzen, und ist auf verschiedenen Plattformen einsetzbar. Unterstützt werden DirectX und OpenGL, Grafikkarten von ATI und NVIDIA und sowohl Windows als auch Linux. Dabei hält sich Brook strikt an die Prinzipien des Stream-Processing-Modells. Das heißt, Funktionen werden als Kernel implementiert und auf jedem Stream-Element identisch ausgeführt und Abhängigkeiten zwischen den Elementen sind nicht erlaubt. Bei derÜbersetzung werden Ker- nels und Streams automatisch auf Fragment-Programme und Texturen abgebildet. Brook und BrookGPU sind aktive Forschungsprojekte an der Stanford-Universität. [BFH+04, OLG+07, Houb, BFH+]

CUDA

NVIDIA CUDA (Compute Unified Device Architecture) ist eine neuartige Hardware- und Softwarearchitektur zur Nutzung der GPU für datenparallele Berechnungen, ohne diese auf die Grafik-API abbilden zu müssen. Das CUDA-Toolkit von NVIDIA beinhaltet einen C- Compiler, eine Zwischenschicht in Assemblersprache sowie ein Laufzeitinterface. Es erlaubt das Programmieren in Standard-C und bietet dazu zahlreiche Bibliotheken mit nützlichen mathematischen Funktionen, wie z. B. aus der linearen Algebra und die Fast-Fourier-Trans- formation. CUDA arbeitet Thread-basiert. Die Thread-Programme unterliegen bei der An- zahl der Instruktionen sowie bei Schleifen und bedingten Sprüngen keinen Beschränkungen mehr und besitzen darüber hinaus uneingeschränkten Zugriff auf den Grafikspeicher. Da- mit sind sogar Gather- und Scatter-Operationen auf der GPU in vollem Umfang möglich.

Des Weiteren verfügen CUDA-Grafikkarten über einen zusätzlichen datenparallelen Cache, der die Kommunikation zwischen den Threads unterstützt und Berechnungen auf der GPU weitere massive Beschleunigungen eröffnet. CUDA kann sowohl unter Windows als auch unter Linux verwendet werden. Ein Nachteil ist jedoch, dass es von der Grafikhardware unterstützt werden muss. Bislang ist dies nur bei der GeForce 8 Serie sowie einigen Teslaund Quadro-Lösungen von NVIDIA der Fall. [Buc, nvib, nvia]

PeakStream

Hinter dem Namen PeakStream verbirgt sich eine kommerzielle Plattform zur Entwicklung von Software für Mehr-Prozessor-Maschinen und GPUs. Programmierer können dazu auf eine API im C/C++-Standard zurückgreifen. Mit PeakStream implementierte Berechnun- gen finden anwendungsfallspezifisch automatisch auf der GPU oder aber der CPU statt. Bei Multi-Core-Prozessoren wird die Arbeit entsprechend auf mehrere CPUs verteilt. Neben ei- nigen Mathematikbibliotheken sind auch Debugging- und Profiling-Tools im Produktumfang enthalten. Dabei unterstützt die Plattform sowohl Windows als auch Linux. [Pap]

Berichten im Internet zufolge [Col07, gol07, Sha07] wurde PeakStream im Juni 2007 von Google aufgekauft. Man darf also gespannt sein, wie sich dieses Projekt weiter entwickelt.

RapidMind

RapidMind Inc. gründete sich im Jahr 2004, um die Weiterentwicklung der Metasprache Sh zu kommerzialisieren. Während Sh primär für die GPU-Programmierung entwickelt wurde, verfolgt RapidMind einen allgemeingültigeren Ansatz, der die Unterstützung weiterer Hard- warearchitekturen wie Cell-Prozessoren erlaubt. Der zur Verfügung stehende Funktionsum- fang hängt vom jeweils verwendeten Backend und der Zielarchitektur ab. Momentan existieren Backends für GLSL, Cell-Prozessoren sowie einen Debugging-Modus, der den Quellcode un- ter Verwendung eines C-Compilers auf dem Host-Prozessor ausführt. Auch die RapidMind- Plattform unterstützt sowohl Windows- als auch Linux-Betriebssysteme.

Die zugrunde liegende Metasprache im C++-Standard stellt mit der Definition von sogenannten Streaming-Shadern ein gutes Werkzeug zur GPGPU-Programmierung bereit. Überdies werden mithilfe eines dynamischen Compilers sämtliche hardwarespezifischen Vorgänge vor dem Programmierer verborgen. Zum Beispiel werden notwendige Texturen beim Datenaustausch zwischen CPU und GPU automatisch generiert und ausgelesen. Jedoch darf auch hier nicht außer Acht gelassen werden, dass der zusätzlicheÜbersetzungsvorgang eine weitere potentielle Fehlerquelle birgt. [McC06, rapa, sh]

Shallows

Shallows ist eine C++-Bibliothek, die das datenabhängige Erzeugen und Auslesen von Texturen, und damit die Kommunikation zwischen GPU und CPU, vereinfacht. Sie ist zusammen mit der OpenGL-API sowohl unter Linux als auch unter Windows mit Visual Studio .NET verwendbar. [shab]

Weitere Ansätze

Als weitere Ansätze GPGPU-spezifischer Programmierung sind vor allem CTM sowie CGiS

hervorzuheben.

”CloseTotheMetal“(oderkurzCTM)wirdinderForschungsabteilungvon

AMD/ATI als Konkurrenz zu NVIDIAs CUDA entwickelt und besteht aus einem SoftwareDevelopement-Kit mit eigener API [Hen, ati06], während CGiS eine an der Universität des Saarlandes konzipierte datenparallele Hochsprache ist, die in Syntax und Semantik C ähnelt [FLCG, OLG+07]. Ferner sind noch zwei Projekte namens Scout und Glift zu erwähnen, die sich ebenfalls dem Thema GPGPU widmen, vgl. hierzu [OLG+07].

3.2.3 Debugging-Tools

Debugging auf GPUs, vgl. [OLG+07], unterlag lange Zeit starken Beschränkungen. Hinzu kam, dass notwendige Funktionen eines guten GPU-Debuggers nicht ausreichend definiert waren. Das Aufkommen der GPGPU-Programmierung machte deutlich, Debugger auf GPUs sollten ähnliche Fähigkeiten wie die traditionellen CPU-Debugger besitzen. Dies beinhaltet Funktionen wie die Überwachung von Variablen, Break-Points und schrittweise Ausführung. Häufig beziehen GPU-Programme Interaktionen des Benutzers in ihren Programmablauf mit ein. Die Debugger müssen solche Programme nicht zwangsweise in Echtzeit ablaufen lassen, dennoch sollte eine gewisse Interaktivität möglich sein. Weitere Anforderungen sind das ein- fache Hinzufügen und Entfernen des Debuggers bezüglich bestehender Anwendungen, den Zustand der GPU so wenig wie möglich zu beeinflussen sowie den Code direkt auf der GPU auszuführen und diese nicht etwa in Software zu emulieren. Zudem sollten die üblichen APIs und herstellerspezifischen Erweiterungen unterstützt werden.

In vielen Fälle kann eine grafische Darstellung der Daten ein besseres Gefühl über die Kor- rektheit der Berechnung vermitteln als die Ausgabe vieler Zahlen. Eine solche Visualisierung kann man als eine Art ”printf“-Debuggingauffassen,beidemdieWertedesInteressesauf den Bildschirm gezeichnet werden. Der ideale GPGPU-Debugger würde eine automatisier- te ”printf“-AusgabezusammenmitSkalierungs-undVerzerrungsfunktionenanbieten,umso auch Werte außerhalb des Darstellungsbereichs (z. B. Fließkommazahlen) abbilden zu können. Darüber hinaus stellt er zu jeder Zeit die exakten Daten als Zahlenwerte zur Verfügung. Es gibt einige Tools zum Debuggen von GPU-Programmen. Ihre wichtigsten Funktionen werden im Folgenden kurz geschildert. Bei fast jedem dieser Systeme fehlt aber mindestens eine der oben diskutierten Funktionen.

Apple OpenGL Profiler

Der Apple OpenGL Profiler [app] ist Debugger und Profiler in einem. Er stammt aus der Entwicklungsabteilung von Apple und ist ausschließlich für Grafikanwendungen mit OpenGL konzipiert. Von ihm werden neben einigen Visualisierungen auch Break-Points, Zugriffe auf den OpenGL-State sowie Stacktraces unterstützt.

gDEBugger

Beim gDEBugger [gra] handelt es sich um ein kommerzielles System von Graphic Remedy zum Debuggen von OpenGL-Anwendungen. Außer dem Setzen von Break-Points gestattet es das Beobachten der OpenGL-Variablen sowie das Modifizieren von Shader-Programmen zur Laufzeit. Darüber hinaus werden Signale der Grafikhardware ausgewertet und zu Profiling- Zwecken verwendet.

[27]

Abbildung in dieser Leseprobe nicht enthalten

ÄhnlichdemgDEBuggeristauchGLIntercept[gli]zumDebuggenvonOpenGL-An- wendungen konzipiert. Sein Funktionsumfang fällt jedoch deutlich geringer aus. Im Gegensatz zum gDEBugger ist GLIntercept nur in der Lage den OpenGL-Status zu protokollieren und Shader-Programme während der Ausführung zu manipulieren. Dafür ist es aber als freies Programm inklusive aller Quelldateien erhältlich.

Microsoft Shader Debugger

Das Shader Debugger Tool [mic] von Microsoft stellt eine DirectX-Erweiterung für das Visual Studio .NET dar. Es ist mit den von CPU-Debuggern gewohnten Funktionen ausgestattet. Dazu zählen schrittweises Debuggen, Anzeigen von Registerinhalten und Gerätestatus sowie das Setzen von Break-Points in nahezu jeder Übersetzungsschicht. Zudem lassen sich C++- und Direct3D-Quellcode simultan debuggen.

Shadesmith Fragment Program Debugger

Der Shadesmith Fragment Program Debugger [shaa] ist ein von Purcell und Sen ent- wickeltes freies System für OpenGL-Anwendungen unter Windows mit automatisiertem

”printf“-Debugging.DaessichineinemfrühenEntwicklungsstadiumbefindet,istderFunk- tionsumfang bisher recht eingeschränkt. So ist diesbezüglich nur das schrittweise Debuggen

hervorzuheben. Auf Funktionen wie dem Setzen von Break-Points, dem Zugriff auf OpenGLVariablen, etc. muss derzeit verzichtet werden. Leider fand die letzte Aktualisierung der zugehörigen Webseite im Jahr [2003] statt, sodass von einer Weiterentwicklung momentan nicht ausgegangen werden kann.

The Image Debugger

Einen anderen Ansatz verfolgt The Image Debugger [Bax] von Baxter. Dieser visualisiert Speicherinhalte im oben genannten

”printf“-Stil.DaernichtspeziellfürShader-Programme entwickelt wurde, muss der Programmierer seine Shader so bauen, dass die Ausgaben in den entsprechenden Ausgabepuffer geschrieben werden. The Image Debugger ist als freie Software für Windows verfügbar.

Weitere Tools

Als kleinere Tools sind noch BuGLe [bug] und GLTrace [haw] zu erwähnen. Während es sich bei BuGLe um einen OpenGL-Debugger für Linux handelt, dient GLTrace lediglich dem Nachvollziehen von OpenGL-Aufrufen.

Des Weiteren ist zu berücksichtigen, dass in größeren GPGPU-Programmiersystemen wie beispielsweise dem CUDA-Toolkit von NVIDIA [nvia] oft eigene Werkzeuge zum Debugging und Profiling enthalten sind.

3.3 Einführung in die RapidMind-Entwicklungsplattform

Für die im Rahmen dieser Arbeit durchgeführte Implementierung des Radiosity-Verfahrens wurde die Entwicklungsplattform von RapidMind verwendet. Als Gründe für diese Wahl sind vor allem die kostenfrei verfügbare Developer-Version, die einfach zu verwendende und übersichtliche API in Standard-C++ sowie der recht gute Funktionsumfang in Bezug auf die GPGPU-Programmierung zu nennen. Hinzu kommt, dass sich die Plattform nicht nur auf einen Grafikkartenhersteller beschränkt und durch das allgemeingültige Programmiermo- dell, mit der Option weitere Hardwarearchitekturen nutzen zu können, das Interesse weckt. Aufgrund seiner Entwicklungsgeschichte ist zudem ein relativ ausgereifter Zustand dieses Programmiersystems zu erwarten.

Die RapidMind-Plattform, vgl. [McC06], ist eine Weiterentwicklung des früheren Systems Sh, welches die Grundlage für das innovative Interface zur Metaprogrammierung bildet und diese Technologie bereits etablierte. Sh wurde an der Universität von Waterloo entwickelt und verfolgte primär das Ziel als eine Art höhere Shading-Language für GPUs zu fungieren. Im Gegensatz dazu wurde die kommerziell orientierte RapidMind-Plattform allgemeingültiger entworfen, um so neben der GPU auch andere Hardwarearchitekturen zu unterstützen. Ihr Programmiermodell fällt deshalb ausgereifter und in sich geschlossener aus als das von Sh. Bemerkbar macht sich dies vor allem an der gestiegenen Bedeutung von Arrays sowie den hinzugekommenen kollektiven Operationen auf den Arrays. Bevor das Programmiermodell mit dem zugehörigen Interface vorgestellt wird, folgen zunächst einige bedeutende Eigenheiten der RapidMind-Plattform:

1. Anwendungen werden in C++ entwickelt, was das Erlernen einer zusätzlichen neuen

Programmiersprache erspart. Das Interface verwendet dabei nur C++-Funktionen des ISO-Standards und ist somit in der Lage mit jedem Standard-C++-Compiler zu arbeiten. Dazu muss lediglich eine Header-Datei eingebunden und auf eine Bibliothek verlinkt werden. Obwohl als Bibliothek implementiert, erlaubt das System aber mehr als die schlichte Anwendung vorgefertigter Funktionen. Die Plattform beinhaltet einen dynamischen Compiler, der Codegenerierung zur Laufzeit ermöglicht. So verhält sich das Ausführungsmodell trotz C++-Interface eher wie FORTRAN und erlaubt es, den vom System generierten Code stark zu optimieren.

2. Neben dem genannten Interface und der dynamischen Compilerkomponente verfügt die

RapidMind-Plattform über eine umfangreiche Laufzeitumgebung. Diese automatisiert Aufgaben wie die Umsetzung von Task-Warteschlangen, Streaming und Transfer von Daten, Synchronisation und Load-Balancing. Sie verwaltet die Ausführung von Tasks auf mehreren Prozessoren sowie den Datenaustausch mit verteiltem Speicher. Da so die Details der parallelen Ausführung verdeckt ablaufen, kann das Interface einfach und übersichtlich bleiben.

3. Das Programmierinterface ist sicher. Das bedeutet, mit der RapidMind-Plattform ge-

schriebene Programme können keine Deadlocks, Read-Write-Hazards oder Synchronisa- tionsprobleme verursachen. Bereits die Sprachstruktur spiegelt die Parallelität wieder und ermutigt so zur Entwicklung und Verwendung effizienter, skalierbarer paralleler Algorithmen.

4. Das Programmiermodell der RapidMind-Plattform ist, im Gegensatz zu Alternativen

wie OpenMP, MPI und Threads, auf eine Vielzahl paralleler Hardwarearchitekturen portierbar. Dies schließt sowohl Vektor- und Streaming-Rechner wie GPUs als auch Maschinen mit verteiltem Speicher wie Cell-Prozessoren ein. Zudem verfügt das System bei Ausführung und Datenhaltung über einen hohen Abstraktionsgrad, der zugleich modular, portabel und effizient ist.

Da sich das Programmierinterface der Plattform sehr nah am datenparallelen Program- miermodell orientiert, ist es günstig beide gemeinsam vorzustellen. Die Basis des Interface bilden drei C++-Haupttypen: Value<N,T>, Array<D,T> und Program. Alle drei Typen sind Container, die ersten beiden für Daten und der letztere für Operationen. Auf den ersten Blick erscheinen der Value- und der Array-Typ nicht besonders interessant, da vermutlich nahezu jeder C++-Programmierer schon einmal n-Tupel und Arrays verwendet hat. Aber der Schein trügt, wurden doch die Semantiken dieser Typen so gewählt, dass ihr Gebrauch intuitiv und leicht erlernbar ist. Zwei Wege führen zur parallelen Berechnung: entweder man wendet ein Program-Objekt auf Array-Objekte an oder aber eine parallele kollektive Operation (z. B. eine Reduktion) auf ein Program-Objekt. Zusammen formen die drei Haupttypen ein Interface, das genug Freiheiten für die oben genannten Eigenschaften der Plattform lässt. Im Mittelpunkt steht dabei der neuartige Program-Typ. Sind seine Operationen einmal festgelegt, können diese direkt vom dynamischen Compiler der Plattform manipuliert und optimiert werden.

3.3.1 Value-Typ

Beim Value<N,T>-Typ handelt es sich um ein N-Tupel, dessen Instanzen N Variablen vom Typ T halten. Dabei kann T einer der numerischen Grundtypen sein. Unterstützt werden neben den üblichen Typen wie beispielsweise Fließkommazahlen einfacher und doppelter Genauigkeit sowie Integer-Werten auch seltenere Typen wie Fließkommazahlen halber Genauigkeit mit 10-Bit-Mantisse und 5-Bit-Exponent. Für Tupel bis zur Länge vier gibt es eine Reihe von Kurzformdeklarationen auf der Basis eines Suffixsystems. Zum Beispiel kann ein 4-Tupel von Fließkommazahlen einfacher Genauigkeit als Value4f spezifiziert werden und ein 3-Tupel vorzeichenloser Byte-Werte als Value3ub. Auch Boolean-Werte können als Komponententyp verwendet werden. Sie nutzen das Suffix bool. Dabei muss immer beachtet werden, dass die Zielarchitektur wie z. B. die GPU möglicherweise nicht alle Datentypen unterstützt und diese deshalb, wenn möglich, vom RapidMind-Compiler automatisch auf einen anderen Datentyp abgebildet werden.

Sämtliche Operatoren wurden für die Value-Typen überladen und arbeiten nun kompo- nentenweise. Dazu zählen neben den Operatoren der Grundrechenarten auch Vergleichs- und logische Operatoren. Trigonometrische, logarithmische und exponentielle Standardfunktionen stehen ebenfalls zur Verfügung. Sie haben, wann immer möglich, dieselben Namen wie in der Mathematikbibliothek von C erhalten und arbeiten nun ebenfalls komponentenweise. Da der Value1f-Typ normalerweise den Typ float ersetzt, lassen sich damit beispielsweise komplexe Zahlen unter Verwendung der entsprechenden Standard-Template-Bibliothek implementieren: std::complex<Value1f>.

Zusätzlich zu den komponentenweisen Operationen sind auf Value-Instanzen das soge- nannte Swizzling sowie die Verwendung von Schreibmasken erlaubt. Das Swizzling ermöglicht die willkürliche Umstellung und Vervielfältigung der Komponenten. Man drückt dies mittels wiederholten Array-Zugriffen und dem

”()“-Operatoraus.DieeinzelnenKomponentensind dabei von Null beginnend durchnummeriert. Präsentiert beispielsweise a eine RGBA-Farbe in Form eines Value<4,float>, so steht der Ausdruck a(2,1,0) für ein 3-Tupel der ersten drei Komponenten in umgekehrter Reihenfolge (also einen BGR-Farbwert). Wünscht man statt- dessen eine dreifache Wiederholung der G-Komponente gefolgt von der A-Komponente, so kann dafür der Ausdruck a(1,1,1,3) geschrieben werden. Erscheint ein solcher Ausdruck auf der linken Seite einer Zuweisung, so spricht man von einer Schreibmaske. Diese ist eigentlich nichts anderes als ein Array von Zeigern auf die jeweiligen Komponenten. In Schreibmasken darf jede der Zielkomponenten nur einmal vorkommen, da ihr sonst mehrere Werte zugewiesen würden.

Des Weiteren können gebräuchliche C++-Konstrukte wie Namensräume, Klassen und Funktionen frei verwendet werden. Zum Beispiel berechnet die folgende Funktion die Reflexionsrichtung zu einem Vektor an einer Ebene im Raum, wobei die Ebene durch ihren Normalenvektor gegeben ist:

Abbildung in dieser Leseprobe nicht enthalten

3.3.2 Array-Typ

Der Array<D,T>-Typ ist ebenfalls ein Datencontainer, allerdings ein mehrdimensionaler. Die Dimensionalität wird durch D angegeben, welches vom Typ Integer ist und die Werte 1, 2 oder 3 annehmen kann. T beschreibt den Typ der Elemente des Arrays und ist momentan auf Instanzen vom Typ Value<N,T> beschränkt. Die Anzahl der Elemente im Array ist dabei nicht Bestandteil des Typs.

Seine Instanzen unterstützen die Operatoren

”[]“und ”()“fürbeliebigeZugriffeauf

ihre Elemente. Dabei dienen Value-Tupel der entsprechenden Dimensionalität als Argumente

beider Operatoren. Während der

”[]“-OperatorInteger-basierteKoordinatenverarbeitet,

die von Null beginnen, sind die Koordinaten des

”()“-OperatorsReal-WerteimIntervall

[[0], [1]] je Dimension. Dies kann zum Beispiel in Verbindung mit Interpolationsverfahren sehr nützlich sein.

Der Zugriff auf Teil-Arrays kann über die Funktionen slice, offset und stride erfol- gen. Zum Festlegen des Verhaltens am Array-Rand verwendet man die Member-Funktion boundary. Diese beinhaltet Modi zum Einschränken der Zugriffe auf eine bestimmte Kante des Arrays sowie zur periodischen Wiederholung von Array-Inhalten oder der Rückgabe ei- nes konstanten Wertes außerhalb der Array-Grenzen. Interpolationsmodi können auf ähnliche Weise festgelegt werden.

Zwei weitere Array-ähnliche Typen sind ArrayReferences und ArrayAccessors. ArrayRefe- rences bieten Zeigersemantik und sind nützlich, wenn beispielsweise das Ziel einer Berechnung im Speicher vom Benutzer gesteuert wird. ArrayAccessors dienen hingegen als Zeiger auf Teil- Arrays. Darüber hinaus werden auch virtuelle Arrays unterstützt. Ein virtuelles Array wird prozedural generiert und belegt keinen Speicher, da seine Werte on-the-fly berechnet werden.

3.3.3 Program-Typ

Der letzte grundlegende Typ ist Program. Ein Objekt dieses Typs speichert eine Abfolge von Operationen und wird durch einen Wechsel in den aufbewahrenden Modus ( ”retainedmode“)

erzeugt. Das Makro BEGIN kennzeichnet einen solchen Wechsel. Normalerweise arbeitet das System im fortlaufenden Modus (

”immediatemode“).DabeiwerdenOperationenaufValue-

Tupeln durchgeführt, als wären sie in einer normalen Matrix-Vektor-Bibliothek spezifiziert: die Berechnung findet auf dem gleichen Prozessor wie das ausführende Programm statt und das numerische Ergebnis wird dem angegebenen Value-Tupel zugewiesen. Durch Erzeugen ei- nes neuen Program-Objekts mittels BEGIN wird der aufbewahrende Modus betreten. Innerhalb dieses Modus werden Operationen auf Value-Tupeln nicht direkt ausgeführt, sondern statt- dessen symbolisch berechnet und im Program-Objekt gespeichert. Der aufbewahrende Modus wird mit dem Makro END wieder verlassen. Damit wird das Program-Objekt geschlossen und als bereit zum Kompilieren markiert. Geschlossene Program-Objekte können berechnet werden. Sie werden mithilfe des dynamischen Compilers kompiliert und können so die Vorteile der Zielhardware genießen.

An dieser Stelle sei nochmals darauf hingewiesen, dass das Ausführungsmodell der RapidMind-Plattform eher FORTRAN als C++ ähnelt. Die objektorientierten Funktionen von C++ können verwendet werden, um die Berechnung zu strukturieren. Im finalen von der Plattform generierten Code erscheinen sie jedoch nicht. C++ liefert zwar die Grundlage zur Codegenerierung, wird aber selbst nicht für die Ausführung verwendet. Hinzu kommt, dass nur Operationen, die Value- und Array-Typen betreffen und im aufbewahrenden Modus gespeichert wurden, hardwarenah berechnet werden können.

Es folgt ein Beispiel, welches die oben definierte reflect-Funktion sowie die vorgefertigte normalize-Funktion aus der RapidMind-Bibliothek verwendet:

Abbildung in dieser Leseprobe nicht enthalten

Dieses Beispiel demonstriert weitere Konzepte. So findet die Deklaration von Ein- und Ausga- bewerten mithilfe der Template-Wrapper In<T> und Out<T> statt. Im Beispiel betrifft dies die Tupel v und r. Bevor ein Program-Objekt ausgeführt wird, werden die als Input markierten Werte mit aktuellen Daten initialisiert. Die resultierenden Werte, der als Output markier- ten Variablen, werden am Ende vom Program-Objekt zurückgegeben. Man beachte, dass die nicht-lokale Variable n innerhalb der Definition des Program-Objekts verwendet wird. Eine solche Variable wird uniform genannt. Bei der späteren Ausführung des Program-Objekts p wird der zu diesem Zeitpunkt gültige Wert von n in die Berechnung eingehen, wie auch immer er dann lautet und auch, wenn die Berechnung auf einem anderen Prozessor oder gar parallel auf mehreren Prozessoren stattfindet. Die RapidMind-Plattform wird dafür Sorge tragen, dass allen Program-Objekten ihre jeweils benötigten Daten zur Verfügung stehen und gegebenenfalls geänderte Werte entsprechend verteilen. Des Weiteren kann man feststellen, dass die normalize-Operation auf eine nicht-lokale Variable angewendet wird. Eine mögli- che automatische Optimierung des dynamischen Compilers könnte sein, einen solchen Aufruf aus der parallelen Berechnung herauszulösen und diesen nur einmal durchzuführen bzw. nur dann, wenn die zugehörige uniforme Variable geändert wird.

Ein einmal erstelltes Program-Objekt wird ausgeführt, indem ihm entweder Value-Tupel oder Array-Objekte übergeben werden. Im Folgenden seien V und R zweidimensionale Arrays (jede andere Dimensionalität wäre auch denkbar, solange sie bei beiden übereinstimmt). Das nachstehende Codefragment wendet die im Program-Objekt p spezifizierten Operationen auf alle Elemente von V an und erzeugt anschließend ein neues Array R der gleichen Größe:

Abbildung in dieser Leseprobe nicht enthalten

3.3 Einführung in die RapidMind-Entwicklungsplattform

Aus diesem Beispiel resultieren eine Million Instanzen der in p spezifizierten Berechnung. Die Initialisierung von V mit Daten könnte sowohl das Ergebnis einer früheren Berechnung sein oder aber auch aus dem Speicher geladen werden. Komplexere Algorithmen benötigen womöglich mehrere solcher Aufrufe. Die Daten können dabei so lange in den Array-Objekten verbleiben, bis die gesamte Berechnung abgeschlossen ist. Es ist auch nicht notwendig für R eine Größe anzugeben, da sich diese automatisch aus der Größe von V ergibt. Während der Berechnung werden alle zuvor in R gespeicherten Daten ersetzt oder gehen verloren. Die Größe eines Arrays ist flexibel und ergibt sich stets aus seinen Werten.

Program-Objekte unterliegen einigen Beschränkungen. Zum Beispiel dürfen sie zwar auf beliebige im fortlaufenden Modus deklarierte Value- und Array-Objekte lesend zugreifen, ihnen ist es aber nicht gestattet, Werte von Variablen zu ändern, die nicht innerhalb ihrer Definition deklariert wurden. Um genauer zu sein, werden solche Modifikationen nicht nach außen sichtbar. Die einzigen Ausgaben eines Program-Objekts sind die explizit mit Out<T> markierten Werte. So fungieren Program-Objekte nach außen wie echte Funktionen. Die aus der Anwendung einer solchen Funktion auf eine Array entstehenden Instanzen einer Berechnung können in beliebiger Reihenfolge, und damit auch parallel, abgearbeitet werden. Da alle resultierenden Werte erst nach der Berechnung freigegeben werden, ist es zudem möglich ein Array gleichzeitig als Input und als Output zu nutzen. Das System verwendet bei Bedarf einen Double-Buffer, um Read-Write-Hazards zu vermeiden.

Ferner ist es möglich, innerhalb des aufbewahrenden Modus Steuerstrukturen wie Schleifen und Verzweigungen zu benutzen und sogar weitere Program-Objekte als Unterprogramme aufzurufen. Dies führt dazu, dass Program-Objekte nicht nach einem reinen SIMD-Modell ausgeführt werden, sondern nach dem allgemeineren SPMD-Modell (vgl. Abschnitt 3.1.4). Herkömmliche C++-Steuerstrukturen können zwar die Codegenerierung von Program-Ob- jekten beeinflussen, haben aber keine direkte Auswirkung auf den Programmablauf während der eigentlichen Berechnung. Zur Steuerung des Ablaufs in Program-Objekten existieren Ma- kros, die entsprechende Steuerstrukturen dynamisch umsetzen. Bei diesen Makros handelt es sich, genauer betrachtet, um Wrapper für einige hardwarespezifische Funktionen, die ein- geführt wurden, um die Syntax für den Programmierer zu vereinfachen. Das folgende Beispiel zeigt ein einfaches Partikelsystem, welches Partikel entlang eines ballistischen Pfades bewegt, bis diese ihre maximale Lebenszeit erreicht haben. Danach werden sie auf Null zurück gesetzt und ihre neue Geschwindigkeit ergibt sich als Hashwert über ihre letzte Position:

Abbildung in dieser Leseprobe nicht enthalten

Zusätzlich zu IF/ELSEIF/ELSE/ENDIF werden auch Schleifen mit FOR/ENDFOR, DO/UNTIL und WHILE/ENDWHILE vom System unterstützt. Ihre Argumente können beliebigen Quellcode enthalten, auch Steuerelemente selbst sind erlaubt. Die InOut<T>-Markierung kennzeichnet Variablen, die sowohl Eingabe- als auch Ausgabewerte sind. Dies bedeutet jedoch nicht, dass die an das Program-Objekt übergebenen Eingabewerte modifiziert werden. Variablen, die eine Eingabe darstellen, werden immer als Werte und nie als Zeiger übergeben. Ein- und AusgabeSlots sind also stets getrennt zu betrachten.

Zudem besitzt das Program-Objekt des obigen Beispiels mehrere Rückgabewerte. Tatsäch- lich geben Program-Objekte immer Instanzen des Bundle-Typs zurück. Dieser stellt einen Container für ein oder mehrere Arrays bzw. Tupel dar. Besteht ein Bundle-Objekt aus nur einem Element, so kann es einem Array bzw. Tupel direkt zugewiesen werden, wobei automa- tisch die entsprechende Konvertierung stattfindet. Program-Objekte, die mehrere Rückgabe- werte besitzen, müssen Bundle-Objekten zugewiesen werden, um ihre Ergebnisse auf mehrere Ziele verteilen zu können. Hierzu vereinfacht die bundle-Funktion das Erstellen temporärer Bundle-Objekte. Wie das folgende Beispiel zeigt, brauchen Bundle-Instanzen nicht ausdrück- lich deklariert werden:

Array<1,Value3f> Vel, Pos; Array<1,Value1f> Time;

bundle(Vel, Pos, Time) = evolve(Vel, Pos, Time);

Das Beispiel berechnet einen Schritt des obigen Partikelsystems und demonstriert dabei die Semantik paralleler Zuweisungen. Auch wenn wie hier die Ausgabe- und Eingabe-Arrays dieselben sind, werden die neuen Werte erst zugewiesen, nachdem alle alten Werte verarbeitet wurden.

3.3.4 Teilweise Berechnung, Programmalgebra und kollektive Operationen

Neben den bereits genannten Funktionen existieren weitere wichtige Besonderheiten bei der Verwendung der RapidMind-Plattform. So können Program-Objekte auch teilweise berechnet werden. Zum Beispiel muss der Ausdruck p(A) nicht alle Elemente von A auf einmal verarbeiten, wobei p vom Typ Program und A ein Array sei. Es besteht also die Möglichkeit, Eingabewerte an ein Program-Objekt zu binden, die eigentliche Berechnung aber verzögert durchzuführen. Die Ausführung eines Program-Objekts wird immer erst dann ausgelöst, wenn es einem Bundle, Array oder Value zugewiesen wird.

Bisher wurde ausschließlich der

”()“-Operatorverwendet,umWerteaneinProgram-Ob-

jekt zu binden. Dies wird feste Bindung genannt. Änderungen an Eingabewerten, die nach dem Eingehen einer festen Bindung stattfinden, haben keinen Einfluss auf die Berechnung, auch wenn diese verzögert wird. Erzeugt man im genannten Beispiel p(A), weist es einem weiteren Program-Objekt q zu und verändert anschließend die Werte von A, so gehen bei einer späteren Ausführung von q trotzdem die zur Zeit des Bindens gültigen Werte ein und nicht die geänderten. Feste Bindungen ermöglichen die auf Eingabewerten basierende Opti- mierung von Program-Objekten. Neben festen Bindungen unterstützt das System auch lose

Bindungen. Diese werden mit dem

”<<“-Operatorgekennzeichnet.DerAusdruckp<<Aist

demnach ähnlich zu p(A), mit dem Unterschied, dass hierÄnderungen an A bei einer späteren Ausführung von p berücksichtigt werden.

An ein Program-Objekt können sowohl Array- als auch Value-Werte gebunden werden. Bei mehreren Arrays verschiedener Größen werden die kleineren Arrays entsprechend ihrer

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 3.3: Berechnungsmuster im Vergleich. Quelle: [McC06].

Randbedingungen aufgefüllt und so an das größte Array angeglichen. Beim Binden von Va- riablen beider Typen werden die Tupel automatisch vervielfältigt, um zu jedem Element

des Arrays zu passen. Lose Bindungen sind invertierbar - der Operator

”>>“kanneinelose

Bindung rückgängig machen. Angenommen ein Program-Objekt p bezieht sich auf eine nichtlokale Variable a, so erzeugt der Ausdruck q = p >> a ein neues Program-Objekt q mit einen zusätzlichen Eingabeparameter. Dieser ersetzt intern alle Referenzen auf a.

Ein weiteres Feature ist das Kombinieren mehrerer Program-Objekte zu einem neuen Program-Objekt. Hierzu existieren zwei Funktionen: die funktionale Zusammensetzung und die Bündelung. Beide Operationen formen zusammen eine über Program-Objekten abge-

schlossene Algebra. Als Operator der funktionalen Zusammensetzung wird

”<<“verwendet.

Seine Anwendung bewirkt, dass alle Ausgabewerte des Objekts zu seiner Rechten die Einga- ben des Objekts zu seiner Linken bilden. p << q erzeugt also eine neues Program-Objekt mit den Eingaben von q und den Ausgaben von p. Dagegen verknüpft die Bündelungsoperation jeweils alle Eingaben und alle Ausgaben der Objekte untereinander. Dies geschieht mithilfe der bundle-Funktion. Sie erzeugt ein neues Program-Objekt, das äquivalent zur Verkettung der Quelltexte ihrer Eingabeprogramme in der angegebenen Reihenfolge ist. Programmalge- brafunktionen sind zum Beispiel nützlich, um die Stufen einer Pipeline in einen einzelnen, effizienteren Durchlauf zu packen.

Außer der parallelen Berechnung der Elemente eines Arrays stehen weitere Berechnungsund Kommunikationsmuster in Form von kollektiven Operationen zur Verfügung. Vollständige Kommunikationsmuster werden z. B. durch die Scatter- und die Gather-Operation bereit gestellt (wenn von der Hardware unterstützt). Ein hierarchisches Berechnungsmuster bietet beispielsweise die Reduktionsoperation. Abbildung 3.3 zeigt diese im Vergleich mit der allgemeinen parallelen Berechnung. Die Kombination aus kollektiven und parallel angewandten Operationen gestaltet das Programmiermodell sehr allgemeingültig und erlaubt so die effiziente Implementierung einer Vielzahl verschiedenster Algorithmen.

Im Folgenden soll die Reduktionsfunktion (reduce) näher betrachtet werden. Sie kombi- niert Elemente eines Arrays durch die Verwendung eines assoziativen Operators mit zwei Eingabewerten und einem Ausgabewert, die alle vom gleichen Typ sind. Zum Beispiel würde eine auf dem Summenoperator basierende Reduktion die Summe aller Elemente eines Arrays berechnen und ein einzelnes Element zurückgeben. Eine allgemeingültige Form der reduce- Funktion verwendet ein Program und ein Array als Argumente. Das Program-Objekt sollte zwei Eingabe- und einen Ausgabewert gleichen Typs haben. Eine solche Reduktion liefert eine einzelne Instanz vom Elementtyp des Arrays als Ergebnis. Ist beispielsweise

Abbildung in dieser Leseprobe nicht enthalten

Die combiner-Funktion darf beliebigen Code enthalten, wobei sich das System aber nicht darum kümmert, ob dieser auch assoziativ operiert. Ist dies nicht der Fall, so hängt das Ergebnis von der zugrunde liegenden Hardware ab. Selbst auf der gleichen Hardware können in verschiedenen Durchläufen unterschiedliche Ergebnisse heraus kommen. Für ein deterministisches Verhalten muss unbedingt ein assoziativer Operator verwendet werden. Einige gebräuchliche Reduktionen sind bereits vorgefertigt: sum, product, min und max werden unterstützt. Hinzu kommen weitere kollektive Operationen. Zum Beispiel ist es möglich, die Anzahl an Elementen in einem Array zu zählen, die eine bestimmte Eigenschaft erfüllen. Für die Zukunft sind noch mehr solcher Operationen geplant.

3.3.5 Parallelisierung eines Beispiels

Um die Verwendung des Systems zu verdeutlichen, folgt nun ein Beispiel, wie einfach sich eine Verschachtelung von Schleifen parallelisieren lässt. Angenommen eine Anwendung beinhaltet folgenden Code:

Abbildung in dieser Leseprobe nicht enthalten

Diese Berechnung kann unter Verwendung der RapidMind-Plattform wie folgt ausgedrückt werden:

3.4 Zusammenfassung

Abbildung in dieser Leseprobe nicht enthalten

Program-Objekte werden dynamisch definiert und kompiliert. Obwohl dies nicht besonders langsam ist, sollten sie dennoch nicht innerhalb einer Schleife platziert oder unnötig wiederholt werden. Das obige Codefragment kann deshalb effizienter gestaltet werden, indem das Program-Objekt als statisch markiert und damit nur beim ersten Aufruf der Funktion func_arrays erstellt wird. Eine alternative Struktur, um Mehrfachdefinitionen von ProgramObjekten zu vermeiden, ist die Verlagerung der Program-Definition in den Konstruktor einer Klasse sowie die Verwendung statischer Instanzen dieser Klasse. Das kompilierte ProgramObjekt kann so als Member-Variable im entsprechenden Objekt der Klasse gespeichert und auf die gewünschten Daten angewendet werden.

Um die durch den dynamischen Compiler der Plattform geleistete Arbeit zu verdeutlichen, ist der zu diesem Beispiel generierte GLSL-Quellcode in Anhang C abgebildet.

3.4 Zusammenfassung

Grafikkarten bieten viel Leistung zu einem vergleichsweise günstigen Preis und sind in nahe- zu jedem handelsüblichen PC zu finden. Hinzu kommt eine steigende Anzahl an Systemen und Tools, welche die Programmierung moderner Grafikhardware stark vereinfachen. Ihre Programmiermodelle sind teilweise so sehr abstrahiert, dass viele grafikspezifische Details im Verborgenen bleiben. Dies macht die Verwendung der GPU auch für allgemeine Berechnungen zugänglich. Aufgrund vieler Beschränkungen (vgl. auch mit Abschnitt 6.2) und der speziellen Hardwarearchitektur sind allerdings einige neue Denkweisen bei der Umsetzung bekannter Algorithmen notwendig. Sich dieser anzunehmen, ist durchaus eine lohnenswerte Investition, sind doch wegen der datenparallelen Berechnung in vielen Fällen enorme Leistungssteigerun- gen gegenüber der traditionellen CPU-Programmierung zu erwarten.

4 Mathematische Hilfsmittel

Zur Umsetzung des Radiosity-Verfahrens sind Kenntnisse über einige mathematische Tech- niken und Algorithmen notwendig. Auf die bei der Implementierung verwendeten Methoden soll deshalb an dieser Stelle genauer eingegangen werden. So beschäftigt sich dieses Kapi- tel zunächst mit dem iterativen Lösen linearer Gleichungssysteme, da diese Verfahren später zur Berechnung der Strahlungswerte genutzt werden. Anschließend wird ein Schnitttest vor- gestellt, der prüft, ob eine Strecke ein Dreieck schneidet. Dieser ist für die Formfaktoren- bestimmung mithilfe der Kreisscheibenapproximation von Bedeutung. Zuletzt wird die zur Verringerung der Anzahl notwendiger Schnitttests verwendete Octree-Datenstruktur näher betrachtet.

4.1 Iteratives Lösen linearer Gleichungssysteme

Laut Goral et al. [GTGB84] kann für die Ermittlung der Strahlungswerte des Radiosity- Gleichungssystems (2.5) jedes Standardverfahren zum Lösen linearer Gleichungssysteme ein- gesetzt werden. Da direkte Verfahren zu aufwendig sind, ihr Aufwand wächst zumeist kubisch mit der Anzahl der Gleichungen, nutzt man stattdessen iterative Verfahren und rechnet bis zum Erreichen der gewünschten Genauigkeit. Goral et al. verwenden hierfür das Gauß-Seidel- Verfahren. Da die Instanzen eines Fragment-Programms aufgrund der in Kapitel 3 genannten Einschränkungen nicht in der Lage sind, untereinander zu kommunizieren, lässt sich der Gauß- Seidel-Algorithmus auf der GPU nicht parallelisieren. An dieser Stelle soll deshalb auf das Jacobi-Verfahren zurückgegriffen werden, welches eine solche Kommunikation nicht benötigt.

Dieser Abschnitt widmet sich folglich dem Gauß-Seidel- und dem Jacobi-Verfahren sowie einer möglichen Parallelisierung des Letzteren. Zuvor soll jedoch auf die klassischen Iterationsverfahren im Allgemeinen eingegangen werden. Bezüglich der verwendeten Literatur sei auf [RR00] verwiesen. Weiterführende Informationen zu den Verfahren, ihren Herleitungen und Konvergenzeigenschaften findet man in [DH93] und [Wer92].

4.1.1 Beschreibung klassischer Iterationsverfahren

Iterative Verfahren zur Lösung eines linearen Gleichungssystems der Form Ax = b mit gege- bener (n×n)-Matrix A ∈ Rn×n und gegebenem Vektor b ∈ Rn berechnen die gesuchte Lösung näherungsweise, indem sie eine Folge von Approximationsvektoren {x(k)}k=1,2,... bilden, die gegen die exakte Lösung x∗ konvergiert. Zur Berechnung eines Approximationsvektors wird jeweils eine Matrix-Vektor-Multiplikation mit einer speziellen Iterationsmatrix durchgeführt. Dabei geht die Matrix A des Gleichungssystems nur über die verwendete Iterationsmatrix in die Berechnung ein. Als ein wesentliches Kriterium zur Bewertung von Iterationsverfahren gilt die Konvergenzgeschwindigkeit.

Klassische Iterationsverfahren beruhen auf einer Zerlegung der gegebenen Matrix A in

(4.1)

Abbildung in dieser Leseprobe nicht enthalten

4 Mathematische Hilfsmittel

mit M, N ∈ Rn×n. Dabei bezeichnet M eine leicht invertierbare, nichtsinguläre Matrix, z. B. eine Diagonalmatrix. Für die gesuchte Lösung x∗ des obigen Gleichungssystems gilt dann

(4.2)

Abbildung in dieser Leseprobe nicht enthalten

woraus sich die Iterationsvorschrift M x(k+[1]) = N x(k) + b für k = 0, 1, . . . bilden lässt. Diese schreibt man typischerweise in der Form

(4.3)

Abbildung in dieser Leseprobe nicht enthalten

wobei C := M−[1]N und d := M−[1]b. Zudem ist die oben genannte Matrix-Vektor-Multi- plikation mit der Iterationsmatrix C in dieser Darstellung gut erkennbar. Gilt unabhängig von der Wahl des Startvektors x([0]) ∈ Rn, dass limk→∞ x(k) = x∗, so heißt das Verfahren konvergent gegen x∗. Liegt Konvergenz vor, so existiert genau ein Vektor x∗ für den x∗ = Cx∗ + d gilt. Durch Subtraktion dieser Gleichung von Gleichung (4.3) sowie Induktion folgt x(k) − x∗ = Ck(x([0]) − x∗), wobei Ck die k-fache Multiplikation der Matrix C bezeichnet. Damit ist die Konvergenz von Iterationsvorschrift (4.3) äquivalent zu

(4.4)

Abbildung in dieser Leseprobe nicht enthalten

4.1.2 Jacobi-Verfahren

Das Jacobi-Verfahren (oder Gesamtschrittverfahren) wird durch die Zerlegung der Matrix A in A = D − L − R mit D, L, R ∈ Rn×n beschrieben. Diesbezüglich ist bei D die Diagonale, bei L die untere und bei R die obere Dreiecksmatrix (jeweils ohne Diagonale) mit den entsprechenden Elementen von A besetzt. Alle anderen Einträge sind jeweils Null. Ein Iterationsschritt nutzt diese Zerlegung in der Form

(4.5)

Abbildung in dieser Leseprobe nicht enthalten

und führt damit zu einer Iterationsmatrix CJ := D−[1](L + R), was

Abbildung in dieser Leseprobe nicht enthalten

(4.6)

entspricht. Während diese Schreibweise in Matrizenform eher für den Nachweis der Konvergenz genutzt wird, ist die folgende Komponentenschreibweise eines Iterationsschritts für die praktische Anwendung günstiger:

(4.7)

Abbildung in dieser Leseprobe nicht enthalten

Bei der Anwendung dieses Verfahrens auf das Radiosity-Gleichungssystem vereinfacht sich ein solcher Berechnungsschritt zu

Abbildung in dieser Leseprobe nicht enthalten da wegen Fii = 0 (vgl. Abschnitt 2.3) alle Diagonalelemente der Matrix gleich Eins sind.

4.1.3 Gauß-Seidel-Verfahren

Das Gauß-Seidel-Verfahren (oder Einzelschrittverfahren) beruht auf der gleichen Zerlegung der Matrix A wie das Jacobi-Verfahren: A = D − L − R. Allerdings nutzt es im Unterschied zu diesem die Form

(4.9)

Abbildung in dieser Leseprobe nicht enthalten

für einen Iterationsschritt. Damit erhält man die Iterationsmatrix CG := (D − L)−[1]R. Für einen Iterationsschritt in Komponentenschreibweise ergibt sich entsprechend

(4.10)

Abbildung in dieser Leseprobe nicht enthalten

Wie beim Jacobi-Verfahren kann auch hier ein solcher Berechnungsschritt bei Anwendung auf das Radiosity-Gleichungssystem vereinfacht werden, da weiterhin aii = 1 für alle Diagonalelemente gilt:

Abbildung in dieser Leseprobe nicht enthalten

Im Gegensatz zum Jacobi-Verfahren gehen in die Berechnung von x(k+[1])i bereits die neu ermittelten Werte x(k+[1])[1] ,...,x(k+[1])i−[1] des aktuellen Iterationsschritts mit ein. Man erhofft sich dadurch eine schnellere Konvergenz, da diese neuen Werte im Allgemeinen genauer als die der vorhergehenden Iteration sind.

4.1.4 Konvergenzkriterium

Häufig wird zum Nachweis der Konvergenz des Jacobi- und des Gauß-Seidel-Verfahrens auf bestimmte Eigenschaften der Matrix A zurückgegriffen, da diese leichter überprüfbar sind. Eine hinreichende Konvergenzbedingung bildet beispielsweise das starke Zeilensummenkriterium. Es besagt, dass eine (n × n)-Matrix A die Ungleichung

(4.12)

Abbildung in dieser Leseprobe nicht enthalten

erfüllt, also A stark diagonaldominant ist. Trifft dies auf die gegebene Matrix des Gleichungssystems zu, so konvergieren sowohl das Jacobi- als auch das Gauß-Seidel-Verfahren. Einen Beweis dieses Zusammenhangs findet man in [Wer92].

Die Matrix des Radiosity-Gleichungssystems (2.5) erfüllt diese Bedingung. Wie bereits weiter oben festgestellt wurde, sind ihre Diagonalelemente gleich Eins. Die Summe der Form- faktoren einer Zeile beträgt nach Abschnitt 2.3 ebenfalls Eins. Zusätzlich wird jeder der Formfaktoren mit einem Reflexionskoeffizienten multipliziert. Dieser liegt im Intervall [0, 1], sodass die Summe über die übrigen Zeilenelemente in der Regel kleiner als Eins ist. Eine Ausnahme liegt vor, wenn alle Reflexionskoeffizienten einer Szene den Wert Eins annehmen, da dann der Betrag des Diagonalelements mit der Summe der übrigen Zeilenelemente über- einstimmt. In diesem Fall kann keine Aussage über die Konvergenz getroffen werden. Für die Praxis hat dies allerdings eine geringe Bedeutung, da solche Szenen sehr selten sind.

Im Allgemeinen führt eine größere Dominanz der Diagonalelemente häufig zu einer besseren Konvergenz. Zudem kann in vielen praktischen Fällen beobachtet werden, dass das GaußSeidel-Verfahren schneller konvergiert als das Jacobi-Verfahren.

4.1.5 Parallele Realisierung des Jacobi-Verfahrens

Gemäß der Berechnungsvorschrift (4.8) für eine Komponente des Jacobi-Verfahrens sind die Berechnungen der einzelnen Elemente des Vektors x(k+[1]) unabhängig voneinander und können somit parallel im Fragment-Shader der GPU durchgeführt werden. Dabei bildet jeder Itera- tionsschritt n Instanzen des verwendeten Fragment-Programms. Weil die Berechnung einer Komponente jeweils den gesamten vorherigen Iterationsvektor x(k) benötigt, ist es günstig, diesen jeder Fragment-Programminstanz als uniforme Variable über den Texturspeicher zur Verfügung zu stellen. Werden darüber hinaus die neu berechneten Werte der Komponenten ebenfalls in den Texturspeicher geschrieben, so lässt sich durch den Austausch von Ziel- und Quelltextur ( ”flip-flop“)desx-VektorsamEndeeinesIterationsschrittsweitererunnötiger Datenaustausch zwischen CPU und GPU vermeiden. Bei einem solchen Vorgehen entspricht ein Durchlauf durch die Grafikpipeline genau einem Iterationsschritt des Jacobi-Verfahrens.

4.2 Ein Schnitttest zwischen Strecke und Dreieck

Die in Abschnitt 2.3.3 erläuterte näherungsweise Berechnung der Formfaktoren mithilfe von Kreisscheiben setzt voraus, dass bekannt ist, ob sich je zwei Flächenelemente gegenseitig sehen oder nicht. Um die Sichtbarkeit zweier Dreiecke festzustellen, muss ihre Verbindungsstrecke mit allen anderen Dreiecken der Szene auf Schnitte untersucht werden. Da die Anzahl der für die gesamte Szene durchzuführenden Strecke-Dreieck-Schnitttests kubisch in der Anzahl der Flächenelemente ist, sollte der verwendete Algorithmus schnell und ressourcenschonend sein.

Einen sehr effizienten und zugleich stabilen Schnitttest zwischen einem Strahl und einem Dreieck bietet der Algorithmus von Möller, vgl. [MT97]. Er kommt in der Praxis häufig beim Raytracing zum Einsatz[9] und lässt sich leicht zu einem Schnitttest zwischen einer Strecke und einem Dreieck anpassen. Es folgt zunächst die Beschreibung einer solchen Adaption, bevor im Anschluss daran einige Vorschläge zur Verbesserung seiner Leistungsfähigkeit diskutiert werden, die sich im Fall dieser speziellen Anwendung eröffnen.

4.2.1 Adaption des Strahl-Dreieck-Schnitttests von Möller

In [MT97] präsentiert Möller gemeinsam mit Trumbore einen Algorithmus, der bestimmt, ob ein Strahl ein Dreieck schneidet oder nicht. Dieser verzichtet im Gegensatz zu vielen anderen solcher Algorithmen auf einen Schnitttest des Strahls mit der Dreiecksebene und spart so die Berechnung bzw. Speicherung der Ebenengleichung. Erreicht wird dies über eine Transformation des Strahlursprungs und anschließendem Wechsel der Basis des dabei entstandenen Vektors. Als Ergebnis entsteht der Vektor (t, u, v)T , wobei t den Abstand vom Ursprung des Strahls zum Schnitt mit der Ebene des Dreiecks beschreibt und (u, v) die Koordinaten innerhalb des Dreiecks darstellen. Wie im Folgenden gezeigt wird, ist für die Anpassung zum Strecke-Dreieck-Schnitttest lediglich die zusätzliche Auswertung des t-Wertes

[9]Um nur drei Arbeiten zu nennen, die im Zusammenhang mit Raytracing den Algorithmus von Möller anführen, sei an dieser Stelle beispielhaft auf [BEL+06, SWW+04, KL04] verwiesen.

4.2 Ein Schnitttest zwischen Strecke und Dreieck

notwendig, wenn man zur Berechnung einen Strahl verwendet, auf dem die zu testende Strecke liegt.

Den Ausgangspunkt im Algorithmus von Möller bildet neben den Eckpunkten V0, V1, V2 des Dreiecks die Gleichung des Strahls R(t) = O + tD mit Ursprung O und Richtung D. Sie lässt sich beim Strecke-Dreieck-Schnitttest leicht mithilfe der Endpunkte C0 und C1 der gegebenen Strecke aufstellen:

(4.13)

Abbildung in dieser Leseprobe nicht enthalten

Im Gegensatz zu Möller, vgl. [MT97], wird hier auf eine Normalisierung des Richtungsvektors verzichtet, da diese nicht unbedingt notwendig ist und sich so die spätere Auswertung des t-Wertes stark vereinfacht. Ein Punkt F (u, v) auf einem Dreieck ist durch

(4.14)

Abbildung in dieser Leseprobe nicht enthalten

festgelegt, wobei (u, v) baryzentrische Koordinaten[10] sind, welche die Bedingungen u ≥ 0, v ≥ 0 und u+v ≤ 1 erfüllen müssen. Zur Berechnung des Schnittpunkts zwischen dem Strahl R(t) und dem Dreieck F (u, v) setzt man beide gleich:

Abbildung in dieser Leseprobe nicht enthalten

(4.15)

Abbildung in dieser Leseprobe nicht enthalten

(4.16)

Abbildung in dieser Leseprobe nicht enthalten

(4.17)

sodass man durch Lösen dieses linearen Gleichungssystems den Abstand t und die baryzentrischen Koordinaten (u, v) erhält. Geometrisch kann man sich dieses Vorgehen als Translation des Dreiecks hin zum Koordinatenursprung und anschließender Transformation in die y-z-Ebene vorstellen, wobei zugleich der Strahl in x-Richtung gedreht wird.

Bezeichnen nun D = C1 − C0, E1 = V1 − V0, E2 = V2 − V0, T = C0 − V0 und |A, B, C| die Determinante einer aus drei Vektoren A, B und C zusammengesetzten Matrix, so gelangt man durch Anwendung der Cramerschen Regel zu folgendem Ergebnis der Gleichung (4.17):

(4.18) weiter zu

Abbildung in dieser Leseprobe nicht enthalten

(4.19)

Abbildung in dieser Leseprobe nicht enthalten

mit P = D × E2 und Q = T × E1 umgeformt werden. Es wird empfohlen die genannten

[10]Baryzentrische Koordinaten beschreiben einen Punkt durch die Koeffizienten einer affinen Kombination, also einer Linearkombination von Punkten, bei der die Summe der Koeffizienten gleich Eins ist.

Abbildung in dieser Leseprobe nicht enthalten

Quelltext 4.1: Pseudocode des angepassten Algorithmus von Möller

Faktoren auch bei der Implementierung des Algorithmus zwischenzuspeichern, um durch ihre Wiederverwendung die Kosten für eine wiederholte Berechnung zu sparen.

Der Quelltext 4.1 zeigt den angepassten Algorithmus in Pseudocode. Dabei wird ein wei- terer Vorteil deutlich: Sobald ein Wert eine notwendige Bedingung nicht erfüllt, kann abge- brochen werden, da folglich kein Schnitt mehr möglich ist. So brauchen beispielsweise v und t erst gar nicht berechnet werden, wenn bereits u außerhalb der geforderten Grenzen liegt. Der Strahl, auf dem die zu testende Strecke liegt, schneidet das Dreieck, sofern u und v die oben genannten Voraussetzungen erfüllen. Damit auch die Strecke selbst das Dreieck schneidet, muss sich die Schnittstelle zwischen ihrem Anfangs- und Endpunkt befinden. Geht man von der Strahlgleichung (4.13) mit nicht normiertem Richtungsvektor aus, muss der Parameter t zwischen Null und Eins liegen, um dem zu entsprechen. Zur Gewährleistung numerischer Stabilität wird zu Beginn die Determinante mit einem kleinen Intervall um Null herum vergli- chen. So können zum Dreieck parallel verlaufende Stahlen eliminiert werden. Weiterführende Informationen findet man in [MT97] und [SF01].

4.2 Ein Schnitttest zwischen Strecke und Dreieck

4.2.2 Vorschläge zur Verbesserung der Leistungsfähigkeit

Betrachtet man den Quelltext 4.1 genauer, so kann man feststellen, dass der Schnitt mit einem Dreieck unter Umständen erst kurz vor Ende des Algorithmus ausgeschlossen wird. Die Idee, den Algorithmus zu beschleunigen, beruht auf der Annahme, dass viele Dreiecke die ersten Bedingungen erfüllen, obwohl die Strecke diese nicht schneidet. Hier kann man ansetzen und mit einer sehr einfach gehaltenen Vorberechnung den Schnitt mit einem Dreieck eventuell schon vorher ausschließen, um so die Kosten der Schnittberechnung zu sparen.

Kann ein Schnitt erst beim Prüfen des Parameters t ausgeschlossen werden, so berechnet der obige Algorithmus 2 Vektorprodukte, 4 Skalarprodukte, 3 Multiplikationen, 1 Division sowie 4 Vektorsubtraktionen, vgl. Quelltext 4.1. Schlüsselt man die Vektoroperationen auf, so ergibt dies einen Gesamtaufwand von 27 Multiplikationen, 1 Division und 18 Additionen/ Subtraktionen (bei einem frühen Abbruch des Algorithmus entsprechend weniger). Ein guter Vorabtest sollte prüfen, ob die Dreiecke in unmittelbarer Umgebung der Strecke liegen, und vergleichsweise kostengünstig zum obigen Algorithmus sein.

Da Schnitttests an Kugeln oft recht einfach zu berechnen sind, könnte man eine oder meh- rere dazu nutzen, um die Strecke in ihnen einzubetten. Ein Kugel-Dreieck-Schnitttest besitzt jedoch eine ähnliche Komplexität wie ein Strecke-Dreieck-Test und kann daher nicht zur Be- schleunigung beitragen. Um die Vorberechnung dennoch einfach zu gestalten, beschränkt man diese stattdessen auf einen Schnitttest zwischen Kugel(n) und der Ebene, in der das Dreieck liegt. Der Vorabtest reduziert sich damit auf die Abstandsbestimmung der Kugelmittelpunkte von der Dreiecksebene.

Vorberechnung mithilfe einer Kugel

Seien A und B die Endpunkte der gegebenen Strecke, V0, V1 und V2 die Eckpunkte des auf Schnitt zu überprüfenden Dreiecks, n = (nx, ny , nz ) ein Vektor senkrecht zur Dreiecksebene und (x, y, z) ein beliebiger Punkt in dieser Ebene (z. B. ein Eckpunkt des Dreiecks). Eine Ku-

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 4.1: Vorberechnung zum Strecke-Dreieck-Schnitttest mithilfe einer Ku- gel. Schneidet die Ebene des Dreiecks (V0, V1, V2) die Kugel mit Mittelpunkt M nicht, so kann auch ein Schnitt mit der Strecke AB ausgeschlossen werden.

der Abstand D des Kugelmittelpunkts zu dieser Ebene damit wie folgt, vergleiche [Wei]:

(4.20)

Abbildung in dieser Leseprobe nicht enthalten

Verwendet man für n den Normalenvektor des Dreiecks (V0, V1, V2), dessen Betrag gleich Eins ist, so vereinfacht sich (4.20) zu

(4.21)

Abbildung in dieser Leseprobe nicht enthalten

Die Bestimmung der Normalenvektoren verursacht hier keinen zusätzlichen Aufwand, weil diese ohnehin bei der Formfaktorberechnung nach Abschnitt 2.3.3 benötigt werden. Um den Schnitt mit einem Dreieck ausschließen zu können, muss D mindestens so groß wie der Radius r der Kugel sein, sodass folgende Bedingung zu einem frühzeitigen Abbruch des Schnitttests führt:

(4.22)

Abbildung in dieser Leseprobe nicht enthalten

Dabei gilt die Relation

”≤“,dadieeinzigenbeidenPunktederKugeloberfläche,diedie

Strecke berühren, die Endpunkte der Strecke selbst sind und damit den beiden Dreiecken angehören, deren Sichtbarkeit bei diesem Test geprüft werden soll. Sowohl Kugelmittelpunkt als auch Radius müssen für jede Strecke nur einmal bestimmt werden und können dann für die Schnitttests mit allen anderen Dreiecken der Szene genutzt werden. Die Kosten ihrer Ermittlung sind damit für die Gesamtberechnung vernachlässigbar klein.

Eine Vorberechnung der hier beschriebenen Form benötigt gemäß ([4].[22]) lediglich [3] Multiplikationen sowie [6] Additionen/Subtraktionen und ist damit im Vergleich zum Algorithmus von Möller, der mindestens [9] Multiplikationen und [7] Additionen/Subtraktionen erfordert (vergleiche Quelltext [4].[1]), recht günstig. Allerdings wird diese Heuristik nur eine Verbesserung bewirken, wenn mit ihrer Hilfe die Schnitte vieler Dreiecke vorab ausgeschlossen werden können, da andernfalls der Mehraufwand höher als die Ersparnis ist. So lassen sich mit hoher Wahrscheinlichkeit Fälle konstruieren, in denen der Algorithmus ohne Vorberechnung besser abschneidet als mit einer solchen Vorberechnung.

Vorberechnung mithilfe von n Kugeln

Stellt man sich vor, die Strecke AB verläuft von einem äußersten Dreieck der Szene zu einem weit entfernten Dreieck an der gegenüberliegenden Außenseite, so wird deutlich, dass die Kugel einen recht großen Raum einnimmt. Je größer dieser eingenommene Raum ist, umso niedriger ist die Wahrscheinlichkeit, dass die anderen Dreiecke die Bedingung ([4].[22]) erfüllen. Um mehr Dreicke ausschließen zu können, kann man den zu überprüfenden Raum um die Strecke verschlanken, indem die Vorberechnung, wie in Abbildung [4].[2] dargestellt, mehrere Kugeln mit geringerem Radius verwendet.

Nutzt man für die Vorberechnung n Kugeln mit n > [1] und platziert diese entsprechend Abbildung [4].[2], so behält die

”≤“-RelationinderBedingung([4].[22])ihreGültigkeit,wennsich die Kugeln etwas durchdringen. Hierzu wählt man den Radius jeder der Kugeln mit

(4.23)

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 4.2: Vorberechnung zum Strecke-Dreieck-Schnitttest mithilfe von drei Kugeln. Schneidet die Ebene, in der das Dreiecks (V0, V1, V2) liegt, keine der drei Kugeln, so kann auch ein Schnitt mit der Strecke AB ausgeschlossen werden.

wobei µ > 0, µ ∈ R ein kleiner Offset ist. Die Mittelpunkte der Kugeln werden so auf der Strecke AB verteilt, dass die beiden Endpunkte der Strecke zur Oberfläche der jeweils äußeren Kugel gehören. Dadurch erfüllen auch Dreiecke in der gleichen Ebene wie die beiden Dreiecke, deren Sichtbarkeit ermittelt wird, die Bedingung (4.22), sodass diese Schnitte ebenfalls vorab ausgeschlossen werden können. Ein Vorteil ergibt sich dabei vor allem in Szenen mit fein unterteilten, ebenen Flächen und Wänden.

Entsprechend (4.22) kostet der Test mit jeder der n Kugeln jeweils 3 Multiplikationen und

6 Additionen/Subtraktionen. Das bedeutet, je mehr Kugeln verwendet werden, umso mehr Dreiecke können zwar beim Vorabtest frühzeitig verworfen werden, desto teurer ist aber auch die Vorberechnung insgesamt. Eine Vorberechnung mithilfe von 10 Kugeln lohnt sich in jedem Fall nicht mehr, da die Kosten dieser Berechnung bereits die eigentlichen Kosten des angepassten Algorithmus von Möller aus Abschnitt 4.2.1 übersteigen.

Die Ergebnisse in Abschnitt 6.1.2 zeigen, dass in vielen Fällen mit nur einer Kugel die besten Resultate erzielt werden und dass ebenso Szenen konstruiert werden können, an denen die hier vorgeschlagenen Vorberechnungen mithilfe von Kugeln zur Verlangsamung der Schnitttests führen.

4.3 Octrees

Wie bereits festgestellt wurde, verlangt die Approximation der Formfaktoren mittels Kreis- scheiben nach der Bestimmung der Sichtbarkeit von je zwei Dreiecken. Die Anzahl der dazu notwendigen Schnittberechnungen liegt in O(n[3]) bei einer aus n Flächenelementen beste- henden Szene, vgl. Abschnitt 4.2. Dies sind bei komplexeren Szenen trotz eines effizienten Schnitttests zu viele. Um von vorn herein weniger Dreiecke in die Schnittberechnung schicken zu müssen, bietet sich die Verwendung eines Raumunterteilungsverfahrens an. Bei einem der- artigen Verfahren wird der Objektraum in Unterräume geteilt, sodass in diesem Fall nur die Objekte in relevanten Unterräumen auf Schnitte überprüft werden müssen.

Als Hauptproblem einer solchen Raumunterteilung gilt die Abwägung zwischen dem hohen Speicherverbrauch und der Genauigkeit. Hier gilt: je feiner eine Szene unterteilt ist, desto weniger Objekte sind für die Schnitttests relevant, desto höher ist jedoch auch der Speicher- verbrauch. Eine Möglichkeit, den benötigten Speicher zu reduzieren, bietet die strukturierte Organisation der Daten.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 4.3: Octree-Darstellung: Zweifach unterteilter kubischer Raum mit dem dazu- gehörigen Octree. Quelle: Nach [Wat02].

Dies geschieht häufig durch die Verwendung eines Octree, vgl. [Wat02]. Dabei handelt es sich um eine hierarchische Datenstruktur, welche die Verteilung der Objekte einer Szene im von der Szene eingenommenen dreidimensionalen Raum beschreibt. Dem liegt die Idee zugrunde, den Raum gleichmäßig in acht gleichgroße Teilräume zu unterteilen. Jeder dieser Teilräume kann wiederum aus acht Unterräumen bestehen. Diese Unterteilung führt man bis zum Erreichen der gewünschten Genauigkeit fort. Dabei brauchen Teilräume, die leer sind oder nur ein Objekt enthalten, nicht weiter unterteilt werden. Jeder Knoten des Octree besitzt somit entweder acht Kinder oder keines. Abbildung 4.3 zeigt einen so unterteilten Raum mit dem dazugehörigen Octree. In der Regel enthalten die Blätter des Baums die Informationen. Sie bestehen häufig aus Zeigern auf eine Datenstruktur für jedes Objekt in diesem Bereich. Da bei einer aus Polygonen bestehenden Szene im Allgemeinen mehrere Polygone in dem durch ein Blatt repräsentierten Segment liegen, stellen in diesem Fall meist Listen von Zeigern die Blätter dar. Daran lässt sich erkennen, dass es sich bei Octrees in der Regel nicht um eine abgeschlossene Darstellungsform handelt, sondern sie oft Teil einer hybriden Struktur sind.

Obwohl die Raumunterteilung mittels Octree für die Berechnung der Sichtbarkeiten ein hohes Einsparpotenzial aufweist, dürfen die Kosten für seine Erstellung nicht außer Acht gelassen werden. Bei der Zerlegung einer Szene in einen Octree handelt es sich durchaus um eine aufwändige Operation, müssen doch die Koordinaten des Begrenzungsrahmens zu jedem Dreieck (also jeweils kleinste und größte x-, y- und z-Koordinate) ermittelt und in entsprechende Einheiten des Octree umgerechnet werden.

Im Allgemeinen sollte die Zerlegung nicht so fein sein, dass die Einsparungen vom höheren Aufwand bei der Auswertung des Octree zunichtegemacht werden. Die Ergebnisse in Abschnitt 6.1.2 zeigen, dass bereits geringe Baumtiefen von vier zu guten Ergebnissen führen können. Einen Nachteil haben Octrees bei Szenen, in denen die Polygone sehr ungleichmäßig verteilt sind. Ist ein Großteil des von der Szene eingenommenen Raums leer und befinden sich einige detailreiche Objekte mit vielen Polygonen auf nur einer Seite des Raums, so führt eine Octree-Zerlegung zur häufigen Unterteilung von leerem Raum.

5 Entwurf und Implementierung

Um die Vor- und Nachteile einer GPU-Implementierung des Radiosity-Verfahrens gegenüber seiner Umsetzung auf der CPU feststellen zu können, wurden beide Herangehensweisen in einem Programm namens Radioso realisiert. Die folgenden Abschnitte widmen sich neben den Besonderheiten der Implementierung auch grundlegenden Designentscheidungen sowie den verwendeten Eingabeformaten und der grafischen Benutzeroberfläche dieser Umsetzung.

5.1 Entwurf

Beim Entwurf von Radioso spielten verschiedene anwendungsspezifische Aspekte eine ent- scheidende Rolle. Beispielsweise sollen die Anwender die Möglichkeit haben, je nach Wunsch und verfügbarer Hardware, zwischen unterschiedlichen Implementierungen einzelner Rechen- schritte des Radiosity-Verfahrens wählen zu können, um so den für sie optimalen Berechnungs- weg zu nutzen oder die unterschiedlichen Implementierungen in ihrer Leistungsfähigkeit zu vergleichen. Zudem ist es aufgrund des hohen Speicherverbrauchs des Radiosity-Verfahrens günstig, wenn alle Implementierungen auf gemeinsame Datenstrukturen zugreifen. Dadurch kann redundante Datenhaltung vermieden werden. Dies gut strukturiert umzusetzen, verlangt nach speziellen Entwurfsmustern.

Ein Entwurfsmuster, das den Wechsel zwischen verschiedenen Implementierungen eines Algorithmus ermöglicht, heißt Strategie, vgl. [GHJV04]. Es erlaubt, das aufrufende Objekt mit einer der zur Verfügung stehenden Varianten zu konfigurieren. Die unterschiedlichen Varianten werden dabei in eigenständigen Klassen oder Klassenhierarchien implementiert. Abbildung 5.1 zeigt die Umsetzung dieses Entwurfsmusters anhand eines Ausschnitts aus dem

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 5.1: Ausschnitt aus dem UML-Klassendiagramm von Radioso.

Klassendiagramm von Radioso. Die drei Varianten FFCalculatorCPU, FFCalculatorGPU und FFCalculatorCPUHemicube zur Formfaktorberechnung bilden jeweils Unterklassen der abstrakten Klasse FFCalculator und implementieren die Methode calculateFormFactors. Das Radiosity-Objekt beinhaltet eine Komponente ffcalculator vom Typ FFCalculator, welche mit einer der drei Unterklassen konfiguriert wird. Die Auslösung der Formfaktorberechnung bleibt damit von der Wahl der Variante unabhängig, wodurch sich relativ einfach neue Varianten hinzufügen lassen. Auf ähnliche Weise wurden die Unterteilung der Flächen, das Lösen des Radiosity-Gleichungssystems und das Rendern realisiert.

Um bei der Umsetzung der in Abschnitt 4.3 beschrieben Octree-Datenstruktur auf eine bestehende Implementierung zurückgreifen zu können, fand außerdem das Entwurfsmuster Adapter Anwendung. Dabei werden in einer Wrapper-Klasse die Schnittstellen der adaptierten Klasse an die von der Anwendung benötigten Schnittstellen angepasst, vgl. [GHJV04].

5.2 Besonderheiten der Implementierung auf der GPU

Trotz Generierung der Shader-Programme durch den dynamischen Compiler der RapidMind- Plattform muss sich der Programmierer um einige grafikhardwarespezifische Einschränkun- gen selbst Gedanken machen. So existieren z. B. Beschränkungen bei der Anzahl möglicher Schleifeniterationen in Shader-Programmen und bei der zulässigen Breite und Höhe von Tex- turen. Im Folgenden werden während der Entwicklung von Radioso aufgetretene Probleme diskutiert und anschließend die vorgeschlagenen Lösungsansätze am Beispiel der Realisierung des Jacobi-Verfahrens aufgezeigt.

5.2.1 Verwendete Entwicklungsumgebung

Zur Entwicklung von Radioso wurde größtenteils das Release 2.0.1 der Entwicklungsplatt- form von RapidMind verwendet. Einige Gründe für diese Wahl wurden bereits in Abschnitt

3.3 genannt. Zudem gestattet das C++-Interface der Plattform die Realisierung der in Abschnitt 5.1 beschriebenen objektorientierten Entwurfsmuster. Ein weiterer Vorteil ist, dass viele grafikhardwarespezifische Details wie der Datenaustausch zwischen CPU und GPU vor dem Programmierer verborgen werden. Einige Beschränkungen heutiger Grafikkarten bleiben dabei dennoch unberücksichtigt. An dieser Stelle ist es notwendig, die zu implementierenden Algorithmen entsprechend anzupassen.

Die Entwicklung fand weiterhin auf zwei unterschiedlichen Grafikkarten statt. Zum Einen auf einer NVIDIA GeForce 7600 GS mit dem Grafiktreiber NVIDIA ForceWare in der Version 93.71, zum Anderen auf einer ATI Radeon X1600 Pro unter dem Treiber ATI Catalyst 7.4. Beide Grafikkarten weisen ein ähnliches Leistungspotenzial auf und unterlie- gen ähnlichen Beschränkungen. Alle übrigen Komponenten der beiden Entwicklungsrechner waren weitestgehend identisch. Nähere Einzelheiten zur verwendeten Hardware können dem Abschnitt 6.1 entnommen werden.

5.2.2 Beschränkte Anzahl an Schleifeniterationen

Der Typ Program des RapidMind-Interface, dessen Objekte vom dynamischen Compiler der Plattform in Shader-Programme überführt werden, unterstützt Schleifen der Formen FOR/ENDFOR, WHILE/ENDWHILE sowie DO/UNTIL. Diese Schlüsselwörter werden in Form von

5.2 Besonderheiten der Implementierung auf der GPU

Makros verwendet, die entsprechende Steuerstrukturen der Grafikhardware kapseln, vergleiche Abschnitt 3.3.3. Sind arrayWideUniform und arrayUniform uniforme Variablen vom Typ Value1ui bzw. Array<2,Value3f>, so addiert das RapidMind-Program des folgenden Beispiels alle Elemente einer Zeile des uniformen zweidimensionalen Arrays:

Abbildung in dieser Leseprobe nicht enthalten

Definiert man ein weiteres Array vom Typ Array<1,Value1ui> und füllt es mit Zeilenindices, so können durch seineÜbergabe an das Program-Objekt row_sum die angegebenen Zeilensum- men unabhängig voneinander und parallel berechnet werden. Man beachte, dass der Ausdruck arrayUniform[Value2ui(j,i)] das Element aus Zeile i und Spalte j des Arrays zurück gibt und nicht anders herum, wie vielleicht erwartet. Leider liefert die Berechnung des daraus er- zeugten Shader-Programms auf der GPU nur korrekte Ergebnisse, solange der Wert von arrayWideUniform nicht größer als 255 ist. Bei mehr als 255 Schleifendurchläufen wird die Berechnung auf beiden Entwicklungsrechnern unabhängig von der Schleifenform abgebrochen, was sich ausschließlich am unerwarteten Ergebnis bemerkbar macht. Eine Fehlermeldung oder Warnung wird nicht ausgegeben und die zugehörige Anwendung läuft normal weiter.

Diese Beschränkung kann durch die Schachtelung von Schleifen gemindert werden. Damit das obige Beispiel auch bei mehr als 255 Iterationen korrekte Ergebnisse liefert, ändert man es beispielsweise wie folgt:

Abbildung in dieser Leseprobe nicht enthalten

5.2.3 Beschränkte Größe von Texturen, Arrays und Eingabedaten

Eine durchaus schwieriger zu handhabende Einschränkung von Grafikhardware ist die be- grenzte Breite und Höhe von Texturen. So führt die Verwendung von zweidimensionalen Arrays, die breiter oder höher als 4096 Elemente sind, auf beiden Entwicklungssystemen zu einem Abbruch der Anwendung mit einem Fehler im GLSL-Backend aufgrund eines ungülti- gen Wertes bei der Texturbehandlung. Eine Abhilfe kann in diesem Fall nur die Verteilung der Daten auf mehrere Arrays schaffen. Da der Array-Typ von RapidMind keine Elemente seines eigenen Typs unterstützt, hat dies, bei von der Eingabe der Anwendung abhängigen Array-Größen, zusätzlich die Verteilung der Berechnung auf mehrere RapidMind-Program- Aufrufe zur Folge.

Des Weiteren werden korrekte Berechnungen auf der GPU nur bei Array-Größen von höchstens 2[22] = 4 194 304 Elementen garantiert. Laut RapidMind-Entwicklerforum [rapb] ergibt sich diese Zahl aus der maximalen Texturgröße abzüglich einiger Elemente, die unter Umständen für Offsets und Schrittweiten bezüglich bestimmter Zugriffsmuster auf die ArrayInhalte verwendet werden. Müssen mehr Daten bewältigt werden, kann auch hier nur die Verteilung auf mehrere Arrays und Program-Aufrufe aushelfen.

Ein weiteres Problem stellt die begrenzte Größe des Grafikspeichers dar, was bei derÜberschreitung aufgrund zu großer Datenmengen zu einem Abbruch der Anwendung mit der Meldung GL_OUT_OF_MEMORY des GLSL-Backend führt. Durch die Verlagerung von Daten aus uniformen Variablen in Eingabedatenströme für Shader-Programme kann der benötigte Speicherplatzbedarf auf der Grafikhardware verringert werden. Einen weiteren Lösungsansatz bietet die Aufteilung der Berechnung, wie sie zuvor beschrieben wurde.

Darüber hinaus wurden bei Berechnungen auf der NVIDIA GeForce 7600 GS weitere Beschränkungen festgestellt, welche nicht näher spezifiziert werden konnten. Diese äußern sich dadurch, dass die Verarbeitung von großen Datenmengen trotz Beachtung der obigen Beschränkungen nicht zu den erwarteten Ergebnissen führt. Verringert man aber die Einga- bedatenmengen mithilfe der oben beschriebenen Aufteilung der Berechnung, so erhält man die erwarteten Ergebnisse.

Der folgende Abschnitt zeigt am Beispiel des Jacobi-Verfahrens, wie sich ein Algorithmus durch die Verteilung der Daten auf mehrere Arrays verändert.

5.2.4 Partitionierung des Jacobi-Verfahrens

Das Jacobi-Verfahren wird aufgrund seiner sehr guten parallelen Eigenschaften zum Lösen des Radiosity-Gleichungssystems (2.5) auf der GPU eingesetzt. Eine Beschreibung des Verfahrens sowie einer möglichen parallelen Realisierung finden sich in Abschnitt 4.1. Das dort geschil- derte parallele Vorgehen ist in Quelltext 5.1 unter Verwendung der RapidMind-Plattform implementiert. Dabei stellen die Radiosity-Matrix M, der Strahlungsvektor B und der Vektor der emittierten Strahlungswerte E sowie die Anzahl der Flächenelemente numPatchesUniform uniforme Variablen dar und sind für alle Instanzen des Program-Objekts jacobi_iteration verfügbar. Den einzigen Eingabewert von jacobi_iteration bildet der Index des zu er- mittelnden Elements von B. Nachdem alle neuen Werte eines Iterationsschritts berechnet wurden, werden die Ausgaben von jacobi_iteration direkt in den Vektor B geschrieben. Dies ermöglicht die beliebige Wiederholung eines solchen Iterationsschritts ohne die Daten von der GPU zur CPU übermitteln zu müssen.

Aufgrund der in Abschnitt 5.2.2 genannten Beschränkung wird eine solche Implementierung nur korrekte Ergebnisse liefern, wenn die Anzahl der Flächenelemente der Szene 255 nicht übersteigt. Durch Schachtelung der FOR-Schleife kann diese Beschränkung gelockert werden. Eine solche Anpassung ermöglicht das Lösen des Gleichungssystems aber nur für Szenen mit bis zu 2048 Flächenelementen, da andernfalls die Matrix M die zulässige Array-Größe überschreitet. Für Szenen mit mehr Flächenstücken ist eine Verteilung der Daten, wie sie in Abschnitt 5.2.3 beschrieben wurde, notwendig. Die hierzu implementierte Lösung zerlegt die Matrix M in gleichgroße quadratische Teilmatrizen der gewünschten maximalen Array-Größe. Fehlende Elemente am rechten und unteren Rand von M sowie am Ende der Vektoren E, B Value1ui numPatchesUniform = numPatches; Array<2,Value3f> M(numPatches, numPatches); Array<1,Value3f> E(numPatches);

Abbildung in dieser Leseprobe nicht enthalten

und tmpB werden gegebenenfalls mit Nullwerten aufgefüllt. Der Quelltext 5.2 zeigt das derart angepasste Jacobi-Verfahren. Daran wird sehr gut deutlich, wie die Komplexität des Quell- codes gegenüber der unpartitionierten Berechnung in Quelltext 5.1 zunimmt. Jeder Aufruf des Program-Objekts jacobi_iteration_segment berechnet nun Zwischenwerte mithilfe der zu diesem Zeitpunkt im Texturspeicher vorhandenen Teilmatrix. Dabei ist zu beachten, auch die Indices entsprechend anzupassen. Da E, B und tmpB vollständig im Texturspeicher liegen, wird der Zeilenindex i doppelt an das Program-Objekt übergeben: einmal für die Zeile in der Teilmatrix und einmal für die zugehörige Zeile in der Gesamtmatrix. Hierfür wird der Typ Value2ui von RapidMind genutzt.

5.3 Szenenbeschreibung

Radioso unterstützt zwei verschiedene Eingabeformate zur Beschreibung der darzustellenden

Szene: MDL (steht für den engl. Begriff

”model“)undRML(RadiosityModelingLanguage).

Bei beiden Formaten handelt es sich um einfache Beschreibungsformen für [3]D-Szenen durch Dreiecksflächen mit einigen zusätzlichen Angaben. Während MDL ein sehr primitives und durchaus übliches Format ist, das anwendungsspezifisch an die benötigten Anforderungen angepasst wird, wurde RML eigens für Radioso entworfen, um das Beschreiben von [3]D- Szenen etwas einfacher zu gestalten.

Damit die Szenenbeschreibungen von Radioso verarbeitet werden können, müssen diese in ASCII-codierten Textdateien mit den Dateiendungen .mdl bzw. .rml vorliegen. In den folgenden Abschnitten werden beide Formate genauer vorgestellt.

5.3.1 MDL

Szenenbeschreibungen im MDL-Format beginnen stets mit einem Integer-Wert, der angibt, aus wie vielen Polygonen die Szene besteht. Darauf folgen für jedes Polygon eine Reihe von Werten nach einem bestimmten Schema. Im Fall von Radioso besteht jedes Polygon aus drei Eckpunkten, einer Farbe (dem diffusen Reflexionskoeffizienten) sowie der vom Polygon emittierten Strahlung. Zur Wahrung der Übersichtlichkeit können zwischen den Werten be- liebige Leerzeichen und Leerzeilen eingefügt werden, sodass jedes Polygon beispielsweise aus

5 Zeilen mit je 3 Fließkommazahlen besteht. Die ersten 3 Zeilen bilden die Eckpunkte des Dreiecks und bestehen jeweils aus x-, y- und z-Koordinate, die 4. und 5. Zeile geben den diffusen Reflexionskoeffizienten und die emittierte Strahlung in Form von RGB-Farbwerten an. Damit die Normalen der Dreiecksflächen korrekt bestimmt werden, sollte die Angabe der Eckpunkte im mathematisch positiven Sinn erfolgen. Das Spezifikationsschema zur hier formulierten Beschreibungssprache ist zusammen mit einem Beispiel in Anhang B zu finden.

5.3.2 RML

MDL besitzt zwei wesentliche Nachteile. Zum Einen müssen für jede Dreiecksfläche sowohl ihre Farbe als auch das von ihr emittierte Licht angegeben werden, obwohl häufig mehrere solcher Flächen ein gemeinsames einfarbiges Objekt bilden, zum Anderen wird die Wieder- verwendung von Objekten erschwert, da für alle Eckpunktkoordinaten die absoluten Werte in der Szene anzugeben sind. RML steht für Radiosity Modeling Language und stellt eine

5.3 Szenenbeschreibung

Abbildung in dieser Leseprobe nicht enthalten

Erweiterung des MDL-Formats dar. Dabei werden mehrere Dreiecksflächen in Flächengruppen zusammen gefasst, sodass der diffuse Reflexionskoeffizient und die emittierte Strahlung je Gruppe nur einmal angegeben werden. Darüber hinaus besteht die Möglichkeit für jede der Flächengruppen Translations-, Rotations- sowie Skalierungsangaben zu machen.

Eine RML-Beschreibung beginnt mit einem Integer-Wert für die Anzahl der vorhandenen Flächengruppen. Jede Flächengruppe besteht aus 5 Zeilen mit je 3 Fließkommawerten für die Angaben von Farbe, emittiertem Licht, Rotation, Translation und Skalierung. Diese werden gefolgt von einem Integer-Wert für die Anzahl der zugehörigen Dreiecksflächen sowie je 9 Fließkommakoordinaten pro Zeile für die Eckpunkte eines jeden Dreiecks. Auch hier dürfen beliebige Leerzeichen und Leerzeilen zur Verbesserung der Übersichtlichkeit verwendet wer- den. Wie bei MDL ist die Angabe der Eckpunkte ebenfalls im mathematisch positiven Sinn vorzunehmen. Dem Anhang B können eine beispielhafte Szenenbeschreibung sowie das RML- Spezifikationsschema entnommen werden.

5.4 Grafische Benutzeroberfläche

Zur Vereinfachung der Bedienung wurde Radioso mit einer grafischen Benutzeroberfläche versehen. Diese ermöglicht neben der einfachen Auswahl der zu verwendenden Algorithmen und dem Tätigen verschiedener Einstellungen auch das Betrachten der erzeugten Szenen. Zudem können berechnete Szenen abgespeichert und geladen werden. Die Umsetzung dieser Bedienoberfläche wurde mithilfe der Qt-Bibliothek realisiert.

Wie die Abbildung 5.2 zeigt, ist die genannte Benutzeroberfläche in drei Bereiche gegliedert: einer OpenGL-Ausgabe links, einiger Bedienelemente rechts sowie einem Textausgabefeld un- ten. Nach dem Öffnen einer Szenenbeschreibungsdatei wählt man anhand der Bedienelemente die gewünschten Einstellungen für die Unterteilung der Flächen, die Formfaktorbestimmung und das Lösen des Radiosity-Gleichungssystems und startet anschließend die Berechnung durch den Button ”StartCalculation“.DasTextausgabefeldinformiertjederzeitüberden aktuellen Stand der Berechnung und gibt die benötigten Rechenzeiten sowie eventuelle Fehler- und Erfolgsmeldungen aus. Nach Abschluss der Berechnung ist die gerenderte Szene im OpenGL-Ausgabebereich zu sehen.

Diese integrierte OpenGL-Ausgabe bietet verschiedene Steueroptionen zum Betrachten einer gerenderten Szene. Um sich durch die [3]D-Szene zu bewegen, verwendet man die Pfeil- tasten zusammen mit der Maus. Die Pfeiltasten ermöglichen dabei das Bewegen nach links, vorn, hinten und rechts. Alternativ können hierzu auch die Tasten A, W, S und D genutzt werden. Durch das Gedrückthalten der linken Maustaste in Kombination mit einer Mausbe- wegung erfolgt eine Rotation um eine Achse durch den Betrachterstandpunkt. Mithilfe der Bedienelemente auf der rechten Seite können überdies das Verschmelzen der Flächenelemente aktiviert und deaktiviert, die Richtung des Up-Vektors zwischen x-, y- und z-Richtung vertauscht und die Geschwindigkeit der Bewegung reguliert werden. Der Button ”Reset“setzt sämtliche Einstellungen der OpenGL-Ausgabe auf ihre ursprünglichen Werte zurück.

[56]

Abbildung 5.2: Die grafische Benutzeroberfläche von Radioso.

Abbildung in dieser Leseprobe nicht enthalten

6 Vergleich von CPU- und GPU- Implementierung

Das primäre Ziel dieser Arbeit ist es, die Vor- und Nachteile einer GPU-Implementierung des Radiosity-Verfahrens gegenüber der herkömmlichen Entwicklung auf der CPU heraus- zuarbeiten. Zu diesem Zweck wurden Leistungsvergleiche durchgeführt, deren Ergebnisse im folgenden Abschnitt diskutiert werden. In einem weiteren Abschnitt werden darüber hinaus Probleme geschildert, die im Zusammenhang mit der GPU-Programmierung aufgetreten sind.

6.1 Ergebnisse der Zeitmessungen

Wie eingangs der Arbeit erwähnt, konnten bereits in der Vergangenheit beachtliche Laufzeitvorteile von GPU-Implementierungen gegenüber herkömmlichen CPU-Berechnungen erzielt werden. Auch das Radiosity-Verfahren wurde unter Einsatz von Shader-Hochsprachen bereits erfolgreich auf der GPU umgesetzt. Im Unterschied zu den genannten bisherigen Arbeiten entstand das im Rahmen dieser Arbeit entwickelte Programm Radioso nicht durch die explizite Programmierung in einer der Shader-Hochsprachen, sondern unter Verwendung eines GPGPU-Werkzeugs (RapidMind-Entwicklungsplattform).

Alle im Folgenden präsentierten Ergebnisse bilden das arithmetische Mittel aus jeweils min- destens fünf Zeitmessungen. Aufgrund der teilweise sehr hohen Laufzeiten von bis zu 21 Stun- den wurde von der Aufnahme von mehr als fünf Messungen Abstand genommen. Sämtliche auf der CPU durchgeführte Vergleichsmessungen fanden auf Intel Pentium 4-Prozessoren mit 3,0 GHz Takt statt. Für die Messungen der GPU-Zeiten standen zwei unterschiedliche Grafikkarten mit ähnlichen Leistungsmerkmalen zur Verfügung: NVIDIA GeForce 7600 GS und ATI Radeon X1600 Pro. Beide verfügen über je 5 Vertex- und 12 Pixel-Shader. Die NVIDIA-Karte besitzt einen Chiptakt von 400 MHz sowie einen 256 MB umfassenden DDR2- Speicher mit einem Speichertakt von 400 MHz (800 MHz effektiv) und erzielt eine Bandbreite von 12,8 GB/s. Im Unterschied dazu verfügt die ATI-Karte zwar über einen Chiptakt von 500 MHz und einen doppelt so großen DDR2-Speicher, jedoch liegt ihr Speichertakt dagegen nur bei 270 MHz (540 MHz effektiv) und die erreichte Bandbreite bei 8,5 GB/s. Die Zeit- messungen fanden unter Verwendung der Grafiktreiber NVIDIA ForceWare 93.71 bzw. ATI Catalyst 7.4 und dem Release 2.0.1 der RadidMind-Plattform satt. Alle Testrechner laufen mit dem Betriebssystem Windows XP und verfügen über 1 GB Arbeitsspeicher. Für eine kompaktere Darstellung werden die beschriebenen Hardwarekomponenten im Folgenden auf die Begriffe CPU bzw. Intel, NVIDIA und ATI reduziert.

Im Vorfeld der Entwicklung von Radioso wurde eine Machbarkeitsstudie bezüglich der Leistungssteigerung mithilfe von RapidMind erzeugter Shader-Programme durchgeführt. Hierzu wurden eine Maximumbestimmung aus k Werten sowie eine Matrixmultiplikation zweier (n × n)-Matrizen sowohl auf der CPU als auch auf der GPU implementiert. Während die gemessenen GPU-Zeiten bei der Maximumbestimmung nicht annähernd an die Laufzeiten der CPU heran kommen, überzeugt die GPU-Implementierung der Matrixmultiplikation und erreicht bei gleichem Algorithmus wie auf der CPU bis zum 30fachen der Geschwindigkeit der Berechnung. Der Grund für dieses Ergebnis liegt in der unterschiedlichen Natur beider Aufgaben. Das erzeugte Shader-Programm zur Maximumbestimmung beinhaltet nur eine Vergleichsoperation, sodass der Rechenaufwand verglichen mit dem Datenaustausch eher ge- ring ausfällt. Dagegen führt das Shader-Programm der Matrixmultiplikation zur Berechnung eines Elements der Ergebnismatrix n Multiplikationen und n Additionen durch.

Zum Laufzeitvergleich des Radiosity-Verfahrens wurden Zeitmessungen an vier verschiede- nen Testszenen mit jeweils unterschiedlich feinen Unterteilungen der Flächenelemente durch- geführt. Die vier Szenen

(L) und

”CornellBox“(CB), ”RadiosGL-Testszene“(R), ”Livingroom“

”Bunny“(B)sindinAnhangAabgebildet.IndennächstenAbschnittenfolgteine Auswertung der aus diesen Zeitmessungen hervorgegangenen Ergebnisse.

6.1.1 Laufzeitvergleich des Radiosity-Verfahrens

Da nicht jede Aufgabe für die Berechnung auf der GPU gleich gut geeignet ist, werden die Messergebnisse der einzelnen Schritte des Radiosity-Verfahrens im Folgenden getrennt betrachtet. Für einen möglichst objektiven Vergleich sind hierzu die jeweils bestmöglichen Laufzeiten von CPU- und GPU-Implementierung bei qualitativ gleichwertigem Renderergebnis in Tabelle 6.1 gegenübergestellt.

Dabei wird deutlich, dass die Unterteilung der Flächenelemente (subdiv) weniger für die GPU geeignet ist, da die auf der CPU gemessene Laufzeit in fast allen Fällen deutlich unter den von den Grafikprozessoren erreichten Zeiten liegt. Eine Ausnahme bildet die ”Bunny“- Testszene (B). Diese besteht bereits vor der Unterteilung aus deutlich mehr Polygonen als die anderen Szenen. Zudem befindet sich in ihr ein Objekt mit sehr vielen kleinen Flächenelemen- ten, die nicht unterteilt werden. Die enormen Zeitunterschiede bei den Szenen mit nur weni- gen Flächenelementen sind auf die Startup-Zeiten der GPU-Berechnungen zurückzuführen. Vor der ersten Ausführung eines Shader-Programms muss dieses zunächst kompiliert werden. Zusätzlich wird vor der ersten Nutzung der GPU durch das GLSL-Backend von RapidMind OpenGL initialisiert. Da es sich bei der Unterteilung der Flächen um den ersten Schritt des Radiosity-Verfahrens handelt, beinhalten die subdiv-Laufzeiten diese OpenGL-Initiali- sierung. Aber auch abzüglich dieser Startup-Zeit erreicht die GPU-Berechnung die CPU-Zeit nur im genannten Ausnahmefall. Ein Grund für die schlechtere Laufzeit ist die Umsetzung des Shader-Programms, welches ein Flächenelement in vier neue Elemente unterteilt, falls der Flächeninhalt ein bestimmte Größe übersteigt. Nach jedem Iterationsschritt müssen die neuen Flächenelemente aus dem Grafikspeicher zurück gelesen werden, was entsprechend viel Zeit kostet. Leider gelang mithilfe der RapidMind-API keine günstigere Realisierung. Im Vergleich der beiden Grafikkarten untereinander lässt sich feststellen, dass die ATI-Karte bei der Unterteilung der Flächen stets schneller als die von NVIDIA ist.

Ein anderes Verhalten wird bei der Berechnung der Sichtbarkeiten (viscalc) deutlich. Mit Laufzeiten, die bis zu [150]-mal schneller als auf der CPU sind (vergleiche Szene ”B-Lo“in

Tabelle [6].[1]), scheint dies die ideale Aufgabe für eine GPU-Berechnung zu sein. Die etwas größeren Differenzen zwischen den GPU-Zeiten von

zu

”R-Lo“zu ”L-Lo“sowievon ”B-Lo“ ”L-Hi“sinddurchdieinAbschnitt[5].[2]genanntenBeschränkungenderGrafikhardware

bedingt. Das entsprechende Programm wurde so entwickelt, dass die Berechnung ab [1001] Flächenelementen auf mehrere Shader-Programmaufrufe verteilt wird. Ab einer Anzahl von [4096] Flächenelementen findet zudem eine zusätzliche Partitionierung aufgrund der begrenzten Texturgröße (vgl. Abschnitt [5].[2].[3]) statt. Erneut erzielt die Grafikkarte von ATI verglichen

6.1 Ergebnisse der Zeitmessungen

Abbildung in dieser Leseprobe nicht enthalten

Tabelle 6.1: Laufzeitvergleich des Radiosity-Verfahrens zwischen Intel Pentium 4

3,0 GHz, NVIDIA GeForce 7600 GS und ATI Radeon X1600 Pro an den Szenen

”CornellBox“(CB), ”RadiosGL-Testszene“(R), ”Livingroom“(L)und ”Bunny“(B).DieSuffixeLoundHigebendenunterschiedlichenUnterteilungs- grad der Flächenelemente an, deren Anzahl der Spalte [2] entnommen werden kann. Spalte [4] zeigt die durchschnittlich zur Unterteilung der Flächenelemente (subdiv) benötigte Zeit in Millisekunden; die Spalten [5], [6], und [7] die Zeiten zur Berechnung der Sichtbarkeiten (viscalc), der Kreisscheibenapproximation bei bekannten Sichtbarkeiten (ffcalc) und zum Lösen des Radiosity-Gleichungssys- tems (radsolv) jeweils in Sekunden. Schließlich gibt total die Gesamtlaufzeit des Verfahrens in Sekunden wieder.

mit der von NVIDIA die besseren Ergebnisse und benötigt teilweise sogar nur die Hälfte an Rechenzeit.

Im Gegensatz dazu finden sich bei der Berechnung des Algorithmus zur Approximation der Formfaktoren mittels Kreisscheiben (ffcalc) die Vorteile bei der GPU von NVIDIA. Liegen die NVIDIA- und ATI-Zeiten anfänglich noch dicht beieinander, so unterscheiden sie sich mit zunehmender Anzahl an Flächenelementen immer deutlicher. Dies ist vermutlich auf das Kompilieren des Shader-Programms zurückzuführen, das bei NVIDIA länger als bei ATI dauert. Möglicherweise führen die damit verbundenen Optimierungen auch zur besseren Laufzeit auf der Grafikkarte von NVIDIA bei den Szenen mit mehr Flächenelementen. Der Vergleich zur benötigten CPU-Zeit endet wiederholt zugunsten der Grafikkarten. Für die Berechnung des Algorithmus benötigt die CPU bis zum 13-fachen der Zeit des NVIDIA- Grafikprozessors. Einzig bei sehr kleinen Szenen kann sie mit den GPUs konkurrieren, da Letztere für das Kompilieren der Shader-Programme die zusätzliche Zeit benötigten.

Ein ähnliches Verhalten findet man beim Lösen des Radiosity-Gleichungssystems (radsolv) vor. Auch hier erzielt die Grafikkarte von NVIDIA die besten Ergebnisse. Zudem scheint die ATI-Karte Probleme mit der Partitionierung der Berechnung zu haben, da die Laufzeiten mit den ersten hinzukommenden Partitionen (also ab 1001 und ab 2001 Flächenelementen) so sehr ansteigen, dass die Berechnung teilweise sogar auf der CPU merklich schneller verläuft. Ein Vorteil der CPU ist die Anwendbarkeit des Gauß-Seidel-Verfahrens. Es wurde festgestellt, dass dieses gegenüber dem Jacobi-Verfahren schneller konvergiert und bereits bei 6 bis 7 Itera- tionen ähnliche Ergebnisse wie das Jacobi-Verfahren bei 10 Iterationen erzielt. Trotz dieser 3 Iterationen weniger konnten die Zeiten der NVIDIA-Berechnung dennoch auf der CPU nicht annähernd erreicht werden. Die NVIDIA-GPU erzielt beispielsweise bei der Szene

”B-Hi“

einen Beschleunigungsfaktor von [27] gegenüber dem Intel-Prozessor (vgl. Tabelle [6].[1]). Bei derselben Testszene reicht es bei der GPU von ATI zumindest noch für eine Beschleunigung um Faktor [6].

Entsprechend der Laufzeiten der einzelnen Berechnungsschritte gelingt mithilfe der Gra- fikkartenprogrammierung gegenüber der herkömmlichen CPU-Programmierung eine enorme Leistungssteigerung in Bezug auf die Gesamtlaufzeit des Radiosity-Verfahrens. So verläuft die gesamte Berechnung gemäß Tabelle [6].[1] auf der GPU mindestens [15]-mal so schnell wie auf der CPU und erreicht im besten Fall einen Beschleunigungsfaktor von mehr als [70]. Insgesamt gesehen weist die Grafikkarte von ATI dabei die beste Gesamtlaufzeit auf. Die Karte von NVIDIA ist nur bei Szenen mit einer Anzahl an Flächenelementen im Intervall von [1001] bis etwa [3000] am schnellsten. In diesem Bereich überwiegen ihre Vorteile beim partitionierten Jacobi-Verfahren die Nachteile bei der Berechnung der Sichtbarkeiten.

6.1.2 Bewertung der Vorschläge zur Beschleunigung

In Kapitel 4 wurden zwei Vorschläge unterbreitet, mithilfe derer sich die Bestimmung der Sichtbarkeit zwischen zwei Flächenelementen möglicherweise beschleunigen lässt. Zum einen wurde in Abschnitt 4.2.2 vorgeschlagen, den zur Sichtbarkeitenbestimmung verwendeten Schnitttest zwischen einer Strecke und einem Dreieck um eine Vorberechnung mit Kugeln zu ergänzen, um dadurch im Mittel Rechenoperationen einzusparen, zum anderen besteht die Möglichkeit mithilfe der in Abschnitt 4.3 beschriebenen Octree-Struktur die Anzahl der notwendigen Schnitttests zu verringern. Beide Ansätze wurden innerhalb von Radioso implementiert und im Zeitverhalten untersucht.

6.1 Ergebnisse der Zeitmessungen

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 6.1: Vergleich der Sichtbarkeitenbestimmung mit Kugelvorberechnun-

gen auf der CPU an grob unterteilten Varianten der Szenen

(CB),

”CornellBox“

Abbildung in dieser Leseprobe nicht enthalten

”RadiosGL-Testszene“(R), ”Livingroom“(L)und ”Bunny“(B).Die Zahlen auf den Balken beziehen sich auf die für die Vorberechnung verwendete

Anzahl an Kugeln. Das jeweils beste Ergebnis ist hervorgehoben.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung [6].[2]: Vergleich der Sichtbarkeitenbestimmung mit Octree-Unterstützung

auf der CPU an grob unterteilten Varianten der Szenen

”CornellBox“(CB),

”RadiosGL-Testszene“(R), ”Livingroom“(L)und ”Bunny“(B).DieBe- schriftungen auf den Balken beziehen sich auf die verwendete Octree-Tiefe,

wobei

”o“fürdieSichtbarkeitenbestimmungohneOctreesteht.Dasjeweils beste Ergebnis ist hervorgehoben.

Umsetzung auf der CPU

Auf der CPU kann durch die Kugelvorberechnung in den meisten Fällen etwas Zeit eingespart werden. Abbildung 6.1 zeigt eine Gegenüberstellung der Sichtbarkeitenbestimmung mithilfe von Vorberechnungen an bis zu sechs Kugeln für alle vier Testszenen. Dabei wird deutlich, dass zwar auch die Verwendung mehrerer Kugeln zur Beschleunigung beiträgt, zumeist jedoch mit nur einer Kugel das beste Ergebnis erzielt wird. Dass dies nicht immer zutrifft, zeigen die Ergebnisse an der

”RadiosGL-Testszene“(R),wokeinederVorberechnungenzueiner

Beschleunigung führt. In Bezug auf die hohe Gesamtlaufzeit der Sichtbarkeitenberechnung fallen die erzielten Beschleunigungswerte mit weniger als Faktor [1],[2] eher gering aus. Eine Steigerung der Komplexität darzustellender Szenen ist damit nicht möglich.

Weitaus bessere Ergebnisse können dagegen mithilfe einer Octree-Struktur erzielt werden. Wie in Abbildung [6].[2] dargestellt, sind damit je nach Aufbau der Szene Beschleunigungen bis zu Faktor [5] möglich. In den meisten Fällen erweist sich dabei ein Octree der Tiefe [4] als vorteilhaft. Nur in der für Octree-Strukturen ungünstigen Testszene

”Bunny“(B),indersich

ein Objekt mit sehr vielen Polygonen in einer Ecke eines größtenteils leeren Raumes befindet, ist es günstiger höhere Octree-Tiefen zu nutzen. Hier wird mit einem Octree der Tiefe [6] das beste Resultat erzielt. Im Vergleich zu den anderen Testszenen fällt die gegenüber einer Berechnung ohne Octree erzielte Beschleunigung jedoch deutlich geringer aus.

Umsetzung auf der GPU

Im Gegensatz zur CPU bewirkt die Umsetzung der Vorberechnung mittels Kugeln auf der GPU eine Verschlechterung der Laufzeit an allen getesteten Szenen. Ein Vergleich der Lauf- zeiten der Sichtbarkeitenbestimmung ohne Beschleunigungsansatz und mit Vorberechnung an einer sowie an zwei Kugeln ist in Tabelle [6].[2] dargestellt. Dementsprechend verschlechtert sich die benötigte Rechenzeit bei Verwendung einer Kugel um etwa [19] % auf der Grafikkarte von NVIDIA und um durchschnittlich [38] % auf der von ATI. Bei Verwendung von zwei Ku- geln beträgt die Verlangsamung sogar durchschnittlich [46] % (auf NVIDIA) bzw. [45] % (auf ATI). Diese Verschlechterung der Laufzeit wird von der durch die Vorberechnung zusätzlich benötigten IF-Bedingung im Shader-Programm verursacht. Wie in Kapitel 3 erwähnt, haben die GPUs der aktuellen Grafikkarten häufig Nachteile beim Verarbeiten von Verzweigungen. Bei den hier verwendeten Modellen stellt es sich deshalb als günstiger heraus, weniger Ver- zweigungen zu verwenden und dafür nicht benötigte Ergebnisse am Ende der Berechnung zu verwerfen. So kann z. B. Rechenzeit eingespart werden, indem der Algorithmus des Strecke- Dreieck-Schnitttests aus Abschnitt 4.2.1 immer alle Parameterwerte berechnet und sämtliche Schnittbedingungen in einer einzigen IF-Bedingung bündelt.

Weitaus schlechter als die Vorberechnung mittels Kugeln wirkt sich die Verwendung eines Octree auf GPU-Laufzeit aus. Zum Vergleich ist das Resultat der Sichtbarkeitenbestimmung mit Octree-Unterstützung bei der jeweils günstigsten Octree-Tiefe ebenfalls in Tabelle 6.2 angegeben. Die erzielten Zeiten liegen dabei fast im Bereich der CPU-Laufzeiten. Ursache für dieses Ergebnis ist hauptsächlich der Art der GPU-Umsetzung zuzuschreiben. So verbleibt die Haltung und Auswertung des Octree komplett auf der CPU und dem Schnitttest auf der GPU werden entsprechend weniger Flächenelemente übergeben. Leider entsteht bei der im- plementierten Variante so viel Datenverkehr zwischen CPU und GPU, dass die Laufzeit der Berechnung ohne Octree um Längen überschritten wird. Einen Ausweg bietet möglicherweise eine Implementierung des kompletten Octree auf der GPU. Aufgrund des relativ hohen Ent-

Abbildung in dieser Leseprobe nicht enthalten

Tabelle 6.2: Vergleich der Vorschläge zur Beschleunigung der Sichtbarkeitenbe- stimmung auf der GPU an den weniger unterteilten Varianten der Szenen

”CornellBox“(CB), ”RadiosGL-Testszene“(R), ”Livingroom“(L)und ”Bun- ny“ (B). Der jeweils verwendete Beschleunigungsansatz ist in Spalte [2] angegeben, wobei die Kugeln für die Vorberechnung mit entsprechender Anzahl an Kugeln beim Strecke-Dreieck-Schnitttest stehen. Für die Berechnung mit OctreeUnterstützung wurde die jeweilig optimale Tiefe verwendet.

wicklungsaufwands war eine solche Umsetzung im Rahmen dieser Arbeit nicht möglich. Mit der Realisierung von Octrees in Grafikhardware haben sich bereits Benson und Davis [BD[02]], Lefebvre et al. [LHN[05]] sowie Kniss et al. [KLS+05] auseinander gesetzt. Aber selbst bei der Verwirklichung des Octree auf der GPU bleibt fraglich, ob sich der zusätzliche Aufwand für Texturzugriffe usw. nicht eher negativ auf die Laufzeit der Sichtbarkeitenbestimmung auswirkt. Zudem könnte die bei höherer Anzahl an Flächenelementen notwendige Partitionierung der Berechnung zu zusätzlichen Problemen führen.

6.1.3 Vergleich von Kreisscheibenapproximation und Halbwürfelverfahren

Bei den bisher dargestellten Ergebnissen erfolgte die Bestimmung der Formfaktoren durch Approximation mittels Kreisscheiben, vgl. Abschnitt 2.3.3. Um einen Vergleich zum häufig verwendeten Halbwürfelverfahren (siehe Abschnitt 2.3.2) herstellen zu können, wurde eine Open Source-Implementierung dieses Verfahrens in Radioso integriert. Diese Implementierung stellt eine reine CPU-Umsetzung dar und kommt ohne die Verwendung von Beschleunigungstechniken durch Grafikhardware aus.

Für einen Laufzeitvergleich der beiden Verfahren wurden die Formfaktoren mithilfe von Halbwürfeln der Pixelbreiten 50, 100, 200, 500 und 1000 bestimmt. Die dafür benötigten Rechenzeiten sind für die Testszene

”Bunny“inAbbildung[6].[3]zusammengestellt.Darinwird

deutlich, dass die Berechnung mit zunehmender Pixelanzahl des Halbwürfels und somit mit zunehmender Genauigkeit der zu berechnenden Formfaktoren stark ansteigt. Liegt die Lauf- zeit bei geringer Pixelbreite des Halbwürfels noch deutlich unter der mittels Kreisscheiben- approximation auf der CPU erreichten Bestzeit, so wird diese mit zunehmender Pixelzahl überschritten. Die benötigte Rechenzeit der GPU-Implementierung der Kreisscheibenappro- ximation bleibt dagegen stets klar unter den Laufzeiten der Halbwürfelmethode. Ein ähnliches Bild ergibt sich für die anderen drei Testszenen. Dort wird die Laufzeit der Kreisscheibenap- proximation auf der CPU aufgrund der besseren Octree-Eigenschaften der Szenen bereits bei

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 6.3: Laufzeitvergleich der Verfahren zur Formfaktorbestimmung. Das Diagramm zeigt die zur Formfaktorbestimmung benötigten Laufzeiten in

Abhängigkeit der Pixelbreite des Halbwürfels an der Szene

”Bunny“([1100]

Flächenelemente nach der Unterteilung). Als Vergleichswerte wurde die jeweils besten Laufzeiten der Approximation der Formfaktoren mittels Kreisscheiben (inklusive Sichtbarkeitenbestimmung) eingetragen.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 6.4: Qualitativer Vergleich der Verfahren zur Formfaktorbestimmung

an der Szene

”CornellBox“.(a)UnterVerwendungdesHalbwürfelverfahrens

entsteht im unteren Drittel der Wände ein dunkler Bereich, hervorgerufen durch die Projektion dieser Flächenelemente auf die Kanten des Halbwürfels.

(b) Bei zu grober Unterteilung der Flächenelemente sind die mit der Kreis- scheibenapproximation ermittelten Formfaktoren zu ungenau, was in diesem Fall nach [100] Iterationen mit dem Jacobi-Verfahren zur Überstrahlung durch den kleineren Quader führt, der RGB-Werte nahe Eins besitzt. (c) Mit nur

[10] Iterationen erscheint die selbe Szene dagegen normal beleuchtet.

geringer Pixelbreite des Halbwürfels überschritten. Die Umsetzung der Kreisscheibenapproximation auf der GPU ist in jedem Fall merklich schneller und benötigt nur etwa ein Siebtel bis ein Zehntel der mit dem 50 Pixel breiten Halbwürfel erreichten Laufzeit.

Ein qualitativer Nachteil entsteht dem Halbwürfelverfahren bei der Projektion von Flächen- elementen auf die Kanten des Halbwürfels. Durch die Verzerrung der Projektion in diesem Bereich können für die entsprechenden Formfaktoren ungenaue Werte ermittelt werden. Die- ser Nachteil tritt durch die Würfelgestalt und der an der Decke zentrierten Lichtquelle bei der

Testszene

”CornellBox“besondersstarkinErscheinung,sieheAbbildung[6].[4]a.Einsolches

Verhalten kann bei der Approximation der Formfaktoren mittels Kreisscheiben ausgeschlos- sen werden. Allerdings besteht hier die Gefahr, dass aufgrund zu ungenauer Näherungswerte die Summe der Formfaktoren zu einem Flächenelement nicht nahe genug bei Eins liegt, sodass die Szene durch zu viele Iterationen beim Lösen des Radiosity-Gleichungssystems möglicher- weise zu dunkel bzw. überstrahlt erscheint (siehe Abbildung [6].[4]b). Die Genauigkeit eines Formfaktors bei der Näherung mittels Kreisscheiben kann durch eine feinere Unterteilung der Szene oder durch die in Abschnitt [2].[3].[3] beschriebene zusätzliche Unterteilung des bestrahlten Flächenelements und der damit verbundenen Zusammensetzung aus mehreren Kreisscheiben erhöht werden.

6.1.4 Zusammenfassung

Zusammenfassend lässt sich feststellen, dass mithilfe der GPU-Programmierung in vielen Fällen hohe Leistungssteigerungen möglich sind. Gut geeignet sind dabei vor allem Aufga- ben, die einen hohen Rechenaufwand verursachen, möglichst wenig Datenaustausch zwischen CPU und GPU verlangen und darüber hinaus so parallelisierbar sind, dass ihre Teilaufgaben unabhängig voneinander berechnet werden können. Außer bei der Unterteilung der Flächen- elemente konnte dies am Radiosity-Verfahren recht erfolgreich umgesetzt werden. So wurden an den untersuchten Testszenen gegenüber der CPU-Implementierung Beschleunigungsfak- toren von bis 70 bei der Gesamtlaufzeit bzw. bis zu 150 bei einzelnen Teilaufgaben erreicht. Dies übertrifft sogar die mit dem Halbwürfelverfahren (auf der CPU) erzielten Laufzeiten bei Weitem.

ÜberdieskonntedurchdieVerwendungeinerneuerenVersiondesNVIDIAForceWare- Treibers (Version 169.21) eine weitere Verbesserung der Laufzeiten auf der GPU festgestellt werden. Diese Verbesserung macht sich mit einem Leistungsanstieg von 3 bis 8 % vor allem an den grob unterteilten Testszenen bemerkbar. An den fein unterteilten Testszenen mit we- sentlich mehr Flächenelementen fällt sie dagegen kaum ins Gewicht. Hier liegt die Vermutung nahe, dass durch die neue Treiberversion primär die zum Kompilieren der Shader-Program- me benötigte Zeit reduziert wird. Eine Ausnahme bildet die Sichtbarkeitenbestimmung mit Octree-Unterstützung, bei der sogar eine Verschlechterung der Laufzeit beobachtet werden kann.

Des Weiteren wurde beobachtet, dass der durch die RapidMind-Plattform generierte GLSL-Quelltext durch redundante Anweisungen und unnötige Variablen nicht immer op- timal zu sein scheint. Es besteht also die Möglichkeit, mit einem verbesserten RapidMind- Release weitere Vorteile in Bezug auf die benötigte Rechenzeit der GPU zu erzielen. Diese Vermutung ließ sich mit dem neueren Release 2.1 der RapidMind-Plattform jedoch nur teil- weise bestätigen. Zwar werden bei der Sichtbarkeitenbestimmung mit Octree-Unterstützung Beschleunigungswerte von 1 bis 4 % erreicht, allerdings verschlechtern sich dafür die Laufzei- ten der anderen Algorithmen marginal. Dabei macht es keinen Unterschied, ob die Berech-nung auf der Grafikkarte von NVIDIA oder der von ATI stattfindet. Im Unterschied zum Release 2.0.1 enthält der generierte GLSL-Quelltext nun zusätzliche Variablen zur Steuerung von Array- bzw. Texturzugriffen. Eventuell tragen zum geänderten Laufzeitverhalten weitere Neuerungen innerhalb der RapidMind-Plattform bei.

In Bezug auf die ursprünglich gemessenen Laufzeiten verringert sich im Zusammenspiel des neueren RapidMind-Release 2.1 mit der neueren Treiberversion 169.21 des NVIDIA ForceWare-Treibers die Berechnungsdauer der Unterteilung der Flächenelemente sowie die der Sichtbarkeitenbestimmung mit Octree-Unterstützung. Bei allen anderen Algorithmen tritt dagegen eine Verschlechterung der Laufzeit ein. Es ist also sehr wahrscheinlich, dass die präsentierten GPU-Laufzeiten im Zuge neuer Grafiktreiber- und RapidMind-Versionen auch in Zukunft leichten Schwankungen unterliegen werden.

Einige Probleme der Programmierung von GPU-Anwendungen wurden bereits in Kapitel 5 diskutiert. Dabei handelt es sich um die beschränkte Anzahl an Schleifeniterationen sowie Beschränkungen in der Größe von Texturen, Arrays und dem Grafikspeicher. Nähere De- tails und hierzu erarbeitete Lösungsansätze können dem Abschnitt 5.2 entnommen werden. Darüber hinaus sind während der Entwicklung von Radioso weiter Probleme aufgetreten, de- ren Ursache vermutlich dem relativ frühen Entwicklungsstadium der GPU-Programmierung zuzuschreiben ist.

So treten in mehreren Versionen des ATI-Grafiktreibers Catalyst unterschiedliche Fehler beim Ausführen der Shader-Programme auf. In den Versionen 7.1 und 7.2 fängt während der Berechnung die Bildschirmausgabe an zu flackern, anschließend startet der Rechner neu. In den Versionen 7.5 bis 7.7 gibt es Probleme bei der Verwaltung der Nutzerrechte. Dabei bricht ein Programm beim Übergang zur Berechnung auf der GPU aus unerklärlichen Gründen ab, wenn es unter einem Benutzerkonto mit eingeschränkten Rechten ausgeführt wird. Das selbe Programm läuft dagegen beim Ausführen mit Administratorrechten ohne Probleme. Hinzu kommt, dass in Version 7.6 einige Berechnungen ungenau durchgeführt werden. Eine Anwendung, die unter den anderen Treiberversionen die erwarteten Ergebnisse liefert, weist unter der Version 7.6 Abweichungen in den Resultaten auf. Dies hat zur Folge, dass das entwickelte GPU-Programm auf der verwendeten ATI-Grafikkarte nur unter Treiberversion 7.4 korrekt arbeitet (höhere Versionen als 7.7 wurden nicht getestet). Für ein marktfähiges Produkt ist dieser Umstand kaum vorstellbar, da die Besitzer einer solchen ATI-Grafikkarte gezwungen werden, eine bestimmte Treiberversion zu installieren.

In Kapitel 3 wurde bereits angedeutet, dass bei der Entwicklung von GPU-Anwendungen jede zusätzliche ÜbersetzungsschichteineweiterepotenzielleFehlerquelledarstellt.Soist auch die verwendete RapidMind-Plattform nicht frei von Fehlern, wie das folgende Beispiel zeigt. Verwendet man in der früheren Version 1.9 der Plattform in einem Program-Objekt mehrere IF-Bedingungen innerhalb einer WHILE-Schleife, führt dies zu einer Endlosschleife im generierten GLSL-Code. Dieser Fehler wurde mit dem Release 2.0.1 der Plattform beseitigt.

Des Weiteren kann bei Ausführung von GPU-Berechnungen beobachtet werden, dass die Grafik- und die Soundausgabe gestört werden und sich die Arbeit des Systems auf die Ausführung des GPU-Task beschränkt. Die Ursache dafür liegt vermutlich beim OpenGLThread, der das System durch aktives Warten auf die GPU-Antwort zu 100 % auslastet. Auch dieser Umstand ist bei kommerziellen Produkten nicht tolerierbar.

7 Resümee und Ausblick

Im Rahmen dieser Arbeit wurde untersucht, inwieweit sich das Radiosity-Verfahren mithilfe programmierbarer Grafikhardware und dem Einsatz von GPGPU-Werkzeugen beschleuni- gen lässt. Für einen Laufzeitvergleich gegenüber herkömmlichen Implementierungsmethoden wurden die einzelnen Schritte des Verfahrens sowohl auf der GPU als auch auf der CPU umgesetzt. Die RapidMind-Entwicklungsplattform diente dabei als Werkzeug für die Pro- grammierung der Grafikhardware. Durch das Verbergen von grafikspezifischen Details und der entwicklerfreundlichen API konnte der Aufwand zur Realisierung der Rechenschritte auf der GPU somit gering gehalten werden. Dennoch zwingen die Einschränkungen der Grafikhard- ware und das spezielle Programmiermodell zu Anpassungen an den verwendeten Algorithmen. Hinzu kommt, dass nicht jeder Algorithmus für eine Implementierung auf der GPU geeignet ist, wie es das Beispiel der Maximumbestimmung und die Unterteilung der Flächenelemente zeigen.

Dass das Radiosity-Verfahren insgesamt gut für eine Implementierung auf der GPU geeignet ist, belegen die Ergebnisse dieser Arbeit mit erzielten Beschleunigungsfaktoren von bis zu 70 gegenüber der Berechnung auf der CPU. In einzelnen Schritten des Verfahrens können die benötigten Rechenzeiten sogar um das 150fache beschleunigt werden. Verglichen mit den Laufzeiten der GPU-Implementierungen von Coombe et al. [CHL04] und Borsdorf et al. [BDS05] ist der hier umgesetzte Algorithmus jedoch deutlich langsamer. Dies ist mit den unterschiedlichen Herangehensweisen zu begründen. So verfolgen sowohl Coombe et al. als auch Borsdorf et al. einen progressiven Ansatz des Radiosity-Verfahrens und nutzen darüber hinaus zur Bestimmung der Sichtbarkeiten eine stereographische Projektion, die wesentlich schneller als der in dieser Arbeit verwendete Ansatz ist. Zudem wird bei beiden das gesamte Verfahren vollständig auf der GPU berechnet. Coombe et al. führen in ihrer Arbeit einen Vergleich ihrer Implementierung zum kommerziellen Programm Lightscape an. Letzteres setzt einen vergleichbaren Algorithmus auf der CPU um und benötigt für die Berechnung an einer Vergleichsszene 96 Sekunden. Die GPU-Implementierung von Coombe et al. ist unter gleichen Kriterien mit benötigten 2,4 Sekunden um ein Vielfaches schneller. Aufgrund des hohen Speicherverbrauchs des Radiosity-Verfahrens stellt der progressive Ansatz für eine praxisrelevante Implementierung die bessere Alternative dar. Insgesamt betrachtet verfügt der im Rahmen dieser Arbeit umgesetzte Algorithmus somit auch auf der verwendeten Hardware über ein deutliches Verbesserungspotenzial.

Ein Nachteil von GPGPU-Werkzeugen ist, dass der Programmierer keinen direkten Ein- fluss auf den generierten Quelltext in einer der Shader-Hochsprachen hat und deshalb bei Optimierungen und Korrektheit auf die Vollkommenheit des verwendeten Werkzeuges ange- wiesen ist. Dies wirkt sich besonders bei Übersetzungsfehlern negativ auf die Entwicklung von GPU-Anwendungen aus. Dass auch die Technik programmierbarer Grafikhardware noch nicht komplett ausgereift ist, zeigen die Fehler und das unterschiedliche Leistungsverhalten verschiedener Grafiktreiber. Demnach ist davon auszugehen, dass sich der erzielte Leistungs- gewinn im Zuge neuer Grafiktreiber und verbesserter Programmierwerkzeuge weiter steigern lässt.

Die nächste Generation an Grafikkarten unterstützt Microsofts Direct3D 10-API (Teil von Microsofts DirectX Multimedia-API). Blythes hat die wichtigsten Änderungen und Merkmale dieser Hard- und Software in [Bly06] zusammengefasst. Die Auswirkungen der neu- en Hardware auf die GPGPU-Programmierung werden aber erst nach und nach wahrzuneh- men sein, wenn die Entwickler ihre Anwendungen auf die neuen Funktionalitäten umstellen und die Leistungsfähigkeit der neuen GPUs entdecken. Für GPGPU-Berechnungen werden nach [OLG+07] vermutlich die folgenden Neuerungen von besonderem Interesse sein:

- DirectX 10 führt eine neue programmierbare Einheit in die Grafikpipeline ein, den Geometrie-Shader, welcher in der Pipeline hinter dem Vertex-Shader platziert ist. Als Eingaben dienen dem Geometrie-Shader vollständige Primitive. Die größten Unterschie- de zwischen ihm und den bisherigen Vertex- und Fragment-Shadern ist ein flexiblerer Speicherzugriff und die Ausgabe beliebig vieler Primitive. Diese Fähigkeit, prozedural neue Elemente zu erzeugen, wird wahrscheinlich nicht nur bei Grafikaufgaben wie der Schattenvolumenberechnung sondern auch bei allgemeinen Aufgaben breiten Nutzen hervorrufen.

- GPGPU-Entwickler bekommen die lange gewünschten flexibleren Operationen auf dem Grafikspeicher. Während Operationen, wie Render-to-Texture und Render-to-Vertex- Array diese Wünsche bereits teilweise erfüllen, verspricht die kommende DirectX 10- Hardware sowohl größere Funktionalität als auch mehr Leistungsfähigkeit auf diesem Gebiet. Eine neue Funktion dieser Hardware wird der Stream-Output sein, welcher dem Geometrie-Shader erlaubt, seine Ausgaben direkt im Grafikspeicher abzulegen.

- Das neue Shader-Model 4.0 vereinheitlicht in Verbindung mit der DirectX 10-Hardwa- re den Grundbefehlssatz programmierbarer Shader-Einheiten. Hatte bisher noch jeder programmierbare Shader seine stufenspezifischen Eigenheiten, so wird in Zukunft für Vertex-, Geometrie- und Fragment-Shader die gleiche Hardware genutzt. Zudem soll eine Reihe von Shader-Beschränkungen, wie Befehlsanzahl, Registerraum und Render- ziel gelockert werden. Auch 32-Bit-Integer-Werte werden nun unterstützt, was sowohl den aktuellen (vor allem bei der Speicheradressierung) als auch neuen GPGPU-An- wendungen wie Kryptografie zu Gute kommt. Darüber hinaus ist zu erwarten, dass die Genauigkeit von Fließkommaoperationen in den kommenden Hardwaregenerationen auf 64 Bit erhöht wird.

Da der Radiosity-Algorithmus die starke Parallelisierung der Grafikpipeline gut ausnutzt, wird im Zuge der stetigen Verbesserungen und Leistungssteigerungen der Grafikkarten auch die Leistungsfähigkeit des Radiosity-Verfahrens weiter zunehmen, sodass in Zukunft mithilfe der GPU-Programmierung sogar eine interaktive Berechnung des Algorithmus möglich er- scheint.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung A.1: Zu den Laufzeitvergleichen herangezogene Testszenen.

A Gerenderte Szenen

B Spezifikationsschemata und Beispiele zu den Szenenbeschreibungsformaten

Das im Rahmen dieser Diplomarbeit implementierte Programm Radioso akzeptiert Szenen- beschreibungen in den Formaten MDL und RML. Nähere Informationen zu beiden Formaten können dem Abschnitt 5.3 entnommen werden. Zur Orientierung beim Erstellen von Szenen hat sich die folgende schemenhafte Darstellung von Szenenbeschreibungssprachen bewährt.[11] Im Anschluss daran finden sich zwei Beispiele zu den hier beschriebenen Formaten.

Abbildung in dieser Leseprobe nicht enthalten

B.1 MDL-Schema

Abbildung in dieser Leseprobe nicht enthalten

B Spezifikationsschemata und Beispiele zu den Szenenbeschreibungsformaten

B.2 RML-Schema

Abbildung in dieser Leseprobe nicht enthalten

B.3 Beispiele zu den Szenenbeschreibungsformaten

Abbildung in dieser Leseprobe nicht enthalten

B.3 Beispiele zu den Szenenbeschreibungsformaten

Abbildung in dieser Leseprobe nicht enthalten

Quelltext B.1: Beipiel einer Szenenbeschreibung in MDL

B Spezifikationsschemata und Beispiele zu den Szenenbeschreibungsformaten

Abbildung in dieser Leseprobe nicht enthalten

Quelltext B.2: Beipiel einer Szenenbeschreibung in RML

C Beispiel eines generierten GLSL-Quelltextes

Der folgende Quelltext zeigt die von der RapidMind-Entwicklungsplattform generierten GLSL-Shader-Programme zur Parallelisierung des Beispiels aus Abschnitt 3.3.5:

Abbildung in dieser Leseprobe nicht enthalten

Quellen

[app] Apple Developer Connection: Graphics & Imaging - OpenGL. http://developer.

apple.com/graphicsimaging/opengl, Abruf: 11. Januar 2008

[ati06] ATI: CTM Guide. Version: 2006. http://ati.amd.com/companyinfo/researcher/

documents/ATI CTM Guide.pdf, Abruf: 10. Januar 2008

[awa] Awaron AG: Produktbeschreibung tucan radiosity. http://www.awaron.de/de/

products/tucan/tucan radiosity.asp, Abruf: 23. November 2007

[Bax] Baxter, W.: Image Debugger Homepage. http://www.billbaxter.com/projects/

imdebug, Abruf: 11. Januar 2008

[BD02] Benson, D. ; Davis, J.: Octree Textures. In: Proceedings of the 29th annual conference on Computer graphics and interactive techniques (2002), S. 785-790

[BDS05] Borsdorf, A. ; Dachsbacher, C. ; Stamminger, M.: Progressive Radiosi- ty auf programmierbarer Grafikhardware. In: Workshop Hardware for Visual Computing, Tübingen (2005)

[BEL+06] Boulos, S. ; Edwards, D. ; Lacewell, J.D. ; Kniss, J. ; Kautz, J. ; Shirley,

P. ; Wald, I.: Interactive Distribution Ray Tracing / SCI Institute, University of Utah. 2006 (UUSCI-2006-022). - Forschungsbericht

[Ber] Berenguier, L.: RealTime Radiosity on GPU. http://berengui.club.fr/rtrad.

html, Abruf: 25. Februar 2008

[BFH+] Buck, I. ; Foley, T. ; Horn, D. ; Sugerman, J. ; Hanrahan, P. ; Hou- ston, M. ; Fatahalian, K.: BrookGPU Homepage. http://graphics.stanford. edu/projects/brookgpu, Abruf: 9. Januar 2008

[BFH+04] Buck, I. ; Foley, T. ; Horn, D. ; Sugerman, J. ; Fatahalian, K. ; Houston,

M. ; Hanrahan, P.: Brook for GPUs: stream computing on graphics hardware. In: ACM Transactions on Graphics (TOG) 23 (2004), Nr. 3, S. 777-786

[ble] Blender Hompage. http://www.blender.org, Abruf: 23. November 2007

[Bly06] Blythe, D.: The Direct3D 10 System. In: ACM Transactions on Graphics

(TOG) 25 (2006), Nr. 3, S. 724-734

[Buc] Buck, I.: SUPERCOMPUTING 2006 GPGPU Tutorial: NVIDIA CUDA.

http://www.gpgpu.org/sc2006/slides/08b.buck.cuda.pdf, Abruf: 10. Januar 2008

[bug] BuGLe: an OpenGL debugging tool. http://www.opengl.org/sdk/tools/BuGLe,

Abruf: 11. Januar 2008

Quellen

[Car] Carlson, W.: History of the Cornell University Program of Computer Graphics. http://design.osu.edu/carlson/history/tree/cornell.html, Abruf: 26. Februar 2008

[CCWG88] Cohen, M.F. ; Chen, S.E. ; Wallace, J.R. ; Greenberg, D.P.: A Progressive Refinement Approach to Fast Radiosity Image Generation. In: Proceedings of the 15th annual conference on Computer graphics and interactive techniques (1988), S. 75-84

[CG85] Cohen, M.F. ; Greenberg, D.P.: The Hemi-Cube: A Radiosity Solution for Complex Environments. In: Proceedings of the 12th annual conference on Com- puter graphics and interactive techniques (1985), S. 31-40

[CH05] Coombe, G. ; Harris, M.J.: Global Illumination Using Progressive Refinement Radiosity. In: GPU Gems II (2005), S. 635-647

[CHH03] Carr, N.A. ; Hall, J.D. ; Hart, J.C.: GPU Algorithms for Radiosity and Sub- surface Scattering. In: Proceedings of the ACM SIGGRAPH/EUROGRAPHICS conference on Graphics hardware (2003), S. 51-59

[CHL04] Coombe, G. ; Harris, M.J. ; Lastra, A.: Radiosity on Graphics Hardware. In: Proceedings of the 2004 conference on Graphics interface (2004), S. 161-168

[Col07] Colehour, C.: Google Blogoscoped: Google Acquired PeakStream. Version: Juni 2007. http://blogoscoped.com/archive/2007-06-06-n38.html, Abruf: 10. Januar 2008

[CW93] Cohen, M.F. ; Wallace, J.R.: Radiosity and Realistic Image Synthesis. Morgan Kaufmann, 1993

[DH93] Deuflhard, P. ; Hohmann, A.: Numerische Mathematik I: Eine algorithmisch orientierte Einführung. de Gruyter, 1993

[Els93] Elsner, N.: Grundlagen der Technischen Thermodynamik 2: Wärmeübertra-gung. Akademie Verlag, 1993

[FDFH90] Foley, J.D. ; Dam, A. van ; Feiner, S.K. ; Hughes, J.F.: Computer Graphics: Principles and Practice. Second Edition. Addison-Wesley, 1990

[FLCG] Fritz, N. ; Lucas, P. ; Cullman, C. ; Gebhard, G.: CGiS Homepage. http:// rw4.cs.uni-sb.de/projects/CGiS/, Abruf: 10. Januar 2008

[GHJV04] Gamma, E. ; Helm, R. ; Johnson, R. ; Vlissides, J.: Entwurfsmuster: Elemente wiederverwendbarer objektorientierter Software. Pearson Education Deutschland, 2004

[gli] GLIntercept Homepage. http://glintercept.nutty.org, Abruf: 11. Januar 2008

[gol07] Golem.de: Google kauft Multi-Core-Spezialisten PeakStream. Version: Juni 2007.

http://www.golem.de/0706/52695.html, Abruf: 10. Januar 2008

[gpga] GPGPU Community Homepage. http://www.gpgpu.org, Abruf: 30. Dezember

Quellen

[gpgb] GPGPU.org Wiki. http://www.gpgpu.org/wiki, Abruf: 12. Januar 2008

[gra] Graphic Remedy Homepage. http://www.gremedy.com, Abruf: 11. Januar 2008

[GTGB84] Goral, C.M. ; Torrance, K.E. ; Greenberg, D.P. ; Battaile, B.: Modeling the Interaction of Light Between Diffuse Surfaces. In: Proceedings of the 11th annual conference on Computer graphics and interactive techniques (1984), S. 213-222

[GWL+05] Goodnight, N. ; Woolley, C. ; Lewin, G. ; Luebke, D. ; Humphreys, G.: A Multigrid Solver for Boundary Value Problems Using Programmable Graphics Hardware. In: International Conference on Computer Graphics and Interactive Techniques (2005)

[har] HardenedCriminal software: RadiosGL Hompage. http://hcsoftware.sourceforge.

net/RadiosGL/RadiosGL.html, Abruf: 23. November 2007

[Hau96] Haußecker, H.: Messung und Simulation von kleinskaligen Austauschvorgän-

gen an der Ozeanoberfläche mittels Thermographie, Ruprecht-Karls-Universität Heidelberg, Diss., 1996

[haw] Hawk Software: GLTrace Programming Utility. http://www.hawksoft.com/

gltrace, Abruf: 11. Januar 2008

[HB05] Kapitel 34: GPU Flow Control Idioms. In: Harris, M. ; Buck, I.: GPU Gems

2. Addison Wesley, 2005, S. 547-555

[Hen] Hensley, J.: SIGGRAPH 2007 GPGPU Course: AMD CTM. http://www.

gpgpu.org/s2007/slides/07-CTM-overview.pdf, Abruf: 10. Januar 2008

[HJ03] Harris, M.J. ; James, G.: Simulation and Animation Using Hardware Acce-

lerated Procedural Textures. In: Proceedings of Game Developers Conference (2003)

[Houa] Houston, M.: SIGGRAPH 2007 GPGPU Course: General-Purpose Computati-on on Graphics Hardware. http://www.gpgpu.org/s2007/slides/01-introduction. pdf, Abruf: 30. Dezember 2007

[Houb] Houston, M.: SIGGRAPH 2007 GPGPU Course: High Level Languages for

GPUs. http://www.gpgpu.org/s2007/slides/05-languages-overview.pdf, Abruf:

9. Januar 2008

[Jac] Jacobson, F.: Autodesk - Discussion Groups - Radiosity Question. http://

discussion.autodesk.com/thread.jspa?threadID=449147, Abruf: 26. Februar 2008

[KF05] Kapitel 30: The GeForce 6 Series GPU Architecture. In: Kilgariff, E. ; Fer-

nando, R.: GPU Gems 2. Addison Wesley, 2005, S. 471-491

[Kil03] Kilgard, M.J.: Cg in Two Pages. Version: Januar 2003. http://xxx.lanl.gov/

ftp/cs/papers/0302/0302013.pdf, Abruf: 12. Januar 2008

[KL04] Karlsson, F. ; Ljungstedt, C.J.: Ray tracing fully implemented on program-

mable graphics hardware, Chalmers University of Technology, Göteborg, Diplomarbeit, 2004

[KLS+05] Kniss, J. ; Lefohn, A. ; Strzodka, R. ; Sengupta, S. ; Owens, J.D.: Oc- tree Textures on Graphics Hardware. In: International conference on Computer graphics and interactive techniques (2005)

[LHN05] Lefebvre, S. ; Hornus, S. ; Neyret, F.: Octree Textures on the GPU. In: GPU Gems II (2005), S. 595-613

[LKHW03] Lefohn, A.E. ; Kniss, J.M. ; Hansen, C.D. ; Whitaker, R.: Interactive De- formation and Visualization of Level Set Surfaces Using Graphics Hardware. In: IEEE Visualization 82 (2003)

[McC06] McCool, M.D.: Data-Parallel Programming on the Cell BE and the GPU Using the RapidMind Development Platform. In: GSPx Multicore Applications Confe- rence, Santa Clara, 2006

[mic] Microsoft Developer Network: Shader Debugger. http://msdn.microsoft.com/archive/default.asp?url=/archive/en-us/directx9 c dec 2004/directx/graphics/ Tools/ShaderDebugger.asp, Abruf: 11. Januar 2008

[MM05] Montrym, J. ; Moreton, H.: The GeForce 6800. In: IEEE Micro 25 (2005), Nr. 2, S. 41-51

[MT97] Möller, T. ; Trumbore, B.: Fast, Minimum Storage Ray-Triangle Intersection. In: Journal of Graphics Tools 2 (1997), Nr. 1, S. 21-28

[nvia] NVIDIA: CUDA Homepage. http://www.nvidia.com/object/cuda home.html,

Abruf: 11. Januar 2008

[nvib] NVIDIA: CUDA Programming Guide 1.1. http://developer.download.nvidia.com/compute/cuda/1 1/NVIDIA CUDA Programming Guide 1.1.pdf, Abruf: 10. Januar 2008

[OLG+07] Owens, J.D. ; Luebke, D. ; Govindaraju, N. ; Harris, M. ; Krüger, J. ; Lefohn, A.E. ; Purcell, T.J.: A Survey of General-Purpose Computation on Graphics Hardware. In: Computer Graphics Forum 26 (2007), Nr. 1, S. 80-113

[ope] OpenGL.org Wiki: Shading languages: GLSL. http://www.opengl.org/wiki/index.php/Shading languages: GLSL, Abruf: 12. Januar 2008

[Pap] Papakipos, M.: SUPERCOMPUTING 2006 GPGPU Tutorial: PeakStream Platform. http://www.gpgpu.org/sc2006/slides/12.papakipos.peakstream.pdf, Ab- ruf: 10. Januar 2008

[PBMH02] Purcell, T.J. ; Buck, I. ; Mark, W.R. ; Hanrahan, P.: Ray Tracing on Programmable Graphics Hardware. In: Proceedings of the 29th annual conference on Computer graphics and interactive techniques (2002), S. 703-712

[PL05] Pohl, R.O. ; Lüders, K.: Pohls Einführung in die Physik 2: Elektrizitätslehre und Optik. Springer, 2005

[POAU00] Peercy, M.S. ; Olano, M. ; Airey, J. ; Ungar, P.J.: Interactive Multi-Pass Programmable Shading. In: Proceedings of the 27th annual conference on Com- puter graphics and interactive techniques (2000), S. 425-432

[pov] POV-Ray Homepage. http://www.povray.org, Abruf: 14. Februar 2008

[Pow] Power, K.: Graphics Notes: Radiosity. http://glasnost.itcarlow.ie/∼powerk/ Graphics/Notes/node13.html, Abruf: 18. Dezember 2007

[rapa] RapidMind Developer Portal. https://developer.rapidmind.net, Abruf: 19. Januar

[rapb] RapidMind Developer Portal - Forums - Development - Definitive array and element sizes. https://developer.rapidmind.net/forums/development/329059068/# 660797769, Abruf: 25. März 2008

[Rau93] Rauber, T.: Algorithmen in der Computergraphik. Teubner, 1993

[RR00] Rauber, T. ; Rünger, G.: Parallele und verteilte Programmierung. Springer, 2000

[SF01] Segura, R.J. ; Feito, F.R.: Algorithms to Test Ray-Triangle Intersection. Comparative Study. In: 9th International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision, Plzen, Czech Republic (2001)

[sh] Sh Homepage. http://libsh.org, Abruf: 12. Januar 2008

[shaa] Shadesmith Homepage. http://graphics.stanford.edu/projects/shadesmith/, Abruf: 11. Januar 2008

[shab] Shallows Homepage. http://shallows.sourceforge.net, Abruf: 10. Januar 2008

[Sha07] Shankland, S.: CNET News.com: Google acquires programming toolmaker PeakStream. Version: Juni 2007. http://www.news.com/2100-1007 3-6188935. html, Abruf: 10. Januar 2008

[SHL91] Siegel, R. ; Howell, J.R. ; Lohrengel, J.: Wärmeübertragung durch Strahlung. Teil 2, Strahlungsaustausch zwischen Oberflächen und in Umhüllungen. Springer-Verlag, 1991

[STM04] Sander, P.V. ; Tatarchuk, N. ; Mitchell, J.L.: Explicit Early-Z Culling for Efficient Fluid Flow Simulation and Rendering / ATI Research Technical Report. 2004. - Forschungsbericht

[Stu07] Stussak, C.: Echtzeit-Raytracing algebraischer Flächen auf der Graphics Processing Unit, Martin-Luther-Universität Halle-Wittenberg, Diplomarbeit, 2007

[SWW+04] Schmittler, J. ; Woop, S. ; Wagner, D. ; Paul, W.J. ; Slusallek, P.: Real- time Ray Tracing of Dynamic Scenes on an FPGA Chip. In: Proceedings of the ACM SIGGRAPH/EUROGRAPHICS conference on Graphics hardware (2004),S. 95-106

[Wag93] Wagner, W.: Wärmeübertragung. Vogel, 1993

[Wat02] Watt, A.: 3D-Computergrafik. Pearson Studium, 2002

[WEH89] Wallace, J.R. ; Elmquist, K.A. ; Haines, E.A.: A Ray Tracing Algorithm for Progressive Radiosity. In: ACM SIGGRAPH Computer Graphics 23 (1989), Nr. 3, S. 315-324

[Wei] Weisstein, E.W.: Point-Plane Distance. http://mathworld.wolfram.com/ Point-PlaneDistance.html, Abruf: 03. Februar 2008. MathWorld-A Wolfram Web Resource

[Wer92] Werner, J.: Numerische Mathematik 1: Lineare und nichtlineare Gleichungs- systeme, Interpolation, numerische Integration. Vieweg, 1992