Część I [1] Część II [2] Część III [3] Część IV
Czym będziemy się zajmować?
Zrezygnujemy teraz z obsługi transmisji i odbicia lustrzanego. Skupimy się wyłącznie na idealnie rozproszonym odbijaniu światła, a nasz program będzie potrafił obliczać odbicie rozproszone światła pochodzącego nie tylko wprost ze źródła światła, ale również światła odbijanego od różnych obiektów. Dzięki temu otrzymamy nowe, ciekawe efekty.
Zmodyfikujemy nasz program (który obsługiwał ścieżki typów E(S|T)*DL) tak, aby obsługiwał wyłącznie ścieżki typów E(D*)L. Poniższa ilustracja prezentuje efekty tej zmiany.
Po lewej stronie fragment obrazu wygenerowanego programem, który pisaliśmy do tej pory. Posługiwaliśmy się w nim tzw. światłem punktowym, przez co mieliśmy efekt cieni o ostrych krawędziach (działo się tak dlatego, że światło punktowe jest albo zasłonięte albo widoczne w całości, stąd gwałtowne przeskoki z pełnego oświetlenia do pełnego cienia). Po prawej stronie fragment obrazu wygenerowanego zmodyfikowaną metodą, w którym rezygnujemy ze światła punktowego na rzecz świateł powierzchinowych. Prowadzi to do miękkich cieni gładko przechodzących od pełnego naświetlenia do pełnego cienia (dzieje się tak dlatego, że światło powierzchniowe może być zasłonięte tylko w części). Ponadto nie musimy już źródłom światła nadawać właściwości ambient (światło tła), gdyż to właśnie ścieżkom typu ED(D+)L odpowiadają za poprawną symulację tego zjawiska. Dodatkowo teraz światło tła będzie bardziej dokładne (zróżnicowane w zleżności od kąta padania) oraz otrzymamy efekt tzw. krwawienia kolorów, w którym powierzchnia barwi się kolorem obiektu sąsiadującego. Modyfikacja, którą będziemy wprowadzać do naszego kodu stworzy algorytm o nazwie śledzenie ścieżek (ang. path tracing lub też monte carlo path tracing).
Zmiany jakich dokonamy w celu uzyskania nowej funkcjonalności będą bardzo bardzo niewielkie - dodamy jedynie dwie krótkie procedury i zmodyfikujemy nieznacznie procedurę .
Losowanie kierunków odbijania
Teraz przyjrzyjmy się jak wygląda scenariusz naszego obliczania oświetlenia. Startujemy z punktu (bedącego kamerą) i trafiamy promieniem w punkt . Losujemy kierunek z półsfery wokół punktu , a następnie śledzimy promień w tym wylosowanym kierunku, trafiamy w punkt , dla którego musimy zrobić to samo co dla punktu , itd. Dla każdego piksela obliczamy wiele w ten sposób losowanych ścieżek światła i obliczamy średnią arytmetyczną uzyskanych wyników.
Kiedy kończy się opisany powyżej proces?
1) Kiedy trafimy w świecącą powierzchnię.
Zakładamy, że powierzchnia emitująca światła nie jest odbijająca.
Powierzchnie świecące na ogół są na tyle jasne, że dodanie właściwości
odbijania nie spowoduje widocznych efektów.
2) Kiedy promień nie trafia w żaden obiekt na scenie.
Nie trafiliśmy w źródło światła, ani w obiekt, od którego można by się odbić,
więc nasza ścieżka niesie zerowe natężenie światła.
3) Kiedy uznamy, że głębokość rekursji jest już zbyt duża.
Teoretycznie może się zdarzyć, że bardzo długo (lub w nawet nieskończoność) nasz promień będzie się
odbijał pomiędzy obiektami sceny. Zauważmy jednak, że im więcej odbić (czyli im większa długość przebytej ścieżki),
tym mniejsze natężenie światła jakie potencjalnie może płynąć po takiej ścieżce. Największy wpływ na jakość
obrazu mają na ogół ścieżki o długościach , oraz , stąd śmiało możemy przyjąć, że całkowicie
ignorujemy ścieżki o długościach np. większych od . Wprowadzimy w ten sposób pewien błąd oszacowania oświetlenia,
ale będzie on na tyle mały, że nie spowoduje widocznych, niepożądanych efektów.
Wyniki renderingu
Słaba zbieżność w oszacowaniu oświetlenia będzie się objawiać ziarnistym szumem. Przyjrzyjmy się kilku przykładom pokazującym jak poprawia się jakość obrazu wraz z rosnącą ilością próbek (losowanych ścieżek światła).
Tutaj obliczenie 100 000 ścieżek na piksel w rozdzielczości 600 x 600 trwało kilkanaście godzin na procesorze o taktowaniu 1.5 GHz. Istnieją rozszerzenia tej metody, które potrafią uzyskać pożądaną jakość obrazu wykorzystujące zaledwie kilkadziesiąt lub kilkaset ścieżek na piksel.
Alternatywna scena
Nasz algorytm daje jeszcze jedną ciekawą możliwość, oprócz miękkich cieni i krwawienia kolorów możemy w łatwy sposób zasymulować oswietlenie pochodzące z atmosfery. Wystarczy umieścić w scenie świecącą sferę o bardzo dużym promieniu i umieścić w jej środku wszystkie inne obiekty sceny.
Co nowego w kodzie?
Na początek potrzebujemy procedury losującej kierunki z półsfery. Procedura daje w wyniku losowy kierunek z półsfery. Nie musisz teraz rozumieć jak ona działa.
1 2 3 4 5 6 7 8 | // losowy wekotr na półsferze o promieniu 1 // rozkład cosinusa inline vec3 sampleHemiSphere() { float Rd = uRand(); float Ad = uRand(); float c = sqrt(1.0 - Rd); return vec3(c * cos(PiMul2*Ad), sqrt(Rd), c * sin(PiMul2*Ad), 1.0f); } |
Kolejną rzeczą jaką będziemy potrzebować jest procedura, która wygeneruje nam kierunek z półsfery zaczepionej w dowolnym punkcie sceny o dowolnie zorientowanym wektorze normalnym do powierzchni. Przyda się do tego procedura , której używalismy w mapowaniu wektorów normalnych.
1 2 3 4 5 6 7 8 | inline vec3 uChangeSpace(vec3::i V, vec3::i Ex, vec3::i Ey, vec3::i Ez) { return vec3( V.x*Ex.x + V.y*Ey.x + V.z*Ez.x, V.x*Ex.y + V.y*Ey.y + V.z*Ez.y, V.x*Ex.z + V.y*Ey.z + V.z*Ez.z, 1.0f ); } |
Spójrzmy na poniższą ilutrację.
Czerwonym promieniem trafiamy w obiekt ( - punkt na jego powierzchni). Z obiektu możemy odczytać wektor normalny do powierzchni. Przy pomocy iloczynów wektorowych obliczamy wektory oraz . Losujemy wektor na półsferze przy pomocy procedury , a następnie mapujemy go na przestrzeń rozpiętą na wektorach , , , otrzymując interesujący nas wektor kierunku odbicia z rozproszeniem (na ilustracji wektor koloru żółtego). Cały scenariusz realizuje poniższa procedura .
1 2 3 4 5 6 7 | // odbicie z rozproszeniem inline ray scatter(ray::i r, vec3::i p, vec3::i N) { vec3 Ex = uCrossProd(r.abackDir, N); vec3 Ez = uCrossProd(Ex, N); vec3 hsV = sampleHemiSphere(); return ray(p, uChangeSpace(hsV, Ex, N, Ez)); } |
W procedurze dokonujemy niewielkiej zmiany. Wcześniej obliczaliśmy trzy składowe - oświetlenie bezpośrednie z punktowego źródła światła, odbicie zwierciadlane i transmisję z refrakcją, a następnie sumowaliśmy ich wyniki. Teraz obsługujemy tylko emisję i odbicie rozproszone.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | // r - śledzony promień // s - scena, w której śledzony jest promień // d - głębokość rekursji color traceRay(ray::i r, scene& s, nat d) { if (d == 0) return black; // koniec rekursji hit h = nearestHit(r, s); // informacje o najbliższym przecięciu if (h.t == undefined) return black; // brak przecięcia material& m = *h.mat; // materiał trafionego obiektu color texMapped = notNull(m.drMap) ? (*m.drMap)(h.uv.u, h.uv.v) : white; if (m.hasNormalMap) // perturbacja wektora normalnego zgodnie z teksturą h.N = uChangeSpace((*m.normalMap)(h.uv.u, h.uv.v), h.Ex, h.N, h.Ez); if (m.isDiffuse) // odbicie rozproszone return m.dr * texMapped * traceRay(scatter(r, h.ip, h.N), s, d-1); else // emisja return m.e; } |
Materiały do artykułu
Zaprogramowaliśmy już wszystkie poprawki w naszym kodzie, więc możesz teraz ściąganąć pełny kod napisanego programu.
Pod niewielkimi zmianami w kodzie kryje się trochę matematyki. Dowiemy się teraz dlaczego to, co zrobiliśmy, działa.
Całkowanie
Niespodzianka, a teraz nauczymy się całkować. Interesować nas będzie geometryczna interpretacja całki, która mówi, że (czyt. całka z funkcji na odcinku ) jest równa polu pod wykresem funkcji , jeśli przyjmuje ona wartości wyłącznie nieujemne. Przyjrzyjmy się kilku przykładom.
Przykład 1
Niech . Zauważmy, że na odcinku pole pod wykresem funkcji jest równe polu prostokąta, którego jeden bok ma długość , a drugi ma długość . A więc:
Przykład 2
Niech . Zauważmy, że na odcinku pole pod wykresem funkcji jest równe połowie pola kwadratu o boku . A więc:
Przykład 3
Niech teraz funkcja , będzie taka jak na poniższej ilustracji. W przypadku takiej funkcji pojawia się problem ("gołym okiem" niezbyt widać ile wynosi pole pod wykresem takiej funkcji). Teraz poradzimy sobie całkowaniem stochastycznym.
Całkowanie stochastyczne
Całkowanie stochastyczne polega na oszacowaniu (przybliżeniu) wartości całki w sposób losowy. Konstruujemy coś, co będzie przybliżało wartość interesującej nas całki (będzięmy to nazywać estymatorem).
Estymator zdefiniowany jest następująco:
Na początek oswoimy się z pojęciem rozkładu. Intuicyjnie rozkład warunkuje to, w których miejscach dziedziny funkcji punkty będą losowane z większym prawdopodobieństwem, a w których z mniejszym (co skutkuje tym, że jak wylosujemy przykładowo 100 punktów, to w jednych obszarach dziedziny mogą powstawać gęstsze skupiska punktów, a w innych punkty będą występować sporadycznie).
Na powyższej ilustracji widać przykład dwóch rozkładów. Wylosowano 16 punktów na odcinku. W rozkładzie jednostajnym punkty rozkładają się na odcinku równomiernie (tj. nie mają wyraźnej tendencji do grupowania się tylko w pewnych miejscach). Natomiast w rozkładzie niejednostajnym punkty mogą sie grupować tak jak na ilustracji na końcach odcinka oraz rzadko występować w okolicy jego środka.
Skupimy sie teraz na rozkładzie jednostajnym. Z rozkładem związane jest pojęcie gęstości. Detale teoretyczne takie jak wyprowadzanie gęstości dla konkretnego rozkładu oraz czym dokładnie jest gęstość wybiegają nieco poza poziom trudności tego artykułu, więc funkcje gęstości dla rozkładów bedą podawane jako fakty, w które należy uwierzyć bez dowodu. Dla rozkładu jednostajnego funkcja gęstości jest równa odwrotności długości odcinka, bądź polu lub objętości zbioru, z którego losujemy punkty. Przykładowo dla odcinka długości 10 funkcja gęstości to , a dla prostokąta o polu 25 funkcja gęstości to .
Przekonajmy się, że to wszystko rzeczywiście działa. Wrócmy do przykładu nr 1 całkowania. Mamy funkcję na pewnym przedziale . Z przykładu nr 1 pamiętamy, że:
Sprawdźmy teraz co się dzie w przypadku przykładu nr 2 całkowania. Mamy funkcję na przedziale . Pamiętamy, że:
Po co właściwie nam to całkowanie?
Wyobraźmy sobie, że funkcja nie jest teraz określona na zbiorze, który jest odcinkiem, tylko na zbiorze kierunków rozpinającym półsferę (tak jak na rysunku poniżej). Będziemy teraz chcieli obliczyć całkę z funkcji określonej na takim zbiorze kierunków, czyli objętość pod wykresem tej funkcji.
Teraz zastanówmy się po co nam taka funkcja, która określona jest na takim półsferycznym zbiorze kierunków. Chcemy symulować krwawienie kolorów oraz miękkie cienie. Będziemy to robić poprzez symulację idealnie rozproszonego odbicia światła. W sytuacji (A) na rysunku poniżej widać przykład idealnie rozproszonego odbicia. Wiązka światła trafia w punkt , a następnie rozprasza się równomiernie we wszystkie kierunki po półsferze umiejscowionej w punkcie . W sytuacji (B) z punktu docieramy do punktu i chcemy obliczyć natężenie światła jakie dociera z punktu do punktu (wniosek z sytuacji (A) jest taki, że nie ma znaczenia jaki jest kierunek od punktu do ). Aby to zrobić, musimy najpierw obliczyć sumaryczne natężenie światła jakie dociera do punktu ze wszystkich kierunków po półsferze wokół punktu (na rysunku jest to ten falisty obszar ponad brzegiem półsfery), dokonując tego własnie przez całkowanie stochastyczne.
Spójrzmy jeszcze raz jak wygląda scenariusz obliczania oświetlenia. Startujemy z punktu (będącego kamerą) i trafiamy promieniem w punkt . Chcemy obliczyć natężenie światła, które dociera z punktu do , a więc musimy oszacować poprzez całkowanie stochastyczne natężenie światła jakie dociera do punktu . Losujemy więc z pewnym rozkładem kierunek z półsfery wokół punktu , a następnie dzielimy natężenie światła przychodzące po wylosowanym kierunku przez gęstość rozkładu z jakim losowaliśmy. Aby obliczyć natężenie światła które przychodzi do punktu z wylosowanego kierunku śledzimy promień, trafiamy w punkt , dla którego musimy zrobić to samo co dla punktu .
Całkowanie w naszym kodzie
Procedura daje w wyniku losowy kierunek z półsfery losując kierunki z rozkładem cosinusa. Poniżej ilustracja pokazująca różnicę pomiędzy rozkładem jednostajnym i cosinusa. W rozkładzie jednostajnym kierunki podczas losowania rozkładają się równomiernie po półsferze, natomiast w rozkładzie cosinusa kierunki bedą występować gęściej wokół wektora normalnego do powierzchni, na której leży punkt, w którym obliczamy oświetlenie.
We wzorze na gęstość w rozkładzie cosinusa to kąt pomiędzy kierunkiem na półsferze a wektorem normalnym do powierzchni obiektu.
1 2 3 4 5 6 7 8 | // losowy wekotr na półsferze o promieniu 1 // rozkład cosinusa inline vec3 sampleHemiSphere() { float Rd = uRand(); float Ad = uRand(); float c = sqrt(1.0 - Rd); return vec3(c * cos(PiMul2*Ad), sqrt(Rd), c * sin(PiMul2*Ad), 1.0f); } |
Wróćmy teraz do procedury .
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | // r - śledzony promień // s - scena, w której śledzony jest promień // d - głębokość rekursji color traceRay(ray::i r, scene& s, nat d) { if (d == 0) return black; // koniec rekursji hit h = nearestHit(r, s); // informacje o najbliższym przecięciu if (h.t == undefined) return black; // brak przecięcia material& m = *h.mat; // materiał trafionego obiektu color texMapped = notNull(m.drMap) ? (*m.drMap)(h.uv.u, h.uv.v) : white; if (m.hasNormalMap) // perturbacja wektora normalnego zgodnie z teksturą h.N = uChangeSpace((*m.normalMap)(h.uv.u, h.uv.v), h.Ex, h.N, h.Ez); if (m.isDiffuse) // odbicie rozproszone return m.dr * texMapped * traceRay(scatter(r, h.ip, h.N), s, d-1); else // emisja return m.e; } |
Tak wygląda całka w przypadku naszego zagadnienia oświetlenia:
Pamiętamy jednak, że prosty estymator był mało zadowalający. Radziliśmy sobie używając wielu prób prostego estymatora, obliczając średnią arytmetyczną uzyskanych wyników. Ten lepszy estymator realizowany jest przez funkcję , która wcześniej służyła jedynie do realizacji antyaliasingu.
Podsumowanie
Nauczyliśmy się jak przy pomocy śledzenia promieni generować zdjęcia wirtualnej sceny. Poznaliśmy zjawiska takie jak odbicia zwierciadlane, odbicia rozproszone, transmisja z refrakcją, miękke cienie, krwawienie kolorów. Dowiedzieliśmy się jak teksturować powierzchnie obiektów. Niedługo pojawią się kolejne części tego cyklu artykułów, w których scalimy wszystkie dotychczas poznane zjawiska tworząc program obsługujący ścieżki typów E (D|S|T)* L. Poznamy zjawisko kaustyk. Dodamy do kamery symulację systemu optycznego, który pozwoli uzyskać efekt głębi ostrości. Dowiemy się trochę więcej o symulacji lokalnych zjawisk związanych z odbijaniem światła. Na sam koniec wprowadzimy obsługę siatek trójkątów i bardziej skomplikowanej geometrii, żeby zobaczyć coś ciekawszego od kilku kulek :)
Odnośniki:
[1] http://informatyka.wroc.pl/node/415
[2] http://informatyka.wroc.pl/node/416
[3] http://informatyka.wroc.pl/node/407
[4] http://informatyka.wroc.pl/upload/raytracing/rt2pack.zip