2011-09-02

Skalowanie - zmianie w kodzie trójkąta

Trójkąt to taki wstęp do implementacji siatki.

W kodzie nic nie zmieniłem. W przypadku trójkąta zamiast transformować promień do współrzędnych lokalnych wszelkich obliczeń dokonujemy dla normalnych wierzchołków i współrzędnych wierzchołków w współrzędnych świata. Dzięki temu unikamy zbędnych transformacji. Korzystamy tutaj z faktu, że wszystkie punkty trójkąta leżą w jednej płaszczyźnie. Płaszczyzna jest innym obiektem dla którego moglibyśmy zrezygnować z transformacji świat-obiekt i obiekt-świat. Swoją drogą płaszczyzna z czasem przestanie być przeze mnie używana, jest to dość nienaturalny obiekt.

Same zaś współrzędne i normalne trójkąta są pamiętane we współrzędnych lokalnych. Każda zmiana w macierzy transformacji (przeskalowanie, translacja, rotacja obiektu) powoduje wyliczenie nowych współrzędnych i normalnych świata. Dzięki temu podczas wielokrotnego zmieniania właściwości obiektu unikamy kumulowania się błędów obliczeniowych.

Żeby było mi łatwiej bezpośrednio manipulować parametrami trójkąta ustawianie ich jest możliwe we współrzędnych świata.

Dodatkowo do właściwości obiektów renderowalnych musiałem wprowadzić nową macierz - macierz transformacji normalnych świat-obiekt, będąca odwrotnością macierzy transformacji normalnych obiekt-świat.

2011-09-01

Skalowanie - zmiany w kodzie płaszczyzny

Zamian wymagał tylko kod intersekcji z promieniem. W przeciwieństwie do sfery kod zwracający normalną w punkcie nie uległ zmiany. Normalna jest dla płaszczyzny wszędzie taka samo i jako taka została wcześniej wyliczona.

Oryginalny kod
public override Intersection GetIntersection(Ray a_ray)
{
    float denom = Normal * a_ray.Dir;

    if (OneSide)
    {
        if (denom > 0)
            return Scene.NoIntersection;
    }

    if (a_ray.PrevHit == this)
        return Scene.NoIntersection;

    if (!denom.IsZero())
    {
        float dist = ((Normal * a_ray.Start) - Distance) / (-denom);

        if (dist <= 0)
            return Scene.NoIntersection;
        else
        {
            bool backHit = (denom > 0);

            if (backHit && OneSide)
                return Scene.NoIntersection;
            else
            {
                return new Intersection()
                {
                    SceneObject = this,
                    SourceRay = a_ray,
                    Dist = dist,
                    Scene = Scene,
                    BackHit = backHit,
                    Pos = a_ray.HitPoint(dist)
                };
            }
        }
    }
    else
        return Scene.NoIntersection;
}
Zmieniony kod:
public override Intersection GetIntersection(Ray a_ray)
{
    Ray ray_local = a_ray.Transform(WorldToLocal);

    bool backHit = false;
    float dist;

    if (ray_local.PrevHit == this)
    {
        if (ray_local.PrevBackHit)
        {
            backHit = true;

            if (OneSide)
                return Scene.NoIntersection;
            else
            {
                float v = -(ray_local.Start * ray_local.Dir);
                float m2 = Radius * Radius - (ray_local.Start.SqrLen - v * v);
                dist = v + (float)Math.Sqrt(m2);
            }
        }
        else
            return Scene.NoIntersection;
    }
    else
    {
        float c2 = ray_local.Start.SqrLen;
        float r2 = Radius * Radius;

        if (c2 > r2)
        {
            float v = -(ray_local.Start * ray_local.Dir);
            if (v <= 0)
                return Scene.NoIntersection;
            else
            {
                float m2 = r2 - (c2 - v * v);
                if (m2 < 0)
                    return Scene.NoIntersection;
                else
                    dist = v - (float)Math.Sqrt(m2);
            }
        }
        else
        {
            backHit = true;

            if (OneSide)
                return Scene.NoIntersection;
            else
            {
                float v = -(ray_local.Start * ray_local.Dir);
                float m2 = r2 - (c2 - v * v);
                dist = v + (float)Math.Sqrt(m2);
            }
        }
    }

    var pos = ray_local.Start + ray_local.Dir * dist;
    pos = LocalToWorld * pos;
    dist = (pos - a_ray.Start).Length;

    return new Intersection()
    {
        SceneObject = this,
        SourceRay = a_ray,
        Dist = dist,
        Scene = Scene,
        BackHit = backHit,
        Pos = pos
    };
}
Zmiany polegały na dodaniu transformacji promienia z współrzędnych świata do lokalnego na początku i policzeniu punktu uderzenia i dystansu w współrzędnych świata na końcu. Dodatkowo w układzie lokalnym płaszczyzna leży (w moim przypadku) w płaszczyźnie XZ, tak więc normalna do wersor Y, zaś dystans płaszczyzny do środka układu współrzędnych to 0. To pozwala nam trochę uprościć kod.

