Comparison of fits to Fresnel

True Fresnel

Unpolarized Fresnel

Form from Cook & Torrance, A Reflectance Model for Computer Graphics, ACM SIGGRAPH 1981
(/. is the Mathematica substitution operator, get the arrow by typing ->)

(1/2)(((((g-c)/(g+c))))^2)(1+((((((g+c)c-1)/((g-c)c+1))))^2))/.g->(Sqrt[(n^2)+(c^2)-1]);

(Subscript[F, 0]) = Value looking at the surface head-on at θ = 0 AKA cos(θ) = c = 1
(use Simplify with a constraint on n so it knows it can convert (Sqrt[(n^2)]) to n)

Simplify[(1/2)(((((g-c)/(g+c))))^2)(1+((((((g+c)c-1)/((g-c)c+1))))^2))/.g->(Sqrt[(n^2)+(c^2)-1])/.c->1,n>0]

(Subscript[F, 0]) for some common materials from https://en.wikipedia.org/wiki/List_of_refractive_indices

IndicesOfRefraction={{name->"Air",n->1},{name->"Water",     n->1.333},{name->"Glass",     n->1.5},{name->"Diamond",n->2.417}};

Panel[TableForm[Table[((n/.IndicesOfRefraction[[j]])/(n/.IndicesOfRefraction[[i]])),{i,1,Length[IndicesOfRefraction]},{j,1,Length[IndicesOfRefraction]}],TableHeadings->{name/.IndicesOfRefraction,name/.IndicesOfRefraction}],"Ratio of indices of refraction"]

Panel[TableForm[Table[(((((η-1)/(η+1))))^2)/.η->((n/.IndicesOfRefraction[[j]])/(n/.IndicesOfRefraction[[i]])),{i,1,Length[IndicesOfRefraction]},{j,1,Length[IndicesOfRefraction]}],TableHeadings->{name/.IndicesOfRefraction,name/.IndicesOfRefraction}],"(Subscript[F, 0]) between materials"]

Conversion from (Subscript[F, 0]) back to index of refraction

Simplify[Solve[(Subscript[F, 0])==(((((n-1)/(n+1))))^2),n]]

First solution is entering material (n/1) vs. exiting (1/n)

{{n->((1+(Sqrt[(Subscript[F, 0])]))/(1-(Sqrt[(Subscript[F, 0])])))},{n->((1-(Sqrt[(Subscript[F, 0])]))/(1+(Sqrt[(Subscript[F, 0])])))}}/.(Subscript[F, 0])->0.04

Fresnel in terms of cosine and (Subscript[F, 0])

Make into a callable function F[c,f]

F[c_,f_]:=(1/2)(((((g-c)/(g+c))))^2)(1+((((((g+c)c-1)/((g-c)c+1))))^2))/.g->(Sqrt[(n^2)+(c^2)-1])/.n->((1+(Sqrt[f]))/(1-(Sqrt[f])))

Graph Fresnel for different values of (Subscript[F, 0])
Manipulate makes an interactive graph with sliders specified as {variable, min, max[, step]} or {{var, initialval}, min, max[, step]})
Click the “+“ next to the slider to see or enter the actual value

Manipulate[Plot[F[c,f],{c,0,1},PlotRange->{0,1}],{{f,0.04},0.01,1}]

Generate some data at the index of refraction of glass

FresnelData=ParallelTable[{c,F[c,0.04]}/.c->RandomReal[],{i,1,1000000}];

Schlick approximation

To get rid of dependence on  (Subscript[f, 0]) , try lerp between (Subscript[f, 0]) and 1: S[c,f] = f + (1-f) * h[c] for h[c] ∈ [0,1] and h[0] = 1, h[1] = 0.

Manipulate[Plot[((F[c,f]-f)/(1-f)),{c,0,1},PlotRange->{0,1}],{{f,0.04},0.02,0.2}]

Looks kind of like a flipped polynomial. What’s the best exponent?

FindFit[FresnelData,0.04+(1-0.04)(((1-c))^p),p,c]

Pretty close to 5, and 5th power is easier to compute in three simple floating point operations as (x^2)*(x^2)*x

S[c_,f_]:=f+(1-f)(((1-c))^5)

Error comparison

Plot Schlick compared to true Fresnel

Manipulate[Plot[Evaluate[Legended[Tooltip[func,label],label]/.{
{func->F[c,f],label->"Fresnel"},
{func->S[c,f],label->"Schlick"}
}],{c,0,1},PlotRange->{0,1}],
{{f,0.04},0.02,0.2}]

Plot difference from true Fresnel

Manipulate[Plot[S[c,f]-F[c,f],{c,0,1},PlotRange->Full],{{f,0.04},0.02,0.2}]

Squared error for all indices of refraction between water and diamond

Plot3D[(((S[c,f]-F[c,f]))^2),{c,0,1},{f,0.02,0.2},AxesLabel->Automatic,PlotRange->Full]

Numeric worst case error is around diamond at about 84º

