BRDF Definitions

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])))
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])
];

Tables Loading

GGX

SetDirectory["D:\\Workspaces\\GodComplex\\Tests\\TestMSBRDF\\Tables"]; 
countX=128;
countY=128;
resultEo = Import["MSBRDF_GGX_E128x128.csv","CSV"];
resultEavg=Import["MSBRDF_GGX_Eavg128.csv","CSV"];
fEo=ListInterpolation[Table[resultEo[[1+countX*(y-1)+(x-1)]][[3]],{x,1,countY},{y,1,countY}],{{0,1},{0,1}}]
fEavg=ListInterpolation[Table[resultEavg[[i]][[2]],{i,1,countY}],{0,1}]
Plot[fEavg[x],{x,0,1},PlotRange->{0,π}, AxesLabel->{"α","(Subscript[E, avg])"}];
Plot3D[fEo[x,y],{x,0,1},{y,0,1},AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}];

BRDFggx[θo_,θi_,ϕi_,α_]:= smith[θo,α]*smith[θi,α]*ggx[change[θo,θi,ϕi],α]
BRDFms[θo_,θi_,α_]:=(((1-fEo[Cos[θo],α])*(1-fEo[Cos[θi],α]))/(π-fEavg[α]))
fWhiteFurnaceSS[α_]:=2*π*NIntegrate[BRDFggx[θo,θi,ϕi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]
fWhiteFurnaceMS[α_]:=2*π*NIntegrate[BRDFms[θo,θi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]
fWhiteFurnaceSS[0.3]
fWhiteFurnaceMS[0.3]
π-fEavg[0.3]
fWhiteFurnaceSS[0.3]+fWhiteFurnaceMS[0.3]

Oren-Nayar

SetDirectory["D:\\Workspaces\\GodComplex\\Tests\\TestMSBRDF\\Tables"]; 
countX=32;
countY=32;
resultOrenEo = Import["MSBRDF_OrenNayar_E32x32.csv","CSV"];
resultOrenEavg = Import["MSBRDF_OrenNayar_Eavg32.csv","CSV"];
fOrenEo=ListInterpolation[Table[resultOrenEo[[1+countX*(y-1)+(x-1)]][[3]],{x,1,countY},{y,1,countY}],{{0,1},{0,1}}]
fOrenEavg=ListInterpolation[Table[resultOrenEavg[[i]][[2]],{i,1,countY}],{0,1}]
Plot[fOrenEavg[x],{x,0,1},PlotRange->{0,π}, AxesLabel->{"α","(Subscript[E, avg])"}];
Plot3D[fOrenEo[x,y],{x,0,1},{y,0,1},AxesLabel->{"cos(θ)", "α","(Subscript[E, i])"}];

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])
];
BRDForen[θo_,θi_,ϕ_,α_]:= oren[θo,0,θi,ϕ,α]
BRDForenms[θo_,θi_,ϕ_,α_]:=(((1-fOrenEo[Cos[θo],α])*(1-fOrenEo[Cos[θi],α]))/(π-fOrenEavg[α]))
fWhiteFurnaceSS[α_]:=2*π*NIntegrate[BRDForen[θo,θi,ϕi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]
fWhiteFurnaceMS[α_]:=2*π*NIntegrate[BRDForenms[θo,θi,ϕi,α]*Cos[θi]*Sin[θi]*Cos[θo]*Sin[θo],{θi,0,(π/2)},{ϕi,-π,π},{θo,0,(π/2)},WorkingPrecision->MachinePrecision]
fWhiteFurnaceSS[0.3]
fWhiteFurnaceMS[0.3]
π-fOrenEavg[0.3]
fWhiteFurnaceSS[0.3]+fWhiteFurnaceMS[0.3]

SH Environment

We start by writing the expression we need to fit for various roughness values:

	(∫_(SuperPlus[Ω])^upper)L((Subscript[ω, i]))  f((Subscript[ω, o]),(Subscript[ω, i])) ((Subscript[ω, i])· n)dω  =(∫_(SuperPlus[Ω])^upper)L((Subscript[ω, i]))  (((1-E((Subscript[ω, o]),α))(1-E((Subscript[ω, i]),α)))/(π-(Subscript[E, avg])(α))) ((Subscript[ω, i])· n)dω
				=(∫_(SuperPlus[Ω])^upper)((Subscript[L, lm]) · (Subscript[Y, lm])((Subscript[ω, i])))  (((1-E((Subscript[ω, o]),α))(1-E((Subscript[ω, i]),α)))/(π-(Subscript[E, avg])(α))) ((Subscript[ω, i])· n)dω
				= (Subscript[L, lm])  ((1-E((Subscript[ω, o]),α))/(π-(Subscript[E, avg])(α)))(∫_(SuperPlus[Ω])^upper) (Subscript[Y, lm])((Subscript[ω, i]))(1-E((Subscript[ω, i]),α))((Subscript[ω, i])· n)dω 
Where (Subscript[Y, lm])((Subscript[ω, i])) is a spherical harmonics coefficient in the direction (Subscript[ω, i]).
We notice that the integral is only dependent on roughness α and we will attempt to compute it for all values of (l,m) and various values of roughness to finally attempt to find a nice analytical fit for the coefficients for each BRDF.

SphericalPlot3D[Abs[Re[SphericalHarmonicY[1,-1,θ,ϕ]]],{θ,0,Pi},{ϕ,0,2Pi}]

GGX Case

fSHBRDFggx= Function[ {i,θ,ϕ,α}, l=Floor[Sqrt[i]];m=i-l*(l+1); Re[SphericalHarmonicY[l,m,θ,ϕ]]*(1-fEo[Cos[θ],α])];
NIntegrate[fSHBRDFggx[2,θ,ϕ,0.5]*Cos[θ]*Sin[θ],{θ,0,(π/2)},{ϕ,-π,π}]

(* Integration method & accuracy test *)
NIntegrate[fSHBRDFggx[1,θ,ϕ,α]*Cos[θ]*Sin[θ],{θ,0,(π/2)},{ϕ,-π,π},Method->"AdaptiveMonteCarlo",MaxRecursion->4, AccuracyGoal->4]

α=1
Table[NIntegrate[fSHBRDFggx[i,θ,ϕ,α]*Cos[θ]*Sin[θ],{θ,0,(π/2)},{ϕ,-π,π},Method->"AdaptiveMonteCarlo",MaxRecursion->4, AccuracyGoal->4],{i,0,8}]

Here I realized only the ZH coefficients yield non-zero values and the rest are useless! No need to bother computing all of them... :)

