Clear[ndf,pdf,cdf]
ndf[θ_,α_]:=((α^2)/(π*(((((Cos[θ])^2)*((α^2)-1)+1))^2)))
(*cdf=Function[{θ,α},Evaluate[Integrate[2*π*ndf[θ,α]*Cos[θ]*Sin[θ],θ]]]*)
cdf[θ_,α_]:=((1-((Cos[θ])^2))/(1+((Cos[θ])^2)*((α^2)-1)))
pdf[θ_,α_]:=ndf[θ,α]*Cos[θ]
pdfCos[cosθ_,α_]:=((α^2)/(π*(((((cosθ)^2)*((α^2)-1)+1))^2)))*cosθ

Compute the difference of the PDFs of 2 neighbor pixels

(*Clear[α]
Integrate[2*π*(ndf[θ,α] - ndf[θ-β,α])*Cos[θ]*Sin[θ],θ]*)

ndf2[μ_,α_]:=((α^2)/(π*((((μ^2)*((α^2)-1)+1))^2)))
test[μ0_,α_]:=cdf[μ0,α]+ NIntegrate[2*π*(ndf2[μ,α] - ndf2[μ-μ0,α])*μ*(Sqrt[1-(μ^2)]),{μ,μ0,0}]

Manipulate[Plot[{test[μ0,α],cdf[μ0,α]},{α,0,1},PlotRange->{0,1}],{μ0,1,0}]

test[θ0_,α_]:=cdf[Cos[θ0],α]+ NIntegrate[2*π*(ndf[θ,α] - ndf[θ-θ0,α])*Cos[θ]*Sin[θ],{θ,θ0,(π/2)}]
Manipulate[Plot[{test[θ,α],cdf[Cos[θ],α]},{α,0,1},PlotRange->{-1,1}],{θ,0.01,(π/2)}]

Compute the correlation between 2 opposed NDFs separated by θ

Here we have 2 lobes aligned along 2 directions separated by an angle θ and we need to find the correlation least correlation between the 2 although there must be a little correlation anyway otherwise we could always choose roughness 0 for maximization.

Test with subtraction of the 2 lobes:

Manipulate[Show[Plot[((pdf[θ,α]*Sin[θ])/(cdf[((θ0)/2),α]))-((pdf[θ0-θ,α]*Sin[θ0-θ])/(cdf[((θ0)/2),α])),{θ,0,θ0},PlotRange->Full],ParametricPlot[{((θ0)/2),t},{t,-100,100},PlotStyle->Red]],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]

Test with multiplication of the 2 lobes (true correlation):

Manipulate[Show[Plot[(((pdf[θ,α]*Sin[θ])*(pdf[θ0-θ,α]*Sin[θ0-θ]))/((cdf[θ0,α])^2)),{θ,0,θ0},PlotRange->Full],ParametricPlot[{((θ0)/2),t},{t,-100,100},PlotStyle->Red]],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]

Showing angular correlation visually:

ColorData[97,"ColorList"]
ColorData[97][1]

defaultBlue=ColorData[97][1];
Manipulate[ParametricPlot[ { 4*θ*{1,Tan[(π/4)-((θ0)/2)]},4*θ*{1,Tan[(π/4)+((θ0)/2)]}, 4*θ*{1,1},
((pdf[Abs[θ-((π/4)-((θ0)/2))],α])/(cdf[((θ0)/2),α]))*{Sin[θ],Cos[θ]},
((pdf[Abs[((π/4)+((θ0)/2))-θ],α])/(cdf[((θ0)/2),α]))*{Sin[θ],Cos[θ]},
(((pdf[Abs[θ-((π/4)-((θ0)/2))],α]*pdf[Abs[((π/4)+((θ0)/2))-θ],α]))/((cdf[((θ0)/2),α])^2))*{Sin[θ],Cos[θ]}
},{θ,0,(π/2)},PlotRange->{{0,1},{0,1}},PlotStyle->{LightGray,LightGray,Orange,defaultBlue,defaultBlue,Red}],{{α,0.1},0,1},{{θ0,(π/2)},0,(π/2)}]

(*Manipulate[Show[Plot[ndf[θ,α]*Cos[θ],{θ,0,(π/2)},PlotRange->Full],ParametricPlot[{((θ0)/2),t},{t,-100,100},PlotStyle->Red]],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]*)

(*Clear[α]
Integrate[(((pdf[θ,α]*Sin[θ])*(pdf[θ0-θ,α]*Sin[θ0-θ]))/((cdf[θ0,α])^2)),θ]*)

(*Manipulate[NIntegrate[2*π*(ndf[θ,α]*Cos[θ]-ndf[θ-θ0,α]*Cos[θ-θ0])*Sin[θ],{θ,0,((θ0)/2)}],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]*)
Manipulate[NIntegrate[(((pdf[θ,α]*Sin[θ])*(pdf[θ0-θ,α]*Sin[θ0-θ]))/((cdf[θ0,α])^2)),{θ,0,((θ0)/2)}],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]

