2011-12-12

Transformata Fouriera

Istnieje wiele transformat, które generalnie można powiedzieć przekształcają funkcję z jednej przestrzeni w inną. Z reguły istnieje także transformata odwrotna. Nie wszystkie funkcje da się także poddać takiej transformacji. Transformata Fouriera to takie przekształcenie, które możemy zinterpretować jako przekształcenie funkcji w dziedzinie czasie w funkcję w dziedzinie częstotliwości.

$f(\omega) = \int_{-\infty}^{\infty} f(x)\ e^{- 2\pi i x t}dt$

Transformata odwrotna:

$f(t) = \int_{-\infty}^{\infty} f(\omega)\ e^{2 \pi i x \omega}d\omega$

W wyniku transformacji funkcji o wartościach rzeczywistych zmiennej rzeczywistej transformata Fouriera jest symetryczna względem zera. Tylko taka symetria po stronie częstotliwości gwarantuje nam otrzymanie funkcji o wartościach rzeczywistych podczas transformacji odwrotnej. Ma to szczególne znaczenie podczas dyskretnej transformaty Fouriera kiedy niedokładności w symetrii są częste.

Po stronie częstotliwości widmo sygnału interpretujemy jako sumę sinusoid o danych częstotliwościach, przesunięciu fazowym i amplitudzie. Ponieważ we wzorach mamy całkę, a nie sumę, musiała by to być suma nieskończona. Przechodząc z całek na sumy przechodzimy na dyskretną transformatę Fouriera. Wynikiem transformaty jest funkcja liczby zespolonej - jej moduł to amplituda dla danej częstotliwości, jej przesunięcie fazowe to przesunięcie fazowe sinusoidy o danej częstotliwości.

Ponieważ z reguły jesteśmy zainteresowaniu amplitudą widma jej wykresy najczęściej podaje się po stronie częstotliwości. Przesunięcia fazowe z reguły pomija się. Zauważmy, że dla różnie poprzesuwanych składowych sinusoid widmo amplitud jest takie same. I drugiej strony przesuniecie funkcji w czasie nie zmienia jej widma. Właściwość ta jest wykorzystywana podczas kompresji obrazów i filmów. Szczególnym przypadkiem dobrania przesunięć fazowych jest funkcja o minimalnym przesunięciu fazowym (jest to taka funkcja w której w dowolnym czasie $<0,t>$ energia sygnału jest maksymalna).

W wyniku tego, że dyskretna transformata Fouriera operuje na dyskretnych wartościach wprowadza ona błędy. Wspomnianą wcześniej asymetrię funkcji rzeczywistej i nierzeczywistą postać jej transformacji odwrotnej. Innym rodzajem błędu jest to, że pojedyncza sinusoida może być reprezentowana poprzez rozmyty pik po stronie częstotliwości. Dlatego po retransformatą warto wziąć jej część rzeczystą (nie moduł).

W dalszej części będziemy mówili tylko o transformacji dyskretnej.

Oto przykład transformaty funkcji sinc:


Teoretycznie jej transformatą powinien być prostokąt. W rzeczywistości ponieważ wartości są dyskretne, a funkcja nie jest nieskończona, transformata nie jest idealnym prostokątem i wykazuje oscylacje na brzegach skoku. Są one tym mniejsze im nasz sinc jest dłuższy.

Dla transformaty odwrotnej: sinc po stronie częstotliwości będzie prostokątem po stronie czasu z takimi samymi oscylacjami na brzegach. Ponieważ nasz sinc jest skończony w czasie brakuje nam widma dla wyższych częstotliwości i one uniemożliwiają nam zreprodukowanie skoku i zlikwidowanie oscylacji - można powiedzieć ich wygaszenie.

Przykład jak może wyglądać transformacja złożenia dwóch sinusoid o f=650 i f=1500. Oczywiście w zależności jak dobierzemy krok próbkowania naszych ciągłych funkcji, ilości próbek czasu i ilości próbek częstotliwości możemy otrzymać wykres bardziej lub mniej przypominający piki.

Warto podkreślić, że dyskretna transformata jest bezwymiarowa podobnie jak wartości wejściowych próbek. Chcąc otrzymać wartości częstotliwości na osi musimy je policzyć.

Skończony sygnał w czasie nieokresowy wymaga do opisu nieskończonego widma częstotliwości. Podobnie jest z nieciągłościami w czasie. I drugą stronę: skończone widmo nie powoli nam na opisanie skończonego w czasie sygnału. Te dwa ograniczenia silnie wiążą się z DFFT, której liczba próbek jest skończona.

Metody resamplowania

Nie chodzi mi tyle o filtry, ale o sposób samego resamplowania.

Sposób pierwszy. Robimy podobnie jak podczas zmiany rozdzielczości obrazka (tak działał resampler po staremu) Tworzymy tablice subpixeli, każdy subpixel przyciągamy do odpowiedniego subpixela. Taką tablicę resamplujemy. Wady to duża ilość pamięci, niedokładność związana z nieuwzględnianiem prawdziwych współrzędnych promienia.

Do metody drugiej i trzeciej musimy użyć innego sposobu pamiętania wartości kolorów promieni.