SchlickMaxError=FindMaximum[{S[c,f]-F[c,f],0≤c≤1,0.02≤f≤0.2},{{c,0.1},{f,0.2}}]

Propagation of error

Error scales as derivative

D[S[c,f],c]

Max for c=0, f=0, where scale factor is 5.

FindMaximum[{Abs[-5 (((1-c))^4) (1-f)],0≤c≤1,0≤f≤0.2},{c,f}]

8-bit texture value tx={r,g,b}±{(Subscript[e, r]),(Subscript[e, g]),(Subscript[e, b])}, where (Subscript[e, r]),(Subscript[e, g]),(Subscript[e, b])≤(2^(-9))
Convert texture to normal in as (OverHat[n])=2*tx-1; Error propagation

D[2*tx-1,tx]

So normal is (OverHat[n])+(Overscript[(Subscript[e, n]), ⟶]), where (Overscript[(Subscript[e, n]), ⟶])≤2*((2^(-9)),(2^(-9)),(2^(-9))) ; where  2*(2^(-9))=(2^(-8)) for 8-bit normal map.
Then the cosine argument c±(Subscript[e, c])=(OverHat[v])•(OverHat[n]) ±(Subscript[e, c])= (OverHat[v])•(OverHat[n]) ±(OverHat[v])• (Overscript[(Subscript[e, n]), ⟶]); (OverHat[v])• (Overscript[(Subscript[e, n]), ⟶])=|(OverHat[v])||(Overscript[(Subscript[e, n]), ⟶])|Cos[θ]. Given (OverHat[v]) on a unit sphere (so |(OverHat[v])|=1), (Overscript[(Subscript[e, n]), ⟶]) in a box with sides (2^(-8)) and max diagonal (Sqrt[3])(2^(-8)), and max Cos[θ] = 1 
The total error is 5(Sqrt[3])(2^(-8))≤0.034for 8-bit normal map

Actual worst case error for Schlick is comparable to this

N[5(Sqrt[3])(2^(-8))]/.SchlickMaxError[[2]]

Spherical Gaussian-inspired Fitting

Spherical-Gaussian-inspired approximation. Not really a Spherical Gaussian, since that’d be (2^(k(c-1))). This is just an exponential of a generic polynomial with the constraint that the function is 1 at c=0.

P[c_,f_,a_,b_]:=f+(1-f)(2^((a c+ b)c))

SG fit to Schlick (matches Sébastien Lagarde blog post from 2012)

Note that Lagarde just fit the non-f term of the equation. Using (((S[c,f] - P[c,f,a,b]))^2) (with variable f) as the error term has an extra (((1-f))^2) factor that doesn’t affect the fitting to a and b, but does prevent numerical integration. Choosing 0.04 (or any value for f) for the fit does not change the a or b that it finds.

SchlickData=ParallelTable[{c,S[c,0.04]}/.c->RandomReal[],{i,1,1000000}];

SGSchlick=FindFit[SchlickData,P[c,0.04,a,b],{a,b},c]

NMaxValue[{(((P[c,0.04,a,b]-F[c,0.04]))^2)/.SGSchlick,0≤c≤1},c]

Fit directly to Fresnel at index of refraction of glass

Timing says how many seconds it took

Timing[SGGlassRMS=FindFit[FresnelData,P[c,0.04,a,b],{a,b},c]]

NMaxValue[{(((P[c,0.04,a,b]-F[c,0.04]))^2)/.SGGlassRMS,0≤c≤1},c]

Minimize max absolute error (this one is a bit slower)

Timing[SGGlassMax=FindFit[FresnelData,P[c,0.04,a,b],{a,b},c,NormFunction->(Norm[#,∞]&),Method->"NMinimize"]]

NMaxValue[{(((P[c,0.04,a,b]-F[c,0.04]))^2)/.SGGlassMax,0≤c≤1},c]

Error plots

Visual comparison

Manipulate[Plot[
Evaluate[Legended[Tooltip[func,label],label]/.{
{func->F[c,f],label->"Fresnel"},
{func->S[c,f],label->"Schlick"},
{func->P[c,f,a,b]/.SGSchlick,label->"Fit to Schlick"},
{func->P[c,f,a,b]/.SGGlassRMS,label->"Fit to Fresnel"},
{func->P[c,f,a,b]/.SGGlassMax,label->"Minimax Fit"}
}],{c,0,1},PlotRange->{0,1}],{{f,0.04},0.02,0.2}]

Absolute error compared to Fresnel

Manipulate[Plot[
Evaluate[Legended[Tooltip[func,label],label]/.{
{func->S[c,f]-F[c,f],label->"Schlick"},
{func->P[c,f,a,b]-F[c,f]/.SGSchlick,label->"Fit to Schlick"},
{func->P[c,f,a,b]-F[c,f]/.SGGlassRMS,label->"Fit to Fresnel"},
{func->P[c,f,a,b]-F[c,f]/.SGGlassMax,label->"Minimax Fit"}
}],{c,0,1},PlotRange->Full],{{f,0.04},0.02,0.2}]