GGX Specular BRDF

(*Clear["Global`*"]*)

ggx[μ_,α_]:=((α^2)/(π((((μ^2)((α^2)-1)+1))^2)))
smith[θ_,α_]:=(1/(Cos[θ]+(Sqrt[(α^2)+(1-(α^2))*((Cos[θ])^2)])))
(* change[] gives cos((Subscript[θ, h])) when provided with (Subscript[ω, i]) and (Subscript[ω, o]) *)
change[θo_,θi_,ϕi_]:=(Cos[θo]+Cos[θi])/(√(2*(1+Cos[θo]*Cos[θi]+Sin[θo]*Sin[θi]*Cos[ϕi])))

Let’s test our formulas:

Plot3D[ggx[Cos[x],y],{x,0,(π/2)},{y,0.001,1}]
Plot3D[smith[x,y],{x,0,(π/2)},{y,0.001,1},PlotRange->{0,2},AxesLabel->{"θ","α","Smith"}]
change[0,(π/2),0]

Manipulate[Plot3D[change[θo,θ,ϕ],{θ,0,(π/2)},{ϕ,-π,π},PlotRange->{0,1},AxesLabel->{"(Subscript[θ, i])","(Subscript[ϕ, i])","(Subscript[θ, h])"}],{{θo,0},0,(π/2)}]

Numerical Integration

Let’s try Mathematica’s numerical integration instead and see if we can go faster:

