Raytracing / Pathtracing / Photon Mapping - Tutorial

XMAMan

Der GeometryTerm

Der GeometryTerm ist ein Faktor (Float-Zahl), um dem die Lichtintensität abgeschwächt wird, wenn zwei Flächen Licht miteinander austauschen. Habe ich zwei Flächen, dessen Mittelpunkte den Abstand d haben und welche beide jeweils ihre eigene Flächenormale N1 und N2 haben, so ist der GeometryTerm wie folgt definiert.

$\mathit{GeomtryTerm}=\frac{\cos (\mathit{w1})\ast \cos (\mathit{w2})}{d\ast d}$

w1 und w2 ist der Winkel zwischen der Flächennormale und der Linie, welche zur jeweils anderen Fläche zeigt.

https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/02_DirectLightingOnly/Images/GeometryTerm.png

Wenn die Flächen genau aufeinander zuzeigen, so ist w1 und w2 0 Grad. Der Cosinus von 0 ist 1. D.h. Es wird mehr Licht zwischen den Flächen übertragen, wenn die Flächen aufeinander zu zeigen, da beide Cosinus-Terme dann 1 sind.

Den Grund für diese Cosinus-Abschwächung versuche ich mit folgenden Bild zu erklären:
https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/02_DirectLightingOnly/Images/CosinusLambda.png

Um so schräger die Fläche das Licht reflektiert oder versendet, um so kürzer wird l1/l2/l3. Wenn ich an den roten Punkten mich befinde und auf die Fläche schaue, so erscheint mir die Fläche immer kleiner, je schräger die Fläche zur mit ist.

Außerdem wird die Übertragene Lichtmenge mit den Quadrat-Abstand abgeschwächt. Um so größer der Abstand, um so größer ist die Abschwächung also.

Den Grund für die Abstandabschwächung versuche ich hiermit zu erklären:
https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/02_DirectLightingOnly/Images/DistanzeFaktor.png

Die gelben Linien stehen für ein jeweils versendetes Photon. F1, F2 und F3 sind Flächen, die alle gleichgroß sind und sich nur im Abstand unterscheiden. F1 erhält 5 Photonen. F2 3 und F3 ein Photon.

Bauen wir den GeometryTerm nun mal in den Raytracer.

var lightPoint = data.LightSource.GetRandomPointOnLightSource(rand);
if (IsVisible(eyePoint.Position, lightPoint.Position))
{
   pixelColor += eyePoint.Color * GeometryTerm(eyePoint, lightPoint);
}

→ Wenn man zwei normalisierte Vektoren multipliziert, dann entspricht die float-Zahl, die dabei rauskommt dem Cosinus von dem Winkel zwischen den beiden Vektoren. Somit ist die Vektor-Punkt-Multiplikation ein guter Weg, um bei Richtungsvektoren den Cosinus-Winkel zu berechnen.

float GeometryTerm(IVertexPoint point1, IVertexPoint point2)
{
    Vector v1To2 = point2.Position - point1.Position;
    float distance = v1To2.Length();
    Vector direction1To2 = v1To2 / distance;
    float lamda1 = Math.Max(Vector.Dot(point1.Normal, direction1To2), 0);
    float lamda2 = Math.Max(Vector.Dot(point2.Normal, -direction1To2), 0);
    return lamda1 * lamda2 / (distance * distance);
}

Nun sieht das Bild schon so aus:
https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/02_DirectLightingOnly/Images/WithGeometryTerm.png

Es sieht nun schon ein wenig besser aus.

Beispielcode (In SimpleRenderer.cs ist das Version 2)

https://github.com/XMAMan/RaytracingTutorials/tree/master/02_DirectLightingOnly/RaytracingTutorials

XMAMan

Diffuse Brdf

Der nächste Schritt wird erst mal nicht sichtbar die Bildqualität verbessern aber er ist trotzdem für später wichtig. Es geht um die Brdf (Bidirektionele Reflektranzverteilungsfunktion).
https://de.wikipedia.org/wiki/Bidirektionale_Reflektanzverteilungsfunktion

Damit wird die Farbe und das Material (Glass, Spiegel, Diffuse Wand, Raues Metall) festgelegt. Hier in unseren Beispiel werde ich erst mal nur diffuse Materialien behandeln, da diese immer das erste Material sind, womit man sich beschäftigt^^

Also was ist eine Brdf? Das ist eine Funktion, welche angibt, wie stark ein Material eine gegebenen Lichtstrahl reflektiert.

Gegeben ist der Richtungsvektor für das einfallende und das ausgehende Licht.
https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/02_DirectLightingOnly/Images/Brdf.png

Die Brdf-Funktion sagt mir nun, wie viel Licht in Richtung des ausgehenden Vektors fließt. Die Lichtmenge wird üblicherweise als Vektor angegeben, wo die einzelnen Komponenten für RGB (Rot Grün Blau) stehen.

Wenn ich jetzt ein weißen Lichtstrahl von der Lichtquelle habe, so könnte der Lichtvektor z.B.
so aussehen (1,1,1). Ich gebe nun der Brdf-Funktion folgende Werte rein:
Welches Licht geht rein? -> (1,1,1)
Richtungsvektor, welcher von der Lichtquelle zum Brdf-Punkt zeigt
Richtungsvektor, welcher vom Brdf-Punkt zu ein beliebigen anderen Punkt (z.B. die Kamera) zeigt
Wenn ich nun eine diffuse Rote Wand habe, dann sieht die Brdf-Funktion wie folgt aus:

Vector RedDiffuseBrdf(Vector inputDirection, Vector outputDirection, Vector inputColor)
{
    return Vector.Mult(inputColor, new Vector(1, 0, 0) / Math.PI);
}

Die Farbe, die reingeht wird mit (1,0,0) / PI gewichtet.
D.h. Das weiß Licht (1,1,1) wird zu (0.32, 0,0)

Man sieht hier, dass die inputDirection und outputDirection egal sind und in der Funktion selber nicht verwendet werden. Das liegt daran, dass eine diffuse Fläche in alle Richtungen gleichstark reflektiert. Wenn ich eine andere Brdf habe, dann könnten diese beiden Angaben auf einmal wichtig werden. Hier bei der diffuse Funktion kann man sie aber auch weglassen.

Der Vektor (1,0,0) steht für die rote Farbe. Aber warum teilt man eigentlich durch PI? Würde es nicht reichen, wenn ich die weiße Farbe nur mit (1,0,0) multipliziere?

Die Antwort auf diese Frage ergibt sich aus der Bedingung, dass eine Brdf nicht mehr Licht reflektieren darf, als wie sie als Input bekommt (Energieerhaltungsregel).

Die Energieerhaltungsregel sieht wie folgt aus (siehe hier)https://de.wikipedia.org/wiki/Bidirektionale_Reflektanzverteilungsfunktion

$\forall w\_i:\int \_{\Omega }f\_r(w\_i,w\_o)\ast \cos \phi\_o\mathit{dw}_o\le 1$

Anmerkung: Das ist ein Lebesgue-Integral (https://de.wikipedia.org/wiki/Lebesgue-Integral). Ω
ist die Menge aller Richtungsvektoren, welche man innerhalb einer Halbkugel hat. 'wo' läuft über jedes Element aus Ω und es wird so die Summe gebildet.

Meiner Meinung nach muss man die Formel aber noch um den Li-Term erweitern, um sie richtig verstehen zu können^^

$\forall w\_i:\int \_{\Omega }f\_r(w\_i,w\_o)\ast L\_i\ast \cos \phi\_o\mathit{dw}\_o\le L_i$

Für jede Inputrichtung wi muss gelten: Die Summe aller ausgesendeten Energie darf nicht größer sein, als die aus wi kommende Lichtmenge Li. Fr ist hier die Brdf-Funktion. Die einkommende Menge Li wird mit dem Brdf-Faktor fr und dem cos-Faktor (stammt aus dem Geometry-Term) gewichtet, um somit Lo (ausgesendete Lichtmenge in Richtung wo) zu berechnen.
Das Integral läuft über alle möglichen wo-Richtungen und bildet darüber die Summe. Wenn also eine Fläche lediglich von einer einzigen Richtung wi Li Lichtenergie bekommt. So sagt mir dieses Integral dann, wie viel Lichtenergie insgesamt die Fläche ausstrahlt. Und diese Menge darf eben nicht mehr sein, als wie sie empfängt.

Wenn wir jetzt nun für fr nun unsere diffuse Brdf sein soll und wir aber nicht wissen, wie sie dann nun aussehen muss so setzen wir erstmal für fr die Zahl 1 ein. D.h. Die Fläche reflektiert so viel Licht wie möglich. Es muss jetzt nur noch geprüft werden, ob die <=Li-Gleichung weiterhin gilt.

Um diese Gleichung überprüfen zu können, muss man den dw-Term verstehen. Der dw-Term ist die türkisfarbene Fläche, welche im nachfolgenden Bild gezeigt wird und mit dΩ hier bezeichnet wird.
https://www.scratchapixel.com/images/upload/shading-intro/shad-light-beam7.png?

Ich ersetze dw durch sin(ɵ)dɵdɸ und fr(wi,wo) durch pd (Diffuse-Konstante, welche angibt, wie viel Prozent reflektiert wird).

$\forall w\_i:\int \_{\theta=0}^{2\pi}\int _{\phi=0}^{\pi/2}p\_d\ast L\_i\ast \cos \phi_o\sin (\phi\_o)\mathit{d\phi}\_o\mathit{d\theta}\_o\le L\_i$

Ich schreibe den Li-Term und pd vor die Integrale (beides Konstanten) und ordne nach ɸ und ɵ

$\forall w\_i:L\_i\ast p\_d\ast \int \_{\theta=0}^{2\pi}\mathit{d\theta}_o\ast \int _{\phi=0}^{\pi/2}\cos \phi_o\sin (\phi\_o)\mathit{d\phi}\_o\le L_i$

Das ɸ-Integral löst sich so:

$\int _{\theta=0}^{2\pi}\mathit{d\theta}=[\theta]_0^{2\pi}=2\pi$

Für das zweite Integral finden wir die Stammfunktion mithilfe von Google^^
https://www.zum.de/Faecher/M/NRW/pm/mathe/stammfkt.htm

$\int _{\phi=0}^{\pi/2}\cos \phi\sin (\phi)\mathit{d\phi}=[\frac 1 2\ast \sin ^2\phi]_0^{\pi/2}=\frac 1 2\ast \sin ^2(\pi/2)-\frac 1 2\ast \sin ^2 0=\frac 1 2-0=\frac 1 2$

Ich setze die beiden aufgelösten Integrale in die Ausgangsformel ein:

$\forall w\_i:L\_i\ast p\_d\ast 2\pi\ast \frac 1 2\le L\_i$

Ich kürze Li weg und multipliziere die Zahlen aus:
$\forall w\_i:p\_d\ast \pi\le 1$

Damit diese Gleichung immer erfüllt ist, darf pd nie größer als 1/П werden.

Will ich nun eine rote diffuse Fläche haben, dann darf die Rot-Komponente maximal 1/П reflektieren. Würde sie mehr reflektieren, dann würde die Lichtmenge in der Szene immer heller werden und wir erhalten ein unrealistisches Bild.

So wird nun der Raytracer um die Brdf erweitert:

static class DiffuseBrdf
{
    public static Vector Evaluate(IntersectionPoint point)
    {
        return point.Color / (float)Math.PI;
    }
}

pixelColor += Vector.Mult(DiffuseBrdf.Evaluate(eyePoint), eyePoint.Color * GeometryTerm(eyePoint, lightPoint));

https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/02_DirectLightingOnly/Images/WithBrdf.png

Das wirkt jetzt erst mal wie ein billige Verbesserung. Einfach nur die Farbe mit den Faktor PI wichten. Das ist nun wirklich kein Hexenwerk aber in ein späteren Abschnitt wird das hier noch wichtig und dann werden wir froh sein, dass wir diese PI-Division haben.

Beispielcode (In SimpleRenderer.cs ist das Version 3)

https://github.com/XMAMan/RaytracingTutorials/tree/master/02_DirectLightingOnly/RaytracingTutorials

XMAMan

Schattenkanten verbessern

Die nächste Verbesserung wird nun das pixelige aussehen von den Schattenkanten verbessern. Dazu muss man lediglich wie folgt vorgehen:

for (int x = 0; x < width; x++)
 for (int y = 0; y < height; y++)
 {
     int sampleCount = 10;
     Vector sum = new Vector(0, 0, 0);
     for (int i = 0; i < sampleCount; i++)
     {
         sum += EstimatePixelColor(x, y, rand);
     }
     sum /= sampleCount;
     image.SetPixel(x, y, VectorToColor(sum));
 }

Anstatt pro Pixel nur einen Farbwert zu berechnen mache ich die vorher beschriebene Berechnung 10 mal. Bei jeden Rechenschritt bestimme ich ein neuen zufälligen Punkt auf der Lichtquelle.

So sieht das Bild dann aus, wenn ich mit verschiedenen Sample-Schritten arbeite:

https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/02_DirectLightingOnly/Images/SoftShadows.png

Beispielcode (In SimpleRenderer.cs ist das Version 4)

https://github.com/XMAMan/RaytracingTutorials/tree/master/02_DirectLightingOnly/RaytracingTutorials

Das ist nun der einfache Raytracer, welcher nur direktes Licht in Betracht zieht. Wenn man sich die Schatten ansieht, dann fällt auf, das diese schon ziemlich dunkel sind. Der Grund dafür ist, dass wenn ein Punkt im Schatten liegt, dann wird die Pixelfarbe fix auf (0,0,0) gesetzt. Wenn ich aber auch indirektes Licht möchte, dann benötige ich ein Algorithmus, welcher globale Beleuchtung genannt wird. Im nächsten Abschnitt werde ich darauf näher eingehen.

XMAMan

Globale Beleuchtung

Bei der globalen Beleuchtung wird berechnet, wie viel Licht zu ein gegebenen Punkt sowohl direkt von der Lichtquelle, als auch indirekt über die Reflektion über die Wände, gesendet wird.

Hier sind einige Flächen und in gelb die Pfade, die ein Photon zurück legt, dargestellt.
https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/03_MonteCarloIntegration/Images/LightPaths.png

Bei direkten Licht besteht ein Lichtpfad entweder aus 2 Punkten (Kamera-Punkt + Punkt auf Lichtquelle) oder aus 3 Punkten (Kamera-Punkt, Lichtquellenpunkt, Punkt in Szene).

Ein indirekter Lichtpfad besteht aus mehr als 3 Punkten.

Wenn ich nun berechnen will, wie viel Licht über so ein Pfad von der Lichtquelle bis zur Kamera fliegt, dann muss ich das Pfadgewicht (PathContribution) für ein gegebenen Lichtpfad berechnen.

Dieses Pfadgewicht ist ein Vektor. Er gibt für die RGB-Farben an, wie viel Licht er jeweils durchlässt. Das Pfadgewicht errechnet sich aus dem Produkt der GeometryTerme (Float-Zahl) und der Brdfs (RGB-Vektor).
https://raw.githubusercontent.com/XMAMan/RaytracingTutorials/master/03_MonteCarloIntegration/Images/PathContribution.png

Wenn ich nun ein Bild mit globaler Beleuchtung berechnen will, so muss ich alle möglich Lichtpfade, die man zwischen der Kamera und der Lichtquelle aufstellen kann, berechnen und für diese jeweils ihren Anteil zur Pixelfarbe hinzuaddieren.

Wenn ich die Summe über alle möglichen Lichtpfade bilde, die man bilden kann, dann weiß ich die Farbe für jedes Pixel und kann somit das Bild ausgeben.

Es gibt dabei folgende Schwierigkeiten. Ein einzelnen Punkt aus dem Lichtpfad ist ja ein 3D-Punkt, dessen Koordinaten über 3 Float-Zahlen angegeben werden. Dieser 3D-Punkt liegt auf ein Dreieck in der Szene. Wenn ich nun alle möglichen Punktpositionen erzeugen will, die es in der jeweiligen Szene gibt, dann sind das ja schier unendlich viele. Die Schleife, welche die Summe bildet, hätten also extrem viele Durchläufe und das würde sehr lange dauern. Außerdem gibt es ja unendlich viel verschiedenen Lichtpfadlängen. Erst erzeuge ich alle möglichen Lichtpfade mit 2 Punkten. Wenn ich die Summiert habe, dann erzeuge ich alle möglichen Pfade mit 3 Punkten. So geht das immer weiter bis ich bei einer unendlichen langen Pfadlänge bin.

Wie kann man diese Aufgabe in endlicher Zeit lösen?

Es gibt zwei Ansätze (die ich kenne):
Erster Ansatz (Radiosity): Ich unterteile die Szene in lauter Dreiecke und/oder Vierecke (Patches genannt) und ein Lichtpfad-Punkt darf immer nur in der Mitte von ein Patch liegen. Somit ist die Menge der möglichen Lichtpfad-Positionen schon stark eingeschränkt. Weiterhin beschränke ich die maximale Pfadlänge auf z.B. 10. Nun gibt es viel weniger Lichtpfade, die ich summiere muss und somit kann ich es schaffen.

Zweiter Ansatz(Monte Carlo Integration): Ich presse den gesamten Beleuchtungsalgorithmus in folgende Formel:

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

f ist hier die PathContribution-Funktion. Xi ist zufällig generierter Lichtpfad (Array von 3D-Punkten). N ist die Anzahl der Samples. Pdf ist die Probabilty density function.

In den folgenden Abschnitten gehe ich näher auf die Monte Carlo Integration ein. Hier in diesen Abschnitt wollte ich erst mal nur klären, was ein Lichtpfad ist und wie man sich das visuell vorstellen kann.

XMAMan

Monte Carlo Integration das Würfelbeispiel

Im vorherigen Abschnitt habe ich erklärt, das man sich den globalen Beleuchtungsalgorithmus als eine Summenfunktion über alle Lichtpfade vorstellen kann. Und wenn wir an Summe denken, da denken wir auch an Integrale^^ In diesen Abschnitt kläre ich nun, wie man mithilfe der Zufallsfunktion ein Integral/Summe lösen kann.

Als Beispiel nehme ich hier die Aufgabe, dass man die Summe der Augenzahlen von ein Würfel berechnen will.

Normalerweise würde man die Summe ja so lösen: Summe = 1 + 2 + 3 + 4 + 5 + 6 = 21

Mit Monte Carlo benutzt man nun folgende Formel zum lösen dieser Summe.
$\frac 1 N\sum _{i=1}^N\frac{f(\mathit{xi})}{\mathit{pdf}(\mathit{xi})}$

Dazu definiere ich:
x = Eine Integerzahl im Bereich von 0 bis 5. Ich generiere x zufällig mithilfe der Rand-Funktion.
f(x) = x + 1 → f ist hier unser Würfel.
pdf(x) = Das ist die Wahrscheinlichkeit, dass ich die Zahl x würfle. D.h. pdf(x) = 1/6
N = Die Sampleanzahl. Eine beliebige Integerzahl, welche ich hier z.B. auf 10 festlege.

//Summe aller Würfelzahlen mit 10 Samples → Ergebnis: sum1 = 21 	sum2 = 25,2
public static string Example1()
{
    Random rand = new Random(0);

    int[] numbers = new int[] { 1, 2, 3, 4, 5, 6 };

    int sum1 = numbers.Sum(); // 21

    //........

    float sum2 = 0;
    int N = 10; //Samplecount
    float pdf = 1.0f / 6;
    for (int i = 0; i < N; i++)
    {
        int x = rand.Next(6); //Erzeuge Zufallszahl im Bereich von  0-5
        sum2 += numbers[x] / pdf;  //numbers entspricht f(x) = x + 1
    }
    sum2 /= N;

    string result = sum1 + " " + sum2;
    return result;
}

Der Grundidee ist also. Erzeuge eine zufällig int-Variable x und addiere dann auf die Summenvariable f(x) / pdf(x).

Bei 10 Samples sagt der Monte Carlo-Algorithmus, dass die Summe der Zahlen von 1-6 25.2 sei. Das ist natürlich Quatsch. Es muss 21 rauskommen. Was ist da also los? Das Ergebnis von diesen numerischen Integrationsprozess ist eine Zufallszahl. Diese Zahl liegt irgendwo in der Nähe vom richtigen Ergebnis. Um so mehr Samples ich nehme, um so genauer ist das Ergebnis. Nehme ich anstatt 10 100 Samples, dann erhalte ich 22.98 und nach 1000 Samples 20.7.
Ich will doch nur 6 Zahlen summieren. Warum muss ich dann auf einmal 1000 Zahlen summieren und mit der rand-Funktion arbeiten, wenn das Ergebnis noch dazu nicht perfekt ist? Nun, bei einfachen Funktionen, wo eine float-Zahl als Input reingeht und eine float-Zahl als Ergebnis kommt , also eine float->float-Abbildung (das gleiche gilt für eine Int->Int-Funktion), da ist Monte Carlo in der Tat eher dürftig. Aber man kann damit auch Integrale von Funktionen lösen, wo z.B. 10 oder 100 float-Zahlen als Input reingehen und 1-3 float-Zahlen als Ergebnis rausgehen.

Genau diesen Fall hat man bei den Lichtpfaden gegeben. Ein Lichtpfad ist eine Menge von 3D-Punkten. Jeder Punkt besteht aus X/Y/Z (float). Habe ich 10 Punkte, dann sind das 30 Float-Zahlen als Input. Die PathContribution-Function errechnet aus diesen 30-Floatzahlen einen RGB-Vektor (3 Floatzahlen). Möchte ich nun die Summe über schier unendlich viele Lichtpfade bilden, dann sind auf einmal 1000 Samples aus ein unendlichen großen Bereich eine schön geringe Untermenge, welche man in endlicher Zeit summieren kann.

Das ist also die Motivation, warum wir uns mit Monte Carlo Integration beschäftigen wollen.

XMAMan

Monte Carlo Integration Cosinus-Integral ohne Importancesampling

In den ersten Würfel-Beispiel habe ich mit Integerzahlen gearbeitet. Ein Lichtpfad besteht aber aus float-Zahlen. Aus diesem Grund möchte ich im nächsten Beispiel nun das Integral von der Funktion cos(x) im Bereich 0 bis PI/2 lösen.

Schauen wir uns die Funktion mal an.

Gesucht ist die Flächeninhalt von der rot schraffierten Fläche.

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

In der Schule würde man das Integral dadurch lösen, indem man die Stammfunktion von cos(x) bildet. Das wäre sin(x). Nun Berechne ich sin(PI/2) – sin(0) = 1.

Das wäre der 'einfache' Weg. Kommen wir nun zum guten alten Monte-Weg.

Diesmal muss ich zufällige x-Werte im Bereich von 0 bis PI/2 generieren. Die Pdf(x) wäre diesmal 1 / (PI/2). Warum ist das so? An die Pdf (Probability density function) gilt die Anforderung, dass die summe aller möglichen pdf-Werte im 0-PI/2-Bereich 1 ergeben muss. Die Pdf gibt an, mit welcher Wahrscheinlichkeit ich eine bestimmte Float-Zahl erzeuge.

Gesucht ist also ein Funktion 'pdf' welche folgende Bedingung erfüllt:

$\int _0^{\pi/2}\mathit{pdf}(x)\mathit{dx}=1$

Jede dieser x-Zahlfallszahlen, welche ich erzeuge sample ich mit gleich hoher Wahrscheinlichkeit. Das Bedeutet, dass die Pdf erst mal eine Konstante sein muss. Eine Konstante Funktion wiederum ist eine Horizontal verlaufende Linie, wenn man diese zeichnen will. Das Integral von einer Konstanten Funktion ist also lediglich der Flächeninhalt von ein Rechteck. Wir wissen über dieses Rechteck, dass der Flächeninhalt 1 sein muss und das eine Seite von diesen Rechteck die Länge PI/2 hat.

Dann stelle ich jetzt mal diese komplizierte Gleichung auf.

$1=a\ast \frac{\mathit{PI}} 2$

a ist gesucht. Pdf(x) = a. Stelle ich nach a um, erhalte ich pdf(x) = 2 / PI

In C# sieht das ganze nun so aus:

Random rand = new Random(0);
float sum = 0;
int sampleCount = 10;
for (int i = 0; i < sampleCount; i++)
{
    float x = (float)(rand.NextDouble() * Math.PI / 2);
    float pdf = 1.0f / (float)(Math.PI / 2);
    sum += (float)Math.Cos(x) / pdf;
}

sum /= sampleCount; 
//sum = 0.8077739

Bei 10 Samples erhalte ich 0.8077 Nehme ich 1000 Samples, erhalte ich 1.01155.

Beispielquelltext (siehe Example2)

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

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^^