Importance sample two Gaussians based on the interpolation weight

Assets/ScriptableRenderLoop/HDRenderPipeline/SceneSettings/SubsurfaceScatteringParameters.cs

+using System;
+
-    [System.Serializable]
+    [Serializable]
     public class SubsurfaceScatteringProfile
     {
         public const int numSamples = 7;
                               EvaluateZeroMeanGaussian(x, variance2), lerpWeight);
         }
 
+        static double RationalApproximation(double t)
+        {
+            // Abramowitz and Stegun formula 26.2.23.
+            // The absolute value of the error should be less than 4.5 e-4.
+            double[] c = {2.515517, 0.802853, 0.010328};
+            double[] d = {1.432788, 0.189269, 0.001308};
+            return t - ((c[2] * t + c[1]) * t + c[0]) / (((d[2] * t + d[1]) * t + d[0]) * t + 1.0);
+        }
+ 
+        // Ref: https://www.johndcook.com/blog/csharp_phi_inverse/
+        static double NormalCDFInverse(double p, double stdDeviation)
+        {
+            double x;
+
+            if (p < 0.5)
+            {
+                // F^-1(p) = - G^-1(p)
+                x = -RationalApproximation(Math.Sqrt(-2.0 * Math.Log(p)));
+            }
+            else
+            {
+                // F^-1(p) = G^-1(1-p)
+                x = RationalApproximation(Math.Sqrt(-2.0*Math.Log(1.0 - p)));
+            }
+
+            return x * stdDeviation;
+        }
+
+        // Ref: http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html
+        double VanDerCorputBase2(uint bits)
+        {
+            bits = (bits << 16) | (bits >> 16);
+            bits = ((bits & 0x00ff00ff) << 8) | ((bits & 0xff00ff00) >> 8);
+            bits = ((bits & 0x0f0f0f0f) << 4) | ((bits & 0xf0f0f0f0) >> 4);
+            bits = ((bits & 0x33333333) << 2) | ((bits & 0xcccccccc) >> 2);
+            bits = ((bits & 0x55555555) << 1) | ((bits & 0xaaaaaaaa) >> 1);
+            return bits * 2.3283064365386963e-10; // 0x100000000
+        }
+
+            if (m_filterKernel == null)
+            {
+                m_filterKernel = new Vector4[numSamples];
+            }
+
             // Our goal is to blur the image using a filter which is represented
             // as a product of a linear combination of two normalized 1D Gaussians
             // as suggested by Jimenez et al. in "Separable Subsurface Scattering".
             // The resulting filter function is a non-Gaussian PDF.
             // It is separable by design, but generally not radially symmmetric.
 
-            // TODO: importance-sample the truncated normal distribution:
-            // https://people.sc.fsu.edu/~jburkardt/cpp_src/truncated_normal/truncated_normal.html
-            // ALternatively, use a quadrature rule for the truncated normal PDF:
-            // https://people.sc.fsu.edu/~jburkardt/cpp_src/truncated_normal_rule/truncated_normal_rule.html
-
-            // For now, we use an ad-hoc approach.
-            // We truncate the distribution at the radius of 3 standard deviations.
-            float maxRadius1 = 3 * Mathf.Sqrt(Mathf.Max(m_filter1Variance.r, m_filter1Variance.g, m_filter1Variance.b));
-            float maxRadius2 = 3 * Mathf.Sqrt(Mathf.Max(m_filter2Variance.r, m_filter2Variance.g, m_filter2Variance.b));
-            float radius     = Mathf.Lerp(maxRadius1, maxRadius2, m_filterLerpWeight);
-
-            // We compute sample positions and weights using Gauss�Legendre quadrature.
-            // Ref: http://keisan.casio.com/exec/system/1329114617
-            float[] unitAbscissae = { 0.0f,        0.40584515f, -0.40584515f, 0.74153118f, -0.74153118f, 0.94910791f, -0.94910791f };
-            float[] unitWeights   = { 0.41795918f, 0.38183005f,  0.38183005f, 0.27970539f,  0.27970539f, 0.12948496f,  0.12948496f };
-
-            if (m_filterKernel == null)
-            {
-                m_filterKernel = new Vector4[numSamples];
-            }
+            // Find the widest Gaussian across 3 color channels.
+            float maxStdDev1 = Mathf.Sqrt(Mathf.Max(m_filter1Variance.r, m_filter1Variance.g, m_filter1Variance.b));
+            float maxStdDev2 = Mathf.Sqrt(Mathf.Max(m_filter2Variance.r, m_filter2Variance.g, m_filter2Variance.b));
-            for (int i = 0; i < numSamples; ++i)
+            // Importance sample two Gaussians based on the interpolation weight.
+            for (uint i = 0; i < numSamples; ++i)
-                // Perform the change of interval: {a, b} = {-radius, radius}.
-                // Ref: https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval
-                float weight   = radius * unitWeights[i];
-                float position = radius * unitAbscissae[i];
-
-                m_filterKernel[i].x = weight * EvaluateGaussianCombination(position, m_filter1Variance.r, m_filter2Variance.r, m_filterLerpWeight);
-                m_filterKernel[i].y = weight * EvaluateGaussianCombination(position, m_filter1Variance.g, m_filter2Variance.g, m_filterLerpWeight);
-                m_filterKernel[i].z = weight * EvaluateGaussianCombination(position, m_filter1Variance.b, m_filter2Variance.b, m_filterLerpWeight);
-                m_filterKernel[i].w = position;
+                double u1 = (i + 0.5) / numSamples;
+                double u2 = VanDerCorputBase2(i + 1);
+                double sd = u1 < m_filterLerpWeight ? maxStdDev1 : maxStdDev2;
+                float pos = (float)NormalCDFInverse(u2, sd);
+                // Since our filter is normalized, f(x) / p(x) = 1.
+                m_filterKernel[i].x = 1.0f / numSamples;
+                m_filterKernel[i].y = 1.0f / numSamples;
+                m_filterKernel[i].z = 1.0f / numSamples;
+                m_filterKernel[i].w = pos;
             }
 
             m_kernelNeedsUpdate = false;