2011-10-26

Błędy numeryczne

Ich występowanie jest szczególnie uciążliwe dla takich obiektów jak sfera i torus. Podejrzewam, że podobne problemy mogą dotyczyć wszystkiego co nie jest płaskie, czyli trójkąta. Patrząc na kod historycznie kod trójkąta (intersekcja, normalna, mapowanie tekstury, mapowanie wypukłości) nie zmienił się ani trochę od samego początku. Nie ma w nim żadnego względnego porównywania względem błędu, a porównanie bezwzględne jest dla wartości z zakresu 0...1.

Z błędami numerycznymi jako takimi można sobie radzić dobierając gdzie trzeba porównania z pożądaną precyzją. Problem zaczyna się kiedy odległości pomiędzy obiektami są bardzo małe albo bardzo duże lub same obiekty są bardzo małe albo bardzo duże. Z pierwszym możemy sobie poradzić normalizując scenę do żądanych rozmiarów np. AABB o wymiarze $\pm 1$. Z drugim możemy sobie poradzić skalując obiekt odpowiednio do współrzędnych lokalnych. W obu przypadkach stabilizujemy wartości na których liczymy i tym samym błąd.

Ciągle pozostaje jednak problem kiedy odległości pomiędzy obiektami różnią się od siebie znacznie. Wtedy po normalizacji sceny odległości pomiędzy nimi mogą stać się mniejsze od precyzji co może wywołać błędy numeryczne. I to chyba będzie waży argument w moim przypadku przenieść obliczenia na double z powrotem. Na razie ten problem zostawiam sobie na później.

Dodatkowo dla sfery i torusa tuż po wyliczeniu wynikowej intersekcji sprawdzamy, czy normalna i promień padający tworzą kąt mniejszy od 90 stopni. Pozwala to nam zapobiegać błędom numerycznym kiedy promień jest prawie prostopadły do powierzchni. Takie błędy lepiej uciąć u źródeł niż dodawać w wielu miejscach kod na ich korekcje.

Kolejna poprawka dotyczyła kontroli argumentów dla funkcji o ograniczonej dziedzinie (pierwiastkowanie, odwrotne funkcje trygonometryczne). Jeśli przekraczają one dziedzinę nieznacznie są do niej zaokrąglane np. -0.000001 do 0 dla pierwiastkowania.

Kolejna poprawka dotyczyła ścisłego kontrolowania współrzędnych UV, co okazało się szczególnie istotnie dla torusa, gdzie błąd może znacznie przekroczyć wartości 0...1.

W ogóle torus to wyjątkowo trudny obiekt. Bardzo niestabilny numerycznie. Wydaje mi się, że w przyszłości wszelkie powierzchnie wyższego rzędu będą służyć tylko do generacji siatki trójkątów. Są one niestabilnie numerycznie i każda nowa rzecz jaką dodawałem do raytracera wymagała jakiś specjalnych działań by wszystko zaczęło dobrze działać.

Kolejna poprawka tym razem na kod sfery dotyczy parametry m2. Generalnie powinien się on zawierać w granicach od 0 do 1. Ale jeśli promień nadchodzi z bardzo daleka niekoniecznie okazuje się to być prawdą. Musimy go z zadaną precyzją przyciąć do wartości 0...1.

I kolejna poprawka tym razem na torus i sferę. Najpierw wyliczamy czy promień uderza w AABB i od razu znajdujemy do niego odległość. Następnie przesuwamy promień na granicę lokalnego AABB wykorzystując tą odległość. Dzięki temu błędy dla promieni nadchodzących z daleka zostały znacznie zredukowane. To także powinno zmniejszyć błędy numeryczne liczenia intersekcji torusa i go ustabilizować pozwalając na jego korygowanie niezależnie jak pada promień.

Jak widać większość błędów kręci się wokół sfery i torusa.

Zmiany wymuszone przez pojawienie się refrakcji

Po pierwsze promień stał się obiektem immutable. Pojawiła się w nim jedna nowa właściwość mówiąca czy jest promieniem odbitym czy załamanym.

Informacja ta jest wykorzystywana przez przez dwie właściwości RefractedSideObject i ReflectedSideObject, których kod gettera korzysta z tej właściwości promienia by wyznaczyć w jaki obiekt wchodzi promień i z jakiego wychodzi. Te zaś służą do wyznaczenia współczynnika pochłaniania światła jak i refrakcji.

