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.

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"))?

Answer 1

It's actually better not to use spherical coordinates (or any angles for that matter) to implement BRDF's, but rather work straight in Cartesian coordinate system and use cosine of the angle between vectors, which is plain dot product between unit vectors as you know. This is both more robust and efficient.

For Oren-Nayar you may think you have to use angles (due to min/max of the angles), but you can simply implement the BRDF straight in Cartesian space: https://fgiesen.wordpress.com/2010/10/21/finish-your-derivations-please

For Torrance-Sparrow or Cook-Torrance microfacet BRDF's you don't need to use spherical coordinates either. In these BRDF's the angle is passed to a trigonometric (usually cosine) function in D/F/G terms & the BRDF denominator, so you can use dot product straight or trigonometric identities without going through spherical coordinates.

Answer 2

You can specify a coordinate system given the normal N and another vector. We'll choose wi. So any vector that has the same direction as wi when projected onto the tangent plane will have an azimuth of 0

First, we project wi on the tangent plane: (assuming wi is already normalized)

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

now, we can do the same with wo:

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

Now, wit and wot lie both on a plane which is orthogonal to N, and tangent to the point of intersection.

We can compute the angle between the two now:

azimuth = arcos ( dot(wit, wot) )

Which is really the azimuth of wot with respect to wit when projected on the tangent plane.

Answer 3

If you know the intersection point, and the point of origin, wouldn't it just be a question of subtracting one from the other so you get the result as if it were from the origin?

If you don't believe the result, and want to get there through the long way, you can also get the rotation transform to get from one point to another via a LookAt matrix, and then decomposing it to get the rotational component. You can also get a quaternion from it if you want.

The results are equal. The proof is a bit long, but not complicated, and is left up to the reader.