Problemy związane z liczbami zmiennoprzecinkowymi

Liczby zmiennoprzecinkowe są przechowywane w systemie binarnym, zaś wyświetlane i zapisywane z reguły w formacie dziesiętnym. Każda konwersja liczby zmiennoprzecinkowej na string i z powrotem może ale nie musi zmienić jej wartości. Głównie chodzi o to, że liczba wymierna w jednym formacie może być niewymierna w drugim. Jeśli mamy zamiar porównywać takie liczby to możemy zawsze przed porównaniem wymusić taką konwersję.

W ostateczności liczby możemy przechowywać w postaci tablicy bajtów (BitConverter), a w tekście w postaci HEX, w ten sposób zachowamy jej bitowy obraz.

Ale nie zawsze. Biorąc losowy ciąg bajtów i robiąc z niego liczbę zmiennoprzecinkową, może być nie poprawna, i w momencie załadowania jej do FPU jej wartość zostanie zmieniona. Tak długo jak to się nie stanie jej wygląd bitowy się nie zmieni. Dotyczy to np. niepoprawnych liczb które podpadają pod nieskończoności.

Porównywanie liczb zmiennoprzecinkowych bezpośrednio nawet jeśli nic z powyższego nie miało miejsca nie jest dobrym pomysłem. Wyniki mogą się różnić.

Spodziewamy się wartości znormalizowanej wektora o długości 1. Możemy dostać coś bliskiego 1, i zawsze z różną precyzją, zależnie czy np. jest to normalna policzona dla sfery, czy trójkąta. Nie pozostaje nam wtedy nic innego jak napisać kod porównujący liczby z pewną precyzją i kontrolować ją w debug assertami. W razie asserta musimy sami zdecydować czy to błąd, czy może brak precyzji i nasza stała precyzji musi się stać bardziej tolerancyjna.

Żeby nie było za łatwo wyniki w debug, release, x86, x64 mogą się od siebie różnic. Przyczyn może być wiele. Główna to taka, że liczby zmiennoprzecinkowe są przechowywane wewnętrznie w FPU w formacie 80-bitowym i taka jest domyślny jego tryb działania w C#. Każdy zapis takiej liczby do pamięci to jej zaokrąglenie. Może ona nastąpić np. na skutek braku optymalizacji pomiędzy wersją debug, a release. A np. pomiędzy x86, a x64 w wyniku innej optymalizacji operacji.

Szczególnie fatalnym przypadkiem takiego zaokrąglania z 80 do 64 bitów jest pojawienie się zera albo nieskończoności.

Liczby zmiennoprzecinkowe mogą występować w dwóch postaciach znormalizowanej i dla bardzo małych liczb w postaci zdenormalizowanej. W moim przypadku operacje na liczbach zdenormalizowanych są 7 razy wolniejsze. Liczby te są bardzo małe i są daleko poza zakresem wszystkich dokładności używanych w raytracerze. Liczby zdenormalizowane służą zasypaniu luki w możliwych do uzyskania liczb zmiennoprzecinkowych w okolicach zera. Istnieje pewna możliwość, że wejdziemy w zakres liczb zdenormalizowanych i znacząco spowolnimy renderowanie. Na razie nie przejmuje się tym.

Odjęcie od siebie dwóch różnych liczb zmiennoprzecinkowych może dać w wyniku zero. Podobnie jak podzielenie jakiejś liczby przez np. 2. Dzieje się tak jeśli brakuje nam ekspomenty na zapisanie wyniku. Podobnie jak w zero możemy w pewnym momencie wejść w nieskończoności. Jednak to wejście w zero jest najmnniej intuicyjne. O błędach tego rodzaju trzeba pamiętać w pętlach. Podczas obliczeń można ich spórbować uniknąć tak przekształcając równanie by zastąpić odejmowanie dodawaniem. Dla raytracingu nie to sensu. Z reguły są to tak małe liczby daleko poza zakresem precyzji z jaką operujemy na obiektach sceny. Klasycznym przykładem podawanym w różnych miejscach jest takie przepisanie wzorów na pierwiastki kwadratowe by wyeliminować odejmowanie jeśli b>0:

$\begin{align*}
x_{1,2} = & \frac{-b-\sqrt{b^2-4ac}}{2ac}=\\& \frac{(-b-\sqrt{b^2-4ac})(-b+\sqrt{b^2-4ac})}{2ac(-b+\sqrt{b^2-4ac})}=\\& \frac{4ac}{2ac(-b+\sqrt{b^2-4ac})}=\\&\frac{2}{-b+\sqrt{b^2-4ac}}
\end{align*}$

Skalowanie - zmiany w kodzie sfery