(*Manipulate[Plot[NIntegrate[2*π*(ndf[θ,α]*Cos[θ]-ndf[θ-θ0,α]*Cos[θ-θ0])*Sin[θ],{θ,0,((θ0)/2)}],{α,0,1},PlotRange->{0,1}],{{θ0,(π/10)},0,(π/2)}]*)
Manipulate[Plot[NIntegrate[(((pdf[θ,α]*Sin[θ])*(pdf[θ0-θ,α]*Sin[θ0-θ]))/((cdf[θ0,α])^2)),{θ,0,((θ0)/2)}],{α,0,1},PlotRange->{0,1}],{{θ0,(π/10)},0,(π/2)}]

Show the correlation

(*Manipulate[Show[Plot[ndf[θ,α]*Cos[θ]*ndf[θ-θ0,α]*Cos[θ-θ0],{θ,0,θ0},PlotRange->{0,Full}],ParametricPlot[{((θ0)/2),t},{t,-100,100},PlotStyle->Red],ParametricPlot[{t,0.05},{t,0,θ0},PlotStyle->Red]],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]*)
(*Manipulate[Show[Plot[(((pdf[θ,α]*Sin[θ])*(pdf[θ0-θ,α]*Sin[θ0-θ]))/((cdf[θ0,α])^2)),{θ,0,θ0},PlotRange->{0,Full}],ParametricPlot[{((θ0)/2),t},{t,-100,100},PlotStyle->Red],ParametricPlot[{t,((((1/π)))^2)},{t,0,θ0},PlotStyle->Red]],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]*)
Manipulate[Show[Plot[((pdf[θ,α]*pdf[θ0-θ,α])/((cdf[((θ0)/2),α])^2)),{θ,0,θ0},PlotRange->{0,Full}],ParametricPlot[{((θ0)/2),t},{t,-100,100},PlotStyle->Red],ParametricPlot[{t,(((((Sqrt[2])/π)))^2)},{t,0,θ0},PlotStyle->Red]],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]

So okay, we can accept a certain amount of correlation that is given by the value at ((Subscript[θ, 0])/2)when α=1 and (Subscript[θ, 0])=(π/2)
Now we can find the roughness for a given (Subscript[θ, 0]) :

	c=(((((ndf(((Subscript[θ, 0])/2),α) cos(((Subscript[θ, 0])/2))sin(((Subscript[θ, 0])/2)))/(cdf(((Subscript[θ, 0])/2),α)))))^2)

We can totally get rid of the squaring and rewrite the same expression for a new c.

By fixing c = ((ndf((π/4),1) cos((π/4))sin((π/4)))/(cdf((π/2),1))) corresponding to the correlation of the worst mip that sees a (π/2) angle (a full face of the cube) and must correspond to the roughness value α=1.

We find:	c=(1/π)

((ndf[(π/4),1]*Cos[(π/4)]*Sin[(π/4)])/(cdf[(π/4),1]))
((pdf[(π/4),1])/(cdf[(π/4),1]))

Solve[((pdf[θ,α]*Sin[θ])/(cdf[θ,α]))==c,α]

Clear[c]
Solve[((pdf[θ,α])/(cdf[θ,α]))==c,α]

(*cmax = (1/π);
c=(cmax/(8π));
sqAlpha[θ_]:=((c*((Cos[θ])^2))/(Sin[θ]-c*((Sin[θ])^2)))
sqAlpha[θ_]:=((-π*cmax* ((Sin[θ])^4))/(π*cmax* ((Cos[θ])^2)-π*cmax* ((Cos[θ])^4)-Cos[θ]*Sin[θ]))*)

cmax=((Sqrt[2])/π);
c=cmax * π;
(*sqAlpha[θ_]:=((-c*((Sin[θ])^4))/(c* ((Cos[θ])^2)-c* ((Cos[θ])^4)-Cos[θ]*Sin[θ])) with sin(θ) *)
sqAlpha[θ_]:=((-c*((Sin[θ])^4))/(c* ((Cos[θ])^2)-c* ((Cos[θ])^4)-Cos[θ]))

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

mipMax=8;
(*mu[mip_]:=1-(1/3)(2^(2(mip-mipMax)))   Solid angles yields θ > (π/2)*)
mu[mip_]:=(1/(Sqrt[1+(2^(2(mip-mipMax)))]))
sqRoughnessFromMip[mip_]:=sqAlpha[ ArcCos[mu[mip]] ]

mu[8]
2*ArcCos[mu[8]]
Table[N[sqRoughnessFromMip[mip]],{mip,0,mipMax}]

Plot[2*ArcCos[mu[mip]],{mip,0,mipMax},PlotRange->{0,(π/2)}]

Plot[{sqRoughnessFromMip[mip],(Sqrt[sqRoughnessFromMip[mip]])},{mip,0,mipMax},PlotRange->{0,1}]

Compute the correlation between 2 opposed NDFs separated by θ

Here we have 2 lobes aligned along 2 directions separated by an angle θ and we need to find the correlation least correlation between the 2 although there must be a little correlation anyway otherwise we could always choose roughness 0 for maximization.

It would look like we’re doing the same thing all over again but the normalization by pdf(0,α) makes much more sense than dividing by the cdf!

Showing angular correlation visually:

