I'm not entirely new to Mathematica but I use it only very sporadically. I'm foremost a C/C++ programmer by profession.

I recently had a few weeks of spare time and launched Mathematica to do some computation on the GGX BRDF model. I wanted to compute the Fresnel Albedo function for all values of roughnesses and incident angles. The Albedo is a directional-hemispherical integral.

I got pretty much everything running the way I wanted but at very small roughness values, I start getting instabilities in the Albedo computations.

To illustrate, I boiled the code with the instability issue down to this:

(* Utility function for creating a unit 3D vector in the upper \
hemisphere from a polar coordinates with radius 1 *)
PolarToNormal[theta_, phi_] := {Sin[theta] Cos[phi], Sin[theta] Sin[phi], Cos[theta]};

(* Utility function for creating a unit 3D vector in the upper hemisphere from a cos *)
CosToNormal[cos_] := {Sqrt[1 - cos^2], 0, cos};

(* Utility function for computing the cosine of the half vector \
between an incident and an exitant vectors *)
cosM[NdotI_, theta_, phi_] := Normalize[
   CosToNormal[NdotI] + PolarToNormal[theta, phi]][[3]];

(* The GGX normal distribution function *)
GGXD[alpha_, NdotI_, theta_, phi_] :=
   alpha^2/(Pi (cosM[NdotI, theta, phi]^2 (alpha^2 - 1) + 1)^2);

(* Plot of the exitant hemispherical integral of GGX D for incident cosines from 0 to 1 *)
   GGXD[0.00005, cosI, theta, phi] Sin[theta] Cos[theta], \
   {phi, -Pi, Pi}, {theta, 0, Pi/2}], {cosI, 0.000000001, 1}]

And this is the plot I get:

And this is the plot I get

How can I debug this?

How can I find which part of the integration process is causing this instability?

Are there any Mathematica parameters such as WorkingPrecision, PrecisionGoal, etc. that I can tweak to fix this?

BTW, Is there any way I can find out what is the current working precision on my installed Mathematica?


Under John Boyd's MAG philosophy (make-a-graph), we should examine the graph of a function when a numerical routine behaves oddly. The support for the integral mainly lies in two spikes:

Block[{cosI = 5/10},
 Plot3D[GGXD[0.00005, cosI, theta, phi] Sin[theta] Cos[theta],
  {phi, -Pi, Pi}, {theta, 0, Pi/2}, PlotPoints -> 25, 
  PlotRange -> All, PlotLabel -> N@ArcCos[cosI]]

Mathematica graphics

These spikes might be missed by the default sampling. We ca add points to the iterators to ensure that they will be found. The occur where the denominator of GGXD[5/100000, cosI, theta, phi] has a minimum. (You can get the denominator with GGXD[5/100000, cosI, theta, phi] /. Abs -> RealAbs // Simplify // Denominator.)

 D[(1 - (399999999 (cosI + Cos[theta])^2)/(800000000 (1 + cosI Cos[theta] + 
         Sqrt[1 - cosI^2] Cos[phi] Sin[theta])))^2 /. phi -> Pi, 
   theta] == 0, theta]

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

{{theta -> -ArcCos[-cosI]}, {theta -> ArcCos[cosI]},
{theta -> -ArcCos[(400000001 cosI - 40000 Sqrt[-1 + cosI^2])/399999999]},
{theta -> ArcCos[(400000001 cosI - 40000 Sqrt[-1 + cosI^2])/399999999]},
{theta -> -ArcCos[(400000001 cosI + 40000 Sqrt[-1 + cosI^2])/399999999]},
{theta -> ArcCos[(400000001 cosI + 40000 Sqrt[-1 + cosI^2])/399999999]}}

The second solution, theta == ArcCos[cosI], is the location of the spike near phi == {-Pi, Pi}. The width in the phi direction of the spikes is quite small, and adding two points near phi == {-Pi, Pi} ensures that the spike will be adequately sampled by NIntegrate.

  GGXD[0.00005, cosI, theta, phi] Sin[theta] Cos[
    theta], {phi, -Pi, -Pi + 10^-3, Pi - 10^-3, Pi}, {theta, 0, 
   ArcCos[cosI], Pi/2}], {cosI, 0.000000001, 1}, MaxRecursion -> 2]

Mathematica graphics

  • $\begingroup$ Thanks Michael. I see what you are doing. Very good debugging technique and I reused it on the same function but on the cos theta dimension. I have a question: The change you suggest will add more samples by overlapping the troublesome region. This will change the probability for the samples in this region. I guess Mathematica will detect the overlap and adjust the samples probabilities accordingly? Is that correct? $\endgroup$ – Yves Poissant Jun 17 '18 at 23:24
  • $\begingroup$ @YvesPoissant Yes NIntegrate computes the integral by adjusting the rule to the density of the sampling in each subregion. (In fact, it almost certainly subdivides further some regions automatically; the process is called "adaptive sampling." It will increase sampling in regions in order to reduce the error in the parts that have the worst error.) $\endgroup$ – Michael E2 Jun 19 '18 at 2:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.