Kod oryginalny:
public override Intersection GetIntersection(Ray a_ray)
{
    bool backHit = false;

    Vector3 c = Pos - a_ray.Start;
    double c2 = c.Length * c.Length;
    double r2 = Radius * Radius;
            
    double dist;

    if (a_ray.PrevHit == this)
    {
        if (a_ray.PrevBackHit)
        {
            backHit = true;

            if (OneSide)
                return Scene.NoIntersection;
            else
            {
                double v = c * a_ray.Dir;
                double m2 = r2 - (c2 - v * v);

                dist = v + Math.Sqrt(m2);
            }
        }
        else
            return Scene.NoIntersection;
    }
    else
    {
        double v = c * a_ray.Dir;
        double m2 = r2 - (c2 - v * v);

        if (c2 > r2)
        {
            if (v <= 0)
                return Scene.NoIntersection;
            else
            {
                if (m2 < 0)
                    return Scene.NoIntersection;
                else
                    dist = v - Math.Sqrt(m2);
            }
        }
        else
        {
            backHit = true;

            if (OneSide)
                return Scene.NoIntersection;
            else
                dist = v + Math.Sqrt(m2);
            }
        }

    return new Intersection()
    {
        SceneObject = this,
        SourceRay = a_ray,
        Dist = dist,
        Scene = Scene,
        BackHit = backHit,
        Pos = a_ray.HitPoint(dist)
    };
}

public override Vector3 GetNormal(Intersection a_intersection)
{
    if (a_intersection.BackHit)
        return (Pos - a_intersection.Pos).Normalized;
    else
        return (a_intersection.Pos - Pos).Normalized;
}
Nowy kod:
public Ray Transform(Matrix4 a_matrix)
{
    return new Ray()
    {
        Dir = (a_matrix.Rotation * this.Dir).Normalized,
        Start = a_matrix * this.Start,
        SourceIntersection = this.SourceIntersection,
        Depth = this.Depth

    };
}

public override Intersection GetIntersection(Ray a_ray)
{
    Ray ray_local = a_ray.Transform(WorldToLocal);
    bool backHit = false;
    float dist;

    if (ray_local.PrevHit == this)
    {
        if (ray_local.PrevBackHit)
        {
            backHit = true;

            if (OneSide)
                return Scene.NoIntersection;
            else
            {
                float v = -(ray_local.Start * ray_local.Dir);
                float m2 = Radius * Radius - (ray_local.Start.SqrLen - v * v);
                dist = v + (float)Math.Sqrt(m2);
            }
        }
        else
            return Scene.NoIntersection;
    }
    else
    {
        float c2 = ray_local.Start.SqrLen;
        float r2 = Radius * Radius;

        if (c2 > r2)
        {
            float v = -(ray_local.Start * ray_local.Dir);
            if (v <= 0)
                return Scene.NoIntersection;
            else
            {
                float m2 = r2 - (c2 - v * v);
                if (m2 < 0)
                    return Scene.NoIntersection;
                else
                    dist = v - (float)Math.Sqrt(m2);
            }
        }
        else
        {
            backHit = true;

            if (OneSide)
                return Scene.NoIntersection;
            else
            {
                float v = -(ray_local.Start * ray_local.Dir);
                float m2 = r2 - (c2 - v * v);
                dist = v + (float)Math.Sqrt(m2);
            }
        }
    }

    var pos = ray_local.Start + ray_local.Dir * dist;
    pos = LocalToWorld * pos;
    dist = (pos - a_ray.Start).Length;

    return new Intersection()
    {
        SceneObject = this,
        SourceRay = a_ray,
        Dist = dist,
        Scene = Scene,
        BackHit = backHit,
        Pos = pos
    };
}

public override Vector3 GetNormal(Intersection a_intersection)
{
    if (a_intersection.BackHit)
        return -(LocalToWorldNormal * a_intersection.UVW).Normalized;
    else
        return (LocalToWorldNormal * a_intersection.UVW).Normalized;
}
Główne zmiany są na początku badania intersekcji kiedy promień jest transformowany do lokalnego układu współrzędnych. Dzięki temu z obliczeń wylatuje Pos. Dodatkowo poprzenosiłem wgłąb wyliczenia niektórych zmiennych do momentu, aż są one potrzebne. Na samym końcu przed zwróceniem wyniku intersekcji wyliczamy miejsce uderzenia i dystans w współrzędnych świata.

W kodzie na wyliczenie normalnej korzystamy od teraz z specjalnej macierzy na transformację normalnych. Zauważmy, że współrzędne UVW są zmapowane na sferę.

Pozostały kod na wyliczanie UVW, UV, tangensów nie ulega zmianie.

Transformacja normalnej