(* Don't evaluate that! Takes FOREVER! *)
(*Ei[θo_,ϕ_,α_]:=Integrate[smith[θi,α]*ggx[change[θo,θi,ϕ],α] * Cos[θi]* Sin[θi],{θi,0,(π/2)}]*)
Ei[0,0,0.5]

Ei[θo_,ϕ_,α_]:=NIntegrate[smith[θi,α]*ggx[change[θo,θi,ϕ],α] * Cos[θi]* Sin[θi],{θi,0,(π/2)}]
θo=0;
ϕ=0;
α=0.5;
Ei[0,0,0.5]

Indeed it IS faster! :D

Technically, the expression we need to evaluate is this:

Eo[θo_,α_]:=smith[θo,α]*Integrate[Integrate[smith[θi,α]*ggx[change[θo,θi,ϕi],α] * Cos[θi]* Sin[θi],{θi,0,(π/2)}],{ϕi,-π,π}]

With numerical integration we get instead

ϕi = -π;
θo=0.1;
α=0.5;
NIntegrate[smith[θi,α]*ggx[change[θo,θi,ϕi],α] * Cos[θi]* Sin[θi],{θi,0,(π/2)}]

θo=0.2;
α=0.005;
smith[θo,α]*NIntegrate[smith[θi,α]*ggx[change[θo,θi,ϕi],α] * Cos[θi]* Sin[θi],{θi,0,(π/2)},{ϕi,-π,π}]
(*smith[θo,α]*NIntegrate[Ei[θo,ϕi,α],{ϕi,-π,π}]*)

Finally, we can build the final expression for irradiance:

Eo[θo_,α_]:=smith[θo,α]*NIntegrate[smith[θi,α]*ggx[change[θo,θi,ϕi],α] * Cos[θi]* Sin[θi],{θi,0,(π/2)},{ϕi,-π,π}]
Eo[0.1,0.5]

α=0.9;
2*π*NIntegrate[smith[θo,α]*smith[θi,α]*ggx[change[θo,θi,ϕi],α] *Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]

Building the Table of Irradiance Values

So let’s build a table of these values!

SetDirectory["D:\\Workspaces\\GodComplex\\Tests\\TestMSBRDF\\Tables"]; 

countX = 128;
countY= 128;

tableFileName="MSBRDF_GGX_E"  <>ToString[countX] <>"x" <> ToString[countY] <>".csv"
tableFileNameAvg="MSBRDF_GGX_Eavg"  <>ToString[countY] <>".csv"

Eo[θo_,α_]:=smith[θo,α]*NIntegrate[smith[θi,α]*ggx[change[θo,θi,ϕi],α] * Cos[θi]* Sin[θi],{θi,0,(π/2)},{ϕi,-π,π}]

(* [2018-07-31] Prefer using the **correlated** Smith masking/shadowing function!!! *)
lambda[θ_,α_]:=(1/2)((Sqrt[1+(((α*Tan[θ]))^2)])-1)
smith2[θl_,θv_,α_]:=(1/(1+lambda[θl,α]+lambda[θv,α]))
BRDFggx[θo_,θi_,ϕi_,α_]:= ((smith2[θo,θi,α]*ggx[change[θo,θi,ϕi],α])/(4*Cos[θo]*Cos[θi]))  (* Here, we restore the geometric term in the denominator *)
Eo[θo_,α_]:=NIntegrate[BRDFggx[θo,θi,ϕi,α] * Cos[θi]* Sin[θi],{θi,0,(π/2)},{ϕi,-π,π}]
tableFileName="MSBRDF_GGX_E"  <>ToString[countX] <>"x" <> ToString[countY] <>"_G2.csv"  (* G2 denotes the correlated Smith term instead *)

buildTable=Function[{i},
x=N[((Mod[i,countX])/(countX-1))];
y=N[((Quotient[i,countX])/(countY-1))];
θo=ArcCos[x];
α=Max[0.01,y];
{x,y,N[Eo[θo,α]]}
];
buildTable[128*5+17]

Function[{i},x=N[((Mod[i,countX])/(countX-1))];y=N[((Quotient[i,countX])/(countY-1))];
θo=ArcCos[x];
α=Max[0.01`,y];
{x,y,N[Eo[θo,α]]}
]

resultEo= Table[buildTable[i],{i,0,countX*countY-1}];
Export[tableFileName,resultEo,"CSV"]

ListPlot3D[resultEo]
Image[Table[Table[1-resultEo[[1+countX*y+x]][[3]],{x,0,countX-1}],{y,0,countY-1}]]

resultEo = Import[tableFileName,"CSV"];
ListPlot3D[resultEo,AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}]

Let’s reorder the table and create an interpolation function from it:

t=Table[resultEo[[1+countX*(y-1)+(x-1)]][[3]],{x,1,countY},{y,1,countY}];
fEo=ListInterpolation[t,{{0,1},{0,1}}]
Plot3D[fEo[x,y],{x,0,1},{y,0,1},AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}]

Computing Average Irradiance

Now let’s integrate the full white furnace irradiance using our table:

	(Subscript[E, avg])(α)=∫(Subscript[E, o])((Subscript[θ, o]),α) cos((Subscript[θ, o])) d(Subscript[θ, o])

resultEavg=Table[{N[(i/(countY-1))],2*π*NIntegrate[fEo[μ,(i/(countY-1))]μ,{μ,0,1}]},{i,0,countY-1}]
Export[tableFileNameAvg,resultEavg,"CSV"]

Export[tableFileNameAvg,resultEavg,"CSV"]

ListLinePlot[resultEavg]

t=Table[resultEavg[[i]][[2]],{i,1,countY}];
fEavg=ListInterpolation[t,{0,1}]
Plot[fEavg[x],{x,0,1},PlotRange->{0,π}, AxesLabel->{"α","(Subscript[E, avg])"}]

Export as binary float arays for C#:

t=Table[resultEo[[i]][[3]],{i,1,countX*countY}];
Export["MSBRDF_GGX_E"<>ToString[countX]<>"x"<>ToString[countY]<>".float",t,"Real32"]
t=Table[resultEavg[[i]][[2]],{i,1,countY}];
Export["MSBRDF_GGX_Eavg"<>ToString[countY]<>".float",t,"Real32"]

Fitting Eavg

Clear[a,b,c,d,x]
fittingMethod="LevenbergMarquardt";
fittingMethod="PrincipalAxis";
fittingMethod="Automatic";

modelEavg=π+0*a+(0+1*b)*x+c*(x^2)+d*(x^3)+e*(x^4);
fitEavg=FindFit[resultEavg,modelEavg,{a,b,c,d,e},x,Method->fittingMethod]
modelEavg/.fitEavg
funcEavg=Function[x,Evaluate[modelEavg/.fitEavg]];
Show[
Plot[fEavg[x],{x,0,1},PlotRange->{0,π},AxesLabel->{"α","(Subscript[E, avg])"}],
Plot[funcEavg[x],{x,0,1},PlotStyle->Red]
]

Let’s check the BRDFs sum to π!

BRDFggx[θo_,θi_,ϕ_,α_]:= smith[θo,α]*smith[θi,α]*ggx[change[θo,θi,ϕi],α]
BRDFms[θo_,θi_,ϕ_,α_]:=(((1-fEo[Cos[θo],α])*(1-fEo[Cos[θi],α]))/(π-fEavg[α]))

whiteFurnaceSS[α_]:=2*π*NIntegrate[BRDFggx[θo,θi,ϕ,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]
whiteFurnaceMS[α_]:=2*π*NIntegrate[BRDFms[θo,θi,ϕ,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]

whiteFurnaceSS[0.3]
whiteFurnaceMS[0.3]

whiteFurnaceResultSS=Table[{α,whiteFurnaceSS[Max[0.01,α]]},{α,0,1,0.05}]
whiteFurnaceResultMS=Table[{α,whiteFurnaceMS[Max[0.01,α]]},{α,0,1,0.05}]

Show[
ListLinePlot[whiteFurnaceResultSS,PlotRange->{0,1.1*π},PlotStyle->Gray,AxesLabel->{"α",""}],
ListLinePlot[whiteFurnaceResultMS,PlotStyle->Blue],
ListLinePlot[whiteFurnaceResultSS+whiteFurnaceResultMS,PlotRange->Full,PlotStyle->Red]
]

Oren-Nayar Diffuse BRDF

Plotting Oren-Nayar Coefficients

Just to make sure the standard deviation σ is in radians! (they keep giving it in degrees, I started to be a bit skeptical for a moment)...

Plot[1-0.5((σ^2)/((σ^2)+0.33)),{σ,0,2}]

Plot[0.45((σ^2)/((σ^2)+0.09)),{σ,0,2}]

oren=Function[ {θo,ϕo,θi,ϕi,α},
σ=(π/2)α;
A = 1-0.5*((σ^2)/((σ^2)+0.33));
B = 0.45*((σ^2)/((σ^2)+0.09));
a=Max[θo,θi];
b=Min[θo,θi];
(1/π)(A+B*Max[0,Cos[ϕo-ϕi]]*Sin[a]*Tan[b])
];
oren[0,0,0.1,0.2,0.2]

Manipulate[Plot3D[oren[θo,0,x,0,α],{x,0,(π/2)},{α,0,1},PlotRange->{0,1},AxesLabel->{"cos(θ)","α","Oren"}],{{θo,(π/2)},0,(π/2)}]

Let’s try integrating the Lambertian case:

α=0.0;
2*π*NIntegrate[oren[θo,0,θi,ϕi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)}]

orenEo[θo_,α_]:=NIntegrate[oren[θo,0,θi,ϕi,α] * Cos[θi]* Sin[θi],{θi,0,(π/2)},{ϕi,-π,π}]
orenEo[1.5,1.0]

SetDirectory["D:\\Workspaces\\GodComplex\\Tests\\TestMSBRDF\\Tables"]; 
countX = 32;
countY= 32;
tableFileName="MSBRDF_OrenNayar_E"  <>ToString[countX] <>"x" <> ToString[countY] <>".csv"
tableFileNameAvg="MSBRDF_OrenNayar_Eavg"  <>ToString[countY] <>".csv"

orenEo[θo_,α_]:=NIntegrate[oren[θo,0,θi,ϕi,α] * Cos[θi]* Sin[θi],{θi,0,(π/2)},{ϕi,-π,π}]

buildTable=Function[{i},
x=N[((Mod[i,countX])/(countX-1))];
y=N[((Quotient[i,countX])/(countY-1))];
θo=ArcCos[x];
α=y;
{x,y,Min[1,orenEo[θo,α]]}
]
buildTable[countX*5+7]

resultOrenEo= Table[buildTable[i],{i,0,countX*countY-1}];
Export[tableFileName,resultOrenEo,"CSV"]

Export[tableFileName,resultOrenEo,"CSV"]

ListPlot3D[resultOrenEo,PlotRange->{0,1},AxesLabel->{"cos((Subscript[θ, o]))","α"}]
Image[Table[Table[1-resultOrenEo[[1+countX*y+x]][[3]],{x,0,countX-1}],{y,0,countY-1}]]

Re-Importing and testing...

resultOrenEo = Import[tableFileName,"CSV"];
ListPlot3D[resultOrenEo,AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}]

Let’s reorder the table and create an interpolation function from it:

t=Table[resultOrenEo[[1+countX*(y-1)+(x-1)]][[3]],{x,1,countY},{y,1,countY}];
fOrenEo=ListInterpolation[t,{{0,1},{0,1}}]
Plot3D[fOrenEo[x,y],{x,0,1},{y,0,1},AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}]

Computing Average Irradiance

Now let’s integrate the full white furnace irradiance using our table:

	(Subscript[E, avg])(α)=∫(Subscript[E, o])((Subscript[θ, o]),α) cos((Subscript[θ, o])) d(Subscript[θ, o])

resultOrenEavg=Table[{N[(i/(countY-1))],2*π*NIntegrate[fOrenEo[μ,(i/(countY-1))]μ,{μ,0,1}]},{i,0,countY-1}]
Export[tableFileNameAvg,resultOrenEavg,"CSV"]

resultOrenEavg=Import[tableFileNameAvg,"CSV"]

ListLinePlot[resultOrenEavg,PlotRange->{0,π}]

t=Table[resultOrenEavg[[i]][[2]],{i,1,countY}];
fOrenEavg=ListInterpolation[t,{0,1}]
Plot[fOrenEavg[x],{x,0,1},PlotRange->{0,π}, AxesLabel->{"α","(Subscript[E, avg])"}]

Export as binary float arays for C#:

t=Table[resultOrenEo[[i]][[3]],{i,1,countX*countY}];
Export["MSBRDF_OrenNayar_E"<>ToString[countX]<>"x"<>ToString[countY]<>".float",t,"Real32"]
t=Table[resultOrenEavg[[i]][[2]],{i,1,countY}];
Export["MSBRDF_OrenNayar_Eavg"<>ToString[countY]<>".float",t,"Real32"]

Fitting Eavg

Clear[a,b,c,d,x]
fittingMethod="LevenbergMarquardt";
fittingMethod="PrincipalAxis";
fittingMethod="Automatic";

modelOrenEavg=π+0*a+(0+0*b)*x+c*(x^2)+d*(x^3)+e*(x^4);
fitOrenEavg=FindFit[resultOrenEavg,modelOrenEavg,{a,b,c,d,e},x,Method->fittingMethod]
modelOrenEavg/.fitOrenEavg
funcOrenEavg=Function[x,Evaluate[modelOrenEavg/.fitOrenEavg]];
Show[
Plot[fOrenEavg[x],{x,0,1},PlotRange->{0,π},AxesLabel->{"α","(Subscript[E, avg])"}],
Plot[funcOrenEavg[x],{x,0,1},PlotStyle->Red]
]

Let’s check the BRDFs sum to π!

BRDForen[θo_,θi_,ϕ_,α_]:= oren[θo,0,θi,ϕ,α]
BRDForenms[θo_,θi_,ϕ_,α_]:=(((1-fOrenEo[Cos[θo],α])*(1-fOrenEo[Cos[θi],α]))/(π-fOrenEavg[α]))

whiteFurnaceSS[α_]:=2*π*NIntegrate[BRDForen[θo,θi,ϕi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]
whiteFurnaceMS[α_]:=2*π*NIntegrate[BRDForenms[θo,θi,ϕi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]

whiteFurnaceSS[0.3]
whiteFurnaceMS[0.3]

whiteFurnaceResultSS=Table[{α,whiteFurnaceSS[α]},{α,0,1,0.05}]
whiteFurnaceResultMS=Table[{α,whiteFurnaceMS[α]},{α,0,1,0.05}]

Show[
ListLinePlot[whiteFurnaceResultSS,PlotRange->{0,1.1*π},PlotStyle->Gray,AxesLabel->{"α",""}],
ListLinePlot[whiteFurnaceResultMS,PlotStyle->Blue],
ListLinePlot[whiteFurnaceResultSS+whiteFurnaceResultMS,PlotRange->Full,PlotStyle->Red]
]

Analytical fitting (not really finished, not that much interesting anyway)

Clear[a,b,c,d,e,f,x,y]
fittingMethod="LevenbergMarquardt";
fittingMethod="PrincipalAxis";
fittingMethod="Automatic";

modelOrenEo=a+b*x+c*y+d*(x^2)+e*(y^2)+f*(x^2)*y+g*x*(y^2)+h*(x^3)+i*(y^3)+j*(x^3)*y+k*x*(y^3)+l*(x^3)*(y^3);
modelOrenEo=(a+b*x+c*(x^2)+d*(x^3)) + (e+f*x+g*(x^2)+h*(x^3))*y+ (i+j*x+k*(x^2)+l*(x^3))*(y^2)+ (m+n*x+o*(x^2)+p*(x^3))*(y^3);
(*fitFuncOren[a_?NumericQ,b_?NumericQ,x_?NumericQ]:=y/.FindRoot[a y+b Log[y]==x,{y,1.}];*)
fitfunc[x_,y_,a_,b_,c_,d_,e_,f_,g_,h_,i_,j_,k_,l_,m_,n_,o_,p_]:=(a+b*x+c*(x^2)+d*(x^3)) + (e+f*x+g*(x^2)+h*(x^3))*y+ (i+j*x+k*(x^2)+l*(x^3))*(y^2)+ (m+n*x+o*(x^2)+p*(x^3))*(y^3)
fitOrenEo=FindFit[resultOrenEo,fitfunc[x,y,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p],{a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p},{x,y},Method->fittingMethod]
Show[
Plot3D[fOrenEo[x,y],{x,0,1},{y,0,1},AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}],
Plot3D[Min[1,modelOrenEo/.fitOrenEo],{x,0,1},{y,0,1},PlotStyle->Red]
]

Clear[x,y,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q]
fitfunc[x_,y_,a_,b_,c_,d_,e_,f_,g_,h_,i_,j_,k_,l_,m_,n_,o_,p_]:=Min[1,(a+b*x+c*(x^2)+d*(x^3)) + (e+f*x+g*(x^2)+h*(x^3))*y+ (i+j*x+k*(x^2)+l*(x^3))*(y^2)+ (m+n*x+o*(x^2)+p*(x^3))*(y^3)]

t=Table[{ ((y-1)/(countX-1)), testResultOrenEo[[1+countX*(y-1)+(x-1)]][[3]]},{x,1,countX},{y,1,countY}];
t[[1]];
ListLinePlot[t[[1]]];
Show[ListLinePlot[t]]

modelOrenEoY=a+b*y+c*(y^2)+d*(y^3)
fitsOrenEoY=Table[FindFit[t[[,modelOrenEoY,{a,b,c,d,e,f,g,h,i},{x,y},Method->fittingMethod]
Show[
Plot3D[fOrenEo[x,y],{x,0,1},{y,0,1},AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}],
Plot3D[Min[1,modelOrenEo/.fitOrenEo],{x,0,1},{y,0,1},PlotStyle->Red]
]

Disney Diffuse BRDF

From Disney’s principled BRDF
https://disney-animation.s3.amazonaws.com/library/s2012_pbs_disney_brdf_notes_v2.pdf

(Subscript[f, d])=(ρ/π)(1-F((Subscript[θ, l])))(1-F((Subscript[θ, v])))
	1-F(θ)=1+((Subscript[F, D90])-1)(((1-cos(θ)))^5)	Schlick Fresnel
	(Subscript[F, D90])=0.5+2·(cos^2)((Subscript[θ, d]))·α
	α = roughness

Clear[α]

furnaceFactor=0.9388180685312748`
disney=Function[ {θo,ϕo,θi,ϕi,α},
ϕ=ϕi-ϕo; (* Only angle difference is relevant *)
NdotL=Cos[θi];
NdotV=Cos[θo];
VdotL=Cos[θi]*Cos[θo]+Sin[θi]*Sin[θo]*Cos[ϕ];
LdotH=((1+VdotL)/(Sqrt[2+2*VdotL])); (* = (Subscript[ω, i])·(((Subscript[ω, i])+(Subscript[ω, o]))/(|(Subscript[ω, i])+(Subscript[ω, o])|))*)
Fd90=0.5+2*LdotH*LdotH*α;
(furnaceFactor/π)*(1+(Fd90-1)*(((1-NdotL))^5))*(1+(Fd90-1)*(((1-NdotV))^5))
];
disney[0,0,0.1,0.2,0.2]

Manipulate[Plot3D[disney[θo,0,x,0,α],{x,0,(π/2)},{α,0,1},PlotRange->{0,1},AxesLabel->{"cos(θ)","α","Disney"}],{{θo,(π/2)},0,(π/2)}]

Let’s try integrating the Lambertian case:

fTest[α_]:= 2*π*NIntegrate[disney[θo,0,θi,ϕi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)}]
fTest[0]
fTest[1]
Plot[fTest[α],{α,0,1},PlotRange->{0,3.5}]

furnaceFactor=(π/(3.3463274290242717`))
furnaceFactor=0.9388180685312748`

disneyEo[θo_,α_]:=NIntegrate[disney[θo,0,θi,ϕi,α] * Cos[θi]* Sin[θi],{θi,0,(π/2)},{ϕi,-π,π}]
disneyEo[1.5,1]

Plot[disneyEo[θ,1],{θ,0,(π/2)},PlotRange->{0,1.4}]

2*π*NIntegrate[disney[θo,0,θi,ϕi,1] * Cos[θi]* Sin[θi]* Cos[θo]* Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)}]

SetDirectory["D:\\Workspaces\\GodComplex\\Tests\\TestMSBRDF\\Tables"]; 
countX = 64;
countY= 64;
tableFileName="MSBRDF_Disney_E"  <>ToString[countX] <>"x" <> ToString[countY] <>".csv"
tableFileNameAvg="MSBRDF_Disney_Eavg"  <>ToString[countY] <>".csv"

buildTable=Function[{i},
x=N[((Mod[i,countX])/(countX-1))];
y=N[((Quotient[i,countX])/(countY-1))];
θo=ArcCos[x];
α=y;
{x,y,disneyEo[θo,α]}
]
buildTable[countX*5+7]

resultDisneyEo= Table[buildTable[i],{i,0,countX*countY-1}];
Export[tableFileName,resultDisneyEo,"CSV"]

Export[tableFileName,resultDisneyEo,"CSV"]

ListPlot3D[resultDisneyEo,PlotRange->{0,2},AxesLabel->{"cos((Subscript[θ, o]))","α"}]
Image[Table[Table[1-resultDisneyEo[[1+countX*y+x]][[3]],{x,0,countX-1}],{y,0,countY-1}]]

Re-Importing and testing...

resultDisneyEo = Import[tableFileName,"CSV"];
ListPlot3D[resultDisneyEo,AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}]

Let’s reorder the table and create an interpolation function from it:

t=Table[resultDisneyEo[[1+countX*(y-1)+(x-1)]][[3]],{x,1,countY},{y,1,countY}];
fDisneyEo=ListInterpolation[t,{{0,1},{0,1}}]
Plot3D[fDisneyEo[x,y],{x,0,1},{y,0,1},AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}]