Ponieważ rezygnujemy z bitmapy subpixeli, która także nawiasem mówiąc byłaby kłopotem podczas resamplownania adaptywnego, które zamierzam dodać, musimy jakoś pamiętać wartości kolorów promieni, które będą potrzebne do zresamplowania próbek które jeszcze nie zostały zrenderowane - z innych prostokątów. Przez prostokąt rozumiem tutaj obszar znajdujący się obok, który nie został jeszcze zrenderowany. W obecnej chwili scena jest dzielona na w miarę równe prostokąty o określonej licznie subpikseli i wszystkie są one renderowane równolegle. Ponieważ w procesie resamplingu wartość subpikseli granicznych może być potrzeba do wyliczenia pikseli w sąsiednim prostokącie musimy je jakoś zapamiętać. Dodatkowo aby uniknąć problemów w wyścigiem zabronione jest renderowanie sąsiednich, wzajemnie na siebie wpływających prostokątów. Możemy na przykład pamiętać je w tablicy hashującej, liście o wielkości równej liczbie pikseli na scenie - każdy element listy byłby null-em, ewentualnie bardzo rzadko tablicą wartości i współrzędnych promieni. Z tym, że podczas zmiany sposobu renderowanie np. z prostokątów na linia po linii mogło by się okazać, że kod wyliczający czy dany piksel będzie jeszcze wykorzystywany mógłby okazać się nie uniwersalny. No i jeśli renderujemy prostokątami to subpiksele brzegowe trzeba pamiętać maksymalnie dla 3 kwadratów - potrzebny byłby jakiś licznik referencji - brzegowe 2, rogowe - 3, albo 2 jeśli są na brzegu sceny. Dosyć kłopotliwe... Dlatego zdecydowałem się na rozwiązanie wymagające więcej pamięci, ale za to bardziej uniwersalne. Resamplowanie działa w 2D. Zmiana rozdzielczości obrazka za pomocą filtru 2D jest wolniejsza niż dwukrotne użycie filtru 1D z pionie i w poziome. Wynik jest ten sam. Kształt filtru powstaje poprzez obrót filtra 1D wokół Y, jeśli tak możemy nieprecyzyjnie powiedzieć. Ogólny wzór na wartość zrekonstruowanego piksela ma postać:

$\displaystyle C(x,y) = \frac{\sum f(x-x_i,y-y_i)c(x_i,y_i) }{\sum f(x-x_i,y-y_i)}$

Wartości wag filtra wyliczone z tego wzory nie wymagają normalizacji, nawet jeśli pole pod filtrem jest różne od 1, a ni tylko jeśli wyliczone wagi obarczone są błędem numerycznym. Normalizacja jest wzór - normalizacja to podzielenie wszystkich wag przez sumę wag, tak by ta suma wynosiła 1.

$C(x,y)$ to wynikowy kolor. $c(x_i,y_i)$ to kolor promienia. $f(x-x_i,y-y_i)$ to wartość filtru. Zgodnie z tym wzorem korzystnie jest pamiętać kolor jako iloraz dwóch liczb, dzięki czemu możemy dodawać do piksela wartości nowych promieni. Znika nam także problem pamiętania wartości wyliczonych promieni na potrzeby sąsiednich prostokątów.

Pomimo normalizacji może dojść do przestrzelenia koloru. Weźmy dwie wagi: 2 i -0.5. Niech oba kolory mają wartość 1, wtedy wartość zrekonstruowana to 1, ale jeśli drugi kolor ma w tym punkcie wartość 0, to otrzymamy 2/1.5. Widać to zwłaszcza podczas Samplowania Jitter 1x1 na Lanczos, który ma przedziały ujemne.

No i teraz dwie pozostałe metody.

Druga to liczenie za każdym razem wagi filtru. W metodzie trzeciej budujemy subtablicę wag o określonej wielkości. Centrujemy ją na pikselu wynikowym. Patrzymy w którą komórkę subtablicy trafił promień i bierzemy z niej wagę. Przy czym subtablica może mieć większą rozdzielczość niż wartość ilości próbek na piksel, dzięki czemu może być dokładniejsza niż metoda pierwsza. Przy odpowiednio dużej wielkości powinna być tak samo efektywna jak metoda druga, a przy tym szybsza.

Metody 2 i 3 w stosunku do 1 generują obraz nieznacznie przesunięty. Przesunięcie jest znacznie mniejsze od 0.5 i ... nie wiem skąd się bierze.

2011-12-11

Zmiany w kodzie

Cały kod obecnie odpowiadający za resamplowanie generuje obrazy nie tak dobrej jakości jakby mógł, gdyż zamiast opierać się na prawdziwej odległości pozycji promienia na płaszczyźnie projekcji od pozycji piksela który chcemy zrekonstruować, obecnie najpierw generowana jest bitmapa subpixeli i ona jest reskalowana. Ta metoda ma też inną wadę, wymaga strasznie dużo pamięci. Tak więc wraz z przepisaniem samplowania (zmianie musi ulec sposób generowania próbek) i resamplowania (sposób rekonstrukcji) zamiast bitmapy subpixeli będziemy tylko pamiętali tyle subpixeli ile potrzebne nam jest do przyszłej rekonstrukcji sąsiednich pixeli.

Niejaki przy okazji pojedyńczy obiekt kamery zostanie zastąpiony kolekcją kamer z jedną aktywną. Każda kamera zostanie połączona z obiektem klasy Film. Klasy Film utworzą swoją własną kolekcję, tak więc wiele kamer będzie mogło korzystać z tego samego filmy. Do klasy Film zostanie zaś przeniesione większość opcji z RenderOptions. W tym wszelkie informacje pozwalające zbudować łańcuch postprocesorów obrazu.

InterestRects zostały przeniesione do Filmu. ... I w trakcie zrezygnowałem z tego. Prawie tego nie używałem. Poza tym dla o ile na początku się starałem by rzeczywiście dobrze definiować obszary zainteresowania to później nie miałem już na to sił - liczba testowych scen to już przeszło tysiąc. Poza tym tak naprawdę to w wyniku zmian w kodzie popsuć się może wszystko i wszędzie nie tylko w obszarach zainteresowań.

