Raytracing / Pathtracing / Photon Mapping - Tutorial

XMAMan

Monte Carlo Integration Cosinus-Integral mit Importancesampling mit Linie

Beim vorherigen Beispiel wurden die Zufallszahlen völlig gleichmäßig im Bereich von 0 bis PI/2 generiert. Man erhält beim Monte Carlo Integrieren schneller gute Ergebnisse, wenn man die Samples so erzeugt, dass sie möglichst die Funktion nachbilden, die man integrieren möchte.

Hier in diesen Bild ist in blau die cos(x)-Funktion dargestellt, welche wir integrieren möchten. In rot ist die pdf(x) aus dem vorherigen Beispiel eingezeichnet.

https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/03_MonteCarloIntegration/Images/ImportaneSamplingLine.png

Würde man es schaffen Zufallszahlen im Bereich 0 bis PI/2 so zu generieren, dass sie von der pdf her so aussehen wie die grüne Linie, dann hätte man nach 10 Samples schon ein viel genaueres Ergebnis. Dummerweise erzeugt die rand-Funktion von C# immer nur gleich verteilte Zufallszahlen.

Was macht man da nun?

Man berechnet die CDF (Cummulative distribution function) und bildet davon die Invers-Funktion.

Die cdf-Funktion ist wie folgt definiert.

$\mathit{cdf}(x)=\int _{-{\infty}}^x\mathit{pdf}(t)\mathit{dt}$

Es werden alle Pdf-Werte bis zu ein gegebenen x-Wert aufsummiert. In unseren Beispiel ist cdf(PI/2) = 1, da unsere Pdf-Funktion nur im Bereich von 0 bis PI/2 definiert ist und die Summe ja 1 ergeben muss.

Das interessante an dieser cdf-Funktion ist folgender Umstand. Als Input bekommt sie unser x. Als Output eine Zahl im Bereich von 0 bis 1.

Hier ist eine CDF in blau dargestellt, welche x-Werte im Bereich von 0 bis PI/2 erwartet und wo dann als maximale Zahl eine 1 zurück gegeben wird.

https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/03_MonteCarloIntegration/Images/CDF.png

Ich habe mit roter Farbe für zwei Beispiele die Beziehung zwischen x und cdf(x) eingezeichnet.

Wenn man nun die Inverse von der CDF bildet, dann erhält man eine Funktion, welche als Input eine Zahl im Bereich von 0 bis 1 erwartet und welche als Output eine Zahl im Bereich von 0 bis PI/2 liefert. Diese Funktion wäre ja genau das was wir brauchen. Als Input liefere ich eine Zahl,welche ich mit rand.NextDouble() erzeuge. Auf die Weise erhalte ich eine Zufallszahl, welche von der Pdf her so ist, die Pdf, welcher in dieser CDF-Funktion drin steckt.

Folgende Schritt müssen nun getan werden, um eine Samplefunktion zu bauen, welche Zufallszahlen mit der Pdf erzeugt, wie die grüne Linie oben aussieht.

1. Funktion aufstellen, welche die grüne Linie beschreibt
2. Pdf-Normalisierungskonstante berechnen.
3. CDF berechnen
4. Inverse von der CDF bilden

Schritt 1: Funktion aufstellen, welche die grüne Linie beschreibt
Die grüne Linie geht durch den Punkt (0, 1) und (PI/2, 0)

$f(x)=y=x\ast \frac{-1}{\mathit{PI}/2}+1$

$f(x)=y=1-\frac 2{\mathit{PI}}\ast x$

Schritt 2: Pdf-Normalisierungskonstante berechnen

Damit die soeben aufgestellte Linien-Funktion als Pdf verwendet werden muss sie der Bedingung genügen, dass das Integral für diese Linie bei 0 .. PI/2 eins ergeben muss.

Das Integral ist die Fläche unter der grünen Linie und somit ein halbes Rechteck mit

PI/2 * 1 / 2 = PI/4.

PI/4 ist ungleich 1. Damit die Pdf also 1 in der Summ ergibt, muss ich die Funktion aus Schritt 1 mit PI/4 dividieren.

Die Normalisierungskonstante pdf_c = PI / 4.

Und somit ist die pdf(x) = f(x) / pdf_c =
$\mathit{pdf}(x)=\frac 4{\mathit{PI}}-\frac 8{\mathit{PI}^2}\ast x$

Schritt 3: CDF berechnen

Um die CDF Funktion zu berechnen, muss ich einfach nur die Stammfunktion von der pdf bilden. Ich integriere nach x und erhalte:
$\mathit{cdf}(x)=\frac 4{\mathit{PI}}\ast x-\frac 4{\mathit{PI}^2}\ast x^2$

Schritt 4: Inverse von der CDF bilden

Jetzt besteht die Aufgabe darin die Gleichung

$y=\frac 4{\mathit{PI}}\ast x-\frac 4{\mathit{PI}^2}\ast x^2$

nach x umzustellen.

Ich wende das Distributivgesetz an und ziehe 4/PI² raus, um somit das x² zu befreien.

$y=\frac{-4}{\mathit{PI}^2}\ast (x^2-\mathit{PI}\ast x)$

Nun wende ich die Quadratische Ergänzung an.
http://www.mathebibel.de/quadratische-ergaenzung
http://www.onlinemathe.de/forum/Umkehrfunktion-von-x2-6x-1-berechnen

Ich addiere erst (-PI/2)² und subtrahiere sofort wieder.

$y=\frac{-4}{\mathit{PI}^2}\ast (x^2-\mathit{PI}\ast x+(\frac{-\mathit{PI}} 2)^2-(\frac{-\mathit{PI}} 2)^2)$

Mit der Binomischen Formel fasse ich x²-\mathit{PI}\ast x+(\frac{-\mathit{PI}} 2)² zusammen

$y=\frac{-4}{\mathit{PI}^2}\ast ((x-\frac{\mathit{PI}} 2)^2-(\frac{-\mathit{PI}} 2)^2)$

Ich multipliziere -4/PI² in die äußere Klammer rein, um diese Wegzumachen

$y=(x-\frac{\mathit{PI}} 2)^2\ast (\frac{-4}{\mathit{PI}^2})+1$

Nun bringe ich den Scheiß, der nichts mit x zu tun hat auf die linke Seite

$\frac{y-1}{\frac{-4}{\mathit{PI}^2}}=(x-\frac{\mathit{PI}} 2)^2$

Wurzel ziehen.

$\sqrt{\frac{y-1}{\frac{-4}{\mathit{PI}^2}}}=x-\frac{\mathit{PI}} 2$

Jetzt nur noch das PI/2 rüber und schon steht x auf einer Seite alleine

$\sqrt{\frac{y-1}{\frac{-4}{\mathit{PI}^2}}}+\frac{\mathit{PI}} 2=x$

Das was wir jetzt hier nach langen Mühen ermittelt haben ist die Samplefunktion.

Nun haben wir zwei Bausteine, womit wir die cos(x)-Integration lösen können:

$\mathit{pdf}(x)=\frac 4{\mathit{PI}}-\frac 8{\mathit{PI}^2}\ast x$

SampleFunktion: Wobei y eine Zufallszahl im Bereich 0 bis 1 ist.