Computing Average Irradiance

Now let’s integrate the full white furnace irradiance using our table:

	(Subscript[E, avg])(α)=∫(Subscript[E, o])((Subscript[θ, o]),α) cos((Subscript[θ, o])) d(Subscript[θ, o])

resultDisneyEavg=Table[{N[(i/(countY-1))],2*π*NIntegrate[fDisneyEo[μ,(i/(countY-1))]μ,{μ,0,1}]},{i,0,countY-1}]
Export[tableFileNameAvg,resultDisneyEavg,"CSV"]

resultDisneyEavg=Import[tableFileNameAvg,"CSV"]

ListLinePlot[resultDisneyEavg,PlotRange->{0,π}]

t=Table[resultDisneyEavg[[i]][[2]],{i,1,countY}];
fDisneyEavg=ListInterpolation[t,{0,1}]
Plot[fDisneyEavg[x],{x,0,1},PlotRange->{0,π}, AxesLabel->{"α","(Subscript[E, avg])"}]

Clear[a,b,c,d,x]
fittingMethod="LevenbergMarquardt";
fittingMethod="PrincipalAxis";
fittingMethod="Automatic";

modelDisneyEavg=π+0*a+(0+0*b)*x+c*(x^2)+d*(x^3)+e*(x^4);
fitDisneyEavg=FindFit[resultOrenEavg,modelDisneyEavg,{a,b,c,d,e},x,Method->fittingMethod]
modelDisneyEavg/.fitDisneyEavg
funcDisneyEavg=Function[x,Evaluate[modelDisneyEavg/.fitOrenEavg]];
Show[
Plot[fDisneyEavg[x],{x,0,1},PlotRange->{0,π},AxesLabel->{"α","(Subscript[E, avg])"}],
Plot[funcDisneyEavg[x],{x,0,1},PlotStyle->Red]
]