SceneElement został zastąpiony przez SceneActor. SceneActor to baza dla wszystkich obiektów mogących występować na scenie - posiadających pozycję. ScemeElement jest od teraz klasą bazową dla SceneActor - posiada nazwę i właściwość Scene. Jest też dolnym ograniczeniem parametru generycznego dla ObjectList.

Klasa Material dziedziczy z SceneElement. Klasa MaterialsList dziedziczy z ObjectsList.

Dodałem do kamery dwie nowe metody UVToPixel i PixelToUV. Przestrzeń UV to zawierający się od 0 do 1 zbiór punktów z których może zostać wypuszczony promień. Dla kamer które wysyłają promienie z prostokątnego obszaru płaszczyzny wartość UV nie różnią niczym w sensie znaczenia od koordynatów tekstury. Co innego kiedy dojdą kamery które wysyłają promienie z wycinka walca (panoramiczne) albo z wycinka sfery (tzw. fish-eye). Procedura Camera.CreateRay wymaga do wygenerowania promienia współrzędnych UV.

Samplery generują próbki zawsze w zakresie 0 do 1. Parametry do generacji do liczba wierszy i kolumn. Ich iloczyn to ilość punktów do wygenerowania. Dla samplerów jednorodnych na każdy prostokąt będzie przypadać jedna próbka, dla samplerów niejednorodnych niekoniecznie.

Postprocesory tworzą łańcuch ustalany w runtime. Parametry dla nich (gamma correction, tone mapping) zniknęły z klasy Film. Łańcuch postprocesorów podłączony jest do klasy Film. Na razie zrezygnowałem z utworzenia osobnej kolekcji łańcuchów postprocesorów. Taki łańcuch postprocesorów jest serializowany do xml bezpośrednio.

Parametrami dla postprocesingu jest jest bitmapa wejściowa i wyjściowa. Z reguły mogą być one takie same. Jeśli tak nie jest to od postprocesora wymaga się by podał jaki ma być jej rozmiar. O ile ma się on zmienić. Wymaga się też od niego by prawidłowo przeliczył na nowy obszar prostokątu. Sama klasa odpowiedzialna za wywołania postprocesorów będzie się starała utrzymać taką bitmapę wyjścia do następnego razu aby nie tworzyć ich za każdym razem. Taka bitmapa wyjściowa będzie mu później wielokrotnie w trakcie generowania poszczególnych obszarów przesyłana mu do uzupełnienia tych obszarów i będzie wykorzystana w dalszym łańcuchu. Jeśli postprocesor nie skaluje obrazu bitmapa wyjściowa będzie równa wejściowej. Posprocesor winien w niej zmodyfikować tylko zadany obszar.

Ponieważ kod resamplujący ulegnie zmianie postanowiłem obecny kod, który jest tak naprawdę kodem na zmniejszanie wielkości obrazka o SampleSize przekształcić w uniwersalny postprocesor pozwalający na skalowanie wyrenderowanej bitmapy.

Problem z mapowaniem

Przy okazji mapowania na torus powierzchni danych wielomianem pojawił się problem, który na razie odłożyłem w przyszłość, ale prędzej czy później trzeba go będzie jakoś rozwiązać. Mianowicie torus dany bezpośrednio mapowany jest w sześcian UVW od 0 do 1. W przypadku powierzchni wyższego rzędu mapowanej na torus nie wiemy nic o geometrii bryły i jest ona mapowana w UVW proporcjonalnie do $r/(R+r)$. Na razie poradziłem sobie z tym rozciągając powierzchnię wielomianową w Y do sześcianu, dzięki czemu dobrze współpracuje z mapperem. A następnie ręcznie ją skalując by przypominała torus. Niestety trzeba to robić ręcznie.

Wyobraźmy sobie teraz siatkę przypominającą spłaszczoną sferę. Chcemy ją zamapować sferycznie. Teraz ponieważ sfera wymaga UVW w pełnym sześcianie musielibyśmy najpierw przeskalować siatkę by zawierała się w sześcianie. A następnie spłaszczyć ją skalowaniem.

Rozwiązanie pierwsze to skalowanie UVW zanim zostanie przekształcone w UV. Z tym, że na razie (nie wiem czy słusznie) rozdzialiłem świat UVW i UV. Nie ma przejścia od Volume do Layer. I raczej nie chciałbym tego robić póki nie będzie to absolutnie potrzebne.

Rozwiązanie drugie to dodanie do mapperów skalowania i translacji, być może nawet macierzy, którą UVW jest transformowane przez mapowaniem.

Problem z nieskończonościami

Użycie nieskończoności (np. jako granice na AABB) powoduje wiele problemów. Niekończoności nie porównują się dobrze ze sobą, wiele standardowych funkcji (Math.Sign) zwraca wyjątek po ich napotkaniu. Zastąpienie ich Double.MaxValue i Double.MinValue jest dużo lepszym pomysłem. Z jednym wyjątkiem liczby takie konwertują się na tekst, ale z powrotem już nie, gdyż otrzymujemy Overflow Exception. Tak więc musimy wprowadzić nasze własne minima i maksima i się ich trzymać: Double.MaxValue/10 i Double.MinValue/10

Powierzchnie dane wielomianem

A dokładniej wielomiany w funkcji trzech zmiennych x, y, z, przeważnie w postaci uwikłanej. Powierzchnie, zwłaszcza te zamknięte, to z reguły wielomiany stopnia 4 lub 6. Przy obliczaniu intersekcji z nimi bardzo ważne jest przybliżenie punktu startu promienia do AABB, by zmniejszyć błędy obliczeniowe, które mogą się pojawić. Dla niektórych powierzchni AABB nie da się określić wprost, gdyż zależy ona od parametrów takiego wielomianu (np. promień dla sfery).