taResultSH=Table[Table[{α,NIntegrate[fSHBRDFggx[i,θ,ϕ,α]*Cos[θ]*Sin[θ],{θ,0,(π/2)},{ϕ,-π,π},Method->"AdaptiveMonteCarlo",MaxRecursion->4, AccuracyGoal->4]},{α,0.05,1,0.05}],{i,0,8}]

We easily notice only the ZH coefficients are non zero, which is even better!

taResultSHGGX={taResultSH[[1]], taResultSH[[3]],taResultSH[[7]]}

Show[Table[ListLinePlot[taResultSHGGX[[i]],PlotRange->{0,0.7},PlotStyle->Hue[(i-1)/3]],{i,1,3}]]

ZHindex={0,2,6};
taResultSHGGX=Table[Table[{α,2*π*NIntegrate[fSHBRDFggx[ZHindex[[i]],θ,0,α]*Cos[θ]*Sin[θ],{θ,0,(π/2)},Method->"AdaptiveMonteCarlo",MaxRecursion->8, AccuracyGoal->6]},{α,0.05,1,0.05}],{i,1,3}];
Show[Table[ListLinePlot[taResultSHGGX[[i]],PlotRange->{0,0.7},PlotStyle->Hue[(i-1)/3],AxesLabel->{"α","(Subscript[E, l])"},PlotLabel->"GGX"],{i,1,3}]
]

Oren-Nayar Case

fSHBRDForen= Function[ {i,θ,ϕ,α}, l=Floor[Sqrt[i]];m=i-l*(l+1); Re[SphericalHarmonicY[l,m,θ,ϕ]]*(1-fOrenEo[Cos[θ],α])];
NIntegrate[fSHBRDForen[2,θ,ϕ,0.5]*Cos[θ]*Sin[θ],{θ,0,(π/2)},{ϕ,-π,π}]
NIntegrate[2*π*fSHBRDForen[2,θ,0,0.5]*Cos[θ]*Sin[θ],{θ,0,(π/2)}]

ZHindex={0,2,6};
taResultSHoren=Table[Table[{α,2*π*NIntegrate[fSHBRDForen[ZHindex[[i]],θ,0,α]*Cos[θ]*Sin[θ],{θ,0,(π/2)},Method->"AdaptiveMonteCarlo",MaxRecursion->8, AccuracyGoal->6]},{α,0.05,1,0.05}],{i,1,3}];
Show[Table[ListLinePlot[taResultSHoren[[i]],PlotRange->{0,0.4},PlotStyle->Hue[(i-1)/3],AxesLabel->{"α","(Subscript[E, l])"},PlotLabel->"Oren-Nayar"],{i,1,3}]
]

Fitting

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

model=a*(1-(b*x+c*(x^2)+d*(x^3)));
model=a*(x^(1/2))+b*(x^(3/2))+c*(x^(5/2))+d*(x^(7/2));
model=a*(x^(1/2))+b*(x^(3/2))+c*(x^(5/2));
fitGGX=Table[FindFit[taResultSHGGX[[i]],model,{a,b,c,d},x,Method->fittingMethod],{i,1,3}]
funcSHGGX=Table[Function[x,Evaluate[model/.fitGGX[[i]]]],{i,1,3}]
Show[
Table[ListLinePlot[taResultSHGGX[[i]],PlotRange->{0,0.7},AxesLabel->{"α","(Subscript[E, 0])"}],{i,1,3}],
Table[Plot[funcSHGGX[[i]][x],{x,0,1},PlotStyle->Red],{i,1,3}]
]

model=0*a+b*x+c*(x^2)+d*(x^3)+e*(x^4);
model=a*(x^(1/2))+b*(x^(3/2))+c*(x^(5/2))+d*(x^(7/2));
fitOren=Table[FindFit[taResultSHoren[[i]],model,{a,b,c,d},x,Method->fittingMethod],{i,1,3}]
funcSHOren=Table[Function[x,Evaluate[model/.fitOren[[i]]]],{i,1,3}]
Show[
Table[ListLinePlot[taResultSHoren[[i]],PlotRange->{0,0.4},AxesLabel->{"α","(Subscript[E, 0])"}],{i,1,3}],
Table[Plot[funcSHOren[[i]][x],{x,0,1},PlotStyle->Red],{i,1,3}]
]

{a->-0.0919559140506979`,b->1.467037714315657`,c->-1.6735448883797404`,d->0.607800523815945`},{a->-0.11366841288600084`,b->1.901273744271233`,c->-2.322322430339633`,d->0.909815621695672`},{a->-0.04124821752212919`,b->1.0933549500536324`,c->-1.4171919237898756`,d->0.5810844359893627`}