Importance sampling of the blinn-phong brdf in pathtracing

Question

tl;dr:

How do you importance sample the blinn-phong-brdf?

Recipe for importance sampling of the phong brdf as far as i understood it (pseudo-code):

Tuple<wi, pdf, brdfValue> samplePhongBRDF(HitInfo hi, Vector3 wo) {

    float sumKD = kd.r + kd.g + kd.b;    // compute diffuse "energy"
    float sumKS = ks.r + ks.g + ks.b;    // compute specular "energy" 
    float sum = sumKD + sumKS;           // compute total sum

    float specProbability = sumKs / sum; // compute probability for specular event
    
    float rn = generateRandomNumber();   // generate random number [0..1]
    Vector3 r = reflect(ray.direction, n);  // compute reflection vector r

    Vector3 wi;
    if rn < specularProbability {
        // sample specular -> generate phong distributed sample
        // around reflection vector r
        Vector3 newDirectionAroundR = phongDistOnHemisphere(...);

        // generate orthonormal basis around r (reflection vector as "normal")
        // -> tangent space
        TangentSpace rts = TangentSpace(r);    // sample around r

        // transform back from tangent space to world coordinates
        wi = rts.toWorldSpace(newDirectionAroundR));
    } else {
        // sample diffuse with cosine-distribution

        // generate orthonormal basis around n (normal)  
        TangentSpace ts(n);

        // generate cosine distributed samples on hemisphere
        Vector3 newDirectionTS = cosineDistOnHemisphere(...);

        // transform back from tangent space to world coordinates
        l = normalize(ts.toWorldSpace(newDirectionTS)); 
    }

    // compute cos theta for pdf
    float cosTheta = dot(l, n);

    // compute diffuse pdf
    float diffusePdfValue = cosTheta / pi;
    
    // compute specular pdf
    float specularPdfValue = 0.5f / pi * (shininess + 1.0f)
        * std::pow(std::max(dot(r, l), 0.0f), shininess);
    
    // do "lerp" linear interpolation between the two 
    float pdf = mix(diffusePdfValue, specularPdfValue, specProbability);

    // compute brdf value for wi (l) and wo (v)
    vec3 attenuation = brdf(wo, wi ...) * cosTheta;

    return Tuple(wi, pdf, brdfValue);
}

Question:

How do i adapt the above code for importance sampling of the phong-brdf to the blinn-phong-brdf according to this importance sampling note?

Thank you very much in advance!

Longer version:

I am following a tutorial on pathtracing (path tracing course) and i am implementing my own pathtracer on the way which is loosely based on the provided source code of the tutorial (source code). At the point of the tutorial where i am now, the author is showing how to importance sample different material brdfs like the lambert diffuse brdf and the phong brdf (Material Importance Sampling).

For the importance sampling of the specular part of the phong brdf a cosine lobe with specular exponent n around the reflected direction (R vector of the phong brdf) has to be sampled according to importance sampling note.

The source code for importance sampling of the phong brdf is the following (slightly modified by me):

class MaterialPhong : public Material
{
public:
    const Texture* k_d;
    const Texture* k_s;
    const Texture* s;

    MaterialPhong(const Texture* k_d, const Texture* k_s, const Texture* s) :
        k_d(k_d), k_s(k_s), s(s)
    {
    }

    // compute phong brdf
    // n  - normal at hitpoint
    // l  - light vector (wi)
    // v  - camera vector (wo)
    // kd - diffuse constant
    // ks - specular constant
    // shininess - phong exponent
    vec3 brdf(const vec3& n, const vec3& l, const vec3& v,
              const vec3& kd, const vec3& ks, float shininess) const
    {
        // compute lambert diffuse
        vec3 diffuse = kd / pi;
        
        // compute phong specular
        float cosRV = std::max(dot(reflect(-l, n), v), 0.0f);
        float normalization = (shininess + 2.0f) / (2.0f * pi);
        vec3 specular = ks * normalization * pow(cosRV, shininess);
        
        // return combined result
        return diffuse + specular;
    }