Wraz z tym musiała się pojawić nowa funkcja FindOuterObject do wyznaczania obiektu zewnętrznego poprzez wypuszczenie promienia. Jest ona także używana do prekalkulowania tłumienia ośrodka w jakim znajduję się źródło światła.

Kolejna dużą rzeczy może nie tyle wymuszoną ale dodaną dla większej realności są półprzeźroczyste cienie. Czyli chodzi nie tyle o zbadanie czy stoi coś na drodze do światła od punktu uderzenia, ale o sprawdzenie czy to coś jest przeźroczyste. Tak długo jak jest cień rzucany przez to źródło światła to nie absolutna czerń. Do jego kalkulacji wykorzystujemy współczynnik pochłaniania refrakcji, który mówi nam jaki procent energii przechodzi przez powierzchnie podczas refrakcji oraz jaki jest współczynnik pochłania światła w danym obiekcie.

Oczywiście przy braku globalnej iluminacji cienie takie są sztuczne - nie są załamane. Ale i tak wygląda to lepiej niż cienie solidne.

No i dodać należało wspomniane wyżej współczynniki pochłaniania energii podczas refrakcji i refleksji. Pochłanianie to nie zależy od drogi promienia. Poza tym w przeciwieństwie do współczynnika pochłaniania jest to warstwa materiału. Można ją np. wykorzystać do symulacji cienkiego foliowego nadruku. Współczynnik pochłaniania jest dany dla materiału i dla obiektu obowiązuje w całej jego objętości w jego wnętrzu.

Każdy materiał od teraz ma współczynnik refrakcji. Jest on dany liczbą. Nie jest warstwą materiału. Trudno sobie wyobrazić jak by to miało działać. Dodatkowo każdy obiekt ma właściwość czy jest domknięty, tylko obiekty domknięte mogą mieć współczynniki refrakcji.

W FindOuterObject musiałem zrezygnować z kodu na pominięcie detekcji obiektu samego ze sobą, by obsłużyć przypadki r>R dla torusa.

Po dodaniu refrakcji poprzednia procedura intersekcji sfery z promieniem okazała się zawodna. Jeśli poprzednie uderzenie nastąpiło w sferę informacja o tym, z której strony nastąpiło uderzenie już nam nie wystarczy do tego by stwierdzić czy uderzenie w ogóle nastąpi. Np. promień uderza z zewnątrz. Ponowna intersekcja jeśli myślimy tylko o promieniu odbitym nie ma sensu. Ale sprawy się mają bardziej skomplikowanie kiedy dochodzi nam promień załamany.

Najważniejszy jest początkowy fragment kodu:
if ((a_source_ray_intersection != null) && 
    a_source_ray_intersection.SceneObject == this)
{
    Debug.Assert((a_ray.RaySurfaceSide == RaySurfaceSide.RefractedSide) ||
                 (a_ray.RaySurfaceSide == RaySurfaceSide.ReflectedSide));
    if (a_source_ray_intersection.BackHit ^ 
       (a_ray.RaySurfaceSide == RaySurfaceSide.ReflectedSide))
    {
        return Scene.NoIntersection;
    }
    else

Od teraz promień niesie z sobą informację, czy jest promieniem odbitym, załamanym albo jakimś innym (np. od kamery). Ta dodatkowa informacja nie pojawiła się tylko by uzupełnić kod sfery. Równie dobrze tych warunków mogłoby nie być. Intersekcję możemy za każdym razem liczyć. Ale niejako wraz z pojawieniem się refrakcji, kod sfery przestał działać, a dodatkowa informacja w promieniu, która pozwala kodowi sfery dalej działać jest wymagana ogólnie do działania refrakcji.

Przejście z powrotem z float na double

Przejście na float nie przyspieszyło obliczeń. Obliczenia na float powodują duże błędy w torusie. Nawet pomijając torus poziom precyzji rzędu $10^{4-}$ wydaje się zbyt mały. Sceny w których różnice w rozmiarach ujawnią te błędy wydają się bardzo możliwe. Dla double poziom precyzji ustalił się na poziomie $10^{-12}$. To już astronomicznie duża liczba. Poza tym liczba błędów w torusie zmniejszyła się. Wprowadzony także błąd na przystawanie trójkątów także mógł zostać wyeliminowany. Float-y ciągle obecne są w tablicach kolorów i wartości by oszczędzać pamięć.