Improve the quality of importance sampling

Assets/ScriptableRenderLoop/HDRenderPipeline/HDRenderPipeline.cs

                 m_ShadowSettings.directionalLightCascades = new Vector3(m_CommonSettings.shadowCascadeSplit0, m_CommonSettings.shadowCascadeSplit1, m_CommonSettings.shadowCascadeSplit2);
                 m_ShadowSettings.maxShadowDistance = m_CommonSettings.shadowMaxDistance;
 
-                sssParameters.profiles[0].filter1Variance  = m_CommonSettings.sssProfileFilter1Variance;
-                sssParameters.profiles[0].filter2Variance  = m_CommonSettings.sssProfileFilter2Variance;
+                sssParameters.profiles[0].filterVariance1  = m_CommonSettings.sssProfileFilterVariance1;
+                sssParameters.profiles[0].filterVariance2  = m_CommonSettings.sssProfileFilterVariance2;
                 sssParameters.profiles[0].filterLerpWeight = m_CommonSettings.sssProfileFilterLerpWeight;
                 sssParameters.bilateralScale               = m_CommonSettings.sssBilateralScale;
             }

Assets/ScriptableRenderLoop/HDRenderPipeline/SceneSettings/CommonSettings.cs

         public float shadowCascadeSplit2 { set { m_ShadowCascadeSplit2 = value; OnValidate(); } get { return m_ShadowCascadeSplit2; } }
 
         // Subsurface scattering
-        [SerializeField] Color m_SssProfileFilter1Variance  = SubsurfaceScatteringProfile.Default.filter1Variance;
-        [SerializeField] Color m_SssProfileFilter2Variance  = SubsurfaceScatteringProfile.Default.filter2Variance;
+        [SerializeField] Color m_SssProfileFilterVariance1  = SubsurfaceScatteringProfile.Default.filterVariance1;
+        [SerializeField] Color m_SssProfileFilterVariance2  = SubsurfaceScatteringProfile.Default.filterVariance2;
-        public Color sssProfileFilter1Variance  { set { m_SssProfileFilter1Variance  = value; OnValidate(); } get { return m_SssProfileFilter1Variance; } }
-        public Color sssProfileFilter2Variance  { set { m_SssProfileFilter2Variance  = value; OnValidate(); } get { return m_SssProfileFilter2Variance; } }
+        public Color sssProfileFilterVariance1  { set { m_SssProfileFilterVariance1  = value; OnValidate(); } get { return m_SssProfileFilterVariance1; } }
+        public Color sssProfileFilterVariance2  { set { m_SssProfileFilterVariance2  = value; OnValidate(); } get { return m_SssProfileFilterVariance2; } }
         public float sssProfileFilterLerpWeight { set { m_SssProfileFilterLerpWeight = value; OnValidate(); } get { return m_SssProfileFilterLerpWeight; } }
         public float sssBilateralScale          { set { m_SssBilateralScale          = value; OnValidate(); } get { return m_SssBilateralScale; } }
 
             m_ShadowCascadeSplit1 = Mathf.Clamp01(m_ShadowCascadeSplit1);
             m_ShadowCascadeSplit2 = Mathf.Clamp01(m_ShadowCascadeSplit2);
 
