Co to jest próbkowanie według ważności? Każdy artykuł, który o nim czytam, wspomina o „PDF”, co to jest?

Z tego, co zbieram, ważność próbkowania jest techniką próbkowania tylko obszarów na półkuli, które mają większe znaczenie niż inne. Idealnie więc powinienem próbować promienie w kierunku źródeł światła, aby zmniejszyć hałas i zwiększyć prędkość. Ponadto, niektóre BRDF przy kątach wypasu nie mają niewielkiej różnicy w obliczeniach, więc użycie ważnego próbkowania, aby tego uniknąć, jest dobre?

Gdybym miał wdrożyć ważność próbkowania dla BRDF Cook-Torrance, jak mogę to zrobić?

Krótka odpowiedź:

Ważność próbkowania jest metodą zmniejszania wariancji w integracji Monte Carlo poprzez wybranie estymatora zbliżonego do kształtu rzeczywistej funkcji.

PDF jest skrótem od funkcji gęstości prawdopodobieństwa . A $pdf(x)$ daje prawdopodobieństwo wygenerowania losowej próbki o wartości $x$ .

Długa odpowiedź:

Na początek przejrzyjmy, czym jest integracja Monte Carlo i jak wygląda matematycznie.

Integracja Monte Carlo to technika szacowania wartości całki. Zwykle jest używany, gdy nie ma zamkniętego rozwiązania całki. To wygląda tak:

$$\int f(x) \, dx \approx \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{pdf(x_{i})}$$

W języku angielskim oznacza to, że można przybliżać całkę, uśredniając kolejne losowe próbki z funkcji. Gdy $N$ staje się duże, przybliżenie zbliża się coraz bardziej do rozwiązania. $pdf(x_{i})$ reprezentuje funkcję gęstości prawdopodobieństwa każdej próbki losowej.

Zróbmy przykład: Oblicz wartość całki $I$ .

$$I = \int_{0}^{2\pi} e^{-x} \sin(x) dx$$

Użyjmy integracji Monte Carlo:

$$I \approx \frac{1}{N} \sum_{i=1}^N \frac{e^{-x} \sin(x_i)}{pdf(x_{i})}$$

Prostym programem pythonowym do obliczenia tego jest:

import random
import math

N = 200000
TwoPi = 2.0 * math.pi

sum = 0.0

for i in range(N):
    x = random.uniform(0, TwoPi)

    fx = math.exp(-x) * math.sin(x)
    pdf = 1 / (TwoPi - 0.0)

    sum += fx / pdf

I = (1 / N) * sum
print(I)

Jeśli uruchomimy program, otrzymamy $I = 0.4986941$

Stosując separację według części, możemy uzyskać dokładne rozwiązanie:

$$I = \frac{1}{2} (1 − e−2π) = 0.4990663$$

Zauważysz, że rozwiązanie Monte Carlo nie jest całkiem poprawne. Jest tak, ponieważ jest to szacunek. To powiedziawszy, gdy przechodzi w nieskończoność, szacunki powinny być coraz bliżej poprawnej odpowiedzi. Już przy niektóre przebiegi są prawie identyczne z poprawną odpowiedzią.$N$$N = 2000$

Uwaga na temat pliku PDF: W tym prostym przykładzie zawsze bierzemy jednolitą losową próbkę. Jednolita losowa próbka oznacza, że ​​każda próbka ma dokładnie takie samo prawdopodobieństwo wyboru. Próbkujemy w zakresie więc$[0, 2\pi]$$pdf(x) = 1 / (2\pi - 0)$

Ważność próbkowania polega na tym, że próbkowanie nie jest jednolite. Zamiast tego staramy się wybierać więcej próbek, które mają duży wpływ na wynik (ważne), a mniej próbek, które tylko nieznacznie przyczyniają się do wyniku (mniej ważne). Stąd nazwa, ważność próbkowania.

