Simulated Multiple-Scattering Fresnel Factor

The idea here is to use my MSBRDF simulator to cast many rays on a tiny surface and measure the total energy that comes out for several scattering orders (up to 6, the 7th order stores the total of the remaining 6+ orders).

I ran the simulation for 16 incidence angles from 0 to 90°, 16 surface roughness values and 20 (Subscript[F, 0]) values from 1 (IOR=∞) down to 0.04 (IOR= 1.5, the classical dielectric value used in many 3D engines).

Theoretically, total energy for order 1 should be roughly the energy of a Cook-Torrance BRDF (because my simulator uses a heightfield with a Beckmann distribution of heights, not a GGX distribution)(using the code from Heitz).
The sum of the energy output by all other orders should be equivalent to our magical multiple-scattering term:

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

The goal then is to check that it is indeed the case.
Then, we must decompose the numerical simulation of orders 2+ to understand what are their general shape (if any can be found, but an approximation will suffice) and what law governs the dependence on (Subscript[F, 0]), typically we expect order 1 to be directly linked to (Subscript[F, 0]) while order 2 should be linked to ((Subscript[F, 0])^2) and subsequent order N to ((Subscript[F, 0])^N).

Eventually, we’re looking for some expression that would look like this:

	(Subscript[E, avg]^j)((Subscript[θ, i]),α,(Subscript[F, 0]))= f((Subscript[θ, i]),α)·g(j)·((h((Subscript[F, 0])))^j)

f((Subscript[θ, i]),α) is some general expression we don’t really care about since the (∑_(j=2)^∞)(Subscript[E, avg]^j)((Subscript[θ, i]),α,(Subscript[F, 0])=1)=π-(Subscript[E, avg])(α)
g(j) is some amplitude function depending on the order j
h((Subscript[F, 0])) is some mapping from (Subscript[F, 0]) to an expression that can be raised to the power of the current order

SetDirectory["D:\\Workspaces\\GodComplex\\Tests\\TestMSBRDF\\Tables"]; 
countTheta = 16;
countRoughness= 16;
countF0=20;
countOrders=7;
taRawTotalReflectances=Table[
Import["TotalReflectance_"  <>ToString[countTheta] <>"x"  <>ToString[countRoughness] <>"x"  <>ToString[countF0] <>"_order"  <>ToString[order] <>".bin", "Real64"],{order,1,countOrders}];

Generate list interpolations

fTotalReflectance=Table[
ListInterpolation[Table[
taRawTotalReflectances[[order]][[countTheta*((roughness-1)+countRoughness*(albedo-1))+theta]],{theta,1,countTheta},{roughness,1,countRoughness},{albedo,1,countF0}],{{0,(π/2)},{1,0},{1,0.04}}],
{order,1,countOrders}];

Manipulate[
Plot3D[fTotalReflectance[[order]][x,α,ρ],{x,0,(π/2)},{α,0,1},PlotRange->{0,(3^(-order+1))}],{{ρ,1},0.04,1},{{order,2},1,countOrders,1}]

Show Total Reflectance Depending on (Subscript[F, 0])

Manipulate[Plot3D[Sum[fTotalReflectance[[order]][x,α,ρ],{order,1,countOrders}],{x,0,(π/2)},{α,0,1},PlotRange->{0,1},AxesLabel->{"θ","α"}],{{ρ,1},0.04,1}]

Manipulate[Plot3D[Sum[fTotalReflectance[[order]][x,α,ρ],{order,2,countOrders}],{x,0,(π/2)},{α,0,1},PlotRange->{0,0.3},AxesLabel->{"θ","α"}],{{ρ,1},0.04,1}]

Visualizing each individual order

Manipulate[
Show[
Table[
Plot3D[fTotalReflectance[[order]][x,α,ρ],{x,0,(π/2)},{α,0,1},
PlotRange->{0,0.4},PlotStyle->Hue[((4*(order-2))/countOrders)]],{order,2,countOrders}]
],
{{ρ,1},0.04,1}]

White Furnace Comparison

