Add clamp to prevent negative metallic Fresnel

[portsmouth]

As discussed in #187, it is useful to quote this formula. (For example, we use it in the multiple scattering compensation for the OpenPBR metal lobe in Arnold).

I put it in a new footnote:

Referenced from the metal section (slightly tweaked to define the $b$ term which appears in the albedo formula, which makes for a clearer presentation anyway closer to the ASM document):

[peterkutz]

Thanks for adding this @portsmouth ! Sorry I had forgotten about it.

This looks good to me.

I had to figure out how you derived this version of the hemispherical albedo formula. Starting with a standard integral of the cosine-weighted F82-tint function over the upper hemisphere in polar coordinates, one can simplify away the phi part and then do a change of variables to convert the integration domain from theta to cos(theta) and ultimately arrive at the form shown in your footnote.

I don't know if any of that needs to be included in the text (or maybe it's already clear if considering the context), but I wanted to understand it at least. The method of deriving the final closed-form solution is also not obvious, and I believe it involves Beta integrals, but presumably most people wouldn't be solving this by hand.

[portsmouth]

I actually just assumed the solution given by Wolfram Alpha in the linked issue was correct.. It's good to verify it by hand though for sure (and can also check numerically of course).

[portsmouth]

I checked numerically and indeed the formula from Wolfram seems to match the numerical results.

Here plotting the albedo versus the F82-tint, for varies $F_0$ values ranging over [0,1] (black to red). Lines are theory, dots are the numerical integral.

However one unexpected thing is that the albedo can actually go negative it seems (!). The boxed point in the lower left for example (with $F_0=0$, and F82-tint zero) has negative albedo: -0.017.

From the formula, if $F_0=0$, then the albedo is negative if $b > 6$. And that can happen, e.g. if the tint is 0 (and $F_0=0$), then $b = 1/[\bar{\mu}(1 - \bar{\mu})] = 8.1666$.

So it seems this F82-tint model can generate negative Fresnel factors. For example with $F_0=0$ and tint 0:

$F_{82}(\mu) = (1 - \mu^5) \left( 1 - \frac{\mu}{\bar{\mu} } \right)$

which is negative for all $\mu \gt \bar{\mu}$ (i.e. all non-grazing angles).

That doesn't seem ideal, and at the least we should stipulate that a clamp has to be applied to avoid that. A clamp introduces a discontinuity through, so we should discuss (maybe there is a nicer modification that avoids this?).

(Applying the clamp would make the nice formula for the average albedo no longer correct though, so we would probably need to do a numerical fit instead).

(I noticed this because it is causing odd results in our testsuite, since it is not that unusual to have a dark $F_0$ and edge-tint, producing a negative Fresnel which creates a massive artifact).

[portsmouth]

Script for the plot, for reference:

import os, sys, math
import numpy as np
import matplotlib.pyplot as plt
import colorsys

def Fresnel_Schlick(mu, F0):
    return F0 + (1.0 - F0)*(1.0 - mu)**5.0

def get_b(F0, F82tint):
    mu_bar = 1.0/7.0
    denom = mu_bar * (1.0 - mu_bar)**6.0
    b = Fresnel_Schlick(mu_bar, F0) * (1.0 - F82tint) / denom
    return b

def Fresnel_F82(mu, F0, F82tint):
    b = get_b(F0, F82tint)
    return Fresnel_Schlick(mu, F0) - b*mu*(1.0 - mu)**6.0

# Compute F82-tint model albedo numerically
def albedo_numerical(F0, F82tint):

    Ntheta = 256
    Integrand = np.empty(Ntheta)
    mu_array = np.linspace(0.0, 1.0, Ntheta)

    # evaluate Integrand for theta_i integral
    for n_theta_i in range(Ntheta):
        mu_i = mu_array[n_theta_i]
        F = Fresnel_F82(mu_i, F0, F82tint)
        Integrand[n_theta_i] = 2.0 * F * mu_i

    # do integral over theta_i
    theta_integral = np.trapz(Integrand, mu_array)
    return theta_integral

# Analytical albedo (from Wolfram)
def albedo_analytical(F0, F82tint):
    b = get_b(F0, F82tint)
    return F0 + (1.0 - F0)/21.0 - b/126.0

N_F0s = 8
F0_values = np.linspace(0.0, 1.0, N_F0s)