W przypadku niektórych kształtów likwidowałem parametry funkcji, gdyż odpowiadały one tak naprawdę za skalowanie. Czasami też przesuwałem powierzchnię (manipulując x, y, z) by znalazła ona się w środku układu współrzędnych oraz by zawierała się w przedziale od -1 do 1.

Na ten moment sprawdziłem renderowanie prawie wszystkich powierzchni zamkniętych podanych na MathWorld. I wystarczy tego. Cel jakim było w ogóle zrenderowanie takich powierzchni został osiągnięty. Na jego potrzeby musiałem napisać całkiem spory kawałek kodu, który bierze wielomian w postaci tekstu, parsuje do tablicy wyrażeń, buduje z niej RPN, z RPN drzewo. Drzewo wyrażeń można zoptymalizować - wymnażanie przez stała, dodawanie, mnożenie wyrazów, uproszczenia związane z 0 i 1. Oprócz tego wyrażenie można rozszerzyć, czyli wymnożyć, wypotęgować co się da - zbudować jedną wielką sumę wyrazów. Dzielenie przez inne wielomiany, czyli funkcje wymierne nie jest obsługiwane, podobnie jak wykładniki potęgowania nie będące liczbami. Dla zmiennych wykładniki te powinny być dodatnimi liczbami całkowitymi. Brakuje kodu na redukcję wyrażeń, czyli czegoś odwrotnego do rozszerzania. Gradient jest liczony wprost z funkcji bez jej rozszerzania. Dla wyliczenia intersekcji podstawiamy do wielomianu równanie kierunkowe linii, rozszerzamy je i grupujemy po stopniu t. Każda taka grupa to współczynnik dla wielomianu intersekcji. Następnie kompilujemy kod który zwraca nam taki wielomian z już wyliczonymi współczynnikami po podaniu mu jako parametrów kierunku i punktu startu promienia. Kompilacja następuje poprzez zbudowanie LINQ Expression. Kod wynikowy jest bardzo szybki.

Szczegółów kodu nie zamierzam podawać. Ale kilka rad po napisaniu czegoś takiego mam. Drzewo wyrażenia powinno być niezmienne. Totalnie niezmienne drzewo będzie wymagało skopiowania całego drzewa przy dowolnej zmianie, a będziemy to robić bardzo często. Najlepiej kod drzewo umieścić w osobnej przestrzeni nazw i poukrywać wszystkie do manipulowania jego strukturą bezpośrednio (rodzic, dzieci) i na zewnątrz wystawić tylko zbiór ustalonych metod. Dlatego najlepiej metody zmieniające drzewo jakoś zaizolować, a na zewnątrz wystawić tylko metody nie modyfikujące drzewa, ale zwracające nowe zmienione drzewo. Co do elementów drzewa to powinno się ono składać tylko z symboli, dodawania, mnożenia, dzielenia, potęgowania. Unarne minusy i odejmowanie powinniśmy zastąpić dodawaniem i mnożeniem przez -1. Dzielenie także jak najszybciej konwertować na mnożenie (obsługujemy dzielenie tylko przez stałe). Przez symbol rozumiemy wyrażenie A*N^B, gdzie A to stała, N to nazwa symbolu, B to jego potęga. A=0 to zero, B=0 to liczby. Dodatkowo warto od początku założyć, że dodawanie i mnożenie może mieć więcej dzieci, dzięki czemu znajdowanie grup typu ABCD+EFGH będzie szybkie i proste. Dzięki temu, że potęgowanie mamy od razu w symbolu nie musimy przeszukiwać drzewa nadmiernie. A stała w symbolu przez którą mnożymy sprawi, że często nie będziemy musieli tworzyć nowych symboli liczb. Na samym początku redukujemy liczby i przenosimy potęgowanie do symboli tak długo jak pozbędziemy się potęgowania i dzielenia. Pozbywamy się także minusów i odejmowania. Spłaszczamy drzewo tak by dodawania i mnożenie nie posiadało innych dzieci - odpowiednio dodawania i mnożenia i w takiej postaci staramy się utrzymywać drzewo.

Sam kod intersekcji nie różni się zbytnio od kodu intersekcji dla torusa. Z tym, że inaczej jest budowany wielomian intersekcji.

Można oczywiście policzyć gradient i współczynniki wielomianu intersekcji na piechotę. Jednakże jest to bardzo pracochłonne. Poza tym mało uniwersalne - mając cały kod przetwarzający możemy implementować dowolne powierzchnie dane wielomianem.

W moim przypadku największe problemy sprawił mi kształt serca: $(x^2+\frac{9}{4} y^2+z^2-1)^3-x^2z^3-\frac{9}{80}y^2z^3$

Serce bez AABB:


Serce z AABB - start promienia jest przybliżany do AABB i dla takiego punktu startu budowany jest wielomian intersekcji.


Serce z AABB - likwidacja miejsc zerowych.


Sprawdziłem jak wyglądają współczynniki i same miejsca zerowego dla promienia odbitego od błędnych punktów. Okazało się, że następuje ponowna detekcja uderzenia w powierzchnie od której się odbijamy. Możemy powiększyć Constants.MINIMIAL_DISTANS = 1e-8 - miejsce od którego szukamy pierwiastków. Ale w tym przypadku ponieważ pierwszy współczynnik wielomianu intersekcji okazał się mniejszy od dokładności Constants.DOUBLE_PRECISION = 1e-12 postanowiłem napisać kod eliminujący takie współczynniki, przesuwający pozostałe współczynniki i zmiejszający stopień wielomianu. Tak mały współczynnik może odpowiadać za miejsce zerowe w granicach 1e-5. Taki zabieg nigdzie nie wyprodukował dodatkowych błędów. Błędy nie znikły całkowicie, ale są już prawie niewidoczne, a na pewno rozmyją przy supersamplowaniu.

