BRDF i współrzędna sferyczna w śledzeniu promieni

Opracowałem ray tracer, który wykorzystuje standardowy model oświetlenia phong / blinn phong. Teraz modyfikuję go, aby obsługiwał renderowanie oparte na fizyce, dlatego wdrażam różne modele BRDF. W tej chwili skupiam się na modelu Oren-Nayar i Torrance-Sparrow. Każda z nich oparta jest na współrzędnych sferycznych używanych do wyrażenia zdarzenia wi i wychodzącego wo kierunku światła.

Moje pytanie brzmi: w którą stronę należy przekonwertować wi i wo ze współrzędnych kartezjańskich na współrzędne sferyczne?

Stosuję standardową formułę podaną tutaj https://en.wikipedia.org/wiki/Spherical_coordinate_system#Coordinate_system_conversions, ale nie jestem pewien, czy robię coś dobrego, ponieważ mój wektor nie jest z ogonem u źródła kartezjański układ współrzędnych, ale są wyśrodkowane na punkcie przecięcia promienia z obiektem.

Tutaj możesz znaleźć moją obecną implementację:

• https://github.com/chicio/Multispectral-Ray-tracing/tree/brdf/RayTracing/RayTracer/Objects/BRDF

• https://github.com/chicio/Multispectral-Ray-tracing/blob/brdf/RayTracing/RayTracer/Math/Vector3D.cpp

Czy ktoś może mi pomóc w wyjaśnieniu prawidłowego sposobu konwersji wektora wi i wo ze współrzędnych kartezjańskich na sferyczne?

AKTUALIZACJA

Tutaj kopiuję odpowiednią część kodu:

obliczanie współrzędnych sferycznych

float Vector3D::sphericalTheta() const {

    float sphericalTheta = acosf(Utils::clamp(y, -1.f, 1.f));

    return sphericalTheta;
}

float Vector3D::sphericalPhi() const {

    float phi = atan2f(z, x);

    return (phi < 0.f) ? phi + 2.f * M_PI : phi;
}

Oren Nayar

OrenNayar::OrenNayar(Spectrum<constant::spectrumSamples> reflectanceSpectrum, float degree) : reflectanceSpectrum{reflectanceSpectrum} {

    float sigma = Utils::degreeToRadian(degree);
    float sigmaPowerTwo = sigma * sigma;

    A = 1.0f - (sigmaPowerTwo / 2.0f * (sigmaPowerTwo + 0.33f));
    B = 0.45f * sigmaPowerTwo / (sigmaPowerTwo + 0.09f);
};

Spectrum<constant::spectrumSamples> OrenNayar::f(const Vector3D& wi, const Vector3D& wo, const Intersection* intersection) const {

    float thetaI = wi.sphericalTheta();
    float phiI = wi.sphericalPhi();

    float thetaO = wo.sphericalTheta();
    float phiO = wo.sphericalPhi();

    float alpha = std::fmaxf(thetaI, thetaO);
    float beta = std::fminf(thetaI, thetaO);

    Spectrum<constant::spectrumSamples> orenNayar = reflectanceSpectrum * constant::inversePi * (A + B * std::fmaxf(0, cosf(phiI - phiO) * sinf(alpha) * tanf(beta)));

    return orenNayar;
}

Torrance-Sparrow

float TorranceSparrow::G(const Vector3D& wi, const Vector3D& wo, const Vector3D& wh, const Intersection* intersection) const {

    Vector3D normal = intersection->normal;
    normal.normalize();

    float normalDotWh = fabsf(normal.dot(wh));
    float normalDotWo = fabsf(normal.dot(wo));
    float normalDotWi = fabsf(normal.dot(wi));
    float woDotWh = fabsf(wo.dot(wh));

    float G = fminf(1.0f, std::fminf((2.0f * normalDotWh * normalDotWo)/woDotWh, (2.0f * normalDotWh * normalDotWi)/woDotWh));

    return G;
}

float TorranceSparrow::D(const Vector3D& wh, const Intersection* intersection) const {

    Vector3D normal = intersection->normal;
    normal.normalize();

    float cosThetaH = fabsf(wh.dot(normal));

    float Dd = (exponent + 2) * constant::inverseTwoPi * powf(cosThetaH, exponent);

    return Dd;
}