Promień uderza w skalowaną nierównomiernie sferę, czyli elipsoidę. Nie chcemy jednak pisać osobnego kodu na obsługę elipsy. Promień przed analizą intersekcji poddajemy transformacji świat-obiekt. W wyniku czego nasz promień uderza w sferę o środku w środku układu współrzędnych. Znajdujemy punkt intersekcji. Transformujemy go z powrotem obiekt-świat. To działa. Teraz chcemy wyliczyć normalną w punkcie uderzenia naszej elipsy. Podobnie transformujemy punkt uderzenia świat-obiekt. Wyliczamy normalną dla sfery. Transformujemy ją obiekt-świat. To nie działa.


  
    
      
    
    
      
    
    
      
    
  
  
  
    
      
        image/svg+xml
        
        
      
    
  
  
    
      
      
    
    
      
      
    
    
      
      
    
    
      
      
    
    
      
      
    
    
    
      
      
    
    
      
      
    
    
      
      
    
    
      
      
    
    
      
      
    
    
    
      
      
    
    
      
      
    
    
      
      
    
    
      
      
    
    
      
      
    
    
  


Na pierwszym rysunku widzimy normalne dla sfery. Na drugim normalne sfery zmapowane obiekt-świat błędnie, za pomocą takiej samej metody jak dla punktów intersekcji. Na trzecim widzimy to co chcielibyśmy uzyskać.

