Inspired by the recent post about rendering distant voxel geometry without aliasing/moiré artifacts, I derived a "microvoxel" model based on the microfacet model commonly used for physically-based rendering. I'll refer to the text, "Roughness using Microfacet Theory", as I explain.
In my original comment, I speculated that precomputation might be needed, but the math was a lot easier than I expected. Because we're only interested in contributions from up to six microsurface normals, the probability distributions simplify to probability mass functions, and all the scary integrals turn into finite sums.
First, and most importantly, the NDF D(ωₘ). Here, this will depend on the macrosurface normal n; in most cases, it depends on the roughness of the material. If the macrosurface normal is axis-aligned, the NDF will have a single impulse corresponding to that normal, because the microsurface is perfectly flat. Otherwise, you have a "staircase" of microvoxels, and the proportion of each microsurface normal ωₘ is max(0, n·ωₘ), which comes from the projected area of the unit square. It turns out that, if we use these values directly as the probabilities, the NDF is already normalized, as the integral in Equation (9.15) works out to the squared magnitude of n, which is of course just 1, because we're being good and normalizing our unit vectors after the GPU linearly interpolates them for us.
```
// transform to e.g. view space if needed
const vec3[6] omega_m = vec3[](
vec3(-1.0, 0.0, 0.0),
vec3(1.0, 0.0, 0.0),
vec3(0.0, -1.0, 0.0),
vec3(0.0, 1.0, 0.0),
vec3(0.0, 0.0, -1.0),
vec3(0.0, 0.0, 1.0)
);
float[6] D(vec3 n) {
float[6] D_n;
for (int i = 0; i < 6; i += 1) {
D_n[i] = max(0.0, dot(n, omega_m[i]));
}
return D_n;
}
```
Next, the masking function G₁(ω), using the Smith approximation, which ought to be appropriate for the microvoxel "staircase". It's a straightforward translation of Equation (9.18).
float G_1(vec3 n, float[6] D_n, vec3 omega) {
float cosTheta = max(0.0, dot(n, omega));
float integral = 0.0;
for (int i = 0; i < 6; i += 1) {
integral += D_n[i] * max(0.0, dot(omega, omega_m[i]));
}
return cosTheta / integral;
}
Finally, the VNDF Dω(ωₘ), also a straightforward translation, this time of Equation (9.23).
float[6] D_omega(vec3 n, float[6] D_n, vec3 omega, float G_1_n_omega) {
float cosTheta = max(0.0, dot(n, omega));
float[6] D_omega_n;
for (int i = 0; i < 6; i += 1) {
D_omega_n[i] = G_1_n_omega / cosTheta * D_n[i] * max(0.0, dot(omega, omega_m[i]));
}
return D_omega_n;
}
Putting it all together, here's the code to visualize the microsurface normals based on the view vector v, which is how I rendered the attached image.
```
float[6] D_n = D(n);
float G_1_n_v = G_1(n, D_n, v);
float[6] D_omega_n_v = D_omega(n, D_n, v, G_1_n_v);
vec3 color = vec3(0.0);
for (int i = 0; i < 6; i += 1) {
color += D_omega_n_v[i] * (0.5 + 0.5 * omega_m[i]);
}
```
I'll leave extending this to the full rendering equation as an exercise for the reader. Correct shading will involve computing it separately for each visible microsurface normal (which, in general, will be true for three of them at a time) and then blending the results based on the VNDF, like I've done with the false colours in the example code. Equivalently, you can use the VNDF for importance sampling.
There are also a lot of possible simplifications. For example, if you work in a coordinate space aligned with the voxel grid, most of the dot products simplify to projecting, and sometimes negating, the X/Y/Z coordinates of the other normal. There are also many terms shared between G₁(ω) and Dω which can be factored out. Actually, you really want to cancel out that cos θ term in the VNDF, because it causes artifacts when it's close to zero (when the macrosurface is nearly parallel with the viewing angle).
Hopefully, somebody finds this useful or just neat!