Let’s check the BRDFs sum to π!

BRDFdisney[θo_,θi_,ϕ_,α_]:= disney[θo,0,θi,ϕ,α]
BRDFdisneyms[θo_,θi_,ϕ_,α_]:=(((1-fDisneyEo[Cos[θo],α])*(1-fDisneyEo[Cos[θi],α]))/(π-fDisneyEavg[α]))

whiteFurnaceSS[α_]:=2*π*NIntegrate[BRDFdisney[θo,θi,ϕi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]
whiteFurnaceMS[α_]:=2*π*NIntegrate[BRDFdisneyms[θo,θi,ϕi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]

whiteFurnaceSS[0.3]
whiteFurnaceMS[0.3]

whiteFurnaceResultSS=Table[{α,whiteFurnaceSS[α]},{α,0,1,0.05}]
whiteFurnaceResultMS=Table[{α,whiteFurnaceMS[α]},{α,0,1,0.05}]

Show[
ListLinePlot[whiteFurnaceResultSS,PlotRange->{0,1.1*π},PlotStyle->Gray,AxesLabel->{"α",""}],
ListLinePlot[whiteFurnaceResultMS,PlotStyle->Blue],
ListLinePlot[whiteFurnaceResultSS+whiteFurnaceResultMS,PlotRange->Full,PlotStyle->Red]
]

Fresnel Formulations

Classical dielectric Fresnel formulation from Walter / Cook&Torrance.
Note: this formulation is perfectly fine for metallic materials as long as the (Subscript[F, 90]) factor is 1, as can be seen at the bottom of this page.

fresnel=Function[{c,F0},
η=((1+(Sqrt[F0]))/(1.0000001-(Sqrt[F0])));
g=(Sqrt[(η^2)-1+(c^2)]);
(1/2)(((((g-c)/(g+c))))^2)(1+(((((c*(g+c)-1)/(c*(g-c)+1))))^2))
];
F0=1.0;
Manipulate[Plot[fresnel[Cos[θ],F0],{θ,0,(π/2)},PlotRange->{0,1.1}],{{F0,1},0,1}]

Half-Vector formulation of Fresnel, before I understand my way of thinking is a load of bullshit... Actually, cos((Subscript[θ, h])) is not the cosine that is used in the Fresnel formulation of BRDFs but we must use cos((Subscript[θ, d])) instead!
That’s because, according to Walter, the perfect mirror condition appears when (Subscript[ω, h])=(Subscript[ω, m]), i.e. the half-vector is perfectly aligned with the microscopic surface’s micro-facet. Then the reflectance of the incoming light due to Fresnel only depends on the angle between the incoming light vector (Subscript[ω, i]) and the micro-facet normal (Subscript[ω, m])=(Subscript[ω, h]) so it’s really dependent on (Subscript[ω, i])·(Subscript[ω, h])=cos((Subscript[θ, d]))...

fresnelH=Function[{θo,θi,F0},
lengthH=(Sqrt[2+2*Cos[θo-θi]]);
cosH=((Cos[θo]+Cos[θi])/lengthH);
fresnel[cosH,F0]
];
Manipulate[Plot[fresnelH[θo,θi,F0],{θi,((-π)/2),(π/2)},PlotRange->{0,1}],{{θo,0},((-π)/2),(π/2)},{{F0,1},0,1}];
Manipulate[Plot[fresnel[Cos[θo],F0]*fresnel[Cos[θi],F0],{θi,((-π)/2),(π/2)},PlotRange->{0,1}],{{θo,0},((-π)/2),(π/2)},{{F0,1},0,1}];

Manipulate[
Show[
Plot[fresnelH[θo,θi,F0],{θi,((-π)/2),(π/2)},PlotRange->{0,1}],
Plot[fresnel[Cos[θo],F0]*fresnel[Cos[θi],F0],{θi,((-π)/2),(π/2)},PlotRange->{0,1},PlotStyle->Red]
]
,{{θo,0},((-π)/2),(π/2)},{{F0,0.8},0,1}];

Artist-Friendly Metal Fresnel formulation, from Ole Gulbrandsen.
• Mostly similar to dielectric Fresnel except at grazing angles (use dielectric Fresnel all the way if you really always have (Subscript[F, 90])=1)
• Fails to properly account for IOR=1 and diverges from dielectric Fresnel for (Subscript[F, 0])<0.2 (you should revert to using dielectric Fresnel in these cases anyway)

fresnelMetal=Function[{c,F0,F90},
r=F0;
g=F90;
c2=(c^2);

(* Compute n and k *)
sqrtR=(Sqrt[r]);
nmin=((1-r)/(1+r));
nmax= ((1+sqrtR)/(1-sqrtR));
n=nmin*(1-g)+nmax*g;

n2=(n^2);
nr=(((n+1))^2)*r-(((n-1))^2);
k2=(nr/(1-r));

(* Compute perpendicular polarized Fresnel *)
numPe=n2+k2-2*n*c+c2;
denPe=n2+k2+2*n*c+c2;
Pe=(numPe/denPe);

(* Compute parallel polarized Fresnel *)
numPa=(n2+k2)*c2-2*n*c+1;
denPa=(n2+k2)*c2+2*n*c+1;
Pa=(numPa/denPa);

0.5*(Pe+Pa)
];
fresnelMetal[0.2,0.8,0.95];

fresnelSchlick=Function[{c,F0},
F0+(1-F0)*(((1-c))^5)
];

Here we show the comparison between dielectric and metallic Fresnel expressions (Schlick is given to show how it differs to more accurate versions):

Manipulate[
Show[
Plot[fresnelSchlick[Cos[θ],F0],{θ,0,(π/2)},PlotRange->{0,1.1},PlotStyle->Green]
,Plot[fresnelMetal[Cos[θ],F0,F90],{θ,0,(π/2)},PlotRange->{0,1}]
,Plot[fresnel[Cos[θ],F0],{θ,0,(π/2)},PlotRange->{0,1.1},PlotStyle->Red]
]
,{{F0,0.0 },0,0.999},{{F90,1},0,1}]