defaultBlue=ColorData[97][1];
(*correl[θ_,α_,θ0_]:=((pdf[Abs[θ-((π/4)-((θ0)/2))],α]*pdf[Abs[((π/4)+((θ0)/2))-θ],α])/((pdf[0,α])^2))
correl[θ_,α_,θ0_]:=((Max[0,Min[pdf[Abs[θ-((π/4)-((θ0)/2))],α],pdf[Abs[((π/4)+((θ0)/2))-θ],α]]])/(pdf[0,α]))
Manipulate[ParametricPlot[ { 4*θ*{1,Tan[(π/4)-((θ0)/2)]},4*θ*{1,Tan[(π/4)+((θ0)/2)]},{Sin[θ],Cos[θ]},4*θ*{1,1},
((pdf[Abs[θ-((π/4)-((θ0)/2))],α])/(pdf[0,α]))*{Sin[θ],Cos[θ]},
((pdf[Abs[((π/4)+((θ0)/2))-θ],α])/(pdf[0,α]))*{Sin[θ],Cos[θ]},
correl[θ,α,θ0]*{Sin[θ],Cos[θ]},
((((1/(Sqrt[2]))))^2)*{Sin[θ],Cos[θ]}
},{θ,0,(π/2)},PlotRange->{{0,1},{0,1}},PlotStyle->{LightGray,LightGray,LightGray,Orange,defaultBlue,defaultBlue,Red,Green}],{{α,0.1},0,1},{{θ0,(π/2)},0,(π/2)}]*)

correl[θ_,α_,θ0_]:=((pdf[θ,α]*pdf[θ0-θ,α])/((pdf[0,α])^2))
correl[θ_,α_,θ0_]:=((Max[0,Min[pdf[θ,α],pdf[θ0-θ,α]]])/(pdf[0,α]))
Manipulate[ParametricPlot[ { 4*θ*{1,Tan[(π/2)-θ0]},{Sin[θ],Cos[θ]},4*θ*{Tan[((θ0)/2)],1},
((pdf[θ,α])/(pdf[0,α]))*{Sin[θ],Cos[θ]},
((pdf[θ0-θ,α])/(pdf[0,α]))*{Sin[θ],Cos[θ]},
correl[θ,α,θ0]*{Sin[θ],Cos[θ]},
((((1/(Sqrt[2]))))^2)*{Sin[θ],Cos[θ]}
},{θ,-(π/2),(π/2)},PlotRange->{{-1,1},{0,1}},PlotStyle->{LightGray,LightGray,Orange,defaultBlue,defaultBlue,Red,Green}],{{α,0.1},0,1},{{θ0,(π/2)},0,(π/2)}]

Plotting angular correlation as a function of angle:

(*Manipulate[Show[Plot[(((pdf[θ,α]*Sin[θ])*(pdf[θ0-θ,α]*Sin[θ0-θ]))/((pdf[0,α])^2)),{θ,0,θ0},PlotRange->Full],ParametricPlot[{((θ0)/2),t},{t,-100,100},PlotStyle->Red]],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]*)
Manipulate[Show[Plot[((pdf[θ,α]*pdf[θ0-θ,α])/((pdf[0,α])^2)),{θ,0,θ0},PlotRange->{0,Full}],ParametricPlot[{((θ0)/2),t},{t,-100,100},PlotStyle->Red],ParametricPlot[{t,((((1/(Sqrt[2]))))^2)},{t,0,θ0},PlotStyle->Red]],{{θ0,(π/10)},0,(π/2)},{{α,0.1},0,1}]

Computing Surface Correlation

Let’s try and integrate the volume coverage of the interpenetration of the 2 lobes, instead of looking at a single value at the half-angle...

(*Manipulate[2π*NIntegrate[pdf[θ,α]*pdf[θ0-θ,α]*Sin[θ],{θ,0,(π/2)}],{{θ0,(π/2)},0,(π/2)},{{α,1},0.0001,1}]
Manipulate[((2π)/(0+1*((pdf[0,α])^2)))*NIntegrate[pdf[θ,α]*pdf[θ0-θ,α]*Sin[θ],{θ,0,(π/2)}],{{θ0,(π/2)},0,(π/2)},{{α,1},0.0001,1}]
Manipulate[((2π)/(1+0*pdf[0,α]))*NIntegrate[Min[pdf[θ,α],pdf[θ0-θ,α]]*Sin[θ ],{θ,0,(π/2)}],{{θ0,(π/2)},0,(π/2)},{{α,1},0,1}]*)

Manipulate[(1/2)*NIntegrate[((correl[θ,α,θ0])^2),{θ,-(π/2),(π/2)}],{{θ0,(π/2)},0,(π/2)},{{α,1},0,1}]
Manipulate[(1/(2*((pdf[0,α])^2)))*NIntegrate[(((pdf[θ,α]))^2),{θ,-(π/2),(π/2)}],{{α,1},0,1}]
(1/2)NIntegrate[((Cos[θ])^2),{θ,-(π/2),(π/2)}]
pdf[0,1]

Finding a Mapping

So okay, we can accept a certain amount of correlation that is given by the value at ((Subscript[θ, 0])/2)when α=1 and (Subscript[θ, 0])=(π/2)
Now we can find the roughness for a given (Subscript[θ, 0]) :

	c=(((((pdf(((Subscript[θ, 0])/2),α) )/(pdf(0,α)))))^2)

We can totally get rid of the squaring and rewrite the same expression for a new c.