$\sqrt{\frac{y-1}{\frac{-4}{\mathit{PI}^2}}}+\frac{\mathit{PI}} 2=x$

In C# sieht das nun so aus:

Random rand = new Random(0);
int sampleCount = 10;

float sum = 0;
for (int i = 0; i < sampleCount; i++)
{
  float x = (float)(Math.Sqrt((rand.NextDouble() - 1) / (-4 / (Math.PI * Math.PI))) + Math.PI / 2);
  float pdf = (float)(4 / Math.PI - 8 / (Math.PI * Math.PI) * x);
  sum += (float)Math.Cos(x) / pdf;
}
sum /= sampleCount;

Nach 10 Samples erhalte ich 1.05.

Es werden nun viel weniger Samples benötigt, um ein gutes Ergebnis zu erhalten.

Beispielquelltext (siehe Example3)

https://github.com/XMAMan/RaytracingTutorials/blob/master/03_MonteCarloIntegration/RaytracingTutorials/MonteCarloIntegration.cs

XMAMan

Monte Carlo Integration Cosinus-Integral mit Importancesampling mit Cosinusfunktion

Das letzte Beispiel mit der grünen Linie war schon heftig, was dieses ganze Mathezeug angeht. Jetzt kommt nochmal ein kleiner Matheteil aber dann habt ihr dieses verfluchte Cosinus-Integral-Problem endlich geschafft.

Diesmal sollen die Zufallswerte direkt mit Cos(x) als Pdf erstellt werden. Das wäre also die perfekte Übereinstimmung zwischen der Funktion, die man integrieren möchte und der pdf. Hier müsste man also am Schnellsten zum Ergebnis kommen.

Es werden wieder die 4 Schritt abgearbeitet.
1. Funktion aufstellen, welche die unnormalisierte Pdf beschreibt
2. Pdf-Normalisierungskonstante berechnen.
3. CDF berechnen
4. Inverse von der CDF bilden

Schritt 1: Cos(x) ist unsere gesuchte funktion.

Schritt 2: Pdf-Normalisierungskonstante berechnen
Das Integral von Cos(x) im Bereich x=0..PI/2 = 1. Das haben wir vorhin bereits raus gefunden. Die Konstante ist also 1 und somit ist die pdf(x) = cos(x)

Schritt 3: CDF berechnen
Die Stammfunktion von cos(x) ist sin(x) somit ist cdf(x) = sin(x)

Schritt 4: Inverse von der CDF bilden um Samplefunktion zu erhalten
Die Inverse von sin(x) ist arcsin(x)

C#-Beispiel:

Random rand = new Random(0);
int sampleCount = 10;

float sum = 0;
for (int i = 0; i < sampleCount; i++)
{
    float x = (float)Math.Asin(rand.NextDouble());
    float pdf = (float)Math.Cos(x);
    sum += (float)Math.Cos(x) / pdf; //Cos(x) / Cos(x) -> Kürzt sich zu 1
}

sum /= sampleCount;

Was sehen wir da? Cos(x) / pdf ergibt immer 1. Somit hat man schon nach einen Sample das richtige Ergebnis. Was heißt das? In dem Moment, wo man die Pdf-Normalisierungskonstante für die Funktion, die man integrieren möchte bereits kennt, kann man sich den restlichen Aufwand sparen. Schließlich ist die Konstante ja genau das gesucht Integral.

In diesen kleinen Ausflug sollte es darum gehen, wie man mit Monte Carlo überhaupt erst mal integriert. Wir haben hier nur float- und Integerzahlen addiert. Im nächsten Abschnitt soll nun auf das Thema eingegangen werden, wie man Lichtpfade aufsummiert.

XMAMan

Monte Carlo Integration und das Path-Integral

Ausgangspunkt ist wie immer unsere geliebte Monte-Formel:

$\frac 1 N\sum _{i=1}^N\frac{f(\mathit{xi})}{\mathit{pdf}(\mathit{xi})}$

Folgene Bausteine benötigen wir nun:

xi = zufällig generiertes Array von 3D-Punkten (Lichtpfad)
f(xi) = PathContributionFunction = Produkt aus allen Brdf- und Geometrytermen und dem Light-Emission-Faktor.
pdf(xi) = Wahrscheinlichkeit(sdichte), ein bestimmten Pfad zu erzeugen

Es wird eine Funktion benötigt, welche ein zufälligen Lichtpfad x erzeugt. So ein Lichtpfad wird beim Raytracing schrittweise Punkt für Punkt aufgebaut. Der erste Punkt auf der Kamera ist gegeben. Diesen muss man nicht zufällig erzeugen. Der nächste Punkt wird mit Hilfe der Bildebene bestimmt, indem ein Strahl durch die Bildebene für ein gegebenes Pixel geschossen wird und dann geschaut wird, welcher Punkt getroffen wird. Somit ist die Erzeugung des zweiten Punktes auch noch machbar. Um nun den nächsten Punkt zu erzeugen, hat man folgende Möglichkeiten:

1. Brdf-Sampling: Ich habe ein Punkt und ein Richtungsvektor, welcher zu diesem Punkt zeigt, gegeben, und bestimme dann eine neue zufällige Richtung und schaue dann, welcher Punkt in diese Richtung liegt.
2. Light-Area-Sampling: Ich entscheide mich dazu ein zufälligen Punkt auf einer zufällig ausgewählten Lichtquelle zu erzeugen und erhalte somit den Endpunkt von dem Lichtpfad.

Man kann also einen Pfad erzeugen, indem ich einfach auf der Kamera starte und dann mit Brdf-Sampling durch den Raum durchwandere. Das wäre der PathTracing-Ansatz. Genau so gut kann man auch auf der Lichtquelle starten und dann einfach sich direkt (ohne Zufall) mit der Kamera verbinden (LightTracing-Ansatz).

Nächste Frage ist nur noch: Wie berechne ich die Pdf für ein gegebenen Pfad?

Die Pdf für den Pfad x mit der Pfadlänge 3 wird wie folgt berechnet:

Pdf(x) = pdfA(x0) * pdfA(x1) * pdfA(x2)

Die Wahrscheinlichkeit, dass ich den Pfad x mit den Punkten x0, x1 und x2 erzeuge ist das Produkt der pdfA-s(A steht für Area) für die einzelnen Punkte. Die pdfA wiederum ist die Wahrscheinlichkeit, dass ich ein 3D-Punkt xi an einer gegebenen Stelle in der Szene erzeuge.

Wenn ich ein Dreieck mit den Flächeninhalt von 5 habe und ich erzeuge zufällig uniform ein Punkt auf diesem Dreieck, so ist die PdfA für diesen Punkt auf diesen Dreieck 1 / 5. Das Integral der PdfA auf dem Dreieck muss ja 1 ergeben. Damit das der Fall ist, muss die PdfA 1 / TriangleSurfaceArea sein. Jedes Dreieck besitzt also seine eigene Samplefunktion, um Zufallspunkte zu erzeugen und jedes Dreieck hat auch seine eigene PdfA.

Wenn ich nun eine Szene mit 100 Dreiecken habe und ich will nun auf einen von diesen 100 Dreiecken einen Zufallspunkt erzeugen, dann könnte eine Samplefunktion z.B. erst Zufällig uniform eins von den Dreiecken auswählen und dann wird auf dem ausgewählten Dreieck ein Zufallspunkt erzeugt. Die PdfA für diesen so erzeugten Punkt ist dann 1/100 * PdfA_von_ausgewählten_Dreieck.