Serce z przesuwaniem punktu startu wzdłuż normalnej.


Serce ciągle ma dwa piksele za ciemne. Nie pozostało mi nic innego jak odsunąć punkt startu od powierzchni wzdłuż normalnej. Wartość przy której znikły ostatnie błędy to Constants.POLYNOMIAL_SURFACE_FROM_MOVE_DELTA = 1e-5. Jest to bardzo dużo.

Przykład powierzchni jabłka z refrakcją:


Ponieważ dysponuję kodem na obiekt torusa liczonego bezpośrednio i teraz mam do dyspozycji torus dany przez wielomian mogę porównać ich jakość i szybkość. Implementacja oparta na wielomianie jest ponad dwa razy wolniejsza. Przede wszystkim największy współczynnik wielomianu, który po uproszczeniu jeden jedynką tutaj jest wyliczany. Prawdopodobnie pozostałe współczynniki też nie są w optymalnej postaci. Coś za coś. Np. dla serca liczba współczynników to 217 po pełnym rozwinięciu. To daje obraz zabawy podczas wyliczania współczynników wielomianu.

Kod obiektu sceny:

public class PolynomialSurface : RenderableObject
{
    protected PolynomialFunction m_polynomial;
    private Func<Vector3, Vector3> m_gradient; 
    private Func<Vector3, Vector3, Polynomial> m_intersection;
    private PolynomialSolver m_solver;

    [YAXNode]
    public string Equation;

    [YAXNode]
    public string Parameter1 = "";

    [YAXNode]
    public string Parameter2 = "";

    [YAXNode]
    public string Parameter3 = "";

    [YAXNode]
    public AABB LocalBoundBox;

    private Vector3 m_aabb_minimum;
    private Vector3 m_aabb_maximum;

    public PolynomialSurface(Vector3 a_right, Vector3 a_up) :
        base(a_right, a_up)
    {
        Name = "Polynomial surface";
        m_aabb_minimum = Vector3.MAXIMUM;
        m_aabb_maximum = Vector3.MINIMUM;
        Closed = false;
        LocalBoundBox = AABB.INFINITY;
        m_uv_mapper = new SphericalUVMapper();

        Update(UpdateFlags.All);
    }

    public override Intersection GetIntersection(
        Intersection a_source_ray_intersection, Ray a_ray)
    {
        Vector3 local_dir;
        Vector3 local_start;
        bool back_hit = false;

        if (!TransformToLocalToBoundBox(a_ray, out local_dir, out local_start))
            return Scene.NoIntersection;

        if ((a_source_ray_intersection != null) && 
            (a_source_ray_intersection.SceneObject == this))
        {
            if (a_ray.RaySurfaceSide == RaySurfaceSide.RefractedSide)
            {
                local_start += ((WorldToLocalNormal * 
                    a_source_ray_intersection.Normal).Normalized) *
                    -Constants.POLYNOMIAL_SURFACE_FROM_MOVE_DELTA;
            }
            else
            {
                local_start += ((WorldToLocalNormal * 
                    a_source_ray_intersection.Normal).Normalized) *
                    Constants.POLYNOMIAL_SURFACE_FROM_MOVE_DELTA;
            }
        }

        Polynomial p = m_intersection(local_start, local_dir);

        double dist = m_solver.Solve(p);

        if (dist == Double.PositiveInfinity)
            return Scene.NoIntersection;

        Vector3 local_pos = local_start + local_dir * dist;

        Vector3 normal = m_gradient(local_pos);

        back_hit = normal * local_dir > 0;

        if (back_hit && OneSide)
            return Scene.NoIntersection;

        Vector3 world_pos = LocalToWorld * local_pos;
        double world_dist = (world_pos - a_ray.Start).Length;

        Intersection intersection = new Intersection()
        {
            PrevIntersection = a_source_ray_intersection, 
            SceneObject = this,
            SourceRay = a_ray,
            LocalPos = local_pos,
            Dist = world_dist,
            Scene = Scene,
            BackHit = back_hit,
            Pos = world_pos
        };

        // Numerical errors.
        if (intersection.Normal * intersection.SourceRay.Dir >= 0)
            return Scene.NoIntersection;

            m_aabb_minimum = Vector3.Minimize(local_pos, m_aabb_minimum);
            m_aabb_maximum = Vector3.Maximize(local_pos, m_aabb_maximum);

        return intersection;
    }

    public override Vector3 GetNormal(Intersection a_intersection)
    {
        if (a_intersection.BackHit)
        {
            return (LocalToWorldNormal * 
                -m_gradient(a_intersection.LocalPos)).Normalized;
        }
        else
        {
            return (LocalToWorldNormal * 
                m_gradient(a_intersection.LocalPos)).Normalized;
        }
    }

    public override Vector2 GetUV(Intersection a_intersection)
    {
        return UVMapper.Map(a_intersection.UVW);  
    }

    public override Vector3 GetUVW(Intersection a_intersection)
    {
        return base.GetUVW(a_intersection) * 0.5 + Vector3.HALF;
    }

    public override string ToString()
    {
        return String.Format("PolynomialSurface: {0}", Name);
    }

