<p>Note that in the importance sampling from the <a href="https://www.mathematik.uni-marburg.de/%7Ethormae/lectures/graphics1/code/ImportanceSampling/importance_sampling_notes.pdf" rel="nofollow noreferrer">importance sampling note provided by you</a> (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 <a href="https://computergraphics.stackexchange.com/questions/13590/importance-sampling-in-a-path-tracer/13591#13591">this answer post</a> for more information on microfacet theory. You can think of this process as sampling the normal to be used for perfect specular reflection.</p>
<p>So the process is shown as follows:</p>
<ul>
<li>We know that Blinn-Phong BRDF (regardless of the normalization term <span class="math-container">$Z$</span> for proper BRDF definition)
<span class="math-container">$$
\text{Blinn-Phong} = (N\cdot H)^\alpha \times Z
$$</span>
Therefore, <span class="math-container">$H$</span> can be importance sampled according to this <span class="math-container">$\cos^\alpha$</span> term. This is actually stated in the importance sampling note:</li>
</ul>
<blockquote>
<p>The Phong BRDF requires sampling a cosine lobe with specular exponent n around the reflected direction <strong>(or half-vector direction in case of BlinnPhong)</strong></p>
</blockquote>
<ul>
<li>We have our <span class="math-container">$H$</span>, so we can actually use this <span class="math-container">$H$</span> to get the reflected direction, then using the <span class="math-container">$w_o$</span> we can get the <span class="math-container">$w_i$</span> with vector <span class="math-container">$H$</span> via the simple law of reflection.</li>
</ul>
<p>Also, I believe you will find this paper useful: <a href="https://www.researchgate.net/profile/Yves-Willems-2/publication/2361953_Using_the_Modied_Phong_Reflectance_Model_for_Physically_Based_Rendering/links/562e17af08ae22b170352271/Using-the-Modied-Phong-Reflectance-Model-for-Physically-Based-Rendering.pdf" rel="nofollow noreferrer">Lafortune &amp; Willems, using the modified phong reflectance model for physically based rendering, 1994</a>. Although old, it is the paper that put forward the <strong>Modified-Phong BRDF model</strong> and what it does is exactly: sample half vector and use the sampled half vector to compute incident direction via reflection.</p>
<p>If you are looking for some specific implementation... I have implemented this (in a simple way) in <a href="https://github.com/Enigmatisms/AdaPT" rel="nofollow noreferrer">my own renderer</a>. You can find the implementation in <a href="https://github.com/Enigmatisms/AdaPT/blob/master/bxdf/brdf.py" rel="nofollow noreferrer">bxdf/brdf.py:line-191 (whole BRDF) and line-209 (sampling the BRDF)</a>, this is implemented in python so it might be easier to read than the implementation of <a href="https://github.com/mitsuba-renderer/mitsuba/blob/master/src/bsdfs/phong.cpp" rel="nofollow noreferrer">mitsuba-Phong-BRDF</a> (and I think mitsuba implements Phong but not Blinn-Phong).</p>
<pre class="lang-python prettyprint-override"><code># 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)
</code></pre>
<p>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.</p>

CC BY-SA 4.0
2023-08-17T08:27:21.203Z
I looked at the code of your implementation of the blinn-phong / phong brdfs in your pathtracer and i am impressed by the results in the image gallery and that you did this with python!

But when i looked correctly at your source code, your implementation of blinn-phong starts on line **164**. Beginning from line 191 on starts your implementation of the modified phong model, which already works in my pathtracer.
Berti Krüger
CC BY-SA 4.0
2023-08-17T08:39:53.903Z
And your code seems to do the same as mine:

- compute the reflection vector vec3 r = reflect(-ray.direction, n) = reflect(wo, n)
- use random number to decide, whether to do diffuse or specular sample
- if diffuse, sample with cosine-distributed sample -> evaluate lambert brdf
- if specular,  generate phong distributed sample ***and use reflection vector r as normal for the orthonormal base / tangent space**, so you **generate a distribution around the reflection vector r** 
- convert to world space -> you get wi
- you have wo (-ray direction) and wi -> you can evaluate the phong brdf
Berti Krüger
CC BY-SA 4.0
2023-08-17T08:50:17.990Z
But for blinn phong, as far as i understand, i would have to generate a phong distribution around the reflection vector h instead of r. Then i could compute wi which is the result i need at the end. But to get h, i need to compute it first. And h = wo + wi. But i don't have wi at this point. I am doing this to get wi at the end.

In your implementation of the importance sampling of blinn-phong you don't use a phong-distribution, but a cosine-distribution which is better than a uniform-distribution. But the importance sampling note mentions using the phong-distribution. How should that work?
Berti Krüger
CC BY-SA 4.0
2023-08-17T09:00:46.877Z
And thanks for the hints on looking at the microfacet theory and brdfs. I already looked at microfacet brdfs and implemented a plastic like material based on the cook-torrance model. But i have a lot of wavefront obj models which use the classic non-physical blinn-phong brdf model and i want to render them as close as possible in my pathtracer to this model. But i have nowhere found a renderer / pathtracer which does importance sampling of blinn-phong. All (mitsuba, tungsten etc.) use modified phong not blinn-phong. Yours seem to be the only one but using cosine-distribution.
Berti Krüger
CC BY-SA 4.0
2023-08-17T09:10:02.023Z
When i looked correctly this is your code for importance sampling blinn-phong:

> def sample_phong(self, it:ti.template(), incid: vec3):
>   local_new_dir, pdf = cosine_hemisphere()
>   ray_out_d, _ = delocalize_rotate(it.n_s, local_new_dir)
>   spec = self.eval_phong(it, incid, ray_out_d)
>   return ray_out_d, spec, pdf

So your are using a cosine-distribution on hemisphere for importance sampling which is clearly better than uniform. But from the notes i linked, better would be phong-distribution which should work too and which you use in your phong code.
Berti Krüger
CC BY-SA 4.0
2023-08-17T09:16:08.457Z
But from the notes i linked, better would be phong-distribution which should work too and which you use in your phong code.

It is a bit strange, that none of the well-known pathtracers (mitsuba, tungsten, renderbitch) has the blinn-phong-brdf but only the phong brdf.

If nothing works, i will do it like you do and use the cosine-distribution for importance sampling of the blinn-phong-brdf. But it should somehow work with blinn-phong too.
Berti Krüger
CC BY-SA 4.0
2023-08-17T14:56:02.013Z
Note that, $w_o$, $H$ and $w_i$ form a linear relationship, i.e once you know any two among these three vectors you can get the remaining one. So all you should do in Blinn-Phong BRDF is to generate direction for half vector, that way you get your reflected incident direction. Also, the code starting from line 164 is actually for Phong model since it does not require us to compute the half vector $H$. The comment in my code should be corrected.
Enigmatisms
CC BY-SA 4.0
2023-08-17T15:00:58.713Z
Yeah, mitsuba and pbrt they all have no Blinn-Phong model implemented. This may come with a reason: microfacet model is a better physically-based model for glossy specular appearance. Since we use Phong/Blinn-Phong to model glossy reflection, yet Phong and Blinn-Phong only differ slightly, and also, the intuition behind Blinn-Phong  model (sampling the half vector as 'microfacet normal') can be covered by the microfacet models. Therefore, Blinn-Phong model might overlap with microfacet models.
Enigmatisms
CC BY-SA 4.0
2023-08-17T15:04:42.367Z
Quote from your second comment: "and use reflection vector r as normal for the orthonormal base / tangent space, so you generate a distribution around the reflection vector r". I didn't use the reflection vector as normal, I will **sample the half vector as normal, transform it to the world space** and then calculate the reflected direction. Since we can, of course, sample $w_i$ first and then we calculate $H$ but this is difficult mathematically. So we do it the other way around, we sample $H$ to get $w_i$. As I said, $w_o$ is known so either $w_i$ or $H$ will get you the other one.
Enigmatisms
CC BY-SA 4.0
2023-08-17T15:05:59.770Z
BTW, the code is not based on python though it is *written in python*. It's based on Taichi lang. Python would be too slow for a renderer.
Enigmatisms
CC BY-SA 4.0
2023-08-16T16:42:52.453Z
2023-08-16T16:42:52.453Z
Enigmatisms
<h1>tl;dr:</h1>
<h2><strong>How do you importance sample the blinn-phong-brdf?</strong></h2>
<p><strong>Recipe for importance sampling of the phong brdf as far as i understood it (pseudo-code):</strong></p>
<pre><code>Tuple&lt;wi, pdf, brdfValue&gt; samplePhongBRDF(HitInfo hi, Vector3 wo) {

    float sumKD = kd.r + kd.g + kd.b;    // compute diffuse &quot;energy&quot;
    float sumKS = ks.r + ks.g + ks.b;    // compute specular &quot;energy&quot; 
    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 &lt; specularProbability {
        // sample specular -&gt; generate phong distributed sample
        // around reflection vector r
        Vector3 newDirectionAroundR = phongDistOnHemisphere(...);

        // generate orthonormal basis around r (reflection vector as &quot;normal&quot;)
        // -&gt; 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 &quot;lerp&quot; 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);
}
</code></pre>
<h2><strong>Question:</strong></h2>
<p><em>How do i adapt the above code for importance sampling of the <strong>phong-brdf</strong> to the <strong>blinn-phong-brdf</strong> according to this <a href="https://www.mathematik.uni-marburg.de/%7Ethormae/lectures/graphics1/code/ImportanceSampling/importance_sampling_notes.pdf" rel="nofollow noreferrer">importance sampling note</a>?</em></p>
<p>Thank you very much in advance!</p>
<h2><strong>Longer version:</strong></h2>
<p>I am following a tutorial on pathtracing (<a href="https://marlam.de/path-tracing/course/" rel="nofollow noreferrer">path tracing course</a>) and i am implementing my own pathtracer on the way which is loosely based on the provided source code of the tutorial (<a href="https://github.com/marlam/path-tracing-mirror/tree/main/src/09" rel="nofollow noreferrer">source code</a>). 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 (<a href="https://marlam.de/path-tracing/path-tracing-07.pdf" rel="nofollow noreferrer">Material Importance Sampling</a>).</p>
<p>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 <a href="https://www.mathematik.uni-marburg.de/%7Ethormae/lectures/graphics1/code/ImportanceSampling/importance_sampling_notes.pdf" rel="nofollow noreferrer">importance sampling note</a>.</p>
<p>The source code for importance sampling of the phong brdf is the following (slightly modified by me):</p>
<pre><code>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&amp; n, const vec3&amp; l, const vec3&amp; v,
              const vec3&amp; kd, const vec3&amp; 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&amp; ray, const HitRecord&amp; hr, Prng&amp; 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-&gt;value(hr.texcoord);
        vec3 ks = k_s-&gt;value(hr.texcoord);
        float shininess = s-&gt;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 &quot;energy&quot; 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) &quot;wavelength&quot;-energies added
        float sumKs = ks.x() + ks.y() + ks.z(); // to compute &quot;total energy&quot; 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() &lt; 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 &quot;normal&quot; for the onb) -&gt; 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 &lt;= 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 &quot;lerp&quot; 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);
    }
};
</code></pre>
<p>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 <a href="https://www.mathematik.uni-marburg.de/%7Ethormae/lectures/graphics1/code/ImportanceSampling/importance_sampling_notes.pdf" rel="nofollow noreferrer">importance sampling note</a>.</p>
<p>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.</p>
<p>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:</p>
<pre><code>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&amp; n, const vec3&amp; l, const vec3&amp; v,
            const vec3&amp; kd, const vec3&amp; 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&amp; ray, const HitRecord&amp; hr, Prng&amp; 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-&gt;value(hr.texcoord);
        vec3 ks = k_s-&gt;value(hr.texcoord);
        float shininess = s-&gt;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 &quot;energy&quot; 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) &quot;wavelength&quot;-energies added
        float sumKs = ks.x() + ks.y() + ks.z(); // to compute &quot;total energy&quot; 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() &lt; 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 &lt;= 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 &quot;lerp&quot; 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);
    }
};
</code></pre>
<p>So does anyone know how to importance sample the blinn-phong brdf in a pathtracer?</p>
<p>Thank you very much in advance.</p>