Weźmy prostą w 2D $\mathbf{f=ax}$, jej pochodna $\mathbf{f'=a}$, jej normalna (nazwijmy ją tak) $\mathbf{f''=-a}$. Teraz tak jak na naszym rysunku skalujemy zwężając dwukrotnie w osi Y. Równanie naszej prostej $\mathbf{f=\frac{1}{2}ax}$, zaś normalnej $\mathbf{f''=-\frac{1}{2}a}$. Załóżmy teraz dwukrotnie rozciągamy w osi X. Równanie naszej prostej, pochodnej i normalnej zmienia się podobnie. Czyli skalując o a skalujemy normalną o $\mathbf{\frac{1}{2}}$. Podobnie moglibyśmy wykazać wychodząc z definicji pochodnej.

Czyli jeśli nasza macierz skalowania to $\mathbf{S}$, to musimy ją zastąpić w macierzy transformacji obiekt-świat macierzą $\mathbf{S^{-1}}$. Translację powinniśmy pominąć, gdyż mamy do czynienia nie z punktem, a z wektorem. Rotację musimy zachować. Jeśli macierz naszej transformacji dla wierzchołków ma postać $\mathbf{M=TRS}$, to dla normalnych będzie miała postać $\mathbf{M=(TRS^{-1})_{ROT}=RS^{-1}}$

Przy okazji zauważmy, że mając macierz obiekt-świat: $\mathbf{M=TRS}$, macierz transformacji normalnej możemy zapisać jako:

$\begin{split}MN &= (TRS^{-1})_{ROT} \\
&= ((T^{-1})^{T}(R^{-1})^{T}(S^{-1})^T)_{ROT} \\
&= (((XRS)^{-1})^T)_{ROT} \\
&= ((M^{-1})^T)_{ROT} \\
&= ((M^{T})^{-1})_{ROT} \end{split}$

W przekształceniach skorzystaliśmy z faktu, że macierz odwrotna macierzy rotacji jest równa jej transpozycji, macierz transponowana macierzy skalowania jest równa jej samej. X to przekształcona macierz transpozycji. Ponieważ bierzemy część rotacji X ulega ona zniesieniu, równie dobrze możemy wpisać dowolną macierz która w części rotacyjnej 3x3 jest macierzą jednostkową, w szczególności T. Zauważmy, że elementy macierzy T nie wpływają na elementy macierzy rotacji 3x3 w macierzy S i R, ani ich iloczynu. Tak więc macierz transpozycji obiekt-świat dla normalnych jest równa odwrotnej transponowanej macierzy obiekt-świat dla wierzchołków, albo jak kto woli odwrotnej transponowanej macierzy obiekt-świat dla wierzchołków. Czasami może nam się to przydać jeśli nie wiemy jak taka macierz została złożona.

Przykład kodu:
private void UpdateTransformationMatrices()
{
    m_local_to_world = Matrix4.CreateTranslation(Pos) *
        new Matrix4(Right, Up, Forward) * Matrix4.CreateScale(Scale);
    m_world_to_local = m_local_to_world.Inverted;
    m_local_to_world_normal = m_world_to_local.Transposed.Rotation;
}

Kontrola oświetlenia sceny

Podczas renderowania sceny może nam wyjść, że jasność danego punktu (a konkretnie jednej ze składowych może przekroczyć dopuszczalny zakres). Mogą nam także wyjść liczby ujemne (w wyniku artefaktów po resamplowaniu, kiedy kernel filtru posiada składowe ujemne). Wartości ujemne powinniśmy zwyczajnie przyciąć do zera.

W tym momencie zakładam, że całe obliczenia dokonywane są na liczbach zmiennoprzecinkowych, gdzie 0 to zero, a 1 to biel.

Niezależnie od podanych niżej metod normalizacji jasności i tak warto podczas konwersji takiego obrazu do bitmapy sprawdzać zakres.

Pierwsze co możemy zrobić to przyciąć problematyczne składowe RGB. Kolory przestrzelonych jasności nie zostaną w tym wypadku zachowane.

Możemy więc cały kolor RGB przeskalować, kolor zostanie zachowany, a punkt zostanie przyciemniony. W tym wypadku zostanie zachwiana dynamika jasności. Jasność przejaśnionego obszaru będzie nieproporcjonalnie jasna w stosunku do otoczenia.

No to może po zrenderowaniu sceny skalujemy wszystkie punkty proporcjonalnie. Wadą jest nienaturalny wygląd sceny. Będzie ona zbyt ciemna, tylko po to by jakiś obszar nie był zbyt jasny.

Możemy zastosować prostą normalizację na całym wyrenderowanym obrazie: $C'=\frac{C}{C+1}$. Zdecydowanie metoda ta powinna dać najlepsze rezultaty.

Oprócz tego jest cała grupa metod do konwersji obrazów HDR. Tą konwersję nazywamy mapowaniem tonów. Algorytmy te są z reguły nieliniowe i starają się one oddać działanie oka i ludzkiego mózgu. Na ich zaimplementowanie przyjdzie czas...

Przykłady poszczególnych metod:

Przycięcie składowych RGB

foreach (var p in a_rect.EnumPixels())
{
    ColorFloat c = a_pixels.GetColor(p.X, p.Y);
    c = new ColorFloat(
        c.R > 1 ? 1 : c.R,
        c.G > 1 ? 1 : c.G,
        c.B > 1 ? 1 : c.B);
    a_pixels.SetColor(p.X, p.Y, c);
}



Skalowanie proporcjonalne składowych RGB

foreach (var p in a_rect.EnumPixels())
{
    ColorFloat c = a_pixels.GetColor(p.X, p.Y);
    float m = Math.Max(Math.Max(c.R, c.G), c.B);
    if (m > 1)
    {
        c = new ColorFloat(c.R / m, c.G / m, c.B / m);
        a_pixels.SetColor(p.X, p.Y, c);
    }
}



Skalowanie proporcjonalne wszystkich pikseli

float m = 0;

foreach (var p in new Rectangle(0, 0,
    a_pixels.Width, a_pixels.Height).EnumPixels())
{
    ColorFloat c = a_pixels.GetColor(p.X, p.Y);
    m = Math.Max(m, c.R);
    m = Math.Max(m, c.G);
    m = Math.Max(m, c.B);
}

if (m > 1)
{
    foreach (var p in new Rectangle(0, 0,
        a_pixels.Width, a_pixels.Height).EnumPixels())
    {
        ColorFloat c = a_pixels.GetColor(p.X, p.Y);
        c = new ColorFloat(c.R / m, c.G / m, c.B / m);
        a_pixels.SetColor(p.X, p.Y, c);
    }
}



Proste normalizacja

foreach (var p in new Rectangle(0, 0,
    a_pixels.Width, a_pixels.Height).EnumPixels())
{
    ColorFloat c = a_pixels.GetColor(p.X, p.Y);

    c = new ColorFloat(
        c.R / (c.R + 1),
        c.G / (c.G + 1),
        c.B / (c.B + 1));

    a_pixels.SetColor(p.X, p.Y, c );
}



2011-08-30

Dlaczego przejście z double na float nie przyspieszyło renderowania

Postanowiłem się temu dokładniej przyjrzeć.

Zanim w ogóle uderzymy w C# przyjrzymy się jak to wygląda w C++. Włączamy release, wszelkie optymalizację, zostawiamy informację o debagowaniu pozwalające w miarę śledzić kod, wyłączamy inlinowanie funkcji.
static float Test1(float a, float b)
{
    return a + 0.5f + b;
}

static float Test2(float a, float b)
{
    return sin(a) + 0.5f + b;
}

static double Test3(double a, double b)
{
    return a + 0.5 + b;
}

static double Test4(double a, double b)
{
    return sin(a) + 0.5 + b;
}

static float Test5(double a, float b)
{
    return a + 0.5 + b;
}

static void Test()
{
    double r = 0;

    {
        Stopwatch sw;
        for (int i=0; i<100000000; i++)
            r += Test1(i, i + 5);
        sw.Stop();
        auto t = sw.GetElapsedMilliseconds();
        printf("float simple: %f\n", t);
    }

    {
        Stopwatch sw;
        for (int i=0; i<100000000; i++)
            r += Test3(i, i + 5);
        sw.Stop();
        auto t = sw.GetElapsedMilliseconds();
        printf("double simple: %f\n", t);
    }

    {
        Stopwatch sw;
        for (int i=0; i<10000000; i++)
            r += Test2(i, i + 5);
        sw.Stop();
        auto t = sw.GetElapsedMilliseconds();
        printf("float sin: %f\n", t);
    }

    {
        Stopwatch sw;
        for (int i=0; i<10000000; i++)
            r += Test4(i, i + 5);
        sw.Stop();
        auto t = sw.GetElapsedMilliseconds();
        printf("double sin: %f\n", t);
    }

    {
        Stopwatch sw;
        for (int i=0; i<1000000000; i++)
            r += Test5(i, i + 5);
        sw.Stop();
        auto t = sw.GetElapsedMilliseconds();
        printf("mixed simple: %f\n", t);
    }

    printf("%f", r);
}
Wyniki x64:
float simple: 312.602579
double simple: 309.030617
float sin: 739.128367
double sin: 1250.607130
mixed simple: 4314.021684
Wyniki x86:
float simple: 682.324569
double simple: 498.409225
float sin: 797.726457
double sin: 745.483504
mixed simple: 12319.519241
Można zgłupieć. Raz jedno jest szybsze, raz drugie. Ale generalnie float nie przyspiesza, jedynie w x64 sin jest liczony szybciej na float, ale to raczej wynika z tego, że ktoś popsuł wersję double.

Skorzystanie z _controlfp( _PC_24, MCW_PC ); (przestawienie FPU w tryb pojedynczej precyzji) nic nie zmienia. Jedynie zwalnia, gdyż biblioteka standardowa przestawia i pewnie przywraca status FPU. Tutaj w sieci znalazłem informacje, że sin nie przyspiesza w pojedynczej precyzji ale takie log, pow, div już tak.

Moja szybka wersja funkcji sin, ale pewnie wymagające wielu założeń:
static float mysin(float a)
{
    __asm
    {
        fld a
        fsin
        fst a
    }
}

static double mysin(double a)
{
    __asm
    {
        fld a
        fsin
        fst a
    }
}
Wyniki wersji przyspieszonej x86:
float simple: 687.831624
double simple: 498.392324
float sin: 554.845588
double sin: 548.098432
mixed simple: 12308.813480
Jak widać nasz sin jest szybszy od tego standardowego. Wersji x64 na SIMD ani FPU nie mam, nie chce mi się z tym bawić. Visual C++ nie wspiera wstawek assemblerowych dla x64.

Po włączeniu inlinowania funkcji wersja przyspieszona jeszcze przyspiesza. Wyniki wersji x64 nie zmieniają się. Natomiast wersja x86 przyspiesza i to bardzo:
float simple: 269.746434
double simple: 198.085651
float sin: 774.067095
double sin: 732.198366
mixed simple: 2793.813390
Raytracer napisany C# jest podobnie wydajny w x64 jak w x86. Zgodnie z testami w C++ float są podobnie wydajne jak double. W C# jest podobnie.
float Test1(float a, float b)
{
    return a + 0.5f + b;
}

float Test2(float a, float b)
{
    return (float)Math.Sin(a) + 0.5f + b;
}

double Test3(double a, double b)
{
    return a + 0.5 + b;
}

double Test4(double a, double b)
{
    return Math.Sin(a) + 0.5 + b;
}

float Test5a(double a, float b)
{
    return (float)a + 0.5f + b;
}

float Test5b(double a, float b)
{
    return (float)(a + 0.5 + b);
}
Wyniki w x86:
float simple: 1112
double simple: 421
float sin: 510
double sin: 439
mixed simple A: 11150
mixed simple B: 9555
Wyniki w x64:
float simple: 620
double simple: 610
float sin: 577
double sin: 573
mixed simple A: 8028
mixed simple B: 9330
Testy typu simple są zdecydowanie wolniejsze, ale testy typu sin są porównywalne szybkością do C++. A mixed jest mixed, jak nazwa wskazuje.

Wniosek float nie jest szybsze od double, ale może być szybsze jeśli często odwołujemy się do pamięci.

Przejście z double na float

Dotychczas wszystkie obliczenia przeprowadzałem na double. Wraz z dodaniem korekcji gamma do sceny okazało się, że tektury najlepiej przechowywać po zdekodowaniu gamma. Z uwagi na błędy trzeba zrezygnować z bajtu jako pojemnika na składowe RGB i przejść na int. Ponieważ w obliczeniach posługuję się liczbami z zakresu 0 do 1, tak że 0 to czerń, a jeden to biel naturalnym pojemnikiem okazał się float. Przechowywanie tekstur na bajtach i ich dekodowanie w trakcie obliczeń okazało się wprowadzać znaczące spowolnienie renderowania.

Przechodząc na float w teksturach, uznałem, że nie będę już ich nigdy przechowywał w bajtach i posłałem w otchłań odpowiednie bajtowe pojemniki. A w zasadzie skonwertowałem je na float. Podobnie jak z pojemnikiem na teksturę postąpiłem z pojemnikiem na wartości (np. na mapę wypukłości).

Zużycie pamięci wzrosło czterokrotnie ale ciągle jest to akceptowalne dla mnie. Być w może w przyszłości trzeba będzie przywrócić bajtowe pojemniki by oszczędzać pamięć kosztem prędkości renderowania.

No niejako przy okazji przeniosłem wszelkie obliczenia z double na float.

Jeśli chodzi o jakość renderowanych scen to jest ona bardzo podobna (plus minus 1/255 głównie w niektórych miejscach). Jeśli chodzi o prędkość to nie zmieniła się ona. I to jest zastanawiające.

Kontrola dokładności obliczeń

Po przejściu z double na float na niektórych scenach pojawiły się dziwne błędy. Nie uśmiechało mi się ich debugować - scena z trzema przeźroczystymi sferami + fresnel + refrakcja + ścianki pudełka z tymi kulkami.

Ale okazało się, że podczas renderowania sceny w debug dostałem asserty. Dotyczyły one braku normalizacji wektorów, długość ustawianych wektorów nie była równa jeden. Oczywiście była bliska jeden, ale większa niż ustalona precyzja. Po zwiększeniu precyzji asserty znikły, a scena zrenderowała się bez tych błędów.

Ta sama precyzja kontroluje bowiem większość porównań.

Wniosek warto gęsto siać assertami w debug. Na pewno proste porównywanie liczb będzie generować błędne sceny. Warto używać tej samej precyzji wszędzie. Obecnie w kodzie posługuje się trzema stałymi kontrolującymi precyzję - do wszystkiego innego, do przylegania trójkątów, do sprawdzania kiedy jeden obiekt leży na tyle blisko drugiego, że można powiedzieć, że nie jest w jego cieniu.

Wszystkie trzy stałe są bardzo bliskie sobie i wydaje się, że w miarę czasu zrównają się. To znaczy obecne są niepoprawne, ale nie widać tych błędów na scenie.

Po wyjściu tego na jaw postanowiłem się przyjrzeć kodowi na intersekcję, gdzie jest wiele porównań, a tylko nieliczne robione są z precyzją. Poprzednio jak kod był na doublach zmieniając sposób porównania tylko pogarszałem sprawę. Ostatecznie zostawiając porównanie z precyzją tylko tam gdzie gdzie okazało się to niezbędne.

Samo porównywanie wygląda tak:
public static bool AlmostEqual(this float a_d1, float a_d2, float a_precision)
{
    float ad1 = Math.Abs(a_d1);
    float ad2 = Math.Abs(a_d2);

    if (ad1 < a_precision)
        a_d1 = 0;
    if (ad2 < a_precision)
        a_d2 = 0;

    if (a_d1 == a_d2)
        return true;

    return (Math.Abs(a_d1 - a_d2) / Math.Max(ad1, ad2)) < (a_precision * 4);
}

public static bool AlmostEqual(this float a_d1, float a_d2)
{
    return AlmostEqual(a_d1, a_d2, Constants.PRECISION);
}

public float SqrLen
{
    get
    {
        return X * X + Y * Y + Z * Z;
    }
}

public bool IsNormalized
{
    get
    {
        return SqrLen.AlmostEqual(1.0f);
    }
}

public Vector3 Dir
{
    get
    {
        return m_dir;
    }
    set
    {
        Debug.Assert(value.IsNormalized);
        m_dir = value;
    }
}

2011-08-28

Korekcja gamma

Pod tym terminem rozumiemy ogół zagadnień związanych z zapisywaniem wartości RGB w sposób nielinearny, ich kodowaniem do postaci nielinearnej oraz ich odtworzeniem na ekranie w postaci liniowej.

Historycznie do zapisu wartości RGB używało się bajtów. Jednakże 8 bitów to zbyt mało by liniowo oddać płynne przejścia jasności. Tak naprawdę potrzeba na to około 12 bitów. Z drugiej strony oko ludzkie jest bardziej czułe na zmiany ciemnych odcieni niż jasnych. Równomiernie rozłożone według ludzkiej percepcji wartości od ciemnego do jasnego w rzeczywistości odpowiadają funkcji potęgowej. Właśnie taką funkcję i jej odwrotność wykorzystuje się do kodowania i dekodowania informacji o jasności każdej składowej RGB.

Funkcja kodująca ma postać: $C'=C^{\gamma}$, gdzie $\gamma > 1$. Funkcja dekodująca ma postać: $C=C'^{\displaystyle{^{1}/_{\gamma}}}$. Mówiąc o kodowaniu i dekodowaniu gamma wartości wyjściowe jak i wyjściowe to liczby w zakresie 0-1, gdzie zero to czerń, zaś jeden to biel.
Z wykresu funkcji widać, że jedna czwarta jasności jest kodowana za pomocą połowy wartości.

Kodowanie i dekodowanie gamma stało się koniecznością wraz z pojawieniem się systemów zdolnych reprodukować kolor 24-bitowy (najpierw z palety). I ponieważ do dziś większość systemów używa do reprodukcji kolorów 24-bitów, one także kodują gammę.

sRGB w postaci 24-bitowej to nie tylko kompresja gamma, to także zapisanie jej w postaci ustalonej skali. Kiedy cały obraz jest ciemny, ale ciągle mamy wysoką czułość lub po prostu go ściemniliśmy i zapisaliśmy w postaci JPG. Informacja o dynamice kolorów ulega zmniejszeniu. To samo się dzieje kiedy zarejestrowane zostały obszary jasne (słońce) i ciemne cienie. Zakodowanie sceny o tak dużej dynamice jasności wymaga specjalnych algorytmów (Tone mapping), by nie utracić wrażenia dynamiki sceny. Scena już tak zapisana do JPG traci informację o początkowej jasności. Dlatego w programach graficznych powinniśmy całe przetwarzanie realizować i zapamiętywać na liczba dużej precyzji (najlepsza wydaje się liczba zmiennoprzecinkowa 32-bitowa). Istnieją także formaty pozwalające na zapis i odczyt takich obrazów. Słowo klucz to HDR (High-dynamic range). Ponieważ pamięć, łącza są coraz szybsze jest szansa, że w przyszłości większość materiałów będzie kodowała jasność na większej ilości bitów niż 8. HDR niekoniecznie oznacza brak kodowania gamma.

Gamma miała jeszcze jedną użyteczną cechę. Monitory CRT, ich część odpowiedzialna za emisję elektronów, także pracuje nieliniowo w funkcji bardzo zbliżonej do dekodującej gamma. W każdym razie kodowanie gamma nie ma związku z CRT, jakby się mogło wydawać, tylko ze zmieszczeniem pełnej informacji o jasności w bajcie. Współcześnie LCD także muszą uwzględnić fakt, że sygnał jest poddany kompresji gamma.

Kodowanie gamma obecnie jest dla nas prawie transparentne, dzięki standardowi jakie dla niego przyjęto. Typową wartością $\gamma$ jest 2.2. Używana była ona w kamerach, telewizji, systemach komputerowych. Wierna reprodukcja kolorów wymaga by wszystkie urządzenia po drodze używały takiej samej funkcji do kodowania gamma i jej odwrotności do dekodowania. Dotyczy to całego procesu przetwarzania który może wyglądać tak: kamera, mpeg, system, karta graficzna, złącze, monitor. Dlatego powstał standard sRGB, opisujący między innymi kodowanie gamma. Jest on współcześnie wykorzystywane tak powszechnie, że myśląc o kodowaniu RGB mamy na myśli sRGB.

Wspomniany wcześniej współczynnik gamma $\gamma=2.2$ został zaadaptowany przez sRGB. Obie metody kodowania są bardzo podobne do siebie ale nie identyczne. sRGB używa innego kodowania by uniknąć problemów numerycznych w okolicach zera i by zwiększyć dokładność numeryczną na niektórych systemach. Różnice są bardzo niewielkie w każdym bądź razie, stosując standardową funkcję potęgową nie popełnilibyśmy zbyt dużego błędu.

Funkcja kodują gamma ma postać:

$C=\begin{cases}
12.92C_{cases} & C_{LINEAR} \leqslant 0.0031308 \\
1.055{C_{LINEAR}}^{\displaystyle\frac{1}{2.4}} & C_{LINEAR} > 0.0031308
\end{cases}$
Dekodująca:

$C_{LINEAR}=\begin{cases}
\frac{C_{SRGB}}{12.92} & C_{SRGB} \leqslant 0.04045 \\
{(\frac{C_{SRGB}+0.055}{1.055}})^{2.4} & C_{SRGB} > 0.04045
\end{cases}$
Funkcje te w porównaniu do standardowej potęgowej o wykładniku 2.2 posiadają takie same wartości w zerze i w jeden. Funkcje sRGB w punktach 0.0031308 i 0.04045 są ciągłe, a ich pochodne są sobie równe.

Implementacja w kodzie:
public static class Gamma
{
    private const float A = 0.055f;
    private const float PHI = 12.92f;
    private const float GAMMA = 2.4f;

    private static float[] s_decode_table = new float[256];

    static Gamma()
    {
        for (int i = 0; i < 256; i++)
            s_decode_table[i] = DecodeGamma(i / 255.0f);
    }

    public static float EncodeGamma(float a_value)
    {
        if (a_value <= 0.0031308f)
            return a_value * PHI;
        else
            return (1 + A) * (float)Math.Pow(a_value, 1 / GAMMA) - A;
    }

    public static float DecodeGamma(byte a_value)
    {
        return s_decode_table[a_value];
    }

    public static float DecodeGamma(float a_value)
    {
        if (a_value <= 0.04045)
            return a_value / PHI;
        else
            return (float)Math.Pow((a_value + A) / (1 + A), GAMMA);
    }
}

Podczas wykonywania operacji arytmetycznych na wartościach kolorów (oświetlenie, resampling, filtrowanie tekstur) musimy pamiętać by te operacje wykonywać na kolorze liniowym. Inaczej do obrazu wynikowego wprowadzimy przekłamania kolorów. Ujawnią się one szczególnie na małych jasnych obszarach (resampling), gdzie ulegną one zwężeniu i podczas raytracingu na dużych obszar jednolicie podświetlonych.