    public override void GetTangents(Intersection a_intersection, 
        out Vector3 a_tangent_x, out Vector3 a_tangent_y)
    {
        if (a_intersection.BackHit)
        {
            a_tangent_x = Vector3.CrossProduct(
                Up, a_intersection.Normal).Normalized;
            a_tangent_y = Vector3.CrossProduct(
                a_tangent_x, a_intersection.Normal).Normalized;
        }
        else
        {
            a_tangent_x = Vector3.CrossProduct(
                a_intersection.Normal, Up).Normalized;
            a_tangent_y = Vector3.CrossProduct(
                a_intersection.Normal, a_tangent_x).Normalized;
        }
    }

    public override Scene Scene
    {
        get
        {
            return base.Scene;
        }
        set
        {
            if (Scene != null)
                Scene.RenderEnd -= RenderEnd;

            base.Scene = value;

            if (Scene != null)
                Scene.RenderEnd += RenderEnd;
        }
    }

    private void RenderEnd(bool a_all)
    {
        if (!a_all)
            return;

        if (LocalBoundBox != AABB.INFINITY)
        {
            if ((!LocalBoundBox.Minimum.IsAlmostRelativeEquals(m_aabb_minimum)) &&
                (!LocalBoundBox.Maximum.IsAlmostRelativeEquals(m_aabb_maximum)))
            {
                LogManager.GetCurrentClassLogger().Error("Bound box for {0}", Name);
                LogManager.GetCurrentClassLogger().Error("Actually: {0}", 
                    LocalBoundBox);
                LogManager.GetCurrentClassLogger().Error("Should be: {0}", 
                    new AABB(m_aabb_minimum, m_aabb_maximum));
            }
        }
    }

    protected override AABB GetLocalBoundBox()
    {
        return LocalBoundBox;
    }

    protected override void RenderStart(RenderStartPhase a_phase)
    {
        base.RenderStart(a_phase);

        if (a_phase == RenderStartPhase.PrepareObjectToRender)
        {
            string equation = ApplyParameters(Equation);
            m_polynomial = new PolynomialFunction(equation);
            m_gradient = m_polynomial.Gradient;
            m_intersection = m_polynomial.Intersection;

            m_solver = new PolynomialSolver(PolynomialSolveMethod.Sturm,
                Scene.RenderOptions.RootFindingMethod);

            m_aabb_minimum = Vector3.MAXIMUM;
            m_aabb_maximum = Vector3.MINIMUM;
        }
    }

    private string ApplyParameters(string a_equation)
    {
        a_equation = ApplyParameter(a_equation, Parameter1);
        a_equation = ApplyParameter(a_equation, Parameter2);
        a_equation = ApplyParameter(a_equation, Parameter3);

        return a_equation;
    }

    private string ApplyParameter(string a_equation, string a_parameter)
    {
        if (String.IsNullOrWhiteSpace(a_parameter))
            return a_equation;

        string[] ar = a_parameter.Split('=');

        if (ar.Length < 2)
            return a_equation;

        a_equation = a_equation.Replace(ar[0], "(" + ar[1].Replace(",", ".") + ")");

        return a_equation;

    }

    public override void ScaleAbsolute(double a_scale)
    {
        Scale *= a_scale;

        base.ScaleAbsolute(a_scale);
    }
}

Twierdzenie Sturma

Twierdzenie Sturma pozwala nam stwierdzić ile rzeczywistych pierwiastków wielomianu $f(x)$ znajduje się w określonym przedziale $\left (a, b \right )$, gdzie $a<b$, $f(a) \neq 0$, $f(b) \neq 0$. Liczba pierwiastków rzeczywistych w takim przedziale to $\sigma(a)-\sigma(b)$, gdzie $\sigma(x)$ to liczba zmian znaków w tzw. sekwencji Sturma. Sekwencje tą budujemy następująco:

$p_0(x)=f(x)$
$p_1(x)=f'(x)$
$p_2(x)=-rem(p_0, p_1)=p_1 q_0 - p_0$
$p_3(x)=-rem(p_1, p_2)=p_2 q_1 - p_1$
...
$p_m(x)=-rem(p_{m-1}, p_{m-2}) = p_{m-1} q_{m-2} - p_{m-2}$

$rem$ oznacza resztę z dzielenia wielomianów.

Przykładowo jeśli obliczeniach otrzymamy: ++--, to $\sigma=1$, jeśli +-+- to $\sigma=2$.