    // importance sampling of the phong brdf
    // ray  - ray which hits the surface
    // hr   - hitrecord contains information of the surface hitpoint (normal, texcoords, etc.)
    // prng - random number generator for sampling
    virtual ScatterRecord scatter(const Ray& ray, const HitRecord& hr, Prng& prng) const override
    {
        // if hit on backside return
        if (hr.backside)
            return ScatterRecord();

        // get diffuse, specular and shininess values from texture
        vec3 kd = k_d->value(hr.texcoord);
        vec3 ks = k_s->value(hr.texcoord);
        float shininess = s->value(hr.texcoord).x();

        // sample the brdf
        vec3 n = hr.normal;                     // n = normal at hitpoint
        vec3 v = -ray.direction;                // v = wo = vector to camera
        vec3 r = reflect(ray.direction, n);     // r = reflection vector for phong

        // use "energy" heuristic to decide wether a diffuse or 
        // specular sample should be taken
        float sumKd = kd.x() + kd.y() + kd.z(); // red (x), green (y) and blue (z) "wavelength"-energies added
        float sumKs = ks.x() + ks.y() + ks.z(); // to compute "total energy" from all wavelengths
        float sum = sumKd + sumKs;              // sum it, so the probabilities add to 1
        
        // compute
        float specularProbability = sumKs / sum;
        specularProbability = clamp(specularProbability, 0.1, 0.9); // keep the probability between 0.1 and 0.9
        
        
        // generate phong distributed sample light vector l (wi) around reflection vector
        vec3 l; // this is the new direction
        
        // sample brdf depending on the random number
        if (prng.in01() < specularProbability) {
            // act specular
            
            // generate phong distributed samples (phi, theta)  
            vec3 newDirectionAroundR = Sampler::phongWeightedOnHemisphere(shininess, prng.in01(), prng.in01());
            
            // generate orthonormal basis around r (reflection vector as "normal" for the onb) -> shading / tangent space
            TangentSpace rts = TangentSpace(r);
            
            // transform back from shading / tangent space to world coordinates
            l = normalize(rts.toWorldSpace(newDirectionAroundR));   
        } else {
            // act diffuse
            
            // generate orthonormal basis around n (normal)  
            TangentSpace ts(n);
            
            // generate cosine distributed samples on hemisphere
            vec3 newDirectionTS = Sampler::cosineWeightedOnHemisphere(prng.in01(), prng.in01());
            
            // transform back from tangent space to world coordinates
            l = normalize(ts.toWorldSpace(newDirectionTS)); 
        }

        // compute cos theta for pdf
        float cosTheta = dot(l, n);
        if (cosTheta <= 0.0f)
            return ScatterRecord();

        // compute diffuse pdf
        float diffusePdfValue = cosTheta / pi;
        
        // compute specular pdf
        float specularPdfValue = 0.5f / pi * (shininess + 1.0f)
            * std::pow(std::max(dot(r, l), 0.0f), shininess);
        
        // do "lerp" linear interpolation between the two depending on the specular probability
        float pdf = mix(diffusePdfValue, specularPdfValue, specularProbability);

        // compute brdf value for wi (l) and wo (v)
        vec3 attenuation = brdf(n, l, v, kd, ks, shininess) * cosTheta;
        
        // return result
        return ScatterRecord(l, pdf, attenuation);
    }
};

Now i want to write an equivalent material class for importance sampling of the blinn-phong brdf. But as i understand for the importance sampling of the specular part of the blinn phong brdf i would need to sample the cosine lobe around the half vector direction h of the specular part of the blinn phong brdf importance sampling note.

But since the half vector is computed by h = (v + l) (half vector = viewer / camera vector plus the light vector) and i don't have the light vector (l = wi), because i do the sampling because i want to get the light vector as a result, i can't compute h and therefore i can't generate samples around h. Because i can't generate samples around h, i can't compute l ... I am going round in circles here.

I tried to generate samples around the normal instead of the half vector using the following code but the result doesn't match the normal sampled blinn-phong version:

class MaterialBlinnPhong : public Material
{
public:
    const Texture* k_d;
    const Texture* k_s;
    const Texture* s;

    MaterialBlinnPhong(const Texture* k_d, const Texture* k_s, const Texture* s) :
        k_d(k_d), k_s(k_s), s(s)
    {
    }

    // compute blinn phong brdf
    // n  - normal at hitpoint
    // l  - light vector (wi)
    // v  - camera vector (wo)
    // kd - diffuse constant
    // ks - specular constant
    // shininess - blinn phong exponent
    vec3 brdf(const vec3& n, const vec3& l, const vec3& v,
            const vec3& kd, const vec3& ks, float shininess) const
    {
        // compute lambert diffuse
        vec3 diffuse = kd / pi;
        
        // compute phong specular
        vec3 h = (v + l).normalize();
        float cosNH = powf(std::max(dot(n, h), 0.0f), shininess);
        
        float normalization = ((shininess + 2.0f) * (shininess + 4.0f)) / (8.0f * M_PI * (powf(2.0, -shininess / 2.0f) + shininess));
        vec3 specular = ks * cosNH * normalization;
        
        // return combined result
        return diffuse + specular;
    }