Jeśli wybierzesz funkcję próbkowania, której pdf bardzo ściśle pasuje do kształtu , możesz znacznie zmniejszyć wariancję, co oznacza, że ​​możesz pobrać mniej próbek. Jeśli jednak wybierzesz funkcję próbkowania, której wartość jest bardzo różna od , możesz zwiększyć wariancję. Zobacz zdjęcie poniżej: Zdjęcie z rozprawy doktorskiej Wojciecha Jarosza Załącznik A$f$$f$

Jednym z przykładów ważnego próbkowania w śledzeniu ścieżki jest sposób wyboru kierunku promienia po uderzeniu w powierzchnię. Jeśli powierzchnia nie jest idealnie błyszcząca (tj. Lustro lub szkło), promień wychodzący może znajdować się w dowolnym miejscu na półkuli.

My mogliśmy równomiernie próbki półkulę, aby wygenerować nowy ray. Możemy jednak wykorzystać fakt, że równanie renderowania zawiera w sobie współczynnik cosinus:

$$L_{\text{o}}(p, \omega_{\text{o}}) = L_{e}(p, \omega_{\text{o}}) + \int_{\Omega} f(p, \omega_{\text{i}}, \omega_{\text{o}}) L_{\text{i}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | d\omega_{\text{i}}$$

W szczególności wiemy, że wszelkie promienie na horyzoncie będą silnie osłabione (w szczególności ). Tak więc promienie generowane w pobliżu horyzontu nie przyczynią się zbytnio do końcowej wartości.$\cos(x)$

Aby temu przeciwdziałać, używamy próbkowania według ważności. Jeśli generujemy promienie zgodnie z półkulą ważoną cosinusem, zapewniamy, że więcej promieni generuje się znacznie powyżej horyzontu, a mniej w pobliżu horyzontu. Obniży to wariancję i zmniejszy hałas.

W twoim przypadku podałeś, że będziesz korzystać z BRDF opartego na mikropacetach Cook-Torrance. Powszechną formą jest:

$$f(p, \omega_{\text{i}}, \omega_{\text{o}}) = \frac{F(\omega_{\text{i}}, h) G(\omega_{\text{i}}, \omega_{\text{o}}, h) D(h)}{4 \cos(\theta_{i}) \cos(\theta_{o})}$$

gdzie

$$F(\omega_{\text{i}}, h) = \text{Fresnel function} \\ G(\omega_{\text{i}}, \omega_{\text{o}}, h) = \text{Geometry Masking and Shadowing function} \\ D(h) = \text{Normal Distribution Function}$$

Blog „Notka faceta z grafiki” ma doskonały opis próbkowania BRDF-ów Cook-Torrance. Odsyłam cię do jego postu na blogu . To powiedziawszy, postaram się stworzyć krótki przegląd poniżej:

NDF jest generalnie dominującą częścią Cook-Torrance BRDF, więc jeśli chcemy uzyskać ważność próbki, powinniśmy spróbować na podstawie NDF.

Cook-Torrance nie określa konkretnego NDF do użycia; mamy swobodę wyboru tego, co pasuje do naszych upodobań. To powiedziawszy, istnieje kilka popularnych NDF:

• GGX
• Beckmann
• Blinn

Każdy NDF ma własną formułę, dlatego każdy musi być próbkowany inaczej. Pokażę tylko ostateczną funkcję próbkowania dla każdego. Jeśli chcesz zobaczyć, jak formuła jest uzyskiwana, zobacz wpis na blogu.

GGX jest zdefiniowany jako:

$$D_{GGX}(m) = \frac{\alpha^2}{\pi((\alpha^2-1) \cos^2(\theta) + 1)^2}$$

Aby pobrać próbkę sferycznego kąta współrzędnych , możemy użyć wzoru:$\theta$

$$\theta = \arccos \left( \sqrt{\frac{\alpha^2}{\xi_1 (\alpha^2 - 1) + 1}} \right)$$

gdzie jest jednolitą zmienną losową.$\xi$

Zakładamy, że NDF jest izotropowy, więc możemy próbkować jednolicie:$\phi$

$$\phi = \xi_{2}$$

Beckmann jest zdefiniowany jako:

$$D_{Beckmann}(m) = \frac{1}{\pi \alpha^2\cos^4(\theta)} e^{-\frac{\tan^2(\theta)}{\alpha^2}}$$

Które można próbkować za pomocą:

$$\theta = \arccos \left(\sqrt{\frac{1}{1 = \alpha^2 \ln(1 - \xi_1)}} \right) \\ \phi = \xi_2$$

Wreszcie, Blinn jest zdefiniowany jako:

$$D_{Blinn}(m) = \frac{\alpha + 2}{2 \pi} (\cos(\theta))^{\alpha}$$

Które można próbkować za pomocą:

$$\theta = \arccos \left(\frac{1}{\xi_{1}^{\alpha + 1}} \right) \\ \phi = \xi_2$$

Realizacja w praktyce

Spójrzmy na podstawowy znacznik ścieżki wstecz:

void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene->Camera.CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);

    // Bounce the ray around the scene
    for (uint bounces = 0; bounces < 10; ++bounces) {
        m_scene->Intersect(ray);

        // The ray missed. Return the background color
        if (ray.geomID == RTC_INVALID_GEOMETRY_ID) {
            color += throughput * float3(0.846f, 0.933f, 0.949f);
            break;
        }

        // We hit an object

        // Fetch the material
        Material *material = m_scene->GetMaterial(ray.geomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene->GetLight(ray.geomID);

        // If we hit a light, add the emmisive light
        if (light != nullptr) {
            color += throughput * light->Le();
        }

        float3 normal = normalize(ray.Ng);
        float3 wo = normalize(-ray.dir);
        float3 surfacePos = ray.org + ray.dir * ray.tfar;

        // Get the new ray direction
        // Choose the direction based on the material
        float3 wi = material->Sample(wo, normal, sampler);
        float pdf = material->Pdf(wi, normal);

        // Accumulate the brdf attenuation
        throughput = throughput * material->Eval(wi, wo, normal) / pdf;


        // Shoot a new ray

        // Set the origin at the intersection point
        ray.org = surfacePos;

        // Reset the other ray properties
        ray.dir = wi;
        ray.tnear = 0.001f;
        ray.tfar = embree::inf;
        ray.geomID = RTC_INVALID_GEOMETRY_ID;
        ray.primID = RTC_INVALID_GEOMETRY_ID;
        ray.instID = RTC_INVALID_GEOMETRY_ID;
        ray.mask = 0xFFFFFFFF;
        ray.time = 0.0f;
    }

    m_scene->Camera.FrameBuffer.SplatPixel(x, y, color);
}

