Troubleshoot half vector sampling from a distribution (cook-torrence, trowbridge-reitz, etc)

Question

I am learning ray tracing and the mathematics behind it. I have a working monte-carlo raytracer I am experimenting on. I have gotten past pure diffuse BRDFS, area lights, acceleration structures, etc, and now I'm working on properly sampling the next ray direction from an anisotropic distribution. I believe I have interpreted the D,G,F factors correctly, and have working tests that verify those function's outputs.

I am having trouble with drawing a "wi" sample from the half vector distribution function. The purpose of the function is to give you a randomized half vector in shading space based on the "incoming" wo ray, which is the eye ray, and a roughness parameter. As I graph out the function outputs, it seems like the half-vector sampling function begins to rotate the spherical coordinates about the wrong axis (I know this makes no sense, but I'm just observing the graphed results).

Knowing that situation, I was hoping some people could post some pointers to good academic articles, open source code examples, or other educational material that could help me investigate and educate myself on the behavior of functions that generate half-vectors from distributions. I have found plenty of the original source papers from Cook-Torrence, etc that explain their algorithms, but I cannot find any that explain in plain language HOW to generate the half vector for the eye ray to reflect about.

Note: I have also extensively explored PBRT and its source code, all of which has been incredibly enlightening, and I have seen their sample_wi function, i even went so far as to directly copy and translate it to webgl for troubleshooting my functions. I know i must be doing or interpreting something wrong - I just can't find enough material out there to help me figure out what.

Here is a visualization of only the half vector normals. They have been put into a range of [0,1] for viewing in the RGB spectrum. I have also boosted the contrast of the image to make the error more easily visible. You can clearly see that instead of randomly orienting themselves about the surface normal, they are instead correlating themselves radially about the camera eye-vector.

UPDATE:
here is a webgl "proof" i put together to isolate and illustrate my problem. What should happen is the color should be fairly blended together, interpolated rather nicely. Instead I get this obvious hard shift in the distribution throughout the 4 quadrants: http://rdtests.ml3ds-test.com/microfacettesting.html

The functions referenced by main function here can all be viewed directly in the source code of the page I linked to:

        vec3 finalColor = vec3(0.);

        // create a changing "wo" vector across screen space (simulates perspective ray change)
        vec2 screenUv = gl_FragCoord.xy / u_resolution;
        // fire rays into positive z (forward) space from "camera"s
        float z = 1.;

        vec3 rd = normalize(vec3(screenUv, z));

        // remap 0,1 to -1,1
        vec3 wo = (rd * 2.) - 1.;
        // wo is outward facing
        wo *= -wo;

        // convert wo to spherical coordinates
        vec3 Nt, Nb;
        // simulate hitting an infinite plane, facing us like a wall
        vec3 normal = vec3(0., 0.,-1.);
        CoordinateSystem(normal, Nt, Nb);

        // put wo into the shading space coordinate system of the "wall" normal
        vec3 shadingWo = WorldToShading(wo, Nt, Nb, normal);

        // sample half vector distribution from trowbridge-reitz
        // make random input a constant to allow error to be more clearly seen
        vec3 wh = Dist_TR_Sample_wh(shadingWo, vec2(.2,.1), .5, .5);

        // Reflect wo about the half vector normal
        vec3 wi = Reflect(shadingWo, wh);

        // put wi back into cartesian coordinate system
        vec3 wiWorld = ShadingToWorld(wi, Nt, Nb, normal);

        // now remap wi to a viewable color: 0,1
        finalColor = normalize((wiWorld + 1.) / 2.);

        gl_FragColor = vec4(finalColor, 1.);

Answer 1

The general idea for sampling half vector based distributions is that you generate $H$ and then compute $w_i$ by reflecting $w_o$ about $H$. This is so $H$ will be the half vector of your $w_i$ and $w_o$ pair. It is standard reflection: $$w_i = -w_o + 2(w_o\cdot H)H $$

How you generate $H$ depends on the specific distribution. Generally, it is done in polar coordinates, with the angle from your distribution's center being picked using some specific function. Then the azimuth will be a uniform distribution, unless your distribution is anisotropic, in which case it gets more complicated.

For example, with Cook-Torrance, it goes something like: $$\tan \theta = m \sqrt{ - \log \xi_1 }$$ $$\phi = \xi_2 2\pi$$

Then, given an orthonormal basis made of vectors $N$, $a$, $b$:

$$H = a\sin\theta\cos\phi + b\sin\theta\sin\phi + N\cos\theta$$

As for problems in your code, replacing this:

vec3 rd = normalize(vec3(screenUv, z));

// remap 0,1 to -1,1
vec3 wo = (rd * 2.) - 1.;
// wo is outward facing
wo *= -wo;

with:

vec3 wo = normalize(vec3(screenUv*2.0-1.0, -z));

Gives me a much nicer plot. Your remapping was incorrect as it was done on a vector which includes z and is no longer in 0,1 for x and y because of the normalization. This explains why the center was not in the center of the frame and the whole thing was leaning. Removing wo *= -wo; (a transformation which makes no sense) and using -z instead fixes what I think you were calling the hard transition in the 4 quadrants. That was caused by squaring the components of the vector.

The remaining glitch in the center appears to be numerical issues. Fixing that will require reworking the code to better handle some cases (likely where you have the 0.0008 constant for a sin and the 0.9999 for a cos).