Die allgemeine Transportgleichung ist unabhängig von der

Die allgemeine Transportgleichung ist unabhängig von der
physikalischen Situation. Soweit wurden lediglich wenige
Annahmen über spontane Kollisionen, monoenergetischen
Transport und vernachlässigbare Interaktionen zwischen Partikeln,
unendlich kleine und zahlreiche Partikel, keine externen Kräfte und
den Rand des Gebietes, mit dem Interaktion stattfindet gemacht.
Wir spezifizieren jetzt die Beschreibung für Photonen um die
Mario Hlawitschka
Digitale Signalverarbeitung
69
Gleichungen auf den Strahlungstransport anzuwenden.
Strahlungstransport ist hierbei synonym mit Photonentransport.
Wir müssen dabei drei Eigenschaften berücksichtigen:
1
Photonen reisen mit konstanter Geschwindigkeit. Dies
korrespondiert zu der Ein-Geschwindigkeits-Annahme.
2
Jedes Photon hat eine Frequenz ν, die die interne Energie
bestimmt. Sie ist abhängig von der Wellenlänge
c = 2.997 924 58 × 108 m·s−1 = νλ . Die interne Energie ergibt
sich also zu
c
E = hν = h [J]
λ
wobei h = 6.626 176 × 10−34 J·s das Planck’sche
Wirkungsquantum ist.
3
Photonentrajektorien sind stark beeinflusst durch Oberflächen;
eine Eigenschaft, die mit Gasmolekülen geteilt wird, jedoch
nicht mit Neutronen
Mario Hlawitschka
Digitale Signalverarbeitung
70
Da wir das System als stationär angenommen haben, ist der Effekt
der Phosphoreszenz ausgeschlossen. Zusätzlich schließen wir
Fluoreszenz aus, das heißt, wir betrachten die Frequenzen getrennt.
Da wir lediglich monoenergetisches Licht betrachten, müssen alle
Gleichungen für alle Wellenlängen separat betrachtet werden.
Wir übersetzen nun die abstrakte Anschauung des Phasenflusses in
das radiometrische Konzept der Radianz.
Mario Hlawitschka
Digitale Signalverarbeitung
71
Radianz ist die Energie P pro projizierter Einheitsfläche pro
Raumwinkel (oder die Energie, die die Oberfläche senkrecht
durchströmt). Die Beziehung zwischen Radianz L, ausgestrahlter
Energie Lev und Quellstrahlung Le zu Phasenfluss Φ ist
L(x, ω) =
d 2 P(x, ω)
dx · dω
=
Lev (x, ω) =
Le (x, ω) =
c
h · Φ(x, ω) [W·m−2 ·sr−1 ]
λ
c
h · qv (x, ω) [W·m−3 ·sr−1 ]
λ
c
h · q(x, ω) [W·m−2 ·sr−1 ]
λ
Dabei ist die Emission ein Phänomen der Materie und trifft nur auf
Volumen zu. Der Volumenstreukern ist über die Phasenfunktion fs
beschrieben durch
k(x, t) ≡ fs (x, t) = σs (x)
φ (t)
4π
[m−1 ·sr−1 ]
wobei t ∈ [−1, 1] den Kosinus der Winkel zwischen den Richtungen
(t)
vor und nach der Streuung darstellt. φ4π
ist normiert, integriert
also zu 1.
Mario Hlawitschka
Digitale Signalverarbeitung
72
Auf dem Rand wählen wir die bidirektionale
Reflexions-Verteilungsfunktion (bidirectional reflectance
distribution function, BRDF) fr
1
kb (ωi , x, ω) ≡ fr (ωi , x, ω) cos θi =
π
n
!
∑ ρj (x)Pr ,j (ωi , x, ω)
j=1
cos θi
B in der Einheit [sr−1 ], wobei der Faktor cos θi = hn(x), ωi von
dem Differenzial
dx · dω = hn(x), ωidxdω
kommt.
Der Faktor cos θ dω wird auch projizierter Raumwinkel genannt.
Mario Hlawitschka
Digitale Signalverarbeitung
73
Die Helmholz-Reziprogität garantiert auch hier
fr (ωi , x, ω) = fr (−ω, x, −ωi )
Mario Hlawitschka
Digitale Signalverarbeitung
74
Gegeben sei (Lev , fs , σs , σa , Le , fr , V ). Das globale
Beleuchtungsproblem ist nun die stationäre Lösung der folgenden
Gleichung für L zu finden:
Z
d
L(x, ω) + ωt (x)L(x, ω) = Lev +
fs (x, ωi · ω)L(x, ωi )dωi
dω
S2
Zfür x ∈ V
L(x, ω) = Le (x, ω) +
2 (x)
S−
fr (ωi , x, ω)L(x, ωi ) cos θi dωi
für x ∈ ∂ V
Mario Hlawitschka
Digitale Signalverarbeitung
(1)
75
BSSRDF
Für einige Simulationen ist es sinnvoll auch bidirektionales
Scattering zu berücksichtigen. Dazu definiert man die bidirectional
”
scattering surface reflectance distribution function“ (BSSRDF) S,
die die Menge an Strahlung charakterisiert, die eine Oberfläche im
Punkt x0 in Richtung ω0 unter der Annahme verlässt, dass
Strahlung im Punkt xi in die Oberfläche aus Richtung ωi eintritt.
L(x0 , ω0 ) =
Z Z
2 (x )
A S−
i
S(x0 , ω0 , xi , ωi )Li (xi , ωi ) cos θi dωi dA(xi )
wobei wir über die gesamte Oberfläche und alle einfallenden
Richtungen integrieren müssen. Im Falle der BRDF ist der Term
nur für x0 = xi ungleich Null (δ -Distribution) und das Integral
kann über den den Wert bei x0 ausgewertet werden.
Mario Hlawitschka
Digitale Signalverarbeitung
76
Strahlen und Strahlschnitte
Bevor wir die Integro-Differenzialgleichung in eine Integralform
überführen, führen wir ein paar Notationen bezüglich der
Sichtbarkeit bezüglich Lichtstrahlen ein.
Wir definieren die Distanzfunktion
d:
R3 × S 2
→
R+
(x, ω) → inf{s > 0 : x + sω ∈ ∂ V }
und den Trefferpunkt (hit point)
h(x, ω) := x + d(x, ω)ω
woraus folgt, dass h(x, ω) ∈ ∂ V und der nächste Punkt in
Strahlrichtung auf der Oberfläche ist.
h ist immer wohldefiniert — zur Not durch Berandung des
Volumens mit einer nicht reflektierenden Oberfläche.
Mario Hlawitschka
Digitale Signalverarbeitung
77
Sichtbarkeit
Eine weitere Operation ist der gegenseitige Sichtbarkeitstest
(
1 falls |y − x| ≤ |x − h(x, ωx,y )|
V (x, y ) :=
0 sonst.
Hierbei müssen weder x noch y auf dem Rand liegen. Der Test
gibt an, ob die beiden Punkte gegenseitig sichtbar sind oder nicht,
x−y
wobei ωx,y = |x−y
| ist und stets gilt V (x, y ) = V (y , x). Die Einheit
von V ist 1 sr.
Mario Hlawitschka
Digitale Signalverarbeitung
78
Die Transformation I
Da die meisten Lösungstechniken für den Lichttransfer auf
numerischer Integration basieren, transformieren wir die
Integro-Differenzialgleichung des Strahlungstransports in eine reine
Integralschreibweise. Dazu definieren wir zuerst die Radianz Lv in
einem Punkt im Volumen basierend auf der Volumenemission und
der Einstreuung, und die Radianz Lb auf dem Rand basierend auf
Emission und Reflexion:
Lv (x, ω) := Lev (x, ω) +
Lb (x, ω) := Le +
Mario Hlawitschka
Z
2 (x)
S−
Z
S2
fs (x, ωi · ω)L(x, ωi )dωi ∀x ∈ V
fr (ωi , x, ω)L(x, ω) cos θ dωi ∀x ∈ ∂ V
Digitale Signalverarbeitung
79
Die Transformation II
Die Grundidee der Transformation ist jetzt, Licht entlang eines
Strahles zu integrieren bis es den Rand erreicht. In dieser
Formulierung wird der Differentialoperator eliminiert und der Wert
auf dem Rand des Volumens mit berücksichtigt. Schlüsselidee ist,
die Richtungsableitung in x umzuschreiben
d
d
d
L(x, ω) = L(x + sω, ω)|s=0 = − L(x, −sω, ω)|s=0
dω
ds
ds
wobei s den Strahl parametrisiert. Vernachlässigt man erst einmal
die Randbedingungen, so kann Gleichung 1 umgeschrieben werden
in
d
L(x − sω, ω) − σt (x − sω)L(x − sω, ω) = −Lv (x − sω, ω) (2)
ds
was für alle Punkte x − sω entlang des Strahls (x, ω) gelten muss.
Der Punkt (x, ω) ist dabei ein fester Punkt im Phasenraum. Man
Mario Hlawitschka
Digitale Signalverarbeitung
80
Die Transformation III
beachte, dass Lv keine explizite Abhängigkeit bzgl. L hat. Die
Abhängigkeit von ω verschwindet durch die Integration. Obwohl L
in dem Einstreuungsterm enthalten ist, trägt es dort nur mit Maß
Null bei. Die Gleichung ist also eine lineare Differenzialgleichung
erster Ordnung in L mit variablen Koeffizienten, was einfach zu
lösen ist bzgl. der Faktoren τ(x, ω, s) angewandt auf obige
Gleichung
d
L(x − sω, ω) − σt (x − sω)L(x − sω, ω) (3)
τ(x, ω, s)
ds
= −τ(x, ω, s)Lv (x − sω, ω)
Fordert man
Mario Hlawitschka
Digitale Signalverarbeitung
81
Die Transformation IV
−τ(x, ω, s)σt (x − sω) =
dτ(x, ω, s)
ds
können wir die linke Seite umschreiben zu
d
[τ(x, ω, s)L(x − sω, ω)] = −τ(x, ω, s)Lv (x − sω, ω)
ds
Um nach τ lösen zu können, setzen wir die Randbedingungen
τ(x, ω, 0) = 1 und τ(x, ω, s) 6= 0 für s > 0 und schreiben
σt (x − sω) = −
1
dτ(x, ω, s)
τ(x, ω, s)
ds
Durch Integration beider Seiten von 0 bis d erhält man
Mario Hlawitschka
Digitale Signalverarbeitung
82
Die Transformation V
Z d
0
= −
σt (x − sω)ds
Z d
0
1
dτ(x, ω, s)
ds
τ(x, ω, s)
ds
= − ln τ(x, ω, s)|d0 = − ln τ(x, ω, d)
weshalb der Faktor in der Integration wie folgt aussieht:
τ(x, ω, d) = e −
Rd
0
τt (x−sω)ds
= e −d0 (x,ω,d)
Man nennt diese Funktion Pfadabsorptionsfunktion entlang des
Strahls (x, ω) von 0 bis d, wobei
d0 (x, ω, d) :=
Z d
0
σt (x − sω)ds
die optische Tiefe ist. Im Falle eines homogenen Mediums, wo σt
Mario Hlawitschka
Digitale Signalverarbeitung
83
Die Transformation VI
konstant ist, wird der Ausdruck unabhängig von Ursprung und
Richtung und vereinfacht sich zu d0 (d) = σt d und somit wird
τ(d) = e −σt d . Wir halten fest, dass τ(0) = 1 und τ(s) 6= 0 wie
gefordert. Integriert man beide Seiten aus Gleichung 3, so erhält
Mario Hlawitschka
Digitale Signalverarbeitung
84
Die Transformation VII
man
Z d
d
0
ds
Z d
[τ(x, ω, s)L(x − sω, ω)]ds
=
−
=
τ(x, ω, s)L(x − sω, ω)|d0
0
τ(x, ω, s)Lv (x − sω, ω)ds
⇔ τ(x, ω, d)L(x − dω, ω) − τ(x, ω, 0)L(x, ω)
=
−
Z d
0
τ(x, ω, s)Lv (x − sω, ω)ds
⇔ L(x, ω) = τ(x, ω, d)L(x − dω, ω)
+
Z d
0
τ(x, ω, s)Lv (x − sω, ω)ds.
Obwohl dieser Ausdruck für alle Werte von d gültig ist, ist er nur
Mario Hlawitschka
Digitale Signalverarbeitung
85
Die Transformation VIII
dann sinnvoll, wenn wir Randbedingungen mit einbeziehen können.
Wir setzen also d = d(x, −ω) als die Distanz von x in Richtung
−ω zum Rand. Dann wird die stationäre Transfergleichung in
Integralform
L(x, ω) = τ(x, ω, d(x, −ω))Lb (h(x, −ω), ω)
+
Z d(x,−ω)
0
τ(x, ω, s)Lv (x − sω, ω)ds
Hierbei sind Lb und Lv Funktionen von L, also ist diese Gleichung
keine geschlossene Form von L. Es ist lediglich eine alternative
Schreibweise der Integro-Differenzialgleichung mit
Berücksichtigung der Randbedingungen. Ähnliche Herleitungen
sind möglich, wenn L zeitabhängig ist oder die Frequenzen
gekoppelt sind. In beiden Fällen wird die Herleitung nach einer
Laplacetransformation ähnlich erfolgen wie oben und man erhält
eine Gleichung ähnlicher Struktur.
Mario Hlawitschka
Digitale Signalverarbeitung
86
Vakuumtransport I
Wählen wir nun V als Vakuum. Daraus ergibt sich, dass weder
Volumenstreuung noch Volumenemission stattfindet, also
σs ≡ 0, σa ≡ 0, Lev ≡ 0.
Dann ist auch σt = 0 und Gleichung 1 wird zu
d
L(x, ω) = 0 für x ∈ V
dω
was bedeutet, dass die Radianz entlang eines Strahls (bzw. jeder
Geraden) im Volumen konstant ist. Folglich ist die Radianz an
jedem inneren Punkt bestimmt durch die Radianz des Randes.
Mit Hilfe der Strahlschnittfunktion kann man diese bestimmen
L(x, ω) = L(h(x, −ω), ω) für x ∈ V , ω ∈ S 2
was die Gleichung wiederum auf die Randbedingungen vereinfacht.
Diese Gleichung heißt dann Radianzgleichung oder Rendering
Mario Hlawitschka
Digitale Signalverarbeitung
87
Vakuumtransport II
Equation:
L(x, ωr ) = Le (x, ωr ) +
Z
2 (x)
S−
fr (ωi , x, ωr )L(h(x, −ωi ), ωi ) cos θi dωi
(4)
wobei jetzt x ∈ ∂ V liegt. In Operatorschreibweise ausgedrückt
ergibt sich
L = Le + Tfr L
wobei Tfr der Integraloperator ist und eine Kurzschreibweise für
das Integral der reflektierten Radianz darstellt. Offensichtlich ist
Mario Hlawitschka
Digitale Signalverarbeitung
88
Vakuumtransport III
die Darstellung der Radianz eine Fredholm’sche
Integralgleichung zweiter Art. Aufgrund der physikalischen
Beschränkungen, dass immer ein Teil der Strahlung in Wärme
umgewandelt wird, gilt immer kTfr k < 1 woraus folgt, dass die
Neumannreihe konvergiert. Beschränken wir weiter
fr = fr (x) =
1
ρ(x)
π
auf den rein isotropen Fall, das heißt rein diffuse Reflexion, erhält
man die kontinuierliche Radiosity Equation
L(x) = Le (x) +
ρ(x)
π
Z
2 (x)
S−
L(h(x, −ωi )) cos θi dωi
(5)
wobei aufgrund der Isotropie des Reflexionskerns die Radianz
ebenfalls isotrop wird.
Mario Hlawitschka
Digitale Signalverarbeitung
89
Die Frage nach einer geeigneten Parametrisierung ist wichtig für
die Entwicklung von Algorithmen. Wir führen deshalb im
Folgenden einen Wechsel der Domäne der Radianz von ∂ V × S 2
nach ∂ V × ∂ V durch, die im folgenden nützlich sein wird.
x −y
y −x
x −y
= L h y,
,
L(x, y ) := L x,
|x − y |
|x − y | |y − x|
Hier bedeutet L(x, y ) die Radianz, die den Punkt x in Richtung
Punkt y verlässt, aber nicht notwendigerweise auch dort ankommt.
Beim Wechsel des Integrationsgebietes muss das Maß auch
Mario Hlawitschka
Digitale Signalverarbeitung
90
entsprechend der Substitutionsregel angepasst werden. Unter
Verwendung des Raumwinkels erhalten wir
dωx =
cos θy
1
dyp =
dy
2
r
r2
wobei x und y die Referenzpunkte des lokalen Koordinatensystems
sind und damit der projizierte Raumwinkel
cos θx dωx =
cos θx cos θy
dy
r2
ergibt. Die Distanz r = |y − x| ist der Abstand der Punkte. Der
Mario Hlawitschka
Digitale Signalverarbeitung
91
symmetrische geometrische Term G ist nun definiert als
G (x, y ) := V (x, y )
cos θx cos θy
= G (y , x)
kx − y k2
und wird als Dichte im Integrator verwendet um
(Tfr , L)(y , ωr ) =
mit ωr =
z−y
ky −zk
Z
2
S−
L(h(y , ωi ), −ωi )fr (ωi , y , ωr ) cos θi dωi
durch
(Tf L)(y , x) =
Z
∂V
L(x, y )fr (x, y , z)G (x, y )dx
zu ersetzen. wobei fr (x, y , z) = fr
Mario Hlawitschka
y −x
kx−y k , y , ωr
Digitale Signalverarbeitung
ist. Die
92
Kosinusfaktoren sorgen dafür, dass x und y jeweils senkrecht auf
der Strecke von x nach y betrachtet werden. Die Division durch
das Quadrat des Abstands sorgt für die r12 -Regel des
Strahlungstransports mit zunehmendem Abstand10 . Die
Sichtbarkeit wird von V überprüft. Damit wird aus der
Mario Hlawitschka
Digitale Signalverarbeitung
93
Radianzgleichung
L(y , z) = Le (y , z) +
Z
∂V
L(x, y )fr (x, y , z)G (x, y )dx
und im Spezialfall der Radiosity erhalten wir
ρ(y )
L(y ) = Le (y ) +
π
Z
∂V
L(x)G (x, y )dx
für das Integrationsgebiet, dem Rand von V = ∂ V .
Aufgrund des r12 -Faktors wird diese Formulierung schwach
singulär genannt, da sie numerische Probleme bei kleinen
Abständen hervorrufen kann, die in der anderen Formulierung über
den Raumwinkel nicht auftreten.
10 Die Strahlung verhält sich proportional zur Kugeloberfläche in der
Entfernung
Mario Hlawitschka
Digitale Signalverarbeitung
94
Wenn Volumenrendering durchgeführt wird, ist die gängige
Annahme, keine Volumenstreuung zu berücksichtigen, also σs ≡ 0.
Dies vereinfacht Lv zu
Lv (x, ω) ≡ Lev (x, ω)
und die dominierende Gleichung in Integralform wird
L(x, ω) = τ(x, ω, d(x, −ω))Lb (h(x, −ω), ω)
+
Z d(x,−ω)
0
τ(x, ω, s)Lev (x − sω, ω)ds
Wenn wir weiter annehmen, dass der Rand nicht reflektierend ist,
d.h., fr ≡ 0, dann gilt
Lb (x, ω) ≡ Le (x, ω)
und somit
L(x, ω) = τ(x, ω, d(x, −ω))Leb (h(x, −ω), ω)
+
Z d(x,−ω)
0
τ(x, ω, s)Lev (x − sω, ω)ds
Weiterhin kann man Oberflächen vollständig vernachlässigen und
Digitale Signalverarbeitung
95
Emissions-Absorptions-Modell
Marioerhält
Hlawitschka
das
Difussionsapproximation I
Basierend auf der Beobachtung, dass in stark streuenden Medien
die Lichtverteilung isotrop wird, ist es ausreichend, die Radianz L
über die ersten beiden Momente anzunähern. In
Kugelflächenfunktionen11 ausgedrückt benötigt man die ersten
Basisfunktionen
1
L(x, ω)
=
l
∑ ∑
cl,m (x)Yl,m (ω)
l=0 m=−l
r
r
r
1
3
3
+ c1,−1 (x)
ωy + c1,0 (x)
ωz + c1,1 (x)
= c0,0 (x)
4π
4π
4π
r
r
1
3
= c0,0 (x)
+
hc1,−1 (x), c1,0 (x), c1,1 (x)), ωi
4π
4π
=: L0 (x) + E (x) · ω
r
wobei L0 die Irradianz (manchmal auch fluence“ genannt) und E
”
die Vektor-Irradianz ist. Die Projektion von L in die Basis der
Mario Hlawitschka
Digitale Signalverarbeitung
96
Difussionsapproximation II
Kugelflächenfunktionen ergibt
L0 (x) =
E (x) =
Z
1
L(x, ω)dω
4π S 2
Z
3
L(x, ω)ωdω.
4π S 2
Projiziert man auch die Volumenlichtquellen Lev (x, ω), so erhält
man die Annäherung
Lev ≈ q0 (x) + q1 (x) · ω.
Da die Kugelflächenfunktionen die Legendre-Polynome12
verwenden, nennt man die Annäherung über die ersten beiden
Polynome auch P1 -Approximation.
11 http://de.wikipedia.org/wiki/Kugelfl\protect\unhbox\voidb@x\
bgroup\U@D1ex{\setbox\z@\hbox{\char127}\[email protected]\advance\
dimen@\ht\z@\fontdimen5\font\dimen@}\accent127\fontdimen5\font\U@
Da\egroupchenfunktionen
12 http://de.wikipedia.org/wiki/Legendre-Polynom
Mario Hlawitschka
Digitale Signalverarbeitung
98