Spectrum<constant::spectrumSamples> TorranceSparrow::f(const Vector3D& wi, const Vector3D& wo, const Intersection* intersection) const {

    Vector3D normal = intersection->normal;
    normal.normalize();

    float thetaI = wi.sphericalTheta();
    float thetaO = wo.sphericalTheta();

    float cosThetaO = fabsf(cosf(thetaO));
    float cosThetaI = fabsf(cosf(thetaI));

    if(cosThetaI == 0 || cosThetaO == 0) {

        return reflectanceSpectrum * 0.0f;
    }

    Vector3D wh = (wi + wo);
    wh.normalize();

    float cosThetaH = wi.dot(wh);

    float F = Fresnel::dieletricFresnel(cosThetaH, refractiveIndex);
    float g = G(wi, wo, wh, intersection);
    float d = D(wh, intersection);

    printf("f %f g %f d %f \n", F, g, d);
    printf("result %f \n", ((d * g * F) / (4.0f * cosThetaI * cosThetaO)));

    Spectrum<constant::spectrumSamples> torranceSparrow = reflectanceSpectrum * ((d * g * F) / (4.0f * cosThetaI * cosThetaO));

    return torranceSparrow;
}

AKTUALIZACJA 2

Po kilku poszukiwaniach znalazłem tę implementację Oren-Nayar BRDF .

W powyższej implementacji uzyskano theta dla wi i wo po prostu wykonując arccos (wo.dotProduct (Normal)) i arccos (wi.dotProduct (Normal)). Wydaje mi się to rozsądne, ponieważ możemy wykorzystać normę punktu przecięcia jako kierunek zenitu dla naszego sferycznego układu współrzędnych i wykonać obliczenia. Obliczenie gamma = cos (phi_wi - phi_wo) dokonuje pewnego rzutu wi i wo na to, co nazywa „przestrzenią styczną”. Zakładając, że wszystko jest prawidłowe w tej implementacji, czy mogę po prostu użyć formuł | Widok - Normalny x (View.dotProduct (Normalny)) | i | Light - Normal x (Light.dotProduct (Normal)) | uzyskać współrzędną phi (zamiast używać arctan („coś”))?

W rzeczywistości lepiej jest nie używać współrzędnych sferycznych (ani żadnych innych kątów) do implementacji BRDF, ale raczej pracować prosto w kartezjańskim układzie współrzędnych i stosować cosinus kąta między wektorami, który, jak wiadomo, jest iloczynem prostych kropek między wektorami jednostek. Jest to zarówno bardziej niezawodne, jak i wydajne.

W przypadku Oren-Nayar możesz myśleć, że musisz użyć kątów (ze względu na min / max kątów), ale możesz po prostu zaimplementować BRDF prosto w przestrzeni kartezjańskiej: https://fgiesen.wordpress.com/2010/10/21 / dokończ-swoje-pochodne-proszę

W przypadku mikropacetek Torrance-Sparrow lub Cook-Torrance BRDF nie musisz także używać współrzędnych sferycznych. W tych BRDF kąt jest przekazywany do funkcji trygonometrycznej (zwykle cosinusowej) w kategoriach D / F / G i mianowniku BRDF, dzięki czemu można używać iloczynu prostego lub tożsamości trygonometrycznej bez przechodzenia przez współrzędne sferyczne.

Możesz określić układ współrzędnych, biorąc pod uwagę normalny N i inny wektor. Wybierzemy wi. Tak więc każdy wektor, który ma ten sam kierunek co wi podczas rzutowania na płaszczyznę styczną, będzie miał azymut 0

Najpierw projektujemy wi na płaszczyźnie stycznej: (zakładając, że wi jest już znormalizowany)

wit = normalize(wi - N * dot(wi, N))

teraz możemy zrobić to samo z wo:

wot = normalize(wo - N * dot(wo, N))

Teraz dowcip i wot leżą zarówno na płaszczyźnie prostopadłej do N, jak i stycznej do punktu przecięcia.

Teraz możemy obliczyć kąt między nimi:

azimuth = arcos ( dot(wit, wot) )

Co tak naprawdę jest azymutem wot w odniesieniu do dowcipu, gdy jest rzutowany na płaszczyznę styczną.

Jeśli znasz punkt przecięcia i punkt początkowy, czy nie będzie to po prostu kwestia odjęcia jednego od drugiego, aby uzyskać wynik tak, jakby był z początku?

Jeśli nie wierzysz w wynik i chcesz przejść długą drogę, możesz również uzyskać transformację obrotu, aby przejść z jednego punktu do drugiego za pomocą macierzy LookAt, a następnie rozłożyć ją, aby uzyskać składnik rotacyjny. Jeśli chcesz, możesz także uzyskać kwaternion.

Wyniki są równe. Dowód jest nieco długi, ale nie skomplikowany i należy do czytelnika.