    // importance sampling of the blinn phong brdf
    // ray  - ray which hits the surface
    // hr   - hitrecord contains information of the surface hitpoint (normal, texcoords, etc.)
    // prng - random number generator for sampling
    virtual ScatterRecord scatter(const Ray& ray, const HitRecord& hr, Prng& prng) const override
    {
        // if hit on backside return
        if (hr.backside)
            return ScatterRecord();

        // get diffuse, specular and shininess values from texture
        vec3 kd = k_d->value(hr.texcoord);
        vec3 ks = k_s->value(hr.texcoord);
        float shininess = s->value(hr.texcoord).x();

        // sample the brdf
        vec3 n = hr.normal;                     // n = normal at hitpoint
        vec3 v = -ray.direction;                // v = wo = vector to camera
        vec3 h;                                 // half vector of blinn phong

        // use "energy" heuristic to decide wether a diffuse or 
        // specular sample should be taken
        float sumKd = kd.x() + kd.y() + kd.z(); // red (x), green (y) and blue (z) "wavelength"-energies added
        float sumKs = ks.x() + ks.y() + ks.z(); // to compute "total energy" from all wavelengths
        float sum = sumKd + sumKs;              // sum it, so the probabilities add to 1
        
        // compute
        float specularProbability = sumKs / sum;
        specularProbability = clamp(specularProbability, 0.1, 0.9); // keep the probability between 0.1 and 0.9
        
        
        // generate phong distributed sample light vector l (wi) around reflection vector
        vec3 l; // this is the new direction
 
        if (prng.in01() < specularProbability) {
            // act specular
            
            // generate phong distributed samples (phi, theta according to  
            vec3 newDirectionAroundH = Sampler::phongWeightedOnHemisphere(shininess, prng.in01(), prng.in01());
            
            // generate orthonormal basis around normal n 
            // (but i would probably need to compute the onb around h here since i want to sample the cosine lobe 
            // around the reflected direction h and not n, but i can't compute h because i need l 
            // to compute h (= (l + v)))
            TangentSpace rts = TangentSpace(n);
            
            // transform back from shading / tangent space to world coordinates
            l = normalize(rts.toWorldSpace(newDirectionAroundH));
            
            // now i can compute h (but it is probably wrong)
            h = (l + v).normalize();    
        } else {
            // act diffuse
            
            // generate orthonormal basis around n (normal)  
            TangentSpace ts(n);
            
            // generate cosine distributed samples on hemisphere
            vec3 newDirectionTS = Sampler::cosineWeightedOnHemisphere(prng.in01(), prng.in01());
            
            // transform back from tangent space to world coordinates
            l = normalize(ts.toWorldSpace(newDirectionTS)); 
        }

        // compute cos theta for pdf
        float cosTheta = dot(l, n);
        if (cosTheta <= 0.0f)
            return ScatterRecord();

        // compute diffuse pdf
        float diffusePdfValue = cosTheta / pi;
        
        // compute specular pdf
        float specularPdfValue = 0.5f / pi * (shininess + 1.0f)
            * std::pow(std::max(dot(h, l), 0.0f), shininess);
        
        // do "lerp" linear interpolation between the two depending on the specular probability
        float pdf = mix(diffusePdfValue, specularPdfValue, specularProbability);

        // compute brdf value for wi (l) and wo (v)
        vec3 attenuation = brdf(n, l, v, kd, ks, shininess) * cosTheta;
        
        // return result
        return ScatterRecord(l, pdf, attenuation);
    }
};

So does anyone know how to importance sample the blinn-phong brdf in a pathtracer?

Thank you very much in advance.

Answer 1

Note that in the importance sampling from the importance sampling note provided by you (equation 48-59), we are not sampling the reflection direction but the half-vector direction. This is somewhat similar to the microfacet theory and you can refer to this answer post for more information on microfacet theory. You can think of this process as sampling the normal to be used for perfect specular reflection.

So the process is shown as follows:

• We know that Blinn-Phong BRDF (regardless of the normalization term $Z$ for proper BRDF definition) $$ \text{Blinn-Phong} = (N\cdot H)^\alpha \times Z $$ Therefore, $H$ can be importance sampled according to this $\cos^\alpha$ term. This is actually stated in the importance sampling note:

The Phong BRDF requires sampling a cosine lobe with specular exponent n around the reflected direction (or half-vector direction in case of BlinnPhong)

• We have our $H$, so we can actually use this $H$ to get the reflected direction, then using the $w_o$ we can get the $w_i$ with vector $H$ via the simple law of reflection.

Also, I believe you will find this paper useful: Lafortune & Willems, using the modified phong reflectance model for physically based rendering, 1994. Although old, it is the paper that put forward the Modified-Phong BRDF model and what it does is exactly: sample half vector and use the sampled half vector to compute incident direction via reflection.

If you are looking for some specific implementation... I have implemented this (in a simple way) in my own renderer. You can find the implementation in bxdf/brdf.py:line-191 (whole BRDF) and line-209 (sampling the BRDF), this is implemented in python so it might be easier to read than the implementation of mitsuba-Phong-BRDF (and I think mitsuba implements Phong but not Blinn-Phong).

# get the local direction for half vector
local_new_dir, pdf = mod_phong_hemisphere(self.mean[2])
# transform the local half vector to the world frame
normal, _ = delocalize_rotate(it.n_s, local_new_dir)
# reflect the incident (actually, outgoing direction in backward ray tracing)
ray_out_d = (-2 * normal * tm.dot(incid, normal) + incid).normalized()
# evaluate throughput
spec = self.eval_mod_phong(it, incid, ray_out_d)

At last I recommend the microfacet theory to you again, since you will find that all we do in microfacet BRDF is changing the distribution of the half vector to more complicated ones (like GGX, Beckmann). I hope I understand your question correctly and can be of help to you. Since it's been a while since I last learned this, feel free to point out anything wrong.