-            m_SssProfileFilter1Variance.r = Mathf.Max(0.05f, m_SssProfileFilter1Variance.r);
-            m_SssProfileFilter1Variance.g = Mathf.Max(0.05f, m_SssProfileFilter1Variance.g);
-            m_SssProfileFilter1Variance.b = Mathf.Max(0.05f, m_SssProfileFilter1Variance.b);
-            m_SssProfileFilter1Variance.a = 0.0f;
-            m_SssProfileFilter2Variance.r = Mathf.Max(0.05f, m_SssProfileFilter2Variance.r);
-            m_SssProfileFilter2Variance.g = Mathf.Max(0.05f, m_SssProfileFilter2Variance.g);
-            m_SssProfileFilter2Variance.b = Mathf.Max(0.05f, m_SssProfileFilter2Variance.b);
-            m_SssProfileFilter2Variance.a = 0.0f;
+            m_SssProfileFilterVariance1.r = Mathf.Max(0.05f, m_SssProfileFilterVariance1.r);
+            m_SssProfileFilterVariance1.g = Mathf.Max(0.05f, m_SssProfileFilterVariance1.g);
+            m_SssProfileFilterVariance1.b = Mathf.Max(0.05f, m_SssProfileFilterVariance1.b);
+            m_SssProfileFilterVariance1.a = 0.0f;
+            m_SssProfileFilterVariance2.r = Mathf.Max(0.05f, m_SssProfileFilterVariance2.r);
+            m_SssProfileFilterVariance2.g = Mathf.Max(0.05f, m_SssProfileFilterVariance2.g);
+            m_SssProfileFilterVariance2.b = Mathf.Max(0.05f, m_SssProfileFilterVariance2.b);
+            m_SssProfileFilterVariance2.a = 0.0f;
             m_SssProfileFilterLerpWeight  = Mathf.Clamp01(m_SssProfileFilterLerpWeight);
             m_SssBilateralScale           = Mathf.Clamp01(m_SssBilateralScale);
 

Assets/ScriptableRenderLoop/HDRenderPipeline/SceneSettings/Editor/CommonSettingsEditor.cs

             public readonly GUIContent[] shadowSplits = new GUIContent[] { new GUIContent("Split 0"), new GUIContent("Split 1"), new GUIContent("Split 2") };
 
             public readonly GUIContent sssCategory                = new GUIContent("Subsurface scattering");
-            public readonly GUIContent sssProfileFilter1Variance  = new GUIContent("SSS profile filter #1 variance", "Determines the shape of the 1st Gaussian filter. Increases the strength of the blur of the corresponding color channel.");
-            public readonly GUIContent sssProfileFilter2Variance  = new GUIContent("SSS profile filter #2 variance", "Determines the shape of the 2nd Gaussian filter. Increases the strength of the blur of the corresponding color channel.");
+            public readonly GUIContent sssProfileFilterVariance1  = new GUIContent("SSS profile filter variance #1", "Determines the shape of the 1st Gaussian filter. Increases the strength of the blur of the corresponding color channel.");
+            public readonly GUIContent sssProfileFilterVariance2  = new GUIContent("SSS profile filter variance #2", "Determines the shape of the 2nd Gaussian filter. Increases the strength of the blur of the corresponding color channel.");
             public readonly GUIContent sssProfileFilterLerpWeight = new GUIContent("SSS profile filter interpolation", "Controls linear interpolation between the two Gaussian filters.");
             public readonly GUIContent sssBilateralScale          = new GUIContent("SSS bilateral filtering scale", "Larger values make the filter more tolerant to depth differences.");
         }
         private SerializedProperty[] m_ShadowCascadeSplits = new SerializedProperty[3];
 
         // Subsurface scattering
-        private SerializedProperty m_SssProfileFilter1Variance;
-        private SerializedProperty m_SssProfileFilter2Variance;
+        private SerializedProperty m_SssProfileFilterVariance1;
+        private SerializedProperty m_SssProfileFilterVariance2;
         private SerializedProperty m_SssProfileFilterLerpWeight;
         private SerializedProperty m_SssBilateralScale;
 
             m_SkyRendererTypeValues.Add(m_SkyRendererTypeValues.Count);
             m_SkyRendererTypes.Add(null);
 