SetDirectory["D:\\Workspaces\\GodComplex\\Tests\\TestMSBRDF\\Tables"]; 
countX=128;
countY=128;
resultEo = Import["MSBRDF_E128x128.csv","CSV"];
resultEavg=Import["MSBRDF_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]

whiteFurnaceMS=Table[{α,π-fEavg[α]},{α,0,1,0.05}]

Here we show that we mustn’t forget that the MS white furnace amounts to the complement of the total energy lost by SS BRDF...

Show[ListPlot[whiteFurnaceMS,PlotRange->Full, AxesLabel->{"α","π-(Subscript[E, avg])"}],
Plot[π-fEavg[x],{x,0,1}]
]

Comparing with Simulation

We computed the total energy bouncing off the surface at various scattering orders for different values of roughness, incidence angle and albedos.
So basically we performed:

	E((Subscript[ω, i]),α,ρ)=(∫_(SuperPlus[Ω])) (Subscript[L, i])((Subscript[ω, i])) (Subscript[f, r])( (Subscript[ω, o]),(Subscript[ω, i]),α,ρ)((Subscript[ω, o])·n)(Subscript[dω, o])
	
With (Subscript[L, i])((Subscript[ω, i]))=1
To compute the equivalent white furnace test we thus need to compute:

	(Subscript[E, total])(α,ρ)=(∫_(SuperPlus[Ω]))E((Subscript[ω, i]),α,ρ)((Subscript[ω, i])·n)(Subscript[dω, i])

taRawSumTotalReflectances=Sum[taRawTotalReflectances[[order]],{order,2,countOrders}]; (* Sum all orders > 2 *)
MSTotalReflectances=ListInterpolation[Table[taRawSumTotalReflectances[[countTheta*((roughness-1)+countRoughness*(F0-1))+theta]],{theta,1,countTheta},{roughness,1,countRoughness},{F0,1,countF0}],{{0,(π/2)},{1,0},{1,0.04}}]
fWhiteFurnaceSimul[α_,ρ_,order_]:=2*π*NIntegrate[fTotalReflectance[[order]][θi,α,ρ]*Cos[θi]*Sin[θi],{θi,0,(π/2)},WorkingPrecision->MachinePrecision];
fSumWhiteFurnaceSimul[α_,ρ_]:=2*π*NIntegrate[MSTotalReflectances[θi,α,ρ]*Cos[θi]*Sin[θi],{θi,0,(π/2)},WorkingPrecision->MachinePrecision];

Plot3D[fSumWhiteFurnaceSimul[α,ρ],{α,0,1},{ρ,0.04,1},PlotRange->{0,1},AxesLabel->{"α","ρ"}]

Here we show that the white furnace tests for the simulated data also complement to π

ρ=1;
Show[
Plot[fWhiteFurnaceSimul[α,ρ,1],{α,0,1},PlotRange->{0,π},AxesLabel->{"α","(Subscript[E, avg])"}]
,Plot[fSumWhiteFurnaceSimul[α,ρ],{α,0,1},PlotRange->Full,PlotStyle->Red]
,Plot[π,{α,0,1},PlotRange->Full,PlotStyle->Green]
]

factor=((π-fEavg[1])/(fSumWhiteFurnaceSimul[1,1]))
Show[
ListLinePlot[whiteFurnaceResult1,AxesLabel->{"α","(Subscript[E, avg])"}]
,Plot[factor*fSumWhiteFurnaceSimul[α,ρ],{α,0,1},PlotStyle->Red]
]

We see here the red curve is not quite the same as the white furnace result given by the theoretical multiple-scattering expression for the GGX BRDF. The global factor 2.1 is quite worrying but otherwise, the difference in shape is certainly due to the GGX distribution of normals that is not the same as the Beckmann distribution I used in the simulation.
Overall, we won’t linger too much on these differences and assume the simulation represents ground-truth data.
We are more interested in the relative behavior of individual orders one to another and the geometric progression that is expected from varying the (Subscript[F, 0]) term.

Finding an Analytical Model

We saw that integrating our order (SuperPlus[2]) irradiance values (the red curve in the above graph) yields π-(Subscript[E, total]) where (Subscript[E, total])(α,ρ) is the total irradiance reflected for the single-scattering Cook-Torrance BRDF (NOT GGX BRDF, because I don’t know how to generate a simulation height field with a Smith distribution of heights) (I should re-read Heitz papers but I remember it’s quite complicated, I’m not even sure he provided the code for that).

We are interested in how each order contributes to that complement.
We assume the roughness-dependent shape is not important (i.e. it’s computed by the analytical (Subscript[E, avg]) expression for the GGX MSBRDF) and concentrate on the maximum values of the white furnace functions for each order, when α=1.
But we are interested in the evolution of that maximum value when varying the albedo!

So let’s create a table of samples of white furnace values for varying albedos, for each order (except last order, which is a compound of all remaining orders and doesn’t interest us):

α=1;
taWhiteFurnace=
Table[fWhiteFurnaceSimul[α,ρ,order],{order,2,countOrders},{ρ,0.04,1,0.02}];

fWhiteFurnace=Table[ListInterpolation[taWhiteFurnace[[order-1]],{0.04,1}],{order,2,countOrders}];
Show[Table[Plot[fWhiteFurnace[[order-1]][ρ],{ρ,0.04,1},PlotRange->Full,PlotStyle->Hue[((order-2)/countOrders)]],{order,2,countOrders}]]
Show[Table[Plot[Log[fWhiteFurnace[[order-1]][ρ]],{ρ,0.04,1},PlotRange->Full,PlotStyle->Hue[((order-2)/countOrders)]],{order,2,countOrders-1}]]

Now let’s try and fit each curve with a polynomial of the form f((Subscript[F, 0]))=(Subscript[a, n])*((Subscript[F, 0])^n) where n is the order:

An=Table[fWhiteFurnace[[order-1]][1],{order,2,countOrders}]
ListPlot[An,PlotRange->{{0,7},{0,0.13}}]

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

exponent=2;
fitter=Function[{order},
model=An[[order-1]]*(x^(a*(order^exponent)));
data=Table[{ρ,fWhiteFurnace[[order-1]][ρ]},{ρ,0.04,1,0.01}];
fit=FindFit[data,model,{{a,1},{b,0},{c,1},{d,0},e},x,Method->fittingMethod]
];
fits=Table[Function[x,Evaluate[An[[order-1]]*(x^(a*(order^exponent)))/.fitter[order]]],{order,2,countOrders-1}];
fitExponents=Table[a/.fitter[order],{order,2,countOrders-1}]
ListPlot[fitExponents,PlotRange->{0,1}]

Great! The exponent factor seems to be (1/2)! ^_^

Kn=Table[a/.fitter[order],{order,2,countOrders-1}]
avgK=((Sum[a/.fitter[order],{order,2,countOrders-1}])/5)

Manipulate[
Show[
Plot[fWhiteFurnace[[order-1]][ρ],{ρ,0.04,1},PlotRange->{0,fWhiteFurnace[[order-1]][1]},AxesLabel->{"(Subscript[F, 0])","(Subscript[E, avg])"}]
,Plot[fits[[order-1]][ρ],{ρ,0.04,1},PlotStyle->Red,PlotRange->Full]
],{order,2,countOrders-1,1}]

So we found an expression of the form:
	f((Subscript[F, 0]))=(Subscript[a, n]) · ((((Sqrt[(Subscript[F, 0])])))^(n^2))

Now we will try to make (Subscript[a, n]) enter this pattern as well by writing:
	(Subscript[a, n])=(Subscript[a, 1])·((((Sqrt[τ])))^(n^2))

So that τ can be factored into (Subscript[F, 0])...

An

model=An[[1]]*a*((((Sqrt[b])))^(x^2));
model=a*((((Sqrt[b])))^(x^2));
data=Table[{order,An[[order-1]]},{order,2,countOrders-1}];
fit=FindFit[data,model,{{a,1},{b,1},{c,1},{d,0},e},x,Method->fittingMethod]
fTau=Function[{x},Evaluate[model/.fit]];
Plot[fTau[x],{x,2,countOrders-1},PlotStyle->Red,PlotRange->Full];
Show[
ListPlot[data,PlotRange->Full]
,Plot[fTau[x],{x,2,countOrders-1},PlotStyle->Red,PlotRange->Full]
]
a2=(a/.fit) (* This can be used to fit experimental data *)
(a/.fit)*((((Sqrt[b/.fit])))^4)
3.8655810610440082`*((((Sqrt[0.4554998477195002`])))^4)
3.865581061044009`
0.45549984771950014`

So we finally get our final expression:

	(Subscript[a, n])((Subscript[F, 0]))= (Subscript[a, 1]) ((Subscript[(OverHat[F]), 0])^(n^2))
	(Subscript[(OverHat[F]), 0])=(Sqrt[τ · (Subscript[F, 0])])
	τ=0.45549984771950014
	
We require that:
	(∑_(n=2)^∞)(Subscript[a, n])(1)=(∑_(n=2)^∞)(Subscript[a, 1]) ((((Sqrt[τ])))^(n^2))=(Subscript[a, 1])·(∑_(n=2)^∞) ((((Sqrt[τ])))^(n^2))=1
	
	(Subscript[a, 1])=(1/((∑_(n=2)^∞) ((((Sqrt[τ])))^(n^2))))

The geometric series is computed using the Jacobi elliptic theta 3 function that converges to:

	(∑_(n=2)^∞) ((((Sqrt[τ])))^(n^2))=0.4168994928925871 (numerically computed below)

And so we get:

	(Subscript[a, 1]) = (1/((∑_(n=2)^∞) ((((Sqrt[τ])))^(n^2))))= 4.193906008161867

τ =0.45549984771950014
Sum[((((Sqrt[τ])))^(n^2)),{n,2,∞}]
a1=(1/(Sum[((((Sqrt[τ])))^(n^2)),{n,2,∞}]))
Plot[Sum[((((Sqrt[τ])))^(n^2)),{n,2,maxN}],{maxN,2,10},PlotRange->Full];

4.193906008161867`
4.193906008161867`

With the “Elliptic Theta” function from Mathematica:
	(Subscript[ϑ, 3])(u,q)=1+2(∑_(n=1)^∞)(q^(n^2))cos(2n u)

If we let u=0 then it simplifies into:
	(Subscript[ϑ, 3])(0,q)=1+2(∑_(n=1)^∞)(q^(n^2))
	
We are interested in:
	(∑_(n=1)^∞)(q^(n^2))= (((Subscript[ϑ, 3])(0,q)-1)/2)
	
Furthermore, we are looking for the series starting from n=2 so:

	(∑_(n=2)^∞)(q^(n^2))=(∑_(n=1)^∞)(q^(n^2))- q= (((Subscript[ϑ, 3])(0,q)-1)/2)-q= (((Subscript[ϑ, 3])(0,q)-1-2q)/2)

(Sqrt[τ])
recA1=((EllipticTheta[3,0, (Sqrt[τ])]-1-2(Sqrt[τ]))/2)
a1=1/recA1 (* = a1 *)
Plot[EllipticTheta[3,0,(Sqrt[ x])],{x,0,τ}]

4.193906008161867`

We check the final equation sums to 1:

a1
Sum[a1*((Sqrt[τ])^(n^2)),{n,2,countOrders}]

Same with the elliptic theta:

fEllipticTheta[F0_]:=a1*(((EllipticTheta[3,0, (Sqrt[τ*F0])]-1-2(Sqrt[τ*F0]))/2))
fEllipticTheta[1]

Plot[Sum[a1*((Sqrt[τ*F0])^(n^2)),{n,2,100}],{F0,0,1}]
Plot[fEllipticTheta[F0],{F0,0,1},PlotRange->Full]

Now we try to fit the elliptic curve with a simple polynomial:

model=0*a+b*x+c*(x^2)+d*(x^3);
data=Table[{F0,fEllipticTheta[F0]},{F0,0,1,0.01}];
fit=FindFit[data,model,{{a,1},{b,1},{c,1},{d,0},e},x,Method->fittingMethod]
fEllipticPoly=Function[{x},Evaluate[model/.fit]];
Plot[fEllipticPoly[x],{x,0,1},PlotStyle->Red,PlotRange->{0,1}];
Show[
Plot[fEllipticTheta[F0],{F0,0,1},PlotRange->{0,1}]
,Plot[fEllipticPoly[F0],{F0,0,1},PlotStyle->Red,PlotRange->Full]
]

Final Expression

Clear[factorMS]
τ=0.4554998477195002
a1=(2/(EllipticTheta[3,0, (Sqrt[τ])]-1-2(Sqrt[τ])))
factorMS[F0_]:=F0*(0.039324148543732715`+F0*(0.6563284963169776` + F0*0.2995566941543759`))
factorMS[F0_]:=F0*(0.04+F0*(0.66+ F0*0.3))
factorMS[1]
Show[
Plot[fEllipticTheta[F0],{F0,0,1},PlotRange->{0,1},AxesLabel->{"(Subscript[F, 0])","(Subscript[F, ms])"}]
,Plot[factorMS[ρ],{ρ,0,1},PlotRange->Full,PlotStyle->Red]
]

Verification

Now we must verify that our expression can be substituted to the actual measurements!
We first verify that white furnace results match for each individual orders (Subscript[f, n])((Subscript[F, 0])) * whiteFurnace[(Subscript[F, 0])=1] = whiteFurnace[ρ]:

Manipulate[
Show[
Plot[fWhiteFurnace[[order-1]][F0],{F0,0.04,1},PlotRange->{0,1*fWhiteFurnace[[order-1]][1]},AxesLabel->{"(Subscript[F, 0])","(Subscript[E, avg])"}]
,Plot[a2*((((Sqrt[τ*F0])))^(order^2)),{F0,0,1},PlotStyle->Red,PlotRange->Full,AxesLabel->{"(Subscript[F, 0])","(Subscript[E, avg])"}]
],{order,2,countOrders-1,1}]

Show[
Table[Plot[fWhiteFurnace[[order-1]][F0],{F0,0.04,1},PlotRange->{0,fWhiteFurnace[[1]][1]},AxesLabel->{"(Subscript[F, 0])","(Subscript[E, avg])"}],{order,2,countOrders-1}]
,Table[Plot[a2*((((Sqrt[τ*F0])))^(order^2)),{F0,0,1},PlotStyle->Red,PlotRange->Full],{order,2,countOrders-1}]
]

Show[
Table[Plot[Log[fWhiteFurnace[[order-1]][F0]],{F0,0.04,1},PlotRange->{0,-20},AxesLabel->{"(Subscript[F, 0])","(Subscript[E, avg])"}],{order,2,countOrders-1}]
,Table[Plot[Log[a2*((((Sqrt[τ*F0])))^(order^2))],{F0,0,1},PlotStyle->Red,PlotRange->Full],{order,2,countOrders-1}]
]

Now we check the sum of all orders:

Show[
Plot[Sum[fWhiteFurnace[[order-1]][ρ],{order,2,countOrders-1}],{ρ,0.04,1},PlotRange->{0,1},AxesLabel->{"ρ","(Subscript[E, avg])"}]
,Plot[factorMS[ρ],{ρ,0,1},PlotStyle->Red]
]

Here we simply check ((Subscript[F, 0])^((n^2)/2)) dependance:

Manipulate[
Show[
Plot3D[fTotalReflectance[[order]][x,α,F0],{x,0,(π/2)},{α,0,1},PlotRange->{0,fTotalReflectance[[order]][0.2,1,1]}]
,Plot3D[((((Sqrt[F0])))^(order^2))*fTotalReflectance[[order]][x,α,1],{x,0,(π/2)},{α,0,1},PlotStyle->Red]
],
{{F0,1},0.04,1},{order,2,countOrders-1,1}]

Simulated Multiple-Scattering Diffuse Reflectance Factor

SetDirectory["D:\\Workspaces\\GodComplex\\Tests\\TestMSBRDF\\Tables"]; 
countTheta = 16;
countRoughness= 16;
countAlbedo=20;
countOrders=7;
taRawTotalDiffuseReflectances=Table[
Import["TotalReflectance_Diffuse_Cosine_"  <>ToString[countTheta] <>"x"  <>ToString[countRoughness] <>"x"  <>ToString[countAlbedo] <>"_order"  <>ToString[order] <>".bin", "Real64"],{order,1,countOrders}];

Generate list interpolations

fTotalReflectance=Table[
ListInterpolation[Table[
taRawTotalDiffuseReflectances[[order]][[countTheta*((roughness-1)+countRoughness*(albedo-1))+theta]],{theta,1,countTheta},{roughness,1,countRoughness},{albedo,1,countAlbedo}],{{0,(π/2)},{1,0},{1,0}}],
{order,1,countOrders}];

Manipulate[
Plot3D[fTotalReflectance[[order]][x,α,ρ],{x,0,(π/2)},{α,0,1},PlotRange->{0,(3.5^(-order+1))}],{{ρ,1},0,1},{order,1,countOrders,1}]

Show Total Reflectance Depending on Albedo

Manipulate[Plot3D[Sum[fTotalReflectance[[order]][x,α,ρ],{order,1,countOrders}],{x,0,(π/2)},{α,0,1},PlotRange->{0,1},AxesLabel->{"θ","α"}],{{ρ,1},0,1}]

Manipulate[Plot3D[Sum[fTotalReflectance[[order]][x,α,ρ],{order,2,countOrders}],{x,0,(π/2)},{α,0,1},PlotRange->{0,0.3},AxesLabel->{"θ","α"}],{{ρ,1},0,1}]

Visualizing each individual order

Manipulate[
Show[
Table[
Plot3D[fTotalReflectance[[order]][x,α,ρ],{x,0,(π/2)},{α,0,1},
PlotRange->{0,0.15},PlotStyle->Hue[((4*(order-2))/countOrders)]],{order,2,countOrders}]
],
{{ρ,1},0,1}]

White Furnace Comparison

We computed the total energy bouncing off the surface at various scattering orders for different values of roughness, incidence angle and albedos.
So basically we performed:

	E((Subscript[ω, i]),α,ρ)=(∫_(SuperPlus[Ω])) (Subscript[L, i])((Subscript[ω, i])) (Subscript[f, r])( (Subscript[ω, o]),(Subscript[ω, i]),α,ρ)((Subscript[ω, o])·n)(Subscript[dω, o])
	
With (Subscript[L, i])((Subscript[ω, i]))=1
To compute the equivalent white furnace test we thus need to compute:

	(Subscript[E, total])(α,ρ)=(∫_(SuperPlus[Ω]))E((Subscript[ω, i]),α,ρ)((Subscript[ω, i])·n)(Subscript[dω, i])

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]

(*whiteFurnaceMS=Table[{α,fWhiteFurnaceMS[Max[0.01,α]]},{α,0,1,0.05}]*)

whiteFurnaceMS=Table[{α,π-fOrenEavg[α]},{α,0,1,0.05}]

Here we show that we mustn’t forget that the MS white furnace amounts to the complement of the total energy lost by SS BRDF...

Show[ListPlot[whiteFurnaceMS,PlotRange->Full, AxesLabel->{"α","π-(Subscript[E, avg])"}],
Plot[π-fOrenEavg[x],{x,0,1}]
]

Comparing with Simulation

We computed the total energy bouncing off the surface at various scattering orders for different values of roughness, incidence angle and albedos.
So basically we performed:

	E((Subscript[ω, i]),α,ρ)=(∫_(SuperPlus[Ω])) (Subscript[L, i])((Subscript[ω, i])) (Subscript[f, r])( (Subscript[ω, o]),(Subscript[ω, i]),α,ρ)((Subscript[ω, o])·n)(Subscript[dω, o])
	
With (Subscript[L, i])((Subscript[ω, i]))=1
To compute the equivalent white furnace test we thus need to compute:

	(Subscript[E, total])(α,ρ)=(∫_(SuperPlus[Ω]))E((Subscript[ω, i]),α,ρ)((Subscript[ω, i])·n)(Subscript[dω, i])

taRawSumTotalReflectances=Sum[taRawTotalDiffuseReflectances[[order]],{order,2,countOrders}]; (* Sum all orders > 2 *)
MSTotalReflectances=ListInterpolation[Table[taRawSumTotalReflectances[[countTheta*((roughness-1)+countRoughness*(ρ-1))+theta]],{theta,1,countTheta},{roughness,1,countRoughness},{ρ,1,countAlbedo}],{{0,(π/2)},{1,0},{1,0}}]
fWhiteFurnaceSimul[α_,ρ_,order_]:=2*π*NIntegrate[fTotalReflectance[[order]][θi,α,ρ]*Cos[θi]*Sin[θi],{θi,0,(π/2)},WorkingPrecision->MachinePrecision];
fSumWhiteFurnaceSimul[α_,ρ_]:=2*π*NIntegrate[MSTotalReflectances[θi,α,ρ]*Cos[θi]*Sin[θi],{θi,0,(π/2)},WorkingPrecision->MachinePrecision];
fWhiteFurnaceSimul[1,1,2]
fSumWhiteFurnaceSimul[1,1]

θ=0;
ρ=1;
Plot[fSumWhiteFurnaceSimul[α,ρ],{α,0,1}]

factor=((1.0940380117056643`)/(fSumWhiteFurnaceSimul[1,1]))

Show[
ListLinePlot[whiteFurnaceMS,AxesLabel->{"α","(Subscript[E, avg])"}]
,Plot[factor*fSumWhiteFurnaceSimul[α,ρ],{α,0,1},PlotStyle->Red]
]

We see here the red curve is not exactly identical to the white furnace result given by the theoretical multiple-scattering expression for the Oren-Nayar BRDF. Maybe it’s due to the many simplifications of the Oren-Nayar model? Maybe there’s a bias in my simulation?
Overall, we won’t linger too much on these differences and assume the simulation represents ground-truth data.

Finding an Analytical Model

We saw that integrating our order (SuperPlus[2]) irradiance values (the green curve in the above graph) yields π-(Subscript[E, total]) where (Subscript[E, total])(α,ρ) is the total irradiance reflected for the single-scattering Oren-Nayar BRDF.

We are interested in how each order contributes to that complement.
We assume the roughness-dependent shape is not important and concentrate on the maximum values of the white furnace functions for each order, when α=1. But we are interested in the evolution of that maximum value when varying the albedo.

So let’s create a table of samples of white furnace values for varying albedos, for each order (except last order, which is a compound of all remaining orders and doesn’t interest us):

α=1;
taWhiteFurnaceDiffuse=
Table[fWhiteFurnaceSimul[α,ρ,order],{order,2,countOrders},{ρ,0,1,0.02}];

fWhiteFurnaceDiffuse=Table[ListInterpolation[taWhiteFurnaceDiffuse[[order-1]],{0,1}],{order,2,countOrders}];
Show[Table[Plot[fWhiteFurnaceDiffuse[[order-1]][ρ],{ρ,0,1},PlotRange->Full,PlotStyle->Hue[((order-2)/countOrders)]],{order,2,countOrders}]]

Show[Table[Plot[Log[fWhiteFurnaceDiffuse[[order-1]][ρ]],{ρ,0,1},PlotRange->Full,PlotStyle->Hue[((order-2)/countOrders)]],{order,2,countOrders-1}]]

Now let’s try and fit each curve with a polynomial of the form f(ρ)=(Subscript[a, n])*(ρ^n) where n is the order:

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

fitter=Function[{order},
model=a*(x^order);
data=Table[{ρ,fWhiteFurnaceDiffuse[[order-1]][ρ]},{ρ,0,1,0.01}];
fit=FindFit[data,model,{a,b,c,d,e},x,Method->fittingMethod];
Function[x,Evaluate[model/.fit]];
a/.fit
];

fitter[3];

values=Table[fitter[order],{order,2,countOrders}]

maxValues=Table[fWhiteFurnaceDiffuse[[order-1]][1],{order,2,countOrders}]

Manipulate[
Show[
Plot[fWhiteFurnaceDiffuse[[order-1]][ρ],{ρ,0,1},PlotRange->{0,fWhiteFurnaceDiffuse[[order-1]][1]},AxesLabel->{"ρ","(Subscript[E, avg])"}]
,Plot[values[[order-1]]*(ρ^order),{ρ,0,1},PlotStyle->Red]
,Plot[maxValues[[order-1]]*(ρ^order),{ρ,0,1},PlotStyle->Green]
],{order,2,countOrders-1,1}]

Since our goal is to make the sum of these values be 1 (since they will guide the factor for the MSBRDF), we compute a normalization factor

limitValue=Sum[whiteFurnaceSimul[1,1,order],{order,2,countOrders}] (* Accurate but won't help *)
limitValue=Sum[maxValues[[i-1]],{i,2,countOrders}] (* We need to use our maximum values instead *)

newValues=(values/limitValue)
Show[
ListPlot[values,PlotRange->{0,1.1*limitValue}]
,ListLinePlot[Table[Sum[maxValues[[i]],{i,1,j}],{j,1,countOrders-1}],PlotRange->{0,1.1}](* Show that the sum converges to 1 *)
]
Sum[maxValues[[i-1]],{i,2,countOrders}]

ListPlot[Log[maxValues],PlotRange->Full]

Now we can see the graph of amplitudes for each order is clearly a geometric progression, as expected.
We’re looking for a function of the type:

	(Subscript[f, n])(ρ)=(Subscript[a, 1]) (((τρ))^n)

We already isolated the (ρ^n) part, now we’re looking for the terms (Subscript[a, 1]) and τ so that (Subscript[a, 1])·(τ^n) follows our table of values computed for ρ=1;

model=a*(b^(x+1));(* We need a geometric progression, x+1 because x=1 is actually order 2! *)
data=maxValues;
fit=FindFit[data,model,{a,b,c,d,e},x,Method->fittingMethod]
ffit=Function[x,Evaluate[model/.fit]];
Show[
ListPlot[maxValues,PlotRange->Full],
Plot[ffit[x],{x,0.1,countOrders},PlotStyle->Red]
]

We have our τ factor:

	(Subscript[f, n])(ρ)= (Subscript[a, 1]) (((τρ))^n)
	τ=0.28430405702379613


Since ρ<=1 then τρ<1 and we have a well-behaved geometric series that can be simplified into the final expression:

	(∑_(n=2)^∞) (((τρ))^n)=((((τρ))^2)/(1 - τρ))
	
The total amplitude for multiply-scattered energy should be:

	(∑_(n=2)^∞)(Subscript[f, n])(1) =1
	
And so:
	(∑_(n=2)^∞) (Subscript[a, 1])(τ^n) = (Subscript[a, 1]) (∑_(n=2)^∞) (τ^n)= (Subscript[a, 1]) ((τ^2)/(1 - τ))=1

And simply:

	(Subscript[a, 1]) =((1 - τ)/(τ^2))

Clear[factorMS]
a1=a/.fit
τ=b/.fit
a2=((1-τ)/(τ^2))
factorMS[ρ_]:=a2*((((τ*ρ))^2)/(1-τ*ρ))
factorMS[1]
Plot[factorMS[ρ],{ρ,0,1},PlotRange->All,AxesLabel->{"ρ","(Subscript[F, ms])(ρ)"}]

Verification

Now we must verify that our expression can be substituted to the actual measurements!
We first verify that white furnace results match for each individual orders (Subscript[f, n])(ρ) * whiteFurnace[ρ=1] = whiteFurnace[ρ]:

Manipulate[
Show[
Plot[fWhiteFurnaceDiffuse[[order-1]][ρ],{ρ,0,1},PlotRange->{0,fWhiteFurnaceDiffuse[[order-1]][1]},AxesLabel->{"ρ","(Subscript[E, avg])"}]
,Plot[a1*(((τ*ρ))^order),{ρ,0,1},PlotStyle->Red]
],{order,2,countOrders-1,1}]

Show[
Table[Plot[fWhiteFurnaceDiffuse[[order-1]][ρ],{ρ,0,1},PlotRange->{0,fWhiteFurnaceDiffuse[[1]][1]},AxesLabel->{"ρ","(Subscript[E, avg])"}],{order,2,countOrders-1}]
,Table[Plot[a1*(((τ*ρ))^order),{ρ,0,1},PlotStyle->Red,PlotRange->Full],{order,2,countOrders-1}]
]

Now we check the sum of all orders:

Show[
Plot[Sum[fWhiteFurnaceDiffuse[[order-1]][ρ],{order,2,countOrders-1}],{ρ,0,1},PlotRange->{0,1},AxesLabel->{"ρ","(Subscript[E, avg])"}]
,Plot[factorMS[ρ],{ρ,0,1},PlotStyle->Red]
]

Here we simply check (ρ^n) dependence:

Manipulate[
Show[
Plot3D[fTotalReflectance[[order]][x,α,albedo],{x,0,(π/2)},{α,0,1},PlotRange->{0,fTotalReflectance[[order]][0.2,1,1]}]
,Plot3D[(albedo^order)*fTotalReflectance[[order]][x,α,1],{x,0,(π/2)},{α,0,1},PlotStyle->Red]
],
{{albedo,1},0,1},{order,2,countOrders-1,1}]