Troubleshoot half vector sampling from a distribution (cook-torrence, trowbridge-reitz, etc)
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).