TO ZNACZY. odbijamy się wokół sceny, gromadząc w miarę upływu czasu kolor i tłumienie światła. Przy każdym odbiciu musimy wybrać nowy kierunek promienia. Jak wspomniano powyżej, mogłyby równomiernie próbki półkulę, aby wygenerować nowy ray. Jednak kod jest mądrzejszy; znaczenie pobiera próbki nowego kierunku na podstawie BRDF. (Uwaga: To jest kierunek wprowadzania, ponieważ jesteśmy znacznikiem ścieżki do tyłu)

// Get the new ray direction
// Choose the direction based on the material
float3 wi = material->Sample(wo, normal, sampler);
float pdf = material->Pdf(wi, normal);

Które można zaimplementować jako:

void LambertBRDF::Sample(float3 outputDirection, float3 normal, UniformSampler *sampler) {
    float rand = sampler->NextFloat();
    float r = std::sqrtf(rand);
    float theta = sampler->NextFloat() * 2.0f * M_PI;

    float x = r * std::cosf(theta);
    float y = r * std::sinf(theta);

    // Project z up to the unit hemisphere
    float z = std::sqrtf(1.0f - x * x - y * y);

    return normalize(TransformToWorld(x, y, z, normal));
}

float3a TransformToWorld(float x, float y, float z, float3a &normal) {
    // Find an axis that is not parallel to normal
    float3a majorAxis;
    if (abs(normal.x) < 0.57735026919f /* 1 / sqrt(3) */) {
        majorAxis = float3a(1, 0, 0);
    } else if (abs(normal.y) < 0.57735026919f /* 1 / sqrt(3) */) {
        majorAxis = float3a(0, 1, 0);
    } else {
        majorAxis = float3a(0, 0, 1);
    }

    // Use majorAxis to create a coordinate system relative to world space
    float3a u = normalize(cross(normal, majorAxis));
    float3a v = cross(normal, u);
    float3a w = normal;


    // Transform from local coordinates to world coordinates
    return u * x +
           v * y +
           w * z;
}