N_tints_numerical = 16
F82tint_values_numerical = np.linspace(0.0, 1.0, N_tints_numerical)
albedos_numerical = np.empty(N_tints_numerical)

N_tints_analytical = 64
F82tint_values_analytical = np.linspace(0.0, 1.0, N_tints_analytical)
albedos_analytical = np.empty(N_tints_analytical)

for n_f0 in range(0, N_F0s):

    F0 = F0_values[n_f0]
    print('Running numerics for F0 = %f' % F0)

    for n_t in range(0, N_tints_numerical):
        F82tint = F82tint_values_numerical[n_t]
        albedos_numerical[n_t]  = albedo_numerical(F0, F82tint)
    for n_t in range(0, N_tints_analytical):
        F82tint = F82tint_values_analytical[n_t]
        albedos_analytical[n_t] = albedo_analytical(F0, F82tint)

    grayscale = F0 / F0_values[-1]
    Clin = colorsys.hsv_to_rgb(1, 0.9, grayscale)
    plt.plot(F82tint_values_numerical, albedos_numerical,  label='F0 = %f' % F0, color=Clin, marker='.', linestyle='none')
    plt.plot(F82tint_values_analytical, albedos_analytical, label='',             color=Clin, marker='', linewidth=1.0, linestyle='solid')

plt.xlabel (r'$F_{82}$ tint')
plt.ylabel (r'$E_\mathrm{avg}$, average albedo')
plt.show()

[portsmouth]

If we add the explicit clamp to keep the F82-tint Fresnel positive, the resulting true albedo (numerical result shown as the dots) is not very different from the original albedo formula stated with a clamp applied (the lines):

There is non-negligible discrepancy only for low $F_0$ and tint. Here is the $F_0=0$ case zoomed in:

The clamped albedo formula has roughly at most a magnitude 0.01 discrepancy with the true albedo (of the clamped F82-tint Fresnel), and only for very dark metals. So probably it would be fine to still quote it, with the caveat that it is not totally accurate for dark metals.

Though, it would be nice to develop a better approximate formula for the true albedo. (Assuming we're not modifying the F82-tint Fresnel, except to impose a clamp which seems unavoidable).

[portsmouth]

I made the clamp explicit in the spec in e0963af (screenshots at the top updated).

[portsmouth]

The specular_weight range was allowed to exceed 1, I recall, for dielectrics, so it must presumably also be allowed for metals. So in fact we need another fix to the spec, to make the range consistent for the metal, with an explicit clamp of the total Fresnel to be less than 1. (We already do this in Arnold, passing the specular_weight along into the conductor Fresnel, where the clamped multiplication is done).

The MaterialX graph cannot implement this currently though (as noted in #240) , since there is no way to implement the weight other than multiplying it into the BSDF, which cannot be clamped. So we need to discuss how/whether to deal with this Fresnel clamping for MaterialX.

EDIT: ah no, we did this already, merged into dev_1.2.. #238

We just need to rebase this onto dev_1.2 (and fix conflicts).

[portsmouth]

Updated to dev_1.2. The new text has a clamp applied to the metallic Fresnel at both the upper and lower limits:

[portsmouth]

So we need to discuss how/whether to deal with this Fresnel clamping for MaterialX.

All we can do for now in the graph, is to clamp the specular_weight into $[0,1]$. (Which will look wrong for specular_weight $> 1$, but at least won't break energy conservation).

This still wouldn't deal with the negative Fresnel though.

[peterkutz]

Nice investigation. The spec changes you've proposed look good to me.

It would be nice if we could avoid the clamping entirely. If I remember correctly, the problem is that the original F82 formula lets you specify the y value of the curve at a particular x value (corresponding to about 82 degrees), but that doesn't mean that the minimum y value coincides with that location. Since the overall curve slopes up from left to right, the minimum y value tends to occur at a lower x value.

Clamping both the curve and the average albedo formula seems like a reasonable solution for handling these edge cases for now.

It might be possible to find an exact closed-form solution for the integral of the clamped function, though it seems possible that it would be be significantly more complex and would start to defeat one of the purposes of using this model in the first place: it's relative simplicity and efficiency. Perhaps we could calculate the x intercepts and then subtract off the integral of the part of the curve between them.

[AdrienHerubel]

There is a consensus that this change fixes unwanted behaviors. This does not require a change in the MaterialX graph and the negative clamp should already be present in the MaterialX implementation.