Empirical confimation of diffuse ray hack

Benjamin Keinert  sent a nice note in regards to an earlier post to quickly get "lambertianish" rays.  He empirically confirmed it is exactly lambertian.  Cool, and thanks Benjamin I always believe such demonstrations more convincing that proofs which can have subtle errors (fine, yes, I am an engineer).  I include his full email with his permission.

I think I found a simple informal "engineer's style proof" that it is indeed lambertian (assuming the sample rejection method yields a uniform distribution on a sphere, which it should).
Instead of using a rejection method I picked the uniform distribution on a sphere using spherical fibonacci point sets [2] and constructed a cosine hemisphere sampling variant of it.
Without the loss of generality it should be sufficient to show that the mapping is lambertian for a single normal  (0,0,1) - given uniformly distributed points on a sphere and rotational invariance.

Sorry, rapid prototyping code, oldschool OpenGL, PHI = (sqrt(5.0)*0.5 + 0.5):

// PDF: 1/(4*PI)
float3 uniformSampleSphereSF(float i, float n) {
    float phi = 2*PI*(i/PHI);
    float cosTheta = 1 - (2*i+1)/n;
    float sinTheta = sqrt(1 - cosTheta*cosTheta);
    return float3(cos(phi)*sinTheta, sin(phi)*sinTheta, cosTheta);
}

// PDF: cos(theta)/PI
float3 cosineSampleHemisphereSF(float i, float n) {
    float phi = 2*PI*(i/PHI);
    float cosTheta = sqrt(1 - (i+0.5)/n);
    float sinTheta = sqrt(1 - cosTheta*cosTheta);
    return float3(cos(phi)*sinTheta, sin(phi)*sinTheta, cosTheta);
}

[...]
void test() {
    [...]
    // Enable additive blending etc.
    [...]
    uint n = 1024;
    glBegin(GL_POINTS);
    for (uint i = 0; i < n; ++i) {
        glColor4f(0,1,0,1); // Green
        float3 p = normalize(uniformSampleSphereSF(i, n) + float3(0,0,1));
        glVertex3fv(&p[0]);

        glColor4f(1,0,0,1); // Red
        float3 q = cosineSampleHemisphereSF(i, n);
        glVertex3fv(&q[0]);
        // Additive blending => Yellow == "good"
    }
    glEnd();
}

This little function results in the attached image (orthogonal projection of cosine distributed points on a hemisphere -> uniformly distributed points on a circle).
With some more effort one can show that normalize(uniformSampleSphereSF(i, n) + float3(0,0,1)) = cosineSampleHemisphereSF(i, n) - instead of using additive blending.


[1] http://psgraphics.blogspot.de/2014/09/random-diffuse-rays.html
[2] Spherical Fibonacci Point Sets for Illumination Integrals, Marques et. al.