float LambertBRDF::Pdf(float3 inputDirection, float3 normal) {
    return dot(inputDirection, normal) * M_1_PI;
}

Po próbkowaniu inputDirection („wi” w kodzie), używamy go do obliczenia wartości BRDF. Następnie dzielimy przez pdf zgodnie ze wzorem Monte Carlo:

// Accumulate the brdf attenuation
throughput = throughput * material->Eval(wi, wo, normal) / pdf;

Gdzie Eval () jest samą funkcją BRDF (Lambert, Blinn-Phong, Cook-Torrance itp.):

float3 LambertBRDF::Eval(float3 inputDirection, float3 outputDirection, float3 normal) const override {
    return m_albedo * M_1_PI * dot(inputDirection, normal);
}

Jeśli masz funkcję 1D i chcesz zintegrować tę funkcję od powiedzmy 0 do 1, jednym ze sposobów przeprowadzenia tej integracji jest pobranie N losowych próbek w zakresie [0, 1], oszacuj dla każdego próbkuj i oblicz średnią próbek. Mówi się jednak, że ta „naiwna” integracja Monte Carlo „powoli zbiega się”, tzn. Potrzeba dużej liczby próbek, aby zbliżyć się do prawdy gruntu, szczególnie jeśli funkcja ma wysokie częstotliwości.$f(x)$$f(x)$

Dzięki ważnemu próbkowaniu, zamiast pobierania N losowych próbek w zakresie [0, 1], pobiera się więcej próbek w „ważnych” regionach które najbardziej przyczyniają się do końcowego wyniku. Ponieważ jednak próbujesz odchylenie w kierunku ważnych obszarów funkcji, próbki te muszą być ważone mniej, aby przeciwdziałać odchyleniu, czyli tam, gdzie pojawia się PDF (funkcja gęstości prawdopodobieństwa). PDF określa prawdopodobieństwo próbki w danej pozycji i służy do obliczenia średniej ważonej próbek poprzez podzielenie każdej próbki przez wartość PDF w każdej pozycji próbki.$f(x)$

W przypadku pobierania próbek o znaczeniu ważności Cook-Torrance powszechną praktyką jest dystrybucja próbek w oparciu o normalną funkcję rozkładu NDF. Jeśli NDF jest już znormalizowany, może służyć bezpośrednio jako PDF, co jest wygodne, ponieważ anuluje to określenie z oceny BRDF. Jedyne, co musisz zrobić, to rozdzielić przykładowe pozycje na podstawie PDF i ocenić BRDF bez terminu NDF, tj. I oblicz średnią wyników próbki pomnożoną przez kąt bryłowy domeny, którą integrujesz (np. dla półkuli).

$$f=\frac{FG}{\pi(n\cdot\omega_i)(n\cdot\omega_o)}$$ $2\pi$

W przypadku NDF należy obliczyć funkcję skumulowanego rozkładu pliku PDF, aby przekonwertować równomiernie rozmieszczoną pozycję próbki na pozycję próbki ważoną w formacie PDF. W przypadku izotropowego NDF upraszcza to funkcję 1D ze względu na symetrię funkcji. Aby uzyskać więcej informacji na temat pochodnej CDF, możesz sprawdzić ten stary artykuł o klejnotach GPU .