Energy-Normalizing Blinn–Phong Specular
Most people typically ignore energy normalization for BRDFs, in favor of energy conservation (that is, energy may be lost (to unaccounted factors like absorption), but at least it is never gained). Most people are satisfied with that because it models the shadowing and masking effects in microfacet theory. Unfortunately, it completely ignores multiple scattering. The correct answer, then, is somewhere in-between the extremely lossy models typically used and an idealized model that loses none.
The Blinn–Phong specular BRDF is widely used in teaching applications for its simplicity, but no energy normalization is known. With much the same motivation, I provided one for the Phong BRDFs, and I want to do the same for the Blinn–Phong too.
This page collects some of my efforts toward deriving a normalization term.
Statement of Energy Normalization
The (scalar variant of the) (modified) Blinn–Phong specular BRDF (at a point \(\vec{x}\) with normal \(\vec{N}\), and with energy-normalization factor \(\alpha\) and specular exponent \(n\)) is:
\[ f_r(\vec{\omega}_i,\vec{\omega}_o) ~:=~ \alpha \cdot (\vec{H} \cdot \vec{N})^n = \alpha \cdot \left( \frac{ \vec{\omega}_i + \vec{\omega}_o }{ \left \| \vec{\omega}_i + \vec{\omega}_o \right \| }\cdot\vec{N} \right)^n \]We want to solve for \(\alpha\).
First, energy-normalization means the following[1]:
\[ \iint_{\Omega_{\vec{N}}} f_r(\vec{\omega}_i,\vec{\omega}_o) \left( \vec{N} \cdot \vec{\omega}_i \right)~d\vec{\omega}_i = 1 \]We can now substitute in and rearrange:
\begin{align*} \iint_{\Omega_{\vec{N}}} \alpha \cdot (\vec{H} \cdot \vec{N})^n \left( \vec{N} \cdot \vec{\omega}_i \right)~d \vec{\omega}_i = 1\\ \iint_{\Omega_{\vec{N}}} \left( \frac{ \vec{\omega}_i + \vec{\omega}_o }{ \left \| \vec{\omega}_i + \vec{\omega}_o \right \| }\cdot\vec{N} \right)^n \left( \vec{N} \cdot \vec{\omega}_i \right)~d\vec{\omega}_i = \frac{1}{\alpha} \end{align*}Substitution for Length
To compute \(\left \| \vec{\omega}_i + \vec{\omega}_o \right \|\) is actually surprisingly simple. We can consider both vectors to lie in a plane, separated by angle \(\gamma\). Without loss of generality, for purposes of computing the length, we can move to 2D: suppose \(\vec{\omega}_o:=\vecinline{0,1}\) and \(\vec{\omega}_i:=\vecinline{\sin(\gamma),\cos(\gamma)}\). Their sum is \(\vecinline{\sin(\gamma), 1+\cos(\gamma) }\) and the length of that sum is (Pythagorean theorem):
\[ \left \| \vec{\omega}_i + \vec{\omega}_o \right \| = \sqrt{ \sin^2(\gamma) + (1+\cos(\gamma))^2 } = \sqrt{ 2 + 2 \cos(\gamma)} = 2 \cos \left(\frac{\gamma}{2}\right) \]Since \(\gamma\) is just \(\arccos(\vec{\omega}_i \cdot \vec{\omega}_o)\), we have:
\[ \left \| \vec{\omega}_i + \vec{\omega}_o \right \| = 2 \cos \left(\frac{\gamma}{2}\right) = 2 \cos \left(\frac{1}{2} \arccos(\vec{\omega}_i \cdot \vec{\omega}_o) \right) = 2 \sqrt{\frac{\vec{\omega}_i \cdot \vec{\omega}_o + 1}{2}} = \sqrt{2} \sqrt{\vec{\omega}_i \cdot \vec{\omega}_o + 1} \]Substitute this into the integral:
\[ \frac{1}{\alpha} = \iint_{\Omega_{\vec{N}}} \left( \frac{ \vec{\omega}_i + \vec{\omega}_o }{ \sqrt{2} \sqrt{\vec{\omega}_i \cdot \vec{\omega}_o + 1} }\cdot\vec{N} \right)^n \left( \vec{N} \cdot \vec{\omega}_i \right)~d \vec{\omega}_i \]Spherical Coordinates
One obvious first step is to (re-)define \(\vec{\omega}_o\), \(\vec{\omega}_i\), and \(\vec{N}\) in spherical coordinates:
\[ \vec{\omega}_o = \begin{bmatrix} \sin\phi_o \\ 0 \\ \cos\phi_o \end{bmatrix} \text{,}\hspace{1cm} \vec{\omega}_i = \begin{bmatrix} \cos(\Delta\theta) \cdot \sin\phi_i \\ \sin(\Delta\theta) \cdot \sin\phi_i \\ \hspace{19mm} \cos\phi_i \end{bmatrix} \text{,}\hspace{1cm} \vec{N} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \] \[ \frac{1}{\alpha} = \int_0^{2\pi} \int_0^{\pi/2} \left( \frac{ \cos\phi_o + \cos\phi_i }{ \sqrt{2} \sqrt{ \sin\phi_o \sin\phi_i \cos(\Delta\theta) + \cos\phi_o \cos\phi_i + 1 } } \right)^n \cos\phi_i \sin\phi_i ~d\phi_i ~d\Delta\theta \]As a shorthand (this gets messy fast), I like to define \(c_i:=\cos\phi_i\), \(s_o:=\sin\phi_o\), etc., giving the following:
\begin{equation} \frac{2^{n/2}}{\alpha} = \int_0^{2\pi} \int_0^{\pi/2} \left( \frac{ c_o + c_i }{\sqrt{ s_o s_i \cos(\Delta\theta) + c_o c_i + 1 } } \right)^n c_i s_i ~d\phi_i ~d\Delta\theta \end{equation}Spherical Coordinates: Integrating \(\Delta\theta\)
Here, we try integrating over \(\Delta\theta\) first.
Some rearranging produces:
\begin{align*} \frac{2^{n/2}}{\alpha} &= \int_0^{\pi/2} \left(\frac{ c_o + c_i }{ \sqrt{ s_o s_i } }\right)^n c_i s_i \int_0^{2\pi} \left( \frac{ 1 }{ \sqrt{\smash[b]{ \cos(\Delta\theta) + \underbrace{(c_o c_i + 1)/(s_o s_i)}_{A_i} }} } \right)^n ~d\Delta\theta ~d\phi_i \\[2em]&= \int_0^{\pi/2} \left(\frac{ c_o + c_i }{ \sqrt{s_o s_i} }\right)^n c_i s_i \int_0^{2\pi} \frac{ d\Delta\theta }{ ( A_i + \cos(\Delta\theta) )^{n/2} } ~d\phi_i \end{align*}The inner integral above is surprisingly difficult[2], even if you assume \(m:=n/2\) is an integer. However, with an immense amount of effort[3], you can get it in terms of the hypergeometric function:
\[ = \int_0^{\pi/2} \left(\frac{ c_o + c_i }{ \sqrt{s_o s_i} }\right)^n \cdot c_i s_i \cdot \left[ \frac{2\pi}{A_i^m} \cdot {_2}F_1\!\left( \frac{m}{2}, \frac{m+1}{2} ; 1 ; \frac{1}{A_i^2} \right) \right] \cdot d\phi_i \]From there, we can expand out its definition, extract constants, bring the integral inside, and simplify:
\begin{align*}&= 2\pi \int_0^{\pi/2} \left(\frac{ c_o + c_i }{ \sqrt{ c_o c_i + 1 } }\right)^n \cdot c_i s_i \cdot {_2}F_1\!\left( \frac{n}{4}, \left(\frac{n}{4}\!+\!\frac{1}{2}\right) ; 1 ; \left(\frac{s_o s_i}{c_o c_i + 1}\right)^{\!2} \right) \cdot d\phi_i \\ &= 2\pi \int_0^{\pi/2} \left(\frac{ c_o + c_i }{ \sqrt{ c_o c_i + 1 } }\right)^n \cdot c_i s_i \cdot \left[ \sum_{k=0}^\infty \frac{(\tfrac{n}{4})_k ~(\tfrac{n}{4}\!+\!\tfrac{1}{2})_k}{(1)_k ~k!} \left(\frac{s_o s_i}{c_o c_i + 1}\right)^{\!2k} \right] \cdot d\phi_i \\ &= 2\pi \int_0^{\pi/2} \left(\frac{ c_o + c_i }{ \sqrt{ c_o c_i + 1 } }\right)^n \cdot c_i s_i \cdot \left[ \sum_{k=0}^\infty \frac{ \Gamma(\tfrac{n}{4}\!+\!k) \Gamma(\tfrac{n}{4}\!+\!\tfrac{1}{2}\!+\!k) \Gamma(1) }{ \Gamma(\tfrac{n}{4}) \Gamma(\tfrac{n}{4}\!+\!\tfrac{1}{2}) \Gamma(1+k) \Gamma(k+1) } \left( \frac{s_o s_i}{c_o c_i + 1} \right)^{\!2k} \right] \cdot d\phi_i \\ &= 2\pi \sum_{k=0}^\infty \left[ \frac{ \Gamma(\tfrac{n}{4}\!+\!k) \Gamma(\tfrac{n}{4}\!+\!\tfrac{1}{2}\!+\!k) }{ \Gamma(\tfrac{n}{4}) \Gamma(\tfrac{n}{4}\!+\!\tfrac{1}{2}) \Gamma^2(1+k) } \int_0^{\pi/2} \left(\frac{ (c_o + c_i)^2 }{ c_o c_i + 1 }\right)^{n/2} \left( \frac{s_o s_i}{c_o c_i + 1} \right)^{\!2k} c_i s_i \cdot d\phi_i \right] \\&= \frac{2\pi}{ \Gamma(\tfrac{n}{4}) \Gamma(\tfrac{n}{4}\!+\!\tfrac{1}{2}) } \sum_{k=0}^\infty \left[ \frac{ \Gamma(\tfrac{n}{4}\!+\!k) \Gamma(\tfrac{n}{4}\!+\!\tfrac{1}{2}\!+\!k) }{ \Gamma^2(1+k) } s_o^{2k} \int_0^{\pi/2} \frac{ (c_o + c_i)^n ~s_i^{2k+1} ~c_i }{ (c_o c_i + 1)^{n/2+2k} } \cdot d\phi_i \right] \end{align*}These results have been validated numerically up to here. In particular, this last equation produces the same value as Equation 1.
The inner integral is difficult to compute. A reasonable first step seems to be a \(u\)-substitution with \(u:=c_i\). From there, by using repeated geometric series expansion, I was able to finally integrate it (again, resulting in a hypergeometric function). However, despite a lot of stuff canceling, I was not able to sum any of the indices I used down again.
Spherical Coordinates: Integrating \(\phi_i\)
Instead applying the Weierstrass substitution to the inner integral in Equation 1 produces a form which may be a useful starting point:
\begin{align*} \int_0^{\pi/2} \left(\cdots\right) ~d\phi_i &= \int_0^1 \frac{4 x (1-x^2)}{(1+x^2)^3} \left(\frac{ (~ (c_o\!-\!1)x^2 + (c_o\!+\!1) ~)^2 \hspace{43mm} }{ (1\!+\!x^2) (~ (1\!-\!c_o)x^2 + 2 s_o \cos(\Delta\theta) x + (c_o\!+\!1) ~) }\right)^{n/2} d x \\ &= \int_0^1 \frac{4 x (1-x^2)}{(1+x^2)^3} \left(\frac{ (c_o\!-\!1)^2 x^4 + 2(c_o^2\!-\!1)x^2 + (c_o\!+\!1)^2 \hspace{46mm} }{ (1\!-\!c_o)x^4 + 2 s_o \cos(\Delta\theta) x^3 + 2 x^2 + 2 s_o \cos(\Delta\theta) x + (c_o\!+\!1) }\right)^{n/2} d x \end{align*}Note that this is the inner integral only, not including the outside pieces. These forms have also been validated numerically.
That \((1+x^2)^{n/2}\) in the denominator very nearly (but doesn't) divide the numerator. That, or just straight polynomial division, could produce a form suitable for series expansion. That doesn't seem to lead anywhere, though.
Vector Approach: Stokes's Theorem
Another idea is to apply various vector integral theorems to solve this.
For example, Stokes's Theorem relates an integral over a surface \(\Sigma\) (with surface normal \(\hat{\vec{n}}\); be sure to note this is different from \(\vec{N}\)) to an integration around the surface's boundary \(\partial\Sigma\):
\[ \oint_{\partial\Sigma} \vec{F} \cdot d l = \iint_{\Sigma} ((\vec{\nabla}\times\vec{F}) \cdot \hat{\vec{n}}) ~d\Sigma \]In our case, the surface normal \(\hat{\vec{n}}=\vec{\omega}_i\). Let's also take \(\vec{\omega}_o=\vecinline{a,0,c}\), \(~~\vec{\omega}_i=\vecinline{x,y,z}\), just so we don't have to write subscripts everywhere. It might be nice if we could transform the integral over the hemisphere (the right-hand-side) into an integral around the circle at the base (the left-hand-side). The trouble is we have to write our integrand in terms of \((\vec{\nabla}\times\vec{F}) \cdot \vec{\omega}_i\).
Let \(\vec{\nabla}\times\vec{F} = \vec{G}\). To compute \(\vec{F}=\vecinline{F_x,F_y,F_z}\), first, we compute it up to a constant, apparently using Poincaré's lemma:
\begin{align*} \bar{F}_x &= \smallint_0^1 \left[~ t~z~G_y( t x, t y, t z ) - t~y~G_z( t x, t y, t z ) ~\right] ~d t\\ \bar{F}_y &= \smallint_0^1 \left[~ t~x~G_z( t x, t y, t z ) - t~z~G_x( t x, t y, t z ) ~\right] ~d t\\ \bar{F}_z &= \smallint_0^1 \left[~ t~y~G_x( t x, t y, t z ) - t~x~G_y( t x, t y, t z ) ~\right] ~d t \end{align*}As a test, I chose \(\vec{G}:=\vec{N}\), so that \(~~(\vec{\nabla}\times\vec{F}) \cdot \vec{\omega}_i ~~=~~ \vec{G}\cdot\vec{\omega}_i ~~=~~ \vec{N}\cdot\vec{\omega}_i~~\) and we'll get the integral of a cosine over the hemisphere, which from elementary methods we know should be \(\pi\).
Since \(\vec{G}=\vec{N}:=\vecinline{0,0,1}\), we easily integrate the above to get \(\bar{\vec{F}}=\vecinline{-y/2,x/2,0}\). We calculate \(\vec{\nabla}\times \bar{\vec{F}} = \vecinline{0,0,1} = 1 \vec{N}\), so the scaling factor is \(1\) and we have our final \(\vec{F}=\vecinline{-y/2,x/2,0}\). We can now apply Stokes's Theorem and integrate this over the circle:
\[ \int_\odot \vec{F} \cdot d l = \int_0^{2\pi} \left\| \vec{F}( \cos\theta, \sin\theta, 0 ) \right\| ~d\theta = \int_0^{2\pi} \sqrt{ \left(\frac{-\sin\theta}{2}\right)^2 + \left(\frac{\cos\theta}{2}\right)^2 } d\theta = \int_0^{2\pi} \frac{d\theta}{2} = \pi \]This is pretty awesome, but the real integrand is a lot more challenging.
As (maybe a less trivial) test, I chose \(\vec{G}:=\vecinline{0,0,z^3}\), so that \(~~(\vec{\nabla}\times\vec{F}) \cdot \vec{\omega}_i ~~=~~ \vec{G}\cdot\vec{\omega}_i ~~=~~ z^4 ~~=~~ (\vec{N}\cdot\vec{\omega}_i)^4~~\). Computing this over the hemisphere is harder, but elementary methods are again capable enough to tell us we should get \(2\pi/(n\!+\!1)\) (where \(n\) is the power), so we're looking for \(2\pi/5\).
Integrating gets \(\bar{\vec{F}}=\vecinline{~~ -y z^3/5,~~ x z^3/5,~~ 0 ~~}\). We calculate \(\vec{\nabla}\times \bar{\vec{F}} = \vecinline{~~ -3x z^2/5,~~ 3y z^2/5,~~ 2z^3/5 ~~}\) which . . . is wrong. We'll come back to why in a moment.
For the real integrand, using the fact that \(\vec{N}:=\vecinline{0,0,1}\):
\[ (\vec{\nabla}\times\vec{F}) \cdot \vec{\omega}_i ~~=~~ \left(\frac{ (\vec{\omega}_o + \vec{\omega}_i) \cdot \vec{N} }{ \sqrt{2}\sqrt{ \vec{\omega}_o\cdot\vec{\omega}_i + 1 } }\right)^n (\vec{N} \cdot \vec{\omega}_i) ~~=~~ \frac{ z~(c+z)^n }{ 2^{n/2} ( a x + c z + 1 )^{n/2} } \]Essentially, to get our \(\vec{G}\), we have to take the inverse-dot-product on the right hand side, and it's not clear how to do that because "inverse-dot-product" isn't really a thing. We could obviously do something dumb like the following (note we can drop the constant because this already needs to be scaled by some unknown factor):
\[ \vec{G}(x,y,z) = \cancel{2^{-n/2}} \vecinline{ 0, 0, \frac{ (c+z)^n }{ ( a x + c z + 1 )^{n/2} } } \]—which will indeed give the integrand when you take \(\vec{G}\cdot\vec{\omega}_i\). However, now we need to integrate it to get \(\bar{\vec{F}}\). The relevant terms are of the form \(\smallint_0^1 t~G_z(tx,ty,tz)~d t\), and this integrates to a mess[4]. After some bookkeeping to get \(\bar{\vec{F}}\), we then find the scaling factor to get \(\vec{F}\). Then we'd have to integrate that over the circle.
I did this, partially. It's challenging, but you can take the mess and substitute in the integration bounds, and it simplifies some. You can also substitute in \(z=0\) because we know we'll be integrating over the circle, and it simplifies a lot (you can interleave the series created from subtracting the higher from the lower bound, and then that factors into a constant and the whole series collapses into a couple gamma-function constants). There is a suspicious \(z^2\) in the denominator we can't sub in, though.
However, it's all for naught. Even doing all the integration numerically won't work. The second test should have clued us in. The reason is that inverting curl assumes implicitly that \(\vec{G}\) is divergence-free. That works for a constant (like in the first test), but not in the second, nor for the real integrand.
Vector Approach: Divergence Theorem
A completely different idea is to use the Divergence Theorem instead of Stokes's Theorem:
\[ \iiint_V \left(\vec{\nabla}\cdot\vec{F}\right) ~d V = \oiint_\Sigma (\vec{F} \cdot \hat{\vec{n}}) ~d\Sigma \]Note that here, the surface must be closed, so to apply to our situation of integrating over just the hemisphere, we need to subtract off (e.g.) the surface integral over the circular base:
\begin{align*} I = \iint_{\Omega_{\vec{N}}} (\vec{F} \cdot \vec{\omega}_i) ~d\vec{\omega}_i &= \iiint_V \left(\vec{\nabla}\cdot\vec{F}\right) ~d V - \iint_\odot (\vec{F} \cdot \vecinline{0,0,-1}) ~d A \\ &= \iiint_V \left(\vec{\nabla}\cdot\vec{F}\right) ~d V + \iint_\odot F_z ~d A \end{align*}Again using \(\vec{F}:=\vec{N}\) with expected result \(I=\pi\) as a test:
\begin{align*} I &= \cancelto{0}{ \iiint_V \left(\vec{\nabla}\cdot\vec{N}\right) ~d V } + \iint_\odot N_z ~d A = \iint_\odot ~d A = \pi \end{align*}Again, this is pretty awesome, but the real integrand is more challenging. We again need to do an inverse-dot-product to get our \(\vec{F}\), then sum the partials for volume and integrate \(\vec{F}\) itself for the circular base. If we use the same inverse-dot-product (and \(\vec{G}\)) as above, then the \(y\) partial is zero, and the integral becomes:
\begin{align*} 2^{n/2} I &= \iiint_V \left( \frac{\partial}{\partial x} + \frac{\partial}{\partial z} \right) G_z ~d V + \iint_\odot G_z ~d A \\ &= \iint_\odot G_z ~d A - \frac{1}{2} n \left[ a \iiint_V \left( \frac{ (c+z)^n }{ (~a x + c z + 1~)^{n/2+1} } \right) ~d V + \iiint_V \left( \frac{ (c+z)^{n-1} }{ (~a x + c z + 1~)^{n/2+1} } ( -2 a x - c z + c^2 - 2 ) \right) ~d V \right] \end{align*}This is pretty nasty. After splitting up the first integrand, as shown, you could maybe integrate first along the axis you just differentiated, but this doesn't seem to make life easier since the innermost integration bounds are messy, and you want them to do as little damage as possible. I wasn't able to integrate more than once for the first term either way.
Another idea is to try to find two functions whose gradients are the split up integrands, then apply the divergence theorem again to go back to (two) surface integrals, which then might be easier to solve. This would also conveniently remove again that extraneous circle integral, which is also nontrivial.