Petr Olšák: Matematika aplikovaná v počítačové grafice
15. května 2009
Matematická kolokvia katedry matematiky FEL ČVUT v Praze
Abstrakt: Přednáška obsahuje stručný úvod do jedné oblasti počítačovégrafiky: do problematiky výpočtů šíření světla v 3D scéně. Formulacematematického modelu a metody řešení: konečné prvky, Monte Carlometody.
1
Formulace problému: fotit bez foťáku, filmovat bez kamery
↓Kompromisy, zjednodušení: paprsek v celé své délce svítí stejně, . . .
↓Matematický model: diferenciální, resp. integrální rovnice
↓Kompromisy, zjednodušení: diskretizace, aproximace plochy, linearizace
↓Soustava lin. rovnic: obří matice soustavy, nepřesné hodnoty
↓Kompromisy, zjednodušení: iterační numerické metody
↓Výsledek (Je ještě k něčemu? Dočkáme se výsledku? poznáme, že byiterační metoda mohla skončit?)
2
Veličiny používané v matematickém modelu
Φin(S, U) . . . světelný tok (flux) [W] . . . světelná energie dopadajícína plochu S ze směru U za jednotku času (resp. v okamžiku t)
Φout(S, U) . . . světelný tok (flux) [W] . . . světelná energie opouštějícíplochu S do směru U za jednotku času (resp. v okamžiku t)
E(x,U) = lim|S|→0,x∈S
Φin(S, U)|S|
. . . irradiance
. . . ozáření bodu na ploše ze směru U , tj. hustota toku Φin na ploše.
Analogicky: B(x,U) = lim|S|→0,x∈S
Φout(S, U)|S|
. . . radiozita
. . . vyzařování bodu do směru U , tj. hustota toku Φout na ploše.
častěji: E(x) = E(x,Ω), B(x) = B(x,Ω).
3
Abstrakce paprsku (radiance)
Lout(x, ω) = lim|S|→0,x∈S
1cos θ · |S|
(lim
|U |→0,ω∈U
Φout(S, U)|U |
)Analogicky Lin(x, ω).
cos θ zaručuje nezávislost radiance L na normále plochy S.
Radiance [Wm−2 st−1] se dá měřit expozimetrem, je úměrná reakci fo-tocitlivého materiálu (filmu), respektive čidla v digitálním fotoaparátu.
Na celé délce paprsku změříme stejné L.
Integrální formulace předchozího:
Φin/out(S, U) =∫
S
∫U
Lin/out(x, ω) cos θ dω ds
4
Odrazivost materiálu
Uvažujme x ∈ S.∫O
p(ω′, x, ω) dω . . . pravděpodobost, že foton, který dopadl ze směru ω′
na bod x se odrazí do směru O, tj. p(ω′, x, ·) je hustota náhodné veličinypopisující tuto prvděpodobnost.
E(x,O′) =∫
O′Lin(x, ω′) cos θ′ dω′ . . . irradiance bodu x ze směru O′.
Jak se tato energie odrazí?∫O
∫O′
Lin(x, ω′) cos θ′p(ω′, x, ω) dω′ dω
. . . množství E(x,O′), která se odrazila do směru O.
5
Odrazivost materiálu
Φout(S, O) obsahující odraženou energii ze všech směrů je rovna dvěmavýrazům:
Φout(S, O) =∫
S
∫O
∫Ω′
Lin(x, ω′) cos θ′p(ω′, x, ω) dω′ dω ds,
Φout(S, O) =∫
S
∫O
Lout(x, ω) cos θ dω ds
Po „zderivováníÿ Φout(S, O) podle S a podle O dostáváme
Lout(x, ω) cos θ =∫Ω′
Lin(x, ω′) cos θ′p(ω′, x, ω) dω′
6
Odrazivost materiálu
Lokální odrazová rovnice:
Lout(x, ω) =∫Ω′
p(ω′, x, ω)cos θ
Lin(x, ω′) cos θ′ dω′
Funkce fr(ω′, x, ω) =
p(ω′, x, ω)cos θ
se nazývá BRDF (Bidirectional re-
flection distribution function), popisuje odrazové vlastnosti materiálu.
Helmholtz reciprocity: fr(ω′, x, ω) = fr(ω, x, ω′)
Zákon zachování energie:
a(x, ω′) =∫Ω
p(ω′, x, ω) dω =∫Ω
fr(ω′, x, ω) cos θ dω ≤ 1
. . . Funkce a se nazývá odrazivost (albedo) materiálu.
7
Zobrazovací rovnice [J. T. Kajiya, 1986]
K odraženému světlu přičteme vlastní vyzářené světlo a dostaneme lo-kální zobrazovací rovnici:
Lout(x, ω) = Le(x, ω) +∫Ω′
fr(ω′, x, ω)Lin(x, ω′) cos θ′ dω′
Substitucí ω′ → y = ray(x, ω′), dω′ → cos θ′′
‖x− y‖2ds a s využitím toho,
že Lin(x, ω′) = Lout(y,−ω′), dostáváme:
Lout(x, ω) = Le(x, ω)+∫
y∈S
fr(ω′, x, ω)Lout(y,−ω′)V (x, y)
cos θ′ cos θ′′
‖x− y‖2ds,
kde ω′ =y − x
‖x− y‖, V (x, y) =
1 když x vidí y0 jindy
Toto je globální zobrazovací rovnice s neznámou Lout, Fredholmova in-tegrální rovnice, matematický model šíření světla ve scéně.
8
Integrál jako lineární operátor
Při značení TL =∫
S
fr(ω′, ·, ·)L(y,−ω′)V (·, y)cos θ
′ cos θ′′
‖· − y‖2ds,
kde Lout = L, přechází zobrazovací rovnice na
L = Le + TL, tj. (I − T )L = Le
Inverzi operátoru (I − T ) lze zapsat jako Neumannovu řadu:
(I − T )−1 =∞∑
k=0
T k, takže L =∞∑
k=0
T kLe.
Život ale není takto jednoduchý. . .
9
Zjednodušení: difúzní povrch
Předpoklad: povrch a zdroje jsou pouze difúzní, tj. Lout(x, ω) = Lout(x)a fr = %(x)/π, nezávisí na směrech ω, ω′. Zobrazovací rovnice přecházína tvar:
Lout(x) = Le(x) +%(x)π
∫S
Lout(y)V (x, y)cos θ′ cos θ′′
‖x− y‖2ds,
Radiozita B(x) =∫ΩLout(x) cos θ dω = Lout(x)
∫Ωcos θ dω = Lout(x)π
B(x) = Be(x) +%(x)π
∫S
B(y)V (x, y)cos θ′ cos θ′′
‖x− y‖2ds,
B(x) = Be(x) + %(x)∫
S
B(y)G(x, y) ds,
když G(x, y) = V (x, y)cos θ′ cos θ′′
π ‖x− y‖2
10
Diskretizace: radiozitní metoda [C. M. Goral & all, 1984]
Hledá se B po částech konstantní na ploškách Ai (⋃
Ai = S). Funkce B
aproximuje řešení rovnice:
B(x) = Be(x) + %(x)∫
S
B(y)G(x, y) dsy
Uvažuji skalární součin 〈f, g〉 =∫
S
f(x)g(x) ds.
Množina M = ci; ci(x) = 1/|Ai| na Ai, jinde nula tvoří ortonormálníbázi. Bi = 〈B, ci〉 jsou souřadnice B. Aplikuji skalární součin 〈., ci〉 napravou stranu rovnice (kde jsem funkci B(y) nahradil B =
∑j Bjcj a
souřadnice Be značím Ei) a dostávám
Bi = Ei + %i1|Ai|
∫Ai
∑j
Bj
∫Aj
G(x, y) dsy dsx = Ei + %i
∑j
BjFi,j ,
při značení Fi,j =1|Ai|
∫Ai
∫Aj
G(x, y) dsx dsy
11
Radiozitní metoda, form faktory
Fi,j =1|Ai|
∫Ai
∫Aj
G(x, y) dsx dsy je tzv. form factor.
Platí:
Fi,j · energie vyzářená z Ai = energie z Ai, která dopadla na Aj , nebo:
Provedeme půmět Aj na hemisféru kolem Ai a tento průmět promítnemena základnu hemisféry. Pak Fi,j je obsah tohoto dvojího průmětu dělenýobsahem základny hemisféry.
Zřejmě: Fi,j ≥ 0,n∑
j=0
Fi,j ≤ 1.
12
Radiozitní metoda, soustava lineárních rovnic
Bi = Ei + %i
∑j
Fi,jBj .
při značení B = (B1, . . . , Bn)T , E = (E1, . . . , En)
T , lze zapsat problémmaticově jako:
KB = E, kde K = I − diag(%)F =
1− %1F1,1 . . . −%1F1,n−%2F2,1 . . . −%2F2,n−%2F2,1 . . . −%2F2,n
. . . . . . . . .−%nFn,1 . . . 1− %nFn,n
Je dáno E (radiozity zdrojů) a K (obludná matice form faktorů), hledáse B (radiozity všech plošek).
13
Radiozitní metoda, numerické metody
Gaussova eliminace nepoužitelná
Iterativní metody: Jacobiho metoda, Gauss-Seidel, Southwell, . . .
Další možnosti vylepšení (progresivní radiozita, hiearchická radiosita)[M. F. Cohen, 1993]
Problémy:
Dělení na plošky je někde zbytečně jemné, jinde hrubé, způsobuje ne-žádoucí artefakty. Metody: za chodu provádět zjemňování dělení podlepotřeby.
Singularita v geometrickém členu.
Náročnost výpočtu form faktorů.
Velikost matice K.
14
Monte Carlo, Integral estimator [J. von Neumann, S. Ulam, 1944]
Nechť X je náhodná veličina s rovnoměrným rozdělením na V ⊂ Rs. Pak
E(f(X)) =∫
V
f(x)p(x) dx =∫
V
f(x)1|V |dx = I · 1
|V |
Podle zákona velkých čísel, jsou-li x1, x2, . . . , xN vzorky náhodné veličinyX a N je dostatečně velké, pak
E(f(X)).=1N
N∑i=1
f(xi)
Takže I =V
N
N∑i=1
f(xi) odhaduje integrál∫
V
f(x) dx.
15
Monte Carlo, odhad chyby
Hledanou hodnotu integrálu∫
V
f(x) dx označím C. Náhodné veličiny
f(X1), f(X2), . . . , f(XN ) budiž nezávislé. Mají stejnou střední hodnotuE(f(Xi)) = C/|V |. Označím σ2 jejich společný rozptyl. Pak podle cent-rální limitní věty je
1N
∑Ni=1 f(Xi)− C
|V |
σ/√
N→ N(0, 1).
Z tabulky kvantilů N(0, 1) plyne například, že pro velká N je číslo 0,997přibližně rovno pravděpodobnosti:
P
(∣∣∣∣∣1N
∑Ni=1 f(Xi)− C
|V |
σ/√
N
∣∣∣∣∣ < 3)= P
(∣∣∣∣∣ |V |N
N∑i=1
f(Xi)− C
∣∣∣∣∣ < 3σ|V |√N
)
takže na hladině významnosti 0,997 je odhad chyby uvedený v kulaté zá-vorce výše zaručen. Odhad chyby obsahuje
√N ve jmenovateli nezávisle
na dimenzi V .
16
Kasické integrační metody, dimensional explosion
Uvažuji f(x) = kx1 a jmu se to integrovati na 〈0, 1〉S obdélníkovou me-todou pomocí N vzorků. V jednom směru potřebuji n = S
√N vzorků.
Odhad chyby:
12
k
n
1nS
nS =k/2S√
N
Odhad chyby ve jmenovateli obsahuje S√
N , což je pro s > 2 horší nežodhad chyby Monte Carlo metody.
Pro jiné metody (lichoběžníková atd.) stačí volit jinak funkci f (stáles variací k) a odhad dopadne stejně.
17
Integrační metody obecně, Quasi Monte Carlo
Integral estimator I =N∑
i=1
f(xi)w(xi)
Odhad chyby (Koksma-Hlawka Inequality):∣∣∣∣∫V
f(x) dx− I
∣∣∣∣ ≤ |V | · v(f) ·D, kde
v(f) je variace f
D je diskrepance bodů xi v objemu V .
To vede k nápadu dělat pseudonáhodný výběr xi: částečně náhodný a čás-tečně řízený k vylepšení diskrepance. To je podstata metod quasi MonteCarlo.
18
Modifikace Monte Carlo: Importance sampling
Nechť náhodná veličina X má rozdělení s hustotou p na V . Pak středníhodnota náhodné veličiny f(X)/p(X) je:
E
(f(X)p(X)
)=∫
V
f(x)p(x)
p(x) dx =∫
V
f(x) dx = C (hledaná hodnota)
E
(f(X)p(X)
).=1N
N∑i=1
f(xi)p(xi)
(estimátor)
Rozptyl estimátoru je E
((f(X)p(X)
− C
)2)=∫
V
(f(x)p(x)
− C
)2p(x) dx
Kdyby p = f/C, máme nulový rozptyl. Když se p bude funkci f/C„podobatÿ, budeme mít „malýÿ rozptyl a tím i odhad chyby bude lepší.Samozřejmý problém: neznáme C. Stačí ale najít funkci q, kterou umímeintegrovat přes V (integrál označím Cq). Dále chceme, aby q se až nanásobek podobala f a volíme p = q/Cq.
19
Metody sledování paprsku
Scénou se vede paprsek od kamery nebo od zdroje světla a nechá se (vzávislosti na variantě metody) náhodně odrazit od překážky a trasuje sedál. Vyhodnocuje se sbíraná/vysílaná energie podél paprsku. Náhodnépokračování paprsku v místě odrazu odpovídá výpočtu integrálu v od-razové rovnici metodou Monte Carlo.
20
Metody sledování paprsku, výhody, nevýhody
Výhody:
Popis povrchu možný různými grafickými primitivy
Není nutné dodatečné zjemňování
Je možné pracovat s libovolnou BRDF
Malé požadavky na paměť
Řešení není vychýleno (unbiased)
Nevýhody:
Pomalá konvergence
Šum
Pohledově závislý výpočet
21
Metody sledování paprsku, varianty
Od kamery (gathering methods):
Ray tracing: klasická metoda, rekurze s odrazem/lomem jen na zrca-dlech/místech lomu, počítá jen přímé osvětlení. [T. Whitted, 1980]
Path tracing: rekurze do stanovené hloubky, resp. ukončení paprsku me-todou „ruská ruletaÿ. [J. T. Kajiya, 1986]
Od světla (shooting metods):
Photon tracing: rekurze odrazem/lomem jen na zrcadlech/místech lomu.(nepoužívá se)
Light tracing: rekurze do stanovené hloubky, resp. ukončení paprskupodle zbytkové energie.
22
Kombinace předchozích metod
Bidirectional path tracing [E. P. Lafortune, I. D. Willems, 1993]
Multiple importance path tracing [E. Veach, 1995]
Photon mapping [H. W. Jensen, 1996]
Další metody:
Precomputed radiance transfer [P. P. Sloan, 2002]
Direct-to-indirect transfer [M. Hašan, 2006]
. . .
23