Twierdzenie Sturma działa tylko dla pierwiastków pojedynczych. Jeśli wielomian posiada podwójne pierwiastki musimy się ich najpierw pozbyć budując wielomian bez nich: $\displaystyle ^P/_{GCD(P, P')}$. Przez dzielenie rozumiemy dzielenie wielomianów. GCD to największy wspólny dzielnik. W moim przypadku sekwencje tą wykorzystujemy do obliczeń związanych z raytracingiem. Tak więc póki co nie musiałem się przejmować przypadkami wielokrotnych pierwiastków. Jedynym jest torus zredukowany do sfery.

Tak więc zakładamy, że nasz wielomian posiada tylko pojedyncze pierwiastki.

Dla wielokrotnych twierdzenie nie będzie działało dobrze. Zmiana w znaku zostanie wykryta zawsze dla pierwiastka nieparzystego. Ale może zostać także stwierdzona dla parzystego. Przynajmniej tak to działa do wielomianów stopnia czwartego - sprawdzone doświadczalnie. Innymi słowy nigdy nie dostaniemy fałszywej liczby pierwiastków ale możemy nie dostać prawdziwej. W przypadku torusa zredukowanego do sfery problemem są raczej błędy numeryczne niż ten błąd.

Sekwencja Sturma posiada następujące właściwości:

$1^\circ \quad p_m(x) \neq 0$ \: dla $x \in (a, b)$
$2^\circ \quad \text{jeśli} \:\: 0<i<m-1$ \text{to} $p_{i-1}(x) = p_{i+1}(x) = 0$ \: \text{dla} \: $x \in (a, b)$
$3^\circ \quad \text{jeśli} \:\: 0<i<m$ to $p_{i-1}(x) \neq 0$ i $p_{i}(x) \neq 0$ \: \text{dla} \: $x \in (a, b)$
$4^\circ \quad$ istnieje takie $\alpha$, że dla $p_0(x)=0$ zachodzi $sign(p_0(x-\alpha)) = sign(p_0(x+\alpha))$

$1^\circ$ Wynika z definicji reszty. Musi ona być co najmniej o stopień mniejsza od dzielnej. Stopień $p_1$ też jest o stopień niższy od $p_0$. Stopień wielomianu maleje wraz z wyrazami sekwencji. Jeśli wielomian jest stopnia zerowego to nasza sekwencja liczy tylko jeden wyraz. Dla wielomianu pierwszego stopnia jego pochodna zawsze będzie różna od zera. Dla wielomianów wyższych stopni, ponieważ pierwiastki są pojedyncze to $p_0(x)$ i $p_1(x)$ nie mają wspólnych miejsc zerowych, reszta $p_2$ musi być różna od zera. I nie ma miejsc zerowych wspólnych z $p_1$ i $p_0$. W konsekwencji kolejne reszty nie mogą być zerowe, ich stopień z braku wspólnych miejsc zerowych będzie się zmniejszał zawsze o jeden do $p_{m-1}=ax-b$, $p_m=a$.

$2^\circ$ Sekwencja musi mieć co najmniej 3 wyrazy. Jeśli $p_1=0$ to z definicji wynika, że $p_2=p_1 q_0 - p_0 = -p_0$. Jeśli $p_i=0$, gdzie $i>1$ to z definicji na ogólny wyraz ciągu: $p_i(x)= p_{i-1} q_{i-2} - p_{i-2}$, wynika, że $p_{i+1}(x)= p_{i} q_{i-1} - p_{i-1}$, dla $p_i=0$ mamy: $p_{i+1}=-p_{i-1}$.

$3^\circ$ Wynika z rozważań z punktu $1^\circ$. Stopień wielomianu w sekwencji zawsze maleje od jeden do stałej i nie mają one wspólnych miejsc zerowych.

$4^\circ$ Wynika z definicji pochodnej. Pochodna wielomianu osiąga zera tylko w ekstremach wielomianu. Pomiędzy nimi, czyli w obszarze występowania miejsc zerowych wielomianu jest monotoniczna. O ile miejsca zerowe nie są wielokrotne, wtedy miejsce zerowe może być jednocześnie ekstremum. Dla pojedynczych miejsc zerowych wielomian zawsze przechodzi przez miejsce zerowe i zmienia znak.

Wszystkie te cztery właściwości pozwalają nam udowodnić poprawność twierdzenia. Udowodnimy, że funkcja $\sigma$ wyznaczająca ilość zmian znaków jest malejąca i maleje o jeden tylko przy przejściu przez miejsce zerowe wielomianu.

Sekwencja składa się z wielomianów, których miejsca zerowe są pojedyncze. Jeśli przechodzimy przez zero w wyrazie sekwencji wyraz ten zmienia znak. Z $3^\circ$ wynika, że taki wyraz otoczony jest przez wyrazy niezerowe. Ale nie jest nigdzie powiedziane, że inne wyrazy nie mogą być zerami.

Przypuśćmy, że miejsce zerowe nastąpiło w jednym lub więcej wyrazach różnych od $p_0$. Zgodnie z $4^\circ$ i $2^\circ$ mamy cztery kombinacje znaków przed miejscem zerowym: -++, --+, +--, ++- i po:
--+, -++, ++-, +--. tak więc liczba znaków w sekwencji nie ulega zmianie.

Przypuśćmy, że miejsce zerowe wystąpiło w wyrazie $p_0$ (może też wystąpić w innych wyrazach sekwencji). Zgodnie z $4^\circ$ $p_1$ nie zmienia znaku w otoczeniu miejsca zerowego. Same wielomian $p_0$ zmienia znak z uwagi na to, że posiada on tylko pojedyncze miejsca zerowe. Jeśli wielomian $p_0$ zmienia znak z + na - to znak pochodnej to -, tak więc możliwa kombinacja znaków przed to: +-, -+, zaś po: --, ++.

Widzimy więc, że tylko podczas przejścia przez miejsca zerowe wielomianu $p_0$ liczba zmian znaków w sekwencji maleje o jeden. Możemy więc powiedzieć, że liczba zmian znaków w przedziale $(a, b)$ równa się $\sigma(a) - \sigma(b)$. Ponieważ zmiana znaków następuje przy przejściu przez miejsce zerowe wielomianu $p_0$ jest to tym samym liczba miejsc zerowych w tym przedziale.

Twierdzenie Sturma wykorzystuje się do znajdowania miejsc zerowych wielomianów. Dla wielomianów stopnia 5 i wyższego nie istnieje rozwiązanie dane wzorami wprost. Poza tym wzory na pierwiastki wielomianu stopnia 3 i 4 są niestabilne numerycznie (dokładność double jest niewystarczająca by policzyć je dokładnie, lub by policzyć je dobrze kiedy są blisko siebie - wtedy często wyliczamy jako pierwiastek ekstremum między nimi). Najpierw metodą np. binarnego poszukiwania wyodrębniamy obszary w których zmiana znaku jest pojedyncza, a następnie jedną z metod aproksymacji (Newton, ModifedRegulaFasli, Bisection) przybliżamy się do miejsca zerowego, korzystając z faktu, że jest pojedyncze, i znak z jednej jego strony jest przeciwny do znaku z drugiej strony.

W przypadku obliczeń do raytracingu, interesuje nas tylko pierwsze rozwiązanie większe od zera, a ściśle od Constants.MINIMUM_DISTANCE.

Kod klasy implementującej sekwencje Sturma:

public struct SturmSequence
{
    private readonly Polynomial[] m_seq;
    public readonly int Length;

    public SturmSequence(Polynomial a_poly)
    {
        a_poly = a_poly.RemoveZeroes();

        m_seq = new Polynomial[a_poly.Order + 1];

        Length = 1;
        m_seq[0] = a_poly;

        if (a_poly.Order == 0)
            return;

        Length++;
        m_seq[1] = a_poly.Differentiate();

        while (m_seq[Length - 1].Order != 0)
        {
            Length++;
            m_seq[Length - 1] = -(m_seq[Length - 3] % m_seq[Length - 2]);
            if (m_seq[Length - 1].Order == 0)
                return;
        }
    }

    public static Polynomial SingleRoots(Polynomial a_poly, 
        double a_precision = Constants.DOUBLE_PRECISION)
    {
        return a_poly / Polynomial.GCD(a_poly, a_poly.Differentiate(), a_precision);
    }

    public Polynomial this[int a_index]
    {
        get
        {
            if (a_index >= Length)
                throw new IndexOutOfRangeException();

            return m_seq[a_index];
        }
    }

    public int NumberOfSignChangesAt(double a_x)
    {
        int sc = 0;

        double c1 = m_seq[0].Evaluate(a_x);

        for (int i = 1; i < Length; i++)
        {
            double c2 = m_seq[i].Evaluate(a_x);

            if ((c1 * c2 < 0) || (c1 == 0.0))
                sc++;

            c1 = c2;
        }

        return sc;
    }

    public int NumberOfSignChangesAtPositiveInfinity()
    {
        int sc = 0;

        double c1 = m_seq[0].LeadCoefficient;

        for (int i = 1; i < Length; i++)
        {
            double c2 = m_seq[i].LeadCoefficient;

            if ((Math.Sign(c1) * Math.Sign(c2) < 0) || ((c1 == 0.0) && (c2 != 0.0)))
                sc++;

            c1 = c2;
        }

        return sc;
    }

    public int NumberOfSignChangesAtNegativeInfinity()
    {
        int sc = 0;

        double c1 = -m_seq[0].LeadCoefficient;
        if (m_seq[0].Order % 2 == 0)
            c1 = -c1;

        for (int i = 1; i < Length; i++)
        {
            double c2 = -m_seq[i].LeadCoefficient;
            if (m_seq[i].Order % 2 == 0)
                c2 = -c2;

            if ((Math.Sign(c1) * Math.Sign(c2) < 0) || ((c1 == 0.0) && (c2 != 0.0)))
                sc++;

            c1 = c2;
        }

        return sc;
    }
}
Kod klasy solvera wykorzystującego powyższą sekwencje znajdujący tylko jedno miejsce zerowe większe od Constants.MINIMUM_DISTANCE większe od zera:
public static class Constants
{
    public const int STURM_ANY_ROOTS_MAXIMUM_ITERATIONS = 32;
    public const int STURM_FIND_ONE_ROOT_AREA_MAXIMUM_ITERATIONS = 800;

    public const double MINIMAL_DISTANT = 1e-8;
}

public struct PolynomialSolverAprox
{
    private readonly SturmSequence m_sturm_sequence;
    private RootFinder m_root_finder;

    public PolynomialSolverAprox(Polynomial a_polynomial, RootFinder a_root_finder)
    {
        m_sturm_sequence = new SturmSequence(a_polynomial);
        m_root_finder = a_root_finder;
    }

    public double FindOneRootInRange(double a_a, double a_b, int a_num_roots_at_a)
    {
        for (int i = 0; 
             i < Constants.STURM_FIND_ONE_ROOT_AREA_MAXIMUM_ITERATIONS; i++)
        {
            double mid = (a_a + a_b) / 2;
            int atmid = m_sturm_sequence.NumberOfSignChangesAt(mid);
            int nroots = a_num_roots_at_a - atmid;

            if (nroots == 1)
            {
                Polynomial p = m_sturm_sequence[0];
                return m_root_finder.FindRoot(a_a, mid, (x) => p.Evaluate(x));
            }
            else if (nroots == 0)
                a_a = mid;
            else
                a_b = mid;
        }

        return a_a;
    }

    public double SolveOne()
    {
        int nroots_at_a = m_sturm_sequence.NumberOfSignChangesAt(
           Constants.MINIMAL_DISTANT);
        int nroots_at_pos_inf = 
            m_sturm_sequence.NumberOfSignChangesAtPositiveInfinity();

        int nroots = nroots_at_a - nroots_at_pos_inf;

        if (nroots == 0)
            return Double.PositiveInfinity;

        double a = Constants.MINIMAL_DISTANT;
        double b = a + 1.5;
        int nroots_at_b = nroots_at_a;

        for (int its = 0; its < Constants.STURM_ANY_ROOTS_MAXIMUM_ITERATIONS; its++)
        {
            nroots_at_b = m_sturm_sequence.NumberOfSignChangesAt(b);

            if (nroots_at_a - nroots_at_b >= 1)
                break;

            b *= 2.0;
        }

        if (nroots_at_a == nroots_at_b)
            return Double.PositiveInfinity;

        return FindOneRootInRange(a, b, nroots_at_a);
    }
}