BRDF and Spherical coordinate in ray tracing

Question

I developed a ray tracer that use standard phong/blinn phong lighting model. Now I'm modifying it to support physically based rendering, so I'm implementing various BRDF models. At the moment I'm focused on Oren-Nayar and Torrance-Sparrow model. Each one of these is based on spherical coordinates used to express incident wi and outgoing wo light direction.

My question is: which way is the right one convert wi and wo from cartesian coordinate to spherical coordinate?

I'm applying the standard formula reported here https://en.wikipedia.org/wiki/Spherical_coordinate_system#Coordinate_system_conversions but I'm not sure I'm doing the right thing, because my vector are not with tail at the origin of the cartesian coordinate system, but are centered on the intersection point of the ray with the object.

Here you can find my current implementation:

• 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

Can anyone help me giving an explanation of the correct way to convert the wi and wo vector from cartesian to spherical coordinate?

UPDATE

I copy here the relevant part of code:

spherical coordinate calculation

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;
}

UPDATE 2

After some search I found this implementation of Oren-Nayar BRDF.

http://content.gpwiki.org/index.php/D3DBook:(Lighting)_Oren-Nayar

In the implementation above theta for wi and wo is obtained simply doing arccos(wo.dotProduct(Normal)) and arccos(wi.dotProduct(Normal)). This seems reasonable to me, as we can use the normal of the intersection point as the zenith direction for our spherical coordinate system and do the calculation. The calculation of gamma = cos(phi_wi - phi_wo) do some sort of projection of wi and wo on what it calls "tangent space". Assuming everything is correct in this implementation, can i just use the formulas |View - Normal x (View.dotProduct(Normal))| and |Light - Normal x (Light.dotProduct(Normal))| to obtain the phi coordinate (instead of using arctan("something"))?

---

N.B.: I'm also reading "Physically based rendering: from theory to implementation", that is the book released to pbrt (http://www.pbrt.org). In this implementation there's some sort of change of coordinate system for the point of intersection (using partial derivates and the parametric coordinate of the surface) (I'm reading it now so what I'm saying could not be accurate). I want to find a straight way (maybe the one above in UPDATE 2 is what I'm searching for, if anyone could confirm it).

Answer 1

The polar coordinate system commonly used in BRDF definitions is set up relative to the surface being shaded, i.e. in tangent space. The $\theta$ angle measures how far you are from the surface normal while $\phi$ is the azimuth around the plane of the surface relative to some reference direction (which doesn't matter unless the BRDF is anisotropic). So if you want to convert to these coordinates, you have to first get your vector in tangent space (with the origin at the intersection point and two of the axes aligned with the surface), then apply the usual Cartesian-to-spherical transformation.

However, normally you can and should evaluate BRDFs without using polar coordinates, trigonometric functions, or angles at all, but just using vector math primitives such as dot products. This is usually more efficient and it's more robust, as you don't have to deal with angle wraparound, factors of pi, out-of-range arguments to inverse trig functions and so on. For instance, you probably know that the cosine of the angle between vectors can be obtained by just dotting the (normalized) vectors. The sine and tangent can be obtained through trig identities from the cosine (i.e. from the dot product).

Fabian Giesen wrote an article on this topic, Finish your derivations, please, that refers to the exact Oren-Nayar article you linked and gives an alternate, trig-free form:

Time to brush up your trigonometric identities. A particularly bad offender can be found here – the relevant section from the simplified shader is this:

float alpha = max( acos( dot( v, n ) ), acos( dot( l, n ) ) );
float beta  = min( acos( dot( v, n ) ), acos( dot( l, n ) ) );
C = sin(alpha) * tan(beta);

Ouch! If you use some trig identities and the fact that acos is monotonically decreasing over its domain, this reduces to:

float vdotn = dot(v, n);
float ldotn = dot(l, n);
C = sqrt((1.0 - vdotn*vdotn) * (1.0 - ldotn*ldotn))
  / max(vdotn, ldotn);

..and suddenly there’s no need to use a lookup texture anymore (and by the way, this has way higher accuracy too).

And note that this formulation only has one transcendental function (sqrt) in place of four (two acos, a sin, and a tan).