<raytracing><pathtracing><brdf><importance-sampling>
Importance sampling of the blinn-phong brdf in pathtracing
Berti Krüger
# Importance sampling of the blinn-phong brdf in pathtracing

(Asked by: username_0 on 2023-08-15)

# 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](https://www.mathematik.uni-marburg.de/%7Ethormae/lectures/graphics1/code/ImportanceSampling/importance_sampling_notes.pdf)?*

Thank you very much in advance!

## **Longer version:**

I am following a tutorial on pathtracing ([path tracing course](https://marlam.de/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](https://github.com/marlam/path-tracing-mirror/tree/main/src/09)). 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](https://marlam.de/path-tracing/path-tracing-07.pdf)).

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](https://www.mathematik.uni-marburg.de/%7Ethormae/lectures/graphics1/code/ImportanceSampling/importance_sampling_notes.pdf).

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](https://www.mathematik.uni-marburg.de/%7Ethormae/lectures/graphics1/code/ImportanceSampling/importance_sampling_notes.pdf).

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.

**Tags:** raytracing, pathtracing, brdf, importance-sampling


---

## 1 answers

### Answer by username_1 on 2023-08-16. Score: 0

Note that in the importance sampling from the [importance sampling note provided by you](https://www.mathematik.uni-marburg.de/%7Ethormae/lectures/graphics1/code/ImportanceSampling/importance_sampling_notes.pdf) (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](https://computergraphics.stackexchange.com/questions/13590/importance-sampling-in-a-path-tracer/13591#13591) 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](https://www.researchgate.net/profile/Yves-Willems-2/publication/2361953_Using_the_Modied_Phong_Reflectance_Model_for_Physically_Based_Rendering/links/562e17af08ae22b170352271/Using-the-Modied-Phong-Reflectance-Model-for-Physically-Based-Rendering.pdf). 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](https://github.com/username_1/AdaPT). You can find the implementation in [bxdf/brdf.py:line-191 (whole BRDF) and line-209 (sampling the BRDF)](https://github.com/username_1/AdaPT/blob/master/bxdf/brdf.py), this is implemented in python so it might be easier to read than the implementation of [mitsuba-Phong-BRDF](https://github.com/mitsuba-renderer/mitsuba/blob/master/src/bsdfs/phong.cpp) (and I think mitsuba implements Phong but not Blinn-Phong).

```python
# 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.

> **Comments:**
> **username_0**: I looked at the code of your implementation of the blinn-phong / phong brdfs in your pathtracer and i am impressed by the results in the image gallery and that you did this with python!  But when i looked correctly at your source code, your implementation of blinn-phong starts on line **164**. Beginning from line 191 on starts your implementation of the modified phong model, which already works in my pathtracer.  
> **username_0**: And your code seems to do the same as mine:  - compute the reflection vector vec3 r = reflect(-ray.direction, n) = reflect(wo, n) - use random number to decide, whether to do diffuse or specular sample - if diffuse, sample with cosine-distributed sample -> evaluate lambert brdf - if specular,  generate phong distributed sample ***and use reflection vector r as normal for the orthonormal base / tangent space**, so you **generate a distribution around the reflection vector r**  - convert to world space -> you get wi - you have wo (-ray direction) and wi -> you can evaluate the phong brdf  
> **username_0**: But for blinn phong, as far as i understand, i would have to generate a phong distribution around the reflection vector h instead of r. Then i could compute wi which is the result i need at the end. But to get h, i need to compute it first. And h = wo + wi. But i don't have wi at this point. I am doing this to get wi at the end.  In your implementation of the importance sampling of blinn-phong you don't use a phong-distribution, but a cosine-distribution which is better than a uniform-distribution. But the importance sampling note mentions using the phong-distribution. How should that work?  
> **username_0**: And thanks for the hints on looking at the microfacet theory and brdfs. I already looked at microfacet brdfs and implemented a plastic like material based on the cook-torrance model. But i have a lot of wavefront obj models which use the classic non-physical blinn-phong brdf model and i want to render them as close as possible in my pathtracer to this model. But i have nowhere found a renderer / pathtracer which does importance sampling of blinn-phong. All (mitsuba, tungsten etc.) use modified phong not blinn-phong. Yours seem to be the only one but using cosine-distribution.  
> **username_0**: When i looked correctly this is your code for importance sampling blinn-phong:  > def sample_phong(self, it:ti.template(), incid: vec3): >   local_new_dir, pdf = cosine_hemisphere() >   ray_out_d, _ = delocalize_rotate(it.n_s, local_new_dir) >   spec = self.eval_phong(it, incid, ray_out_d) >   return ray_out_d, spec, pdf  So your are using a cosine-distribution on hemisphere for importance sampling which is clearly better than uniform. But from the notes i linked, better would be phong-distribution which should work too and which you use in your phong code.  
> **username_0**: But from the notes i linked, better would be phong-distribution which should work too and which you use in your phong code.  It is a bit strange, that none of the well-known pathtracers (mitsuba, tungsten, renderbitch) has the blinn-phong-brdf but only the phong brdf.  If nothing works, i will do it like you do and use the cosine-distribution for importance sampling of the blinn-phong-brdf. But it should somehow work with blinn-phong too.  
> **username_1**: Note that, $w_o$, $H$ and $w_i$ form a linear relationship, i.e once you know any two among these three vectors you can get the remaining one. So all you should do in Blinn-Phong BRDF is to generate direction for half vector, that way you get your reflected incident direction. Also, the code starting from line 164 is actually for Phong model since it does not require us to compute the half vector $H$. The comment in my code should be corrected.  
> **username_1**: Yeah, mitsuba and pbrt they all have no Blinn-Phong model implemented. This may come with a reason: microfacet model is a better physically-based model for glossy specular appearance. Since we use Phong/Blinn-Phong to model glossy reflection, yet Phong and Blinn-Phong only differ slightly, and also, the intuition behind Blinn-Phong  model (sampling the half vector as 'microfacet normal') can be covered by the microfacet models. Therefore, Blinn-Phong model might overlap with microfacet models.  
> **username_1**: Quote from your second comment: "and use reflection vector r as normal for the orthonormal base / tangent space, so you generate a distribution around the reflection vector r". I didn't use the reflection vector as normal, I will **sample the half vector as normal, transform it to the world space** and then calculate the reflected direction. Since we can, of course, sample $w_i$ first and then we calculate $H$ but this is difficult mathematically. So we do it the other way around, we sample $H$ to get $w_i$. As I said, $w_o$ is known so either $w_i$ or $H$ will get you the other one.  
> **username_1**: BTW, the code is not based on python though it is *written in python*. It's based on Taichi lang. Python would be too slow for a renderer.  

---