By fixing c = ((pdf((π/4),1) )/(pdf(0,1))) corresponding to the correlation of the worst mip that sees a (π/2) angle (a full face of the cube) and must correspond to the roughness value α=1.

We find:	c=(1/(Sqrt[2]))

((ndf[(π/4),1]*Cos[(π/4)])/(pdf[0,1]))
((pdf[(π/4),1])/(pdf[0,1]))

Clear[c,θ,α]
(*Solve[((pdf[θ,α])/(pdf[0,α]))==c,α]*)
Solve[(((α^4)*b)/((((b^2)((α^2)-1)+1))^2))==c,α]

Résultat impossible à inverser

(*Solve[(((α^4)*b)/((((b^2)((α^2)-1)+1))^2))==c,b]*)

c=(1/(Sqrt[2]));
sqAlpha:=Function[θ,
μ=Cos[θ];
(((μ^4)*c-(μ^2)*c-(Sqrt[μ*(((1-(μ^2)))^2)*c]))/((μ^4)*c-μ))
]
sqAlphaMu[μ_]:=(((μ^4)*c-(μ^2)*c-(Sqrt[μ*(((1-(μ^2)))^2)*c]))/((μ^4)*c-μ))

Plot[{(Sqrt[sqAlpha[(θ/2)]]),sqAlpha[(θ/2)]},{θ,0,(π/2)},PlotRange->{0,1},AxesLabel->{"θ","α"}]
Plot[{(Sqrt[sqAlphaMu[μ]]),sqAlphaMu[μ]},{μ,(1/(Sqrt[2])),1},PlotRange->{{0.70,1},{0,1}},AxesLabel->{"cos(θ)","α"}]

mipMax=8;
mu[mip_]:=1-(1/3)(2^(2(mip-mipMax)))   (*Solid angles yields θ > (π/2)*)
mu2[mip_]:=(1/(Sqrt[1+(2^(2(mip-mipMax)))]))
sqRoughnessFromMip[mip_]:=sqAlpha[ArcCos[mu[mip]] ]
sqRoughnessFromMip2[mip_]:=sqAlpha[ArcCos[mu2[mip]] ]
mu[8]
N[2*ArcCos[mu[8]]]
Table[N[(Sqrt[sqRoughnessFromMip[mip]])],{mip,0,mipMax}]
Plot[{sqAlpha[(θ/2)],(Sqrt[sqAlpha[(θ/2)]])},{θ,0,(π/2)},PlotRange->{0,1}]
Plot[{2*ArcCos[mu[mip]],2*ArcCos[mu2[mip]]},{mip,0,mipMax},PlotRange->{0,1.1*(π/2)}]
Plot[{(Sqrt[sqRoughnessFromMip[mip]]),(Sqrt[sqRoughnessFromMip2[mip]])},{mip,0,mipMax},PlotRange->{0,1}]

Manipulate[ParametricPlot[ { 4*θ*{1,Tan[(π/4)-((θ0)/2)]},4*θ*{1,Tan[(π/4)+((θ0)/2)]},{Sin[θ],Cos[θ]},4*θ*{1,1},
((pdf[Abs[θ-((π/4)-((θ0)/2))],(Sqrt[sqAlpha[((θ0)/2)]])])/(pdf[0,(Sqrt[sqAlpha[((θ0)/2)]])]))*{Sin[θ],Cos[θ]},
((pdf[Abs[((π/4)+((θ0)/2))-θ],(Sqrt[sqAlpha[((θ0)/2)]])])/(pdf[0,(Sqrt[sqAlpha[((θ0)/2)]])]))*{Sin[θ],Cos[θ]},
(((pdf[Abs[θ-((π/4)-((θ0)/2))],(Sqrt[sqAlpha[((θ0)/2)]])]*pdf[Abs[((π/4)+((θ0)/2))-θ],(Sqrt[sqAlpha[((θ0)/2)]])]))/((pdf[0,(Sqrt[sqAlpha[((θ0)/2)]])])^2))*{Sin[θ],Cos[θ]},
((((1/(Sqrt[2]))))^2)*{Sin[θ],Cos[θ]}
},{θ,0,(π/2)},PlotRange->{{0,1},{0,1}},PlotStyle->{LightGray,LightGray,LightGray,Orange,defaultBlue,defaultBlue,Red,Green}],{{θ0,(π/2)},0.001,(π/2)}]

Fitting

The expression to get μ from roughness is a nightmare to invert so let’s do a curve fitting instead...

data=Table[{ N[sqAlphaMu[μ]],N[μ]},{μ,(1/(Sqrt[2])),1,((1-(1/(Sqrt[2])))/100)}]
ListLinePlot[data,PlotRange->Full]

Clear[a,b,c,d,e,f]
model=((a+b*x+c*(x^2))/(d+e*x+f*(x^2)))
model=((a+b*x+c*(x^2))/(a+e*x+f*(x^2)))

f0=Function[x,Evaluate[model/.FindFit[data,model,{a,b,c,d,e,f},x]]]
Show[ListLinePlot[data,PlotRange->{(1/(Sqrt[2])),1}],Plot[f0[x],{x,0,1},PlotStyle->Red]]
(Sqrt[2])f0[1]

3D Lobes Inter-Penetration Modeling

Here we’re attempting to port the same reflection to the 3D realm. The idea is measuring how much 2 adjacent lobes interpenetrate each other depending on the separation angle (i.e. the lobe’s aperture angle) and the roughness of the lobe.
Once again we “normalize a lobe” by dividing by its largest amplitude at θ=0 and we’re looking to compute the integral:

	I=(∫_(-(Subscript[θ, 0]))^(Subscript[θ, 0]))(∫_0^(π/2))(((((D(θ,α) D( θ', α ) cos(θ)  cos(θ'))/((D( 0, α ))^2))))^3) sin(θ)dθdϕ

Here, cos(θ')= sin( θ ) cos( ϕ ) sin( (Subscript[θ, 0])) + cos( θ ) cos( (Subscript[θ, 0]) ) and (Subscript[θ, 0]) is the cone’s aperture angle.

cosAngleChange[θ_,ϕ_,γ_]:=Sin[θ]*Cos[ϕ]*Sin[γ]+Cos[θ]*Cos[γ]
penelobe[θ_,ϕ_,γ_,α_]:=pdf[θ,α]*pdfCos[cosAngleChange[θ,ϕ,γ],α]

γ=0.1
α=0.05
NIntegrate[((penelobe[θ,ϕ,γ,α])^3)*Sin[θ],{θ,0,(π/2)-0.01},{ϕ,-(π/2),(π/2)}]

lobeSeparationAngle[θ0_]:=ArcCos[((Cos[θ0]-((Cos[θ0])^2))/(1-((Cos[θ0])^2)))]
lobesCount[θ0_]:=((2π)/(lobeSeparationAngle[Max[0.0001,θ0]]))
lobeSeparationAngle[0.0001]
Plot[lobesCount[θ0],{θ0,0,(π/2)},PlotRange->{0,6}]

Here it’s interesting to realize the limit converges to 6 possible neighbors, I would have expected an unlimited amount of neighbors as the lobes get thinner and thinner...
But after all, if you zoom in on things, even though they’re getting very small, maybe you always end up with 6 neighbor lobes no matter what.

The other extremity of 4 lobes when the angle is (π/2)is totally expected, on the other hand: each face of the cube has exactly 4 neighbor faces, thus 4 neighbor lobes.

Manipulate[(1+0*lobesCount[θ0])*NIntegrate[(((((penelobe[θ,ϕ,θ0,α])/((pdf[0,α])^2))))^3)*Sin[θ],{θ,0,(π/2)-0.01},{ϕ,-π,π}],{{θ0,(π/10)},0.001,(π/2)},{{α,0.9},0.001,1}]

Manipulate[(1+0*lobesCount[θ0])*NIntegrate[(((((penelobe[θ,ϕ,θ0,α])/((pdf[0,α])^2))))^3)*Sin[θ],{θ,0,(π/2)-0.01},{ϕ,-π,π}],{{θ0,(π/10)},0.001,(π/2)},{{α,0.9},0.001,1}]

Manipulate[NIntegrate[((penelobe[θ,ϕ,θ0,α])/((pdf[0,α])^2))*Sin[θ],{θ,0,(π/2)-0.01},{ϕ,-(π/2),(π/2)}],{{θ0,(π/10)},0.001,(π/2)},{{α,0.9},0.001,1}]

Manipulate[lobesCount[θ0]*NIntegrate[penelobe[θ,ϕ,θ0,α]*Sin[θ],{θ,0,(π/2)-0.01},{ϕ,-(π/2),(π/2)}],{{θ0,(π/10)},0.001,(π/2)},{{α,0.9},0.0001,1}]

Manipulate[NIntegrate[((pdf[θ,α])/(pdf[0,α]))*Sin[θ],{θ,0,(π/2)-0.01},{ϕ,-(π/2),(π/2)}],{{α,0.9},0.001,1}]

Manipulate[((lobesCount[θ0])/2)*NIntegrate[pdf[θ,α]*Sin[θ],{θ,0,(π/2)-0.01},{ϕ,-(π/2),(π/2)}],{{θ0,(π/10)},0.001,(π/2)},{{α,0.9},0.0001,1}]

Manipulate[ParametricPlot3D[{
((pdf[θ,α])/(pdf[0,α]))*{Sin[θ]*Cos[ϕ], Sin[θ]*Sin[ϕ], Cos[θ]},
((pdfCos[cosAngleChange[θ,ϕ,γ],α])/(pdf[0,α]))*{Sin[θ]*Cos[ϕ], Sin[θ]*Sin[ϕ], Cos[θ]},
((Min[pdf[θ,α],pdfCos[cosAngleChange[θ,ϕ,γ],α]])/(pdf[0,α]))*{Sin[θ]*Cos[ϕ], Sin[θ]*Sin[ϕ], Cos[θ]}
},{θ,0,(π/2)},{ϕ,-π,π},PlotRange->{{-1,1},{-1,1},{0,1}},AxesLabel->{"x","y","z"},PlotStyle->{Opacity[0.125],Opacity[0.125],Directive[Red,Opacity[1]]}],{{α,1},0.0001,1},{{γ,(π/2)},0.0001,(π/2)}]

penelobeMin[θ_,ϕ_,γ_,α_]:=Min[pdf[θ,α],pdfCos[cosAngleChange[θ,ϕ,γ],α]]
Manipulate[NIntegrate[penelobeMin[θ,ϕ,θ0,α]*Sin[θ],{θ,0,(π/2)},{ϕ,-(π/2),(π/2)}, PrecisionGoal->2,MaxRecursion->4,WorkingPrecision->$MachinePrecision],{{θ0,(π/2)},0.001,(π/2)},{{α,1},0.0001,1}]

SetDirectory["D:\\Workspaces\\Patapom.com\\Web\\blog\\patapom\\docs\\BRDF"]; 
countX=128;
countY=64;

(*ta=ParallelTable[{N[θ0],N[α],
NIntegrate[penelobeMin[θ,ϕ,θ0,α]*Sin[θ],{θ,0,(π/2)},{ϕ,-(π/2),(π/2)}, PrecisionGoal->2,MaxRecursion->4,WorkingPrecision->$MachinePrecision]
}
,{θ0,0,(π/2),((π/2)/countX)}
,{α,0.001,1,((1-0.001)/countY)}]
Export["VolumeInterPenetration_"<>ToString[countX]<>"x"<>ToString[countY]<>".csv",ta,"CSV"];*)

ta=Import["VolumeInterPenetration_"<>ToString[countX]<>"x"<>ToString[countY]<>".csv","CSV"];
ta=ToExpression[ta]

ListPlot3D[Flatten[ta,1],PlotRange->Full]

ta2=Table[ta[[x]][[y]] * {1,1,((lobesCount[ta[[x]][[y]][[1]]])/(lobesCount[(π/2)]))},{x,1,countX+1},{y,1,countY+1}];
ta3=Table[ta[[x]][[y]] [[3]],{x,1,countX+1},{y,1,countY+1}]; (* This table doesn't account for lobes increase with angle *)
ta3=Table[((lobesCount[ta[[x]][[y]][[1]]])/(lobesCount[(π/2)]))*ta[[x]][[y]] [[3]],{x,1,countX+1},{y,1,countY+1}];  (* This table will increase the amount of lobes from 4 at (π/2) to 6 at 0° *)
roughnessFunction=ListInterpolation[ta3,{{0,(π/2)},{0,1}}]
ta[[countX+1]][[countY+1]]
ref=ta2[[countX+1]][[countY+1]][[3]]

Show[(*ListPlot3D[Flatten[ta2,1],PlotRange->Full,AxesLabel->{"θ","α","Penetration"}],*)
Plot3D[roughnessFunction[θ,α],{θ,0,(π/2)},{α,0,1}],
Plot3D[ref,{θ,0,(π/2)},{α,0,1},PlotStyle->Darker[Green,0.1]]]

Clear[α]
NSolve[roughnessFunction[0.2,α]==ref,α]

roughnessMap=Table[{N[θ0],α/.NSolve[roughnessFunction[Max[0.0001,θ0],α]==ref,α][[1]]},{θ0,0,(π/2),((π/2)/200)}]

ListLinePlot[roughnessMap,PlotRange->{0,1}]

ta3[[129]]
ListInterpolation[ta3[[129]],{0,1}]

taRoughnessFunctions=Table[ListInterpolation[ta3[[θ]],{0,1}],{θ,1,countX+1}];

Plot[Table[taRoughnessFunctions[[x]][α],{x,1,countX+1}],{α,0,1},PlotRange->Full]

roughnessMap2=Table[{N[((θ0-1)/countX)*(π/4)],α/.NSolve[taRoughnessFunctions[[θ0]][α]==ref,α][[1]]},{θ0,2,countX+1}]
roughnessMapCos2=Table[{N[Cos[((θ0-1)/countX)*(π/4)]],α/.NSolve[taRoughnessFunctions[[θ0]][α]==ref,α][[1]]},{θ0,2,countX+1}]
ListLinePlot[roughnessMap2]

roughnessMap3={};
For[i=1,i<=countX,i++,If[roughnessMap2[[i]][[2]]≥0,AppendTo[roughnessMap3,roughnessMap2[[i]]]]]
ListLinePlot[roughnessMap3]

Let’s fit this...

Clear[a,b,c,d,e,f,x,model]
model=0+b*x+c*(x^2)+d*(x^3)+e*(x^4)+f*(x^5);
model=((0+b*x+c*(x^2))/(1+e*x+0*f*(x^2)));
FindFit[roughnessMap3,model,{a,b,c,d,e,f},x]
f2=Function[x,Evaluate[model/.FindFit[roughnessMap3,model,{a,b,c,d,e,f},x]]];
Show[ListLinePlot[roughnessMap3,PlotRange->{0,1}],Plot[f2[x],{x,0,(π/2)},PlotRange->{0,1},PlotStyle->Red]]

Sum[(((f2[roughnessMap3[[i]][[1]]] - roughnessMap3[[i]][[2]]))^2),{i,1,Length[roughnessMap3]}]

Let’s plot and compare with the 2D model

c=(1/(Sqrt[2]));
sqAlpha:=Function[θ,
μ=Cos[θ];
(((μ^4)*c-(μ^2)*c-(Sqrt[μ*(((1-(μ^2)))^2)*c]))/((μ^4)*c-μ))
]
sqAlphaMu[μ_]:=(((μ^4)*c-(μ^2)*c-(Sqrt[μ*(((1-(μ^2)))^2)*c]))/((μ^4)*c-μ))
sqRoughnessFromMip[mip_]:=sqAlphaMu[mu[mip]] 
sqRoughnessFromMip3D[mip_]:=((f2[mu[mip]])^2)
sqRoughnessFromMip3D[mip_]:=((f2[ArcCos[mu[mip]]])^2)
(*Plot[{sqAlpha[(θ/2)],(Sqrt[sqAlpha[(θ/2)]])},{θ,0,(π/2)},PlotRange->{0,1}]*)
Plot[{((sqAlpha[(θ/2)])^(1/2)),f2[(θ/2)]},{θ,0,(π/2)},PlotRange->{0,1},AxesLabel->{"θ","α"}]
Plot[{(Sqrt[sqRoughnessFromMip[mip]]),(Sqrt[sqRoughnessFromMip3D[mip]])},{mip,0,mipMax},PlotRange->{0,1},AxesLabel->{"mip","α"}]

Let’s map the inverse function now:

data=Table[{ N[f2[ArcCos[μ]]],N[μ]},{μ,(1/(Sqrt[2])),1,((1-(1/(Sqrt[2])))/100)}];
ListLinePlot[data,PlotRange->Full]

Clear[a,b,c,d,e,f]
model=((a+b*x+c*(x^2))/(d+e*x+f*(x^2)))
model=((a+b*x+c*(x^2))/(a+e*x+f*(x^2)))
f3=Function[x,Evaluate[model/.FindFit[data,model,{a,b,c,d,e,f},x]]]
Show[ListLinePlot[data,PlotRange->{(1/(Sqrt[2])),1}],Plot[f3[x],{x,0,1},PlotStyle->Red]]
(Sqrt[2])f3[1]

Let’s plot and compare with the 2D model:

Plot[{f0[μ],f3[μ]},{μ,0,1},AxesLabel->{"μ","α"}]

3D Lobes Inter-Penetration Modeling

We’re starting over again but this time measuring the ratio of volume that get interpenetrated depending on the angle between lobes.
We know a volume element is given by dV=(r^2)sin(θ) dr dθ dϕ and the volume integral is:

	I=(∫_(-π)^π)(∫_0^(π/2))(∫_0^(r(θ,ϕ,α)))(r^2) sin(θ)dθdϕ=(∫_(-π)^π)(∫_0^(π/2)) (1/3)((r(θ,ϕ,α))^3) sin(θ)dθdϕ

In one case, we’re looking for the entire volume so r(θ,ϕ,α)=r(θ )=((D(θ,α))/π) cos(θ) the BRDF’s normal distribution function times the cosine of the angle.
In the second case we’re looking for the intersection of 2 volumes and r(θ,ϕ,α)=r(θ )=min( D(θ,α) cos(θ), D(θ',α) cos(θ') ).

Here, cos(θ')= sin( θ ) cos( ϕ ) sin( (Subscript[θ, 0])) + cos( θ ) cos( (Subscript[θ, 0]) ) and (Subscript[θ, 0]) is the cone’s aperture angle.

ParametricPlot[Cos[θ]*{Sin[θ],Cos[θ]},{θ,-(π/2),(π/2)}]

cosAngleChange[θ_,ϕ_,γ_]:=Sin[θ]*Cos[ϕ]*Sin[γ]+Cos[θ]*Cos[γ]
penelobe[θ_,ϕ_,γ_,α_]:=π*Max[0,Min[pdf[θ,α],pdfCos[cosAngleChange[θ,ϕ,γ],α]]]
Manipulate[((2π)/3)*NIntegrate[((Cos[θ])^3)*Sin[θ],{θ,0,(π/2)}],{{α,1},0.001,1}]
Manipulate[((2π)/3)*NIntegrate[(((π*pdf[θ,α]))^3)*Sin[θ],{θ,0,(π/2)}],{{α,1},0.001,1}]
Manipulate[(1/3)NIntegrate[((penelobe[θ,ϕ,θ0,α])^3)*Sin[θ],{θ,0,(π/2)},{ϕ,-π,π}],{{α,1},0.001,1},{{θ0,(π/2)},0,(π/2)}]

SetDirectory["D:\\Workspaces\\Patapom.com\\Web\\blog\\patapom\\docs\\BRDF"]; 
countX=128;
countY=64;
fileName="VolumeInterPenetrationProper_"<>ToString[countX]<>"x"<>ToString[countY]<>".csv"

ta=ParallelTable[{N[θ0],N[α],
(1/3)NIntegrate[((penelobe[θ,ϕ,θ0,α])^3)*Sin[θ],{θ,0,(π/2)},{ϕ,-π,π}, PrecisionGoal->2,MaxRecursion->4,WorkingPrecision->$MachinePrecision]
}
,{θ0,0,(π/2),((π/2)/countX)}
,{α,0.001,1,((1-0.001)/countY)}]
Export[fileName,ta,"CSV"];

ta2=Import[fileName,"CSV"];
ta2=ToExpression[ta2]

Let’s create a normalized table by dividing by the total lobe volume

ta2[[129]][[65]]

lobeVolume[α_]:=((2π)/3)*NIntegrate[(((π*pdf[θ,α]))^3)*Sin[θ],{θ,0,(π/2)}]
taLobeVolumes=Table[lobeVolume[0.001+((1-0.001)/countY)*y],{y,0,countY}]

taNorm=Table[ta2[[x]][[y]]*{1,1,(1/(taLobeVolumes[[y]]))},{x,1,countX+1},{y,1,countY+1}]
ref=taNorm[[countX+1]][[countY+1]][[3]]

Show[ListPlot3D[Flatten[taNorm,1],PlotRange->{0,1},AxesLabel->{"θ","α","Lobe Penetration Ratio"}],Plot3D[ref,{x,0,(π/2)},{y,0,1},PlotStyle->Directive[Green,Opacity[0.5]]]]

ListInterpolation[Table[taNorm[[1]][[y]][[3]],{y,1,countY+1}],{0,1}]

Clear[α]
taRoughnessFunctions=Table[ListInterpolation[Table[taNorm[[θ0]][[y]][[3]],{y,1,countY+1}],{0,1}],{θ0,1,countX+1}];
Plot[taRoughnessFunctions[[64]][α],{α,0,1}]
ref
NSolve[taRoughnessFunctions[[64]][α]==ref,α]

Clear[α]
roughnessMap=Table[{N[((θ0-1)/countX)*(π/4)],α/.NSolve[taRoughnessFunctions[[θ0]][α]==ref,α][[1]]},{θ0,2,countX+1}];
roughnessMapCos=Table[{N[Cos[((θ0-1)/countX)*(π/4)]],α/.NSolve[taRoughnessFunctions[[θ0]][α]==ref,α][[1]]},{θ0,2,countX+1}];
ListLinePlot[roughnessMap]

roughnessMap=Delete[roughnessMap,100];
roughnessMapCos=Delete[roughnessMapCos,100];
ListLinePlot[roughnessMap]

Let’s fit this...

Clear[a,b,c,d,e,f,x,model]
model=0+b*x+c*(x^2)+d*(x^3)+e*(x^4)+f*(x^5);
model=((0+b*x+c*(x^2))/(1+e*x+0*f*(x^2)));
FindFit[roughnessMap,model,{a,b,c,d,e,f},x]
f2=Function[x,Evaluate[model/.FindFit[roughnessMap,model,{a,b,c,d,e,f},x]]];
Show[ListLinePlot[roughnessMap,PlotRange->{0,1}],Plot[f2[x],{x,0,(π/2)},PlotRange->{0,1},PlotStyle->Red]]

Sum[(((f2[roughnessMap[[i]][[1]]] - roughnessMap[[i]][[2]]))^2),{i,1,Length[roughnessMap]}]

Let’s plot and compare with the 2D model

c=(1/(Sqrt[2]));
sqAlpha:=Function[θ,
μ=Cos[θ];
(((μ^4)*c-(μ^2)*c-(Sqrt[μ*(((1-(μ^2)))^2)*c]))/((μ^4)*c-μ))
]
sqAlphaMu[μ_]:=(((μ^4)*c-(μ^2)*c-(Sqrt[μ*(((1-(μ^2)))^2)*c]))/((μ^4)*c-μ))
sqRoughnessFromMip[mip_]:=sqAlphaMu[mu[mip]] 
sqRoughnessFromMip3D[mip_]:=((f2[mu[mip]])^2)
sqRoughnessFromMip3D[mip_]:=((f2[ArcCos[mu[mip]]])^2)
(*Plot[{sqAlpha[(θ/2)],(Sqrt[sqAlpha[(θ/2)]])},{θ,0,(π/2)},PlotRange->{0,1}]*)
Plot[{((sqAlpha[(θ/2)])^(1/2)),f2[(θ/2)]},{θ,0,(π/2)},PlotRange->{0,1},AxesLabel->{"θ","α"}]
Plot[{(Sqrt[sqRoughnessFromMip[mip]]),(Sqrt[sqRoughnessFromMip3D[mip]])},{mip,0,mipMax},PlotRange->{0,1},AxesLabel->{"mip","α"}]

Let’s map the inverse function now:

data=Table[{ N[f2[ArcCos[μ]]],N[μ]},{μ,(1/(Sqrt[2])),1,((1-(1/(Sqrt[2])))/100)}];
ListLinePlot[data,PlotRange->Full]

Clear[a,b,c,d,e,f]
model=((a+b*x+c*(x^2))/(d+e*x+f*(x^2)))
model=((a+b*x+c*(x^2))/(a+e*x+f*(x^2)))
f3=Function[x,Evaluate[model/.FindFit[data,model,{a,b,c,d,e,f},x]]]
Sum[(((f3[data[[i]][[1]]] - data[[i]][[2]]))^2),{i,1,Length[data]}]
Show[ListLinePlot[data,PlotRange->{(1/(Sqrt[2])),1}],Plot[f3[x],{x,0,1},PlotStyle->Red]]
(Sqrt[2])f3[1]

Let’s plot and compare with the 2D model:

Plot[{f0[μ],f3[μ]},{μ,0,1},AxesLabel->{"μ","α"}]

Show a plot of mip level as a function of roughness