Wenn ich ein Zufallspunkt erzeuge, indem ich eine zufällige Richtung von ein gegebenen Punkt erzeuge und dann mit der Szene-Schnittpunkt-Funktion diesen neuen Punkt ermittle, dann sample ich ja diesmal kein Dreieck sondern eine Brdf. Die Wahrscheinlichkeit für eine Richtung nennt man pdfW. Diese muss erst mit folgender Formel in eine PdfA umgewandelt werden, um somit die Path-PDF bestimmen zu können:

Gegeben ist P1 und ich habe die Richtung R gesampelt. Es wird nun P2 mithilfe der Strahl-Szene-Schnittpunktabfrage bestimmt, indem von P1 in Richtung R ein Strahl geschossen wird und geschaut wird, welchen Punkt P2 er in der Szene trifft. Die PdfA von P2 kann ich nun mit

pdfA(P2) = pdfW(R) * (-R) * N / d²

pdfw = Wahrscheinlichkeitsdichte von der Brdf-Samplefunktion
N = Normale bei P2
d = Abstand zwischen P1 und P2.

ermitteln. Diese geheime Umrechnungsformel habe ich zum ersten mal bei Eric Veach gesehen:

http://graphics.stanford.edu/papers/veach_thesis/thesis.pdf
Seite 254

https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/04_PathPdf/Images/pdfWToPdfA_ConverstionFaktor.png

Die Pfad-Pdf muss also als Produkt von allen PdfA's angegeben werden. Die PdfA-Funktion bekommt einen 3D-Punkt (Punkt auf ein Dreieck aus der Szene) als Input und gibt eine float-Zahl zurück, welche sagt, wie hoch die Wahrscheinlichkeit ist, dass genau an dieser Stelle ein 3D-Pfad-Punkt zufällig erzeugt wird.

Die Frage die jetzt noch offen ist, ist, wie funktioniert Brdf-Sampling und all die Sachen, die hier in diesen Abschnitt theoretisch erklärt wurden an mehreren Praxisbeispielen.

XMAMan

Brdf-Sampling

In diesen Abschnitt wollte ich eigentlich darüber schreiben, wie man bei ein diffusen Material ein Richtungsvektor R erzeugt, welcher die Funktion R*Normale importancesampelt. Mir ist es leider nicht gelungen die Herleitung dafür zu finden. Dank SmallVCM weiß ich aber, dass die Funktion so aussehen muss:

Vector SampleDirection(Vector normale, Random rand)
{
    float r1 = (float)(2 * Math.PI * rand.NextDouble());
    float r2 = (float)rand.NextDouble();

    float r2s = (float)Math.Sqrt(1 - r2);

    Vector w = normale,
           u = Vector.Cross((Math.Abs(w.X) > 0.1f ? new Vector(0, 1, 0) : new Vector(1, 0, 0)), w).Normalize(),
           v = Vector.Cross(w, u);

    Vector d = ((u * (float)Math.Cos(r1) * r2s + v * (float)Math.Sin(r1) * r2s + w * (float)Math.Sqrt(r2))).Normalize();

    return d;
}

Die dazugehörte PdfW-Funktion:

float PdfW(Vector normal, Vector sampledDirection)
{
    return Math.Max(1e-6f, Math.Max(0, Vector.Dot(normal, sampledDirection)) / (float)Math.PI);
}

Es wird mit u, v und normale ein Koordinatensystem aufgespannt. Die erste Zufallszahl r1 geht gleichmäßig von 0 bis 2*PI. Hier ist kein Importancesampling nötig, da dieser Anteil des gesampelten Richtungsvektor kein Einfluss auf R*N hat. Die Zahl r2 ist dann die Cosinusgewichtete Importance-Sampling-Zufallszahl. Es klingt jetzt vielleicht dürftig von der Erklärung aber vielleicht findet sich ja jemand, der diese Herleitung findet. Diese Funktion wird von den nachfolgenden Beispielen auch verwendet. Um das Nachfolgende zu verstehen, ist die Herleitung aber nicht nötig.

XMAMan

Pathtracing

Kommen wir nun endlich zu unseren ersten Algorithmus, womit globale Beleuchtung möglich ist. Mit Pathtracing werden zufällige Lichtpfade erzeugt, welche bei der Kamera starten, dann zufällig durch den Raum mit Brdf-Sampling wandern, und somit versuchen die Lichtquelle zu treffen. Ist dies geglückt, so hat man eine Menge von 3D-Punkten, wo der erste Punkt der Kamera-Punkt ist, dann liegen N Punkte irgendwo in der Szene und der letzte Punkt dann auf einer Lichtquelle.

Schauen wir uns im Detail mal ein Lichtpfad mit 3 Punkten an, welcher mit Pathtracing erstellt wurde.

Punkt P0 wurde ohne Sampling erzeugt. Es wurde einfach der Kamera-Punkt genommen. Dann wurde ein Pixel vorgegeben, für das wir die Farbe bestimmen wollen. Auf diese Weise erhalten wir den ersten Richtungsvektor R0, welche von P0 in Richtung der Szene zeigt. Mit dem Strahl-Szene-Schnittpunkttest erhalten wir P1. Nun wird mit Brdf-Sampling beim Punkt P1 eine neue Richtung R1 erzeugt. Der Schnittpunkttest liefert nun P2, welcher auf der Lichtquelle liegt.

So sieht der Pfad nun aus:
https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/05_Pathtracing/Images/3PointPath.png

Wenn wir nun somit diese 3 Punkt erzeugt haben, so müssen nun zwei Dinge als nächstes mithilfe dieser 3 Punkte berechnet werden: Die PathContribution und die Pfad-PDF. Mit diesen beiden Angaben können wir dann ein Schätzwert für die Pixelfarbe erzeugen und wenn man den Durchschnitt von vielen von solchen Pixelschätzwerten bildet, dann weiß man die Pixelfarbe. Somit bauen wir uns unser fertiges Bild zusammen.

Die PathContribution (Vektor) berechnet sich aus:

PathContribution = GeometryTerm(P0,P1) * Brdf(P1) * GeometryTerm(P1, P2)

Dabei ist

$\mathit{GeometryTerm}(\mathit{P0},\mathit{P1})=\frac{(\mathit{N0}\ast \mathit{R0})\ast ((-\mathit{R0})\ast \mathit{N1})}{\mathit{d_1}^2}$

und

$\mathit{GeometryTerm}(\mathit{P1},\mathit{P1})=\frac{(\mathit{N1}\ast \mathit{R1})\ast ((-\mathit{R1})\ast \mathit{N2})}{\mathit{d_2}^2}$

Wir gehen von einer diffusen Brdf aus:

$\mathit{Brdf}(\mathit{P1})=\frac{\mathit{Farbe}\mathit{von}\mathit{P1}}{\mathit{PI}}$

Beim erzeugen dieses Pfades gab es zwei zufällige Dinge:
1. Um R0 zu bestimmen wurde zufällig für ein gegebenen Pixel auf der Bildebene ein Punkt erzeugt. Für dieses Richtungssampling von R0 haben wir die PdfW(R0)
2. R1 wurde zufällig mit Brdf-Sampling bestimmt. Auch hier hat die Brdf-Sample-Funktion eine PdfW(R1)
Nehmen wir mal an, ich schieße immer durch die Pixelmitte einen Strahl. Dann ist PdfW(R0) = 1. Da P1 ein diffuser Punkt ist, gilt hier $\mathit{PdfW}(\mathit{R1})=\frac{\mathit{N1}\ast \mathit{R1}}{\mathit{PI}}$

Die Pfad-PDF muss immer als PdfA angegeben werden. Also

Path-PDF = PdfA(P0) * PdfA(P1) * PdfA(P2)

PdfA(P0) = 1 (Kamerapunkt wird nicht gesampelt)

Die PdfA vom P1 erhalten wir mit dem PdfW-to-PdfA-Conversionfaktor.

$\mathit{PdfA}(\mathit{P1})=\frac{\mathit{N1}\ast (-\mathit{R0})}{d_1^2}\ast \mathit{PdfW}(\mathit{R0})$

Das gleiche gilt für die PdfA von P2

$\mathit{PdfA}(\mathit{P2})=\frac{\mathit{N2}\ast (-\mathit{R1})}{d_2^2}\ast \mathit{PdfW}(\mathit{R1})$

Jetzt haben wir die PathContribution und die Path-Pdf erst mal einzeln so hingeschrieben. Wenn man sich aber die Monte-Carlo-Formel ansieht:

$\frac 1 N\sum _{i=1}^N\frac{\mathit{PathContribution}(\mathit{xi})}{\mathit{PathPdfA}(\mathit{xi})}$

dann interessiert uns ja nur das Verhältnis zwischen diesen beiden Termen. Wenn man das hinschreibt, so sieht man, dass sich ein Teil vom GeometryTerm und aus dem PdfW-to-PdfA-Conversionfaktor gegenseitig wegkürzt.

$\frac{\mathit{GeometryTerm}(\mathit{P0},\mathit{P1})\ast \mathit{Brdf}(\mathit{P1})\ast \mathit{GeometryTerm}(\mathit{P1},\mathit{P2})}{\mathit{PdfA}(\mathit{P0})\ast \mathit{PdfA}(\mathit{P1})\ast \mathit{PdfA}(\mathit{P2})}$

Ich teile den Bruch in vier Unterbereiche auf:

$\frac 1{\mathit{PdfA}(\mathit{P0})}\ast \frac{\mathit{GeometryTerm}(\mathit{P0},\mathit{P1})}{\mathit{PdfA}(\mathit{P1})}\ast \frac{\mathit{Brdf}(\mathit{P1})} 1\ast \frac{\mathit{GeometryTerm}(\mathit{P1},\mathit{P2})}{\mathit{PdfA}(\mathit{P2})}$

Unterbereich 2:

$\frac{\mathit{GeometryTerm}(\mathit{P0},\mathit{P1})}{\mathit{PdfA}(\mathit{P1})}=\frac{\frac{(\mathit{N0}\ast \mathit{R0})\ast ((-\mathit{R0})\ast \mathit{N1})}{d\_1^2}}{\frac{\mathit{N1}\ast (-\mathit{R0})}{d\_1^2}\ast \mathit{PdfW}(\mathit{R0})}=\frac{\mathit{N0}\ast \mathit{R0}}{\mathit{PdfW}(\mathit{R0})}$

Unterbereich 4:

$\frac{\mathit{GeometryTerm}(\mathit{P1},\mathit{P2})}{\mathit{PdfA}(\mathit{P2})}=\frac{\frac{(\mathit{N1}\ast \mathit{R1})\ast ((-\mathit{R1})\ast \mathit{N2})}{\mathit{d2}^2}}{\frac{\mathit{N2}\ast (-\mathit{R1})}{d_2^2}\ast \mathit{PdfW}(\mathit{R1})}=\frac{\mathit{N1}\ast \mathit{R1}}{\mathit{PdfW}(\mathit{R1})}$

Wenn ich Unterbereich 1 und 2 zusammenfasse und für PdfA(P0) = 1 und PdfW(R0) 1 einsetze, dann erhalte ich:

N0 * R0

Unterbereich 3 und 4 fasse ich auch zusammen:
$\frac{\mathit{Brdf}(\mathit{P1})\ast (\mathit{N1}\ast \mathit{R1})}{\mathit{PdfW}(\mathit{R1})}$

Was sehe ich da? Alle Angaben die dort stehen, sind zum Zeitpunkt des Brdf-Samplings beim Punkt P1 vorhanden. D.h. Ich kenne die Position von P1. Ich habe per Brdf-Sampling R1 ermittelt aber P2 wurde noch nicht über die Strahl-Szene-Schnittpunktabfrage ermittelt. Ich kann den Brdf-Wert und den Cosinus-Wert von N1 zu R1 = (N1 * R1) ausrechnen. Und erhalte somit ein RGB-Vektor, welche sich aus den 3 Angaben:
Brdf = RGB-Vektor
N1 * R1 = float-Zahl
PdfW(R1) = float-Zahl

errechnet.

Würde ich jetzt noch einen weiteren Pfad-Punkt dranhängen. So hätte ich beim Punkt P2 und den dort erfolgten Brdf-Sampling dann wieder alle drei benötigten Angaben, um den Gesamt-Pfadgwicht damit weiter auszurechnen.

Folgender C#-Code berechnet für ein einzelnen Pixel ein Farbschätzwert mit Pathtracing

Vector EstimateColorWithPathtracing(int x, int y, Random rand)
{
    Ray primaryRay = data.Camera.GetPrimaryRay(x, y);
    float cosAtCamera = Vector.Dot(primaryRay.Direction, this.data.Camera.Forward);
    Vector pathWeight = new Vector(1, 1, 1) * cosAtCamera;

    Ray ray = primaryRay;

    int maxPathLength = 5;
    for (int i = 0; i < maxPathLength; i++)
    {
        var point = this.intersectionFinder.GetIntersectionPoint(ray);
        if (point == null)
        {
            return new Vector(0, 0, 0);
        }

        if (point.IsLocatedOnLightSource && Vector.Dot(point.Normal, -ray.Direction) > 0)
        {
            pathWeight = Vector.Mult(pathWeight, point.Emission);
            return pathWeight;
        }

        //Abbruch mit Russia Rollete
        float continuationPdf = Math.Min(1, point.Color.Max());
        if (rand.NextDouble() >= continuationPdf)
            return new Vector(0, 0, 0); // Absorbation

        Vector newDirection = DiffuseBrdf.SampleDirection(point.Normal, rand);
        pathWeight = Vector.Mult(pathWeight * Vector.Dot(point.Normal, newDirection) / DiffuseBrdf.PdfW(point.Normal, newDirection), DiffuseBrdf.Evaluate(point)) / continuationPdf;
        ray = new Ray(point.Position, newDirection);
    }

    return new Vector(0, 0, 0);
}

Hier sieht man deutlich, dass zum ausrechnen des Pfadgewichts kein Geometryterm oder bereits der nächste Pfad-Punkt erforderlich ist.

Außerdem wird die Pfad-Wanderung hier zufällig mit etwas, was sich russisches Rollete nennt, beendet.

float continuationPdf = Math.Min(1, point.Color.Max());
if (rand.NextDouble() >= continuationPdf)
   return new Vector(0, 0, 0); // Absorbation

Nehmen wir mal an eine rote Fläche soll 50% seiner erhaltenen Photonen reflektieren und die anderen 50% soll es absorbieren (In Wärmeenergie umwandeln).

Dann hat continuationPdf den Wert 0.5. Wenn nun die rand-Funktion (erzeugt Zahlen zufällig im Bereich von 0 bis 1) eine Zahl größer 0.5 erzeugt, dann wird 0 zurück geben. Wenn ein Wert kleiner 0.5 erzeugt, dann wird z.B. (1,0,0) zurück gegeben. Rufe ich 20 mal diese Funktion auf, dann erhalte ich im Mittel 10 mal den Wert 0 und 10 mal den Wert (1,0,0). Die Summe davon ist (10,0,0). Da wir aber das Endergebnis durch die Sampleanzahl (=20) teilen müssen, erhalten wir somit als Farbe für ein Pixel (0.5 ; 0 ; 0). Also ein mittelhelles rot.

Die continuationPdf ist eine zufällige Zahl, welche den Gesamt-Lichtpfad beeinflusst. Jede zufällig erzeugte Zahl, sei es nun eine Zufallsrichtung oder ein zufälliger Punkt auf einer Lichtquelle muss sich in der Pfad-Pdf widerspiegeln. Aus diesen Grund teile ich das Pfadgewicht durch die continuationPdf so wie die PdfW ja auch.

Wenn ich 1000 Samples pro Pixel nehme, dann bekomme ich folgendes Ergebnis:
https://github.com/XMAMan/RaytracingTutorials/blob/master/05_Pathtracing/Images/1000Samples.jpg?raw=true

Das Bild ist jetzt vom Beleuchtungsalgorithmus her jetzt realistischer aber dafür sieht es sehr verrauscht aus. Ich könnte jetzt anstatt 1000 Samples auch 10 oder 100 mal so viele Samples nehmen. Dann hätte man dann ein rauschfreies Bild. Aber das würde dann auch entsprechend länger dauern. Aus diesem Grund wollen wir uns nun Algorithmen zuwenden, welche viel schneller ein gutes Bild erzeugen.

Beispielquellcode für den Pathtracer: https://github.com/XMAMan/RaytracingTutorials/tree/master/05_Pathtracing/RaytracingTutorials

XMAMan

Photonmap

Beim Pathtracing wurden Lichtpfade dadurch erstellt, indem man bei der Kamera startet und dann zufällig durch den Raum wandert, und darauf hofft, dass man die Lichtquelle trifft. Um so kleiner die Lichtquelle ist oder um so weiter weg von der Kamera sie ist, um so unwahrscheinlicher ist es, dass die Lichtquelle somit erreicht wird. Das Photonmapping versucht ein anderen Ansatz. Ein Photon startet auf ein zufällig auswählten Punkt auf der Lichtquelle und wandert dann zufällig durch den Raum. Diesmal ist es aber egal, ob man die Kamera trifft oder nicht. Jedes mal, wenn man auf eine Wand trifft, dann wird die Position/Farbe von dem Photon, was hier zufällig wandert, gespeichert und die Zufallswanderung geht weiter. Wenn ich mit einer Rekursionstiefe von 10 arbeite, dann kann somit ein einzelnes Photon 10 Wand-Photon-Treffer-Ereignisse generieren.
Es werden sehr viele Photonen ausgesendet und jedes Wand-Auftreffen wird gespeichert. Auf diese Weise erhalte ich eine Liste von Wand-Auftreff-Ereignissen. Diese Liste nennt sich Photonmap.

Nachfolgend möchte ich das Erstellen der Photonmap am C#-Beispiel erläutern.

Die Klasse Photon speichert das Auftreffen eines Photons gegen eine Wand. Wichtig ist jetzt erstmals nur, das du siehst, hier gibt es die Klasse Photon, welche eine Position und eine PathWeight-Property hat. Die anderen Sachen werden später wichtig und können jetzt aber erst mal überlesen werden.

class Photon : IPoint
{
    public float this[int key]
    {
        get
        {
            return this.Position[key];
        }
    }

    public Vector Position { get; private set; }
    public Vector Normal { get; private set; }
    public Vector PathWeight { get; private set; }
    public Vector DirectionToThisPoint { get; private set; }

    public Photon(Vector position, Vector normal, Vector pathWeight, Vector directionToThisPoint)
    {
        this.Position = position;
        this.Normal = normal;
        this.PathWeight = pathWeight;
        this.DirectionToThisPoint = directionToThisPoint;
    }
}

Die nächste Methode lässt ein Photon zufällig durch den Raum wandern und gibt eine Liste von Photon-Objekten zurück. Hier erst mal im ganzen. Wir gehen sie dann einzeln durch.

List<Photon> TracePhoton(Random rand)
{
    var lightPoint = this.lightSource.GetRandomPointOnLightSource(rand);
    Vector direction = DiffuseBrdf.SampleDirection(lightPoint.Normal, rand);
    float pdfW = DiffuseBrdf.PdfW(lightPoint.Normal, direction);
    float lambda = Vector.Dot(lightPoint.Normal, direction);
    Vector pathWeight = lightPoint.Color * lambda / (lightPoint.PdfA * pdfW);
    Ray ray = new Ray(lightPoint.Position, direction);

    int maxPathLength = 5;

    List<Photon> photons = new List<Photon>();

    for (int i = 0; i < maxPathLength; i++)
    {
        var point = this.intersectionFinder.GetIntersectionPoint(ray);
        if (point == null) return photons;

        photons.Add(new Photon(point.Position, point.Normal, pathWeight / this.photonCount, ray.Direction));

        //Abbruch mit Russia Rollete
        float continuationPdf = Math.Min(1, point.Color.Max());
        if (rand.NextDouble() >= continuationPdf)
           return photons; // Absorbation

        Vector newDirection = DiffuseBrdf.SampleDirection(point.Normal, rand);
        pathWeight = Vector.Mult(pathWeight * Vector.Dot(point.Normal, newDirection) / DiffuseBrdf.PdfW(point.Normal, newDirection), DiffuseBrdf.Evaluate(point)) / continuationPdf;
        ray = new Ray(point.Position, newDirection);
    }

    return photons;
}

Zuerst wird Zufallspunkt auf einer Lichtquelle erzeugt:

var lightPoint = this.lightSource.GetRandomPointOnLightSource(rand);

Eine Lichtquelle besteht aus einer Liste von Dreiecken. Ich wähle zufällig eins davon aus. Dann bestimme ich auf dem Dreieck ein Zufallspunkt. Die PdfA von diesen Punkt berechnet sich dann über triangleSelectionPdf * trianglePdfA. Ich frage zuerst, wie hoch ist die Wahrscheinlichkeit, dass ich genau das Dreieck 'triangle' aus der this.triangles-Liste auswähle? Danach: Wie hoch ist die Wahrscheinlichkeit(sdichte), dass ich genau diesen Punkt auf dem Dreieck ausgewählt habe? Die Multiplikation bedeutet, es muss erst das erste zufällie Ereigniss eintreten UND danach dann das andere. Multiplikation bedeutet also eine logische UND-Verknüpfung von diesen beiden Zufallsdingen.

public LightSourcePoint GetRandomPointOnLightSource(Random rand)
{
    var triangle = this.triangles[rand.Next(this.triangles.Length)];
    Vector position = triangle.GetRandomPointOnSurface(rand);

    float triangleSelectionPdf = 1.0f / this.triangles.Length;
    float trianglePdfA = 1.0f / triangle.Area;

    return new LightSourcePoint()
    {
        Position = position,
        PdfA = triangleSelectionPdf * trianglePdfA,
        Color = triangle.Color * this.EmissionPerArea,
        Normal = triangle.Normal
    };
}

Ich wähle eine zufällige Richtung aus, in die das Photon fliegen soll:

Vector direction = DiffuseBrdf.SampleDirection(lightPoint.Normal, rand);
float pdfW = DiffuseBrdf.PdfW(lightPoint.Normal, direction);

Ich berechne das Pfadgewicht indem ich die Farbe von der Lichtquelle nehme und ähnlich wie beim Pathtracing gilt hier: Mal Cosinus-Wert durch alle zufällig getroffenen Entscheidungen.

float lambda = Vector.Dot(lightPoint.Normal, direction);
Vector pathWeight = lightPoint.Color * lambda / (lightPoint.PdfA * pdfW);
Ray ray = new Ray(lightPoint.Position, direction);

Der restliche Teil von der Funktion ist dann wie das Pathtracing aufgebaut. Wenn ein Photon-Objekt gespeichert wird, dann wird nicht einfach nur das aktuelle Pfadgewicht genommen sondern es wird noch durch die Anzahl der Photonen, die ausgesendet werden dividiert.

photons.Add(new Photon(point.Position, point.Normal, pathWeight / this.photonCount, ray.Direction));

Würde ich das nicht machen, dann würde das dazu führen, dass die Photonmap um so heller erscheint, um so mehr Photonen ich aussende. Wenn ich aber die Gesamthelligkeit unabhängig von der Photonmap-Auflösung will, dann muss ich durch die Anzahl dividieren.

Die Klasse Photonmap, welche nun die ganz viele Photonen versendet und speichert, sieht wie folgt aus:

class Photonmap
{
    private IntersectionFinder intersectionFinder;
    private LightSource lightSource;
    private int photonCount;
    private KdTree kdTree;

    public Photonmap(IntersectionFinder intersectionFinder, LightSource lightSource, int photonCount)
    {
        this.intersectionFinder = intersectionFinder;
        this.lightSource = lightSource;
        this.photonCount = photonCount;

        Random rand = new Random(0);

        List<Photon> photons = new List<Photon>();
        for (int i = 0; i < photonCount; i++)
        {
            photons.AddRange(TracePhoton(rand));
        }

        this.kdTree = new KdTree(photons.Cast<IPoint>().ToArray());
    }

    private List<Photon> TracePhoton(Random rand)
    {
	 …
    }
}

Es werden 'photonCount' Photonen von der Lichtquelle los gesendet und jeder Photon-Wand-Treffer in der Lise 'photons' gespeichert. Diese 'photons'-Liste wird dann in ein Kd-Baum speichert, da damit dann die Abfrage der Photonmap schneller geht.

Bis jetzt haben wir uns nur mit der Erstellung der Photonmap beschäftigt. Um damit nun die Farbe für ein einzelnes Pixel zu berechnen, gehen wir ähnlich wie beim einfachen Raytracer so vor, dass wir einfach nur ein Primärstrahl losschicken. An der Stelle, wo der Strahl dann auf ein Punkt in der Szene trifft, da ermitteln wir dann alle Photonen im Umkreis.

https://github.com/XMAMan/RaytracingTutorials/blob/master/06_Photonmapping/Images/PhotonmapDirekt.png?raw=true

Vector GetPixelColorFromPhotonmap(int x, int y, Photonmap photonmap)
{
    Ray primaryRay = data.Camera.GetPrimaryRay(x, y);
    float cosAtCamera = Vector.Dot(primaryRay.Direction, this.data.Camera.Forward);
    Vector pathWeight = new Vector(1, 1, 1) * cosAtCamera;

    var eyePoint = this.intersectionFinder.GetIntersectionPoint(primaryRay);

    Vector pixelColor = new Vector(0, 0, 0);

    if (eyePoint != null)
    {
        float searchRadius = 0.00634257728f * 2;
        var photons = photonmap.SearchPhotons(eyePoint, searchRadius);
        float kernelFunction = 1.0f / (searchRadius * searchRadius * (float)Math.PI);
        foreach (var photon in photons)
        {
           pixelColor += Vector.Mult(pathWeight, Vector.Mult(DiffuseBrdf.Evaluate(eyePoint), photon.PathWeight) * kernelFunction);
        }
    }

    return pixelColor;
}

Ich bilde die Summe über alle Photonen, die ich finde und dividiere den Photon-Farbwert durch den Flächeninhalt des 2D-Suchkreises. Habe ich ein großen Suchradius, dann bekomme ich viele Photonen und dividiere diese große Summe dann durch ein großen Kreis-Flächinhalt. Ist das Suchradius ganz klein, finde ich auch nur wenige Photonen und die Summe ist somit auch geringer. Diese geringe Summe dividiere ich dann aber auch durch eine kleinere Zahl. Auf die Weise wird dafür gesorgt, dass die Helligkeit des Pixels immer ungefähr gleichhell ist. Egal ob ich mit ein großen oder kleinen Suchradius arbeite.

Das Bild sieht nun so aus.

https://github.com/XMAMan/RaytracingTutorials/blob/master/06_Photonmapping/Images/Photonmap.png?raw=true

Beispielquelltext: https://github.com/XMAMan/RaytracingTutorials/tree/master/06_Photonmapping/Photonmap-Source/RaytracingTutorials

Der Renderzeit ist jetzt nur noch ein Bruchteil von dem Pathtracing-Bild aber das Ergebnis ist leider nur eine Näherung. Dadurch, dass die Auflösung der Photonmap begrenzt ist, ist auch die Bildqualität nur begrenzt. Um nun endlich mal ein gutes Ergebnis mit so einer Photonmap zu erhalten, ist ein weiterer Schritt erforderlich: Final Gathering. Was das ist erkläre ich im nächsten Abschnitt.

XMAMan

Photonmapping

Beim Final Gathering wird erst ein Primärstrahl von der Kamera gesendet. Bei der ersten diffusen Fläche, wo der Primärstrahl dann auftritt, wird per Brdf-Sampling ein weiterer sogenannter Final Gather Ray versendet. Dort wo der Final Gather Ray dann auftrifft wird letztendlich erst die Photonmap ausgelesen.
https://github.com/XMAMan/RaytracingTutorials/blob/master/06_Photonmapping/Images/FinalGathering.png?raw=true

Das Final Gathering berechnet für den Punkt P1 (siehe Bild) nur den indirekten Lichtanteil. Das direkte Licht wird wie beim einfachen Raytracer aus unserem Anfangsbeispiel dadurch berechnet, indem ich von P1 aus ein Schattenstrahl zu ein zufällig ausgewählten Punkt P2 auf einer Lichtquelle versende und somit den Lichtpfad (P0, P1, P2) erzeuge.

Die Funktion, welche die Farbe für ein Pixel bestimmt, ist nun also die Summe von direkten Licht und indirekten Licht.

Vector EstimateColorWithPhotonmapping(int x, int y, Random rand, Photonmap photonmap)
{
    Ray primaryRay = data.Camera.GetPrimaryRay(x, y);
    float cosAtCamera = Vector.Dot(primaryRay.Direction, this.data.Camera.Forward);
    Vector pathWeight = new Vector(1, 1, 1) * cosAtCamera;

    var eyePoint = this.intersectionFinder.GetIntersectionPoint(primaryRay);

    Vector pixelColor = new Vector(0, 0, 0);

    if (eyePoint != null)
    {
        if (eyePoint.IsLocatedOnLightSource)
        {
            return eyePoint.Emission;
        }

        pathWeight = Vector.Mult(pathWeight, DiffuseBrdf.Evaluate(eyePoint));

        //Direct Lighting
        var lightPoint = data.LightSource.GetRandomPointOnLightSource(rand);
        if (IsVisible(eyePoint.Position, lightPoint.Position))
        {
            pixelColor += Vector.Mult(pathWeight, lightPoint.Color * GeometryTerm(eyePoint, lightPoint) / lightPoint.PdfA);
        }

         //Final Gathering
         Vector newDirection = DiffuseBrdf.SampleDirection(eyePoint.Normal, rand);
         pathWeight = pathWeight * Vector.Dot(eyePoint.Normal, newDirection) / DiffuseBrdf.PdfW(eyePoint.Normal, newDirection);
         var finalGatherPoint = this.intersectionFinder.GetIntersectionPoint(new Ray(eyePoint.Position, newDirection));

         if (finalGatherPoint != null)
         {
            float searchRadius = 0.00634257728f * 2;
            var photons = photonmap.SearchPhotons(finalGatherPoint, searchRadius);
            float kernelFunction = 1.0f / (searchRadius * searchRadius * (float)Math.PI);
            foreach (var photon in photons)
            {
                 pixelColor += Vector.Mult(pathWeight, Vector.Mult(DiffuseBrdf.Evaluate(finalGatherPoint), photon.PathWeight) * kernelFunction);
            }
        }
    }

    return pixelColor;
}

Der Suchradius ist in diesen Beispiel fest im Quellcode hinterlegt. Er hängt von zwei Dingen ab. Welche Szene möchte ich rendern und wie viel Photonen sende ich aus? Ein Beispiel, wie man diesen Radius automatisch berechnen kann geht wie folgt. Erst werden N Photonen ausgesendet, um somit die Photonmap zu erstellen. Dann wird an 50 zufällig ausgewählten Pixeln ein Primärstrahl verschossen und geschaut ob er ein Dreieck in der Szene trifft. Wenn ja, wird dort dann nach den 5 nächsten Photonen gesucht (K-Nearest-Neighbor-Search). Das geht mit ein KD-Baum. Ein Beispiel dafür findest du hier https://github.com/XMAMan/KNearestNeighborSearch.
Ich schaue also, welchen Radius muss ich verwenden, um genau 5 Photonen zu bekommen. Für alle zufällig ausgewählten Pixel mache ich das. Dann nehme ich einfach den Mittelwert oder Median von den ganzen Radien, und das ist dann der Suchradius für meine Photonmap.

Nach nur 10 Samples erhalten wir dann folgendes Bild. Obwohl es viel weniger Samples als der Pathtracer verwendet, sie die Qualität viel besser aus. Außerdem ist es auf Grund der geringeren Sampelanzahl auf viel schneller als der Pathtracer.

https://github.com/XMAMan/RaytracingTutorials/blob/master/06_Photonmapping/Images/Vergleich.png?raw=true

Eine Lichtkomponente fehlt allerdings noch bei dem Photonmapping-Beispiel hier: Kaustiken. Um diese sehen zu können, benötigen wir aber neben ein diffusen Material ein spekulars Material. Damit ist Glas, Metall oder ein Spiegel gemeint. Also ein Material was Licht beim reflektieren nur wenig streut. Wenn man so ein Material hat, dann kann man beim verfolgen eines Photons nun schauen, auf was für ein Material es nach verlassen der Lichtquelle zuerst trifft. Wenn es auf eine diffuse Fläche zuerst auftrifft, dann ist das ein diffuses Photon. Trifft es zuerst auf ein spekulares Material auf, ist es ein Kaustik-Photon. Alle Kaustik-Photonen werden in einer eigenen Kaustik-Photonmap gespeichert. Will ich nun Kaustiken mit anzeigen, so muss ich beim Punkt P1 einfach nur aus der Kaustik-Map die Daten lesen und diese als dritte Komponente auf die Pixelfarbe drauf addieren. Ein Pixel setzt sich dann aus direktem Licht, indirekten Licht über Final Gathering und der Kaustik-Map zusammen.

Ich habe in diesen Beispiel auf Kaustiken verzichtet, da das Tutorial nicht zu umfangreich werden sollte. Ich denke, wenn jemand wirklich so weit gekommen ist, dass er das Photonmapping wie hier beschrieben umgesetzt hat, dann ist er auch in der Lage Kaustiken hinzubekommen.

Beispielquelltext fürs Photonmapping: https://github.com/XMAMan/RaytracingTutorials/tree/master/06_Photonmapping/Photonmapping-Source/RaytracingTutorials

XMAMan

Schlusswort

In diesen Tutorial habe ich mich hauptsächlich auf das Photonmapping und das dazu benötigte Vorwissen konzentriert. Ich habe dabei einige Dinge weggelassen, die man bereits ab den einfachen Raytracer implementieren kann. Ich halte es spätestens ab jetzt sinnvoll, das der Leser versucht, nun bei sein eigenen Raytracer diese Dinge zu implementieren.

-Beschleunigungsstruktur für die Strahl-Dreiecksliste-Schnittpunktabfrage
→ Siehe kd tree ingo wald
→ Siehe Bounding interval hierarchy
-Weitere Materialien: Glas, Spiegel
-Texturmapping
-Antialiasing Tent filter siehe http://www.kevinbeason.com/smallpt/ (Zeile 86/87)
-Tiefenunschärfe
-Flat-Shadding/Gouraud-Shadding
-Mit Threadprogrammierung das Rendern beschleunigen
-Den Raytracer so bauen, das er verschiedene Anzeigemodi hat (OpenGL/Pathtracer,Photonmapping)

Ich warne vor folgenden Taschenspielertricks:
-Fireflys (Einzelne helle Pixel) dadurch zu entfernen, indem man all die Farbwerte aus der Monte-Carlo-Summe raus nimmt, welche zu krasse Ausreißer haben. Spätestens wenn man anfängt sich mit Kaustiken zu beschäftigen, führt dieser Trick dazu, dass die Kaustiken falsch aussehen.
-Helligkeitsfaktoren, welche je nach Renderverfahren unterschiedlich Werte haben. Wenn ich ein Bild mit ein Pathtracer gerendert habe, dann muss das Bild mit Photonmapping / Photonmap genau so hell/dunkel sein wie der Pathtracer. Ist das nicht der Fall, dann hat man ein Fehler gemacht.

Folgendes Bild ist von mein eigenen Raytracer. Es zeigt, dass mit diffusen Material und Glas schon schöne Effekte möglich sind:
https://i.imgur.com/77QCQN9.jpg

Es wäre nicht schlecht, wenn Leute, die wirklich all das hier durchgearbeitet haben und ein paar von den hier aufgezählten Erweiterungsideen umgesetzt haben, dann hier auch mal ein Bild zeigen könnten. Als Motivation für die anderen^^

rapso

Habs noch nich durchgelesen, aber danke dass du dir die ganze muehe gemacht hast so ein detailiertes tutorial niederzuschreiben! habs gepinnt.

XMAMan

Aufgrund des Wunsches eines einzelnen Herren habe ich jetzt noch Lighttracing hinzugefügt. Ich habe wieder die Cornell-Box-Szene genommen und dabei versucht so minimal wie möglich alles zu halten. Kein Threading/Keine KD-Bäume.

https://github.com/XMAMan/RaytracingTutorials/tree/master/07_Lighttracing

https://github.com/XMAMan/RaytracingTutorials/blob/master/07_Lighttracing/RaytracingTutorials/Lighttracing.cs

https://github.com/XMAMan/RaytracingTutorials/blob/master/07_Lighttracing/Images/Lighttracing.png

Hier ist der Grundgedanke, dass man von der Kamera lauter Photonen aussendet und jedes mal, wenn ein Photon auf eine diffuse Fläche trifft, dann verbindet man den Punkt mit der Kamera und erstellt somit eine Menge von Pfadpunkten, welche von der Camera bis zur Lichtquelle reichen.

XMAMan

Da ich heute mal etwas Zeit habe, habe ich jetzt noch ein Pathtracer mit Next Event Estimation hinzugefügt. Dieser ist meiner Meinung nach der effektivste Algorithmus, wenn man neben dem Bildrauschen auch noch die Quellcode-Anzahl mit in Betracht zieht. Obwohl der Algorithmus nur so klein ist, erzeugt er trotzdem weitaus bessere Ergebnisse als der normale Pathtracer.

Grundidee ist hier: Man erzeugt zufällige Pfade zwischen Camera und Lichtquelle, indem man einerseits das normale Pathtracing verwendet und anderseits verwendet man LightSource-Sampling, welches auch schon beim Berechnen des Direkten Lichtes beim Photonmapping verwendet wurde. Nur dieses mal nimmt man das LightSourcesampling nicht nur für den ersten Punkt, den die Kamera sieht sondern für alle Punkte, den der Pathtracer erzeugt.

Ich habe also zwei Sampler zur Berechnung eines Farbwertes für ein Pixel: Brdf-Sampling und Lightsource-Sampling. Jeder Sampler erzeugt ein zufälligen Schätzwert für die Lichtmenge, die ein Pixel empfängt.

Ich bezeichne die beiden Farbwerte, welche mit die beiden Sampler für ein Pixel erzeugen jetzt mal mit f1 und f2. Die endgültige Pixel-Farbe könnte ich nun z.B: so berechnen:

f = (f1 + f2)/2.

Das wäre einfach nur der Durchschnittswert von Beiden. Dieses Idee würde zwar funktionieren aber sie erzeugt unnötig hohes Bildrauschen.

Eine bessere Idee ist der Einsatz von Multiple Importance Sampling. Dort gibt erzeugt man ein MIS-Faktor, welcher als Gewicht dient. Dieser Gewichtsfaktor ist eine Zahl zwischen 0 und 1 und gibt an, zu wie viel Prozent man das ein oder das andere Sampling-Verfahren in die Summe einfließen läßt. Unsere Formel sieht nun so aus:

f = (1-MIS) * f1 + MIS * f2

Dieser MIS-Faktor wird dadurch berechnet, indem man die Wahrscheinlichkeitsdichte im Bezug zum Flächenmaß (Pdf with Respect to Surface Area = PdfA) verwendet. Dieses PdfA ist eine Funktion, welche für einen konkret gesampelten Pfad P angibt, wie hoch die Wahrscheinlichkeit ist, dass genau dieser Pfad erzeugt wurde. Ich verwende hier Brdf- und LightSource-Sampling. Also habe ich für jedes dieser beiden Sampling-Verfahren auch eine zugehörige PdfA-Funktion:

Brdf-PdfA(P)
LightSampling-PdfA(P)

Diese Funktionen bekommen ein Pfad (Menge von Punkten) als Input und sagen dann: Wie groß ist die Wahrscheinlichkeit, dass ich den Pfad P mit Brdf/Light-Sampling erzeuge?

Ich verwende nun also zwei Sampler und jeder dieser Sampler erzeugt mir jeweils ein Pfad. Brdf-Sampling erzeugt mir den Pfad P1 und LightSourceSampling den Pfad P2.

Der Farbwert(Radiance) für die Pfade wird PathContribution genannt und wird über die PathContribution-Funktion berechnet. Also Contribut(P1) und Contribut(P2).

Jeder dieser beiden PathContribution-Werte muss nun mit den MIS-Factor gewichtet werden:

PixelColor = MIS1 * Contribut(P1) / Brdf-PdfA(P1) + MIS2 * Contribut(P2) / LightSampling-PdfA(P2)

Ich berechne MIS1 nun über folgende Formel, wenn ich davon ausgehe, dass P1 mit BrdfSampling erzeugt wurde:

MIS1 = Brdf-PdfA(P1) / (Brdf-PdfA(P1) + LightSampling-PdfA(P1))

Bei MIS2 gehe ich davon aus, dass P2 mit LightSampling erzeugt wurde.

MIS2 = LightSampling-PdfA(P2) / (Brdf-PdfA(P2) + LightSampling-PdfA(P2))

Die Idee von Multiple Importance Sampling ist also, dass man den Sampling-Verfahren mehr Gewicht verleiht, welche eine größere PdfA haben da großer PdfA-Wert wenig Rauschen bedeutet.

Diese Formel (PixelColor=...), wo die Mis-Gewichtet Summe von mehreren Verfahren gebildet wird, nennt sich Secondary-Estimator.

Der Primary-Estimator ist dann die Durchschnittswert von N Secondary-Estimator-Werten. Dieser Wert ist es dann, welcher am Ende in die Bilddatei geschrieben wird.

Der Quelltext für dieses Pathtracing, wo Brdf-Sampling mit Lightsource-Sampling via Multipe Importance Sampling (MIS) vereint, gezeigt wird, habe ich hier:

https://github.com/XMAMan/RaytracingTutorials/blob/master/08_PathtracingNextEventEstimation/RaytracingTutorials/PathtracingNEE.cs

Hinweis: Wenn die PathContribution-Funktion den Wert 0 liefert, dann mache ich mir garnicht erst die Mühe das zugewhörige MIS-Gewicht zu bilden.

Ausgabebild:
https://github.com/XMAMan/RaytracingTutorials/blob/master/08_PathtracingNextEventEstimation/Images/PathtracingWithNextEventEstimation.png

Um das ganze besser zu verstehen, wäre ein Blick in die Dissertation von Eric Veach hilfreich. Er setzt aber vorraus, dass man den Begriff 'Wahrscheinlichkeitsdichte im Bezug zu ein Maß' kennt. Ich denke, dass wäre mal ein eigenes Thema, was ich mal gesondert erklären müsste, wenn Bedarf besteht^^