-            m_SssProfileFilter1Variance  = serializedObject.FindProperty("m_SssProfileFilter1Variance");
-            m_SssProfileFilter2Variance  = serializedObject.FindProperty("m_SssProfileFilter2Variance");
+            m_SssProfileFilterVariance1  = serializedObject.FindProperty("m_SssProfileFilterVariance1");
+            m_SssProfileFilterVariance2  = serializedObject.FindProperty("m_SssProfileFilterVariance2");
             m_SssProfileFilterLerpWeight = serializedObject.FindProperty("m_SssProfileFilterLerpWeight");
             m_SssBilateralScale          = serializedObject.FindProperty("m_SssBilateralScale");
         }
         {
             EditorGUILayout.LabelField(styles.sssCategory);
             EditorGUI.indentLevel++;
-            EditorGUILayout.PropertyField(m_SssProfileFilter1Variance,  styles.sssProfileFilter1Variance);
-            EditorGUILayout.PropertyField(m_SssProfileFilter2Variance,  styles.sssProfileFilter2Variance);
+            EditorGUILayout.PropertyField(m_SssProfileFilterVariance1,  styles.sssProfileFilterVariance1);
+            EditorGUILayout.PropertyField(m_SssProfileFilterVariance2,  styles.sssProfileFilterVariance2);
             EditorGUILayout.PropertyField(m_SssProfileFilterLerpWeight, styles.sssProfileFilterLerpWeight);
             EditorGUILayout.PropertyField(m_SssBilateralScale,          styles.sssBilateralScale);
             EditorGUI.indentLevel--;

Assets/ScriptableRenderLoop/HDRenderPipeline/SceneSettings/SubsurfaceScatteringParameters.cs

     {
         public const int numSamples = 7;
 
-        Color     m_filter1Variance;
-        Color     m_filter2Variance;
+        Color     m_filterVariance1;
+        Color     m_filterVariance2;
         float     m_filterLerpWeight;
         Vector4[] m_filterKernel;
         bool      m_kernelNeedsUpdate;
-        public Color filter1Variance
+        public Color filterVariance1
-            get { return m_filter1Variance; }
-            set { if (m_filter1Variance != value) { m_filter1Variance = value; m_kernelNeedsUpdate = true; } }
+            get { return m_filterVariance1; }
+            set { if (m_filterVariance1 != value) { m_filterVariance1 = value; m_kernelNeedsUpdate = true; } }
-        public Color filter2Variance
+        public Color filterVariance2
-            get { return m_filter2Variance; }
-            set { if (m_filter2Variance != value) { m_filter2Variance = value; m_kernelNeedsUpdate = true; } }
+            get { return m_filterVariance2; }
+            set { if (m_filterVariance2 != value) { m_filterVariance2 = value; m_kernelNeedsUpdate = true; } }
         }
 
         public float filterLerpWeight
             get
             {
                 SubsurfaceScatteringProfile profile = new SubsurfaceScatteringProfile();
-                profile.filter1Variance  = new Color(0.3f, 0.3f, 0.3f, 0.0f);
-                profile.filter2Variance  = new Color(1.0f, 1.0f, 1.0f, 0.0f);
+                profile.filterVariance1  = new Color(0.3f, 0.3f, 0.3f, 0.0f);
+                profile.filterVariance2  = new Color(1.0f, 1.0f, 1.0f, 0.0f);
                 profile.filterLerpWeight = 0.5f;
                 profile.ComputeKernel();
                 return profile;
-        static float EvaluateZeroMeanGaussian(float x, float variance)
+        static float Gaussian(float x, float variance)
-        static float EvaluateGaussianCombination(float x, float variance1, float variance2, float lerpWeight)
+        static float GaussianCombination(float x, float variance1, float variance2, float lerpWeight)
-            return Mathf.Lerp(EvaluateZeroMeanGaussian(x, variance1),
-                              EvaluateZeroMeanGaussian(x, variance2), lerpWeight);
+            return Mathf.Lerp(Gaussian(x, variance1), Gaussian(x, variance2), lerpWeight);
-        static double RationalApproximation(double t)
+        static float RationalApproximation(float t)
-            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);
+            float[] c = {2.515517f, 0.802853f, 0.010328f};
+            float[] d = {1.432788f, 0.189269f, 0.001308f};
+            return t - ((c[2] * t + c[1]) * t + c[0]) / (((d[2] * t + d[1]) * t + d[0]) * t + 1.0f);
-        static double NormalCDFInverse(double p, double stdDeviation)
+        static float NormalCdfInverse(float p, float stdDev)
-            double x;
+            float x;
-                x = -RationalApproximation(Math.Sqrt(-2.0 * Math.Log(p)));
+                x = -RationalApproximation(Mathf.Sqrt(-2.0f * Mathf.Log(p)));
-                x = RationalApproximation(Math.Sqrt(-2.0*Math.Log(1.0 - p)));
+                x = RationalApproximation(Mathf.Sqrt(-2.0f * Mathf.Log(1.0f - p)));
-            return x * stdDeviation;
+            return x * stdDev;
+        }
+
+        // Ref: https://en.wikipedia.org/wiki/Halton_sequence
+        static float VanDerCorput(uint b, uint i)
+        {
+            float r = 0;
+            float f = 1;
+            
+            while (i > 0) 
+            {
+                f = f / b;
+                r = r + f * (i % b);
+                i = i / b;
+            }
+
+            return r;
-        double VanDerCorputBase2(uint bits)
+        static float VanDerCorputBase2(uint i)
-            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
+            i = (i << 16) | (i >> 16);
+            i = ((i & 0x00ff00ff) << 8) | ((i & 0xff00ff00) >> 8);
+            i = ((i & 0x0f0f0f0f) << 4) | ((i & 0xf0f0f0f0) >> 4);
+            i = ((i & 0x33333333) << 2) | ((i & 0xcccccccc) >> 2);
+            i = ((i & 0x55555555) << 1) | ((i & 0xaaaaaaaa) >> 1);
+
+            return i * (1.0f / 4294967296);
         }
 
         void ComputeKernel()
             // It is separable by design, but generally not radially symmmetric.
 
             // 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));
+            float maxVariance1 = Mathf.Max(m_filterVariance1.r, m_filterVariance1.g, m_filterVariance1.b);
+            float maxVariance2 = Mathf.Max(m_filterVariance2.r, m_filterVariance2.g, m_filterVariance2.b);
-            for (uint i = 0; i < numSamples; ++i)
+            float sd = Mathf.Lerp(Mathf.Sqrt(maxVariance1), Mathf.Sqrt(maxVariance2), m_filterLerpWeight);
+
+            Vector3 weightSum = new Vector3(0, 0, 0); 
+
+            for (uint i = 0; i < numSamples; i++)
-                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;
+                float u   = VanDerCorputBase2(i + 1);
+                float pos = NormalCdfInverse(u, sd);
+                float pdf = Gaussian(pos, sd * sd);
+
+                Vector3 val;
+                val.x = GaussianCombination(pos, m_filterVariance1.r, m_filterVariance2.r, m_filterLerpWeight);
+                val.y = GaussianCombination(pos, m_filterVariance1.g, m_filterVariance2.g, m_filterLerpWeight);
+                val.z = GaussianCombination(pos, m_filterVariance1.b, m_filterVariance2.b, m_filterLerpWeight);
+
+                m_filterKernel[i].x = val.x / (pdf * numSamples);
+                m_filterKernel[i].y = val.y / (pdf * numSamples);
+                m_filterKernel[i].z = val.z / (pdf * numSamples);
+
+                weightSum.x += m_filterKernel[i].x;
+                weightSum.y += m_filterKernel[i].y;
+                weightSum.z += m_filterKernel[i].z;
+            }
+
+            // Renormalize the weights to conserve energy.
+            for (uint i = 0; i < numSamples; i++)
+            {
+                m_filterKernel[i].x *= 1.0f / weightSum.x;
+                m_filterKernel[i].y *= 1.0f / weightSum.y;
+                m_filterKernel[i].z *= 1.0f / weightSum.z;
             }
 
             m_kernelNeedsUpdate = false;