DGI initial state from baked (#99)

* A prototype of DGI radiance initialization from baked SH.

* Disable ZHWindow in ProbePropagationInitialize for now as it gives too smooth and dark radiance to axes.

* Minor fixes in propagation and a cleanup.

* Partner/crazyhunter 10.7.0 custom.1 dgi initial state from baked sh to from basis improvement (#107)

Improve ambient dice basis coefficients, improve basis normalization, and improve methods for projection from and to SH, choosing the right method for the context.

Dynamic GI: Ambient Dice to ZH: Switched over to using a raw ZH fit instead of deringed zh fit. It results in significantly lower error when compared to a ground truth monte carlo projection, which means slightly punchier + higher contrast + less washed out GI. And the ringing concern seems to not be an issue due to the way light is propagated. Across a number of test scenes no ringing is present.

Co-authored-by: Nicholas Brancaccio <[email protected]>

Author: AndrewSaraevUnity

com.unity.render-pipelines.high-definition/Runtime/Lighting/ProbeVolume/DynamicGI/AmbientDice.hlsl

@@ -54,6 +54,56 @@ float AmbientDiceAndClampedCosineProductIntegral(const in AmbientDice a, const i
     return res;
 }
 
+// Intentionally disabled deringing for zonal harmonics. In this context, the deringing results in loss of sharpness (over convolution)
+// without much gain - from our dynamic GI data, ringing does not occur in my tests, so does not need to be protected against after all.
+// #define ZONAL_HARMONICS_FROM_AMBIENT_DICE_DERING_ENABLED
+
+#if !defined(ZONAL_HARMONICS_FROM_AMBIENT_DICE_DERING_ENABLED)
+// https://www.desmos.com/calculator/abirkxdgko
+// Simple Least Squares Fit - raw data generated by directly projecting the ambient dice lobe to a zonal harmonic via integration.
+// Likely contains ringing for some sharpness levels.
+float ComputeZonalHarmonicC0FromAmbientDiceSharpness(float sharpness)
+{
+    // sharpness abs is to simply make the compiler happy.
+    return pow(abs(sharpness) * 1.62301 + 1.59682, -0.993255) * 2.83522 + -0.001;
+}
+
+float ComputeZonalHarmonicC1FromAmbientDiceSharpness(float sharpness)
+{
+    return exp2(-3.56165 * pow(abs(sharpness * -7.42401 + -13.5339), -0.998517)) * -9.1717 + 9.1715;
+}
+
+float ComputeZonalHarmonicC2FromAmbientDiceSharpness(float sharpness)
+{
+    float lhs = 0.239989 + 0.42846 * sharpness + -0.202951 * sharpness * sharpness + 0.0303908 * sharpness * sharpness * sharpness;
+    float rhs = exp2(-7.07205 * pow(abs(sharpness * 0.852362 + -0.469403), -0.898339)) * -0.557094 + 0.539681;
+    return sharpness > 2.33 ? rhs : lhs;
+}
+
+
+// https://www.desmos.com/calculator/1ajnhdbg6j
+// Simple Least Squares Fit - raw data generated by directly projecting the ambient dice lobe to a zonal harmonic via integration.
+// Likely contains ringing for some sharpness levels.
+float ComputeZonalHarmonicC0FromAmbientDiceWrappedSharpness(float sharpness)
+{
+    // sharpness abs is to simply make the compiler happy.
+    return pow(abs(sharpness) * 0.816074 + 0.809515, -0.99663) * 2.87866 + -0.001;
+}
+
+float ComputeZonalHarmonicC1FromAmbientDiceWrappedSharpness(float sharpness)
+{
+    return exp2(-7.01231 * pow(abs(sharpness * 1.36435 + -1.3829), -0.786179)) * -1.14392 + 1.04683;
+}
+
+float ComputeZonalHarmonicC2FromAmbientDiceWrappedSharpness(float sharpness)
+{
+    float lhs = -0.438588 + 0.542959 * sharpness + -0.112098 * sharpness * sharpness + 0.00800693 * sharpness * sharpness * sharpness;
+    float rhs = exp2(-6.39474 * pow(abs(sharpness * 0.488674 + -2.37692), -0.559034)) * -0.816382 + 0.473311;
+    return sharpness > 5.0 ? rhs : lhs;
+}
+
+
+#else
 // https://www.desmos.com/calculator/umjtgtzmk8
 // Fit to pre-deringed data.
 // The dering constraint is evaluated post diffuse brdf convolution.
@@ -75,6 +125,7 @@ float ComputeZonalHarmonicC2FromAmbientDiceSharpness(float sharpness)
     float rhs = exp2(-1.44747 * pow(abs(sharpness * 0.644014 + -0.188877), -0.94422)) * -0.970862 + 0.967661;
     return sharpness > 2.33 ? rhs : lhs;
 }
+#endif
 
 float3 ComputeZonalHarmonicFromAmbientDiceSharpness(float sharpness)
 {
@@ -85,4 +136,13 @@ float3 ComputeZonalHarmonicFromAmbientDiceSharpness(float sharpness)
     );
 }
 
+float3 ComputeZonalHarmonicFromAmbientDiceWrappedSharpness(float sharpness)
+{
+    return float3(
+        ComputeZonalHarmonicC0FromAmbientDiceWrappedSharpness(sharpness),
+        ComputeZonalHarmonicC1FromAmbientDiceWrappedSharpness(sharpness),
+        ComputeZonalHarmonicC2FromAmbientDiceWrappedSharpness(sharpness)
+    );
+}
+
 #endif

com.unity.render-pipelines.high-definition/Runtime/Lighting/ProbeVolume/DynamicGI/ProbePropagationAxes.compute

@@ -87,7 +87,7 @@ int ProbeCoordinateToIndex(uint3 probeCoordinate)
 
 float3 ProbeCoordinatesToWorldPosition(float3 probeCoordinates, float3x3 probeVolumeLtw)
 {
-    float3 localPosition = ((probeCoordinates * _ProbeVolumeDGIResolutionInverse) * 2.0 - 1.0) * _ProbeVolumeDGIBoundsExtents;
+    float3 localPosition = (((probeCoordinates + 0.5) * _ProbeVolumeDGIResolutionInverse) * 2.0 - 1.0) * _ProbeVolumeDGIBoundsExtents;
     return mul(localPosition, probeVolumeLtw) + GetAbsolutePositionWS(_ProbeVolumeDGIBoundsCenter);
 }
 
@@ -221,6 +221,77 @@ float3 ComputeNeighborBoundaryEdgeWorldPosition(float3 probeWorldPosition, int a
     return probeWorldPosition + neighborDirectionWS * neighborAxisLength * _ProbeVolumeDGIMaxNeighborDistance;
 }
 
+float3 PropagateLightFromNeighborSHWithSamplePeaks(SHIncomingIrradiance shIncomingIrradiance, NeighborAxisLookup neighborAxisLookup, float3x3 probeVolumeLtw, float hitWeight)
+{
+    float3 basisAxisHitDirectionWS = mul(neighborAxisLookup.neighborDirection, probeVolumeLtw);
+    return ProbeVolumeEvaluateIncomingRadiance(shIncomingIrradiance, basisAxisHitDirectionWS) * hitWeight;
+}
+
+float3 PropagateLightFromNeighborSHWithZHFit(SHIncomingIrradiance shIncomingIrradiance, NeighborAxisLookup neighborAxisLookup, float3x3 probeVolumeLtw, float hitWeight, BasisAxisHit basisAxisNeighborHit)
+{
+    float3 basisAxisHitDirectionWS = mul(neighborAxisLookup.neighborDirection, probeVolumeLtw);
+    ZHWindow basisAxisHitZHWindow = ComputeZHWindowFromBasisAxisHit(basisAxisNeighborHit);
+    SHIncomingIrradianceConvolveZHWindowWithoutDeltaFunction(shIncomingIrradiance, basisAxisHitZHWindow);
+    return ProbeVolumeEvaluateIncomingRadiance(shIncomingIrradiance, basisAxisHitDirectionWS) * hitWeight;
+}
+
+float3 PropagateLightFromNeighborSHWithMonteCarloNaiveProjection(SHIncomingIrradiance shIncomingIrradiance, NeighborAxisLookup neighborAxisLookup, float3x3 probeVolumeLtw, float hitWeight, BasisAxisHit basisAxisHit, BasisAxisHit basisAxisNeighborHit)
+{
+    BasisAxisHit basisAxisNeighborHitWS = basisAxisNeighborHit;
+    basisAxisNeighborHitWS.mean = mul(basisAxisNeighborHitWS.mean, probeVolumeLtw);
+
+    const float SOLID_ANGLE_SPHERE = 4.0 * PI;
+    const float GOLDEN_ANGLE = PI * (3.0 - sqrt(5.0));
+    const float GOLDEN_ANGLE_HALF = GOLDEN_ANGLE * 0.5;
+    const int SAMPLE_COUNT = 512;
+    const float SAMPLE_COUNT_INVERSE = 1.0 / (float)SAMPLE_COUNT;
+    
+    float3 incomingHitRadiance = 0.0;
+    for (int sampleIndex = 0; sampleIndex < SAMPLE_COUNT; ++sampleIndex)
+    {
+        float3 sampleDirectionWS;
+        {
+            float offset = (float)sampleIndex + 0.5;
+            float theta = offset * GOLDEN_ANGLE + GOLDEN_ANGLE_HALF;
+            float z = -(offset * SAMPLE_COUNT_INVERSE * 2.0 - 1.0);
+            float r = sqrt(1.0 - z * z);
+
+            float3 sampleDirectionBasisOS = float3(
+                r * cos(theta),
+                r * sin(theta),
+                z
+            );               
+
+            sampleDirectionWS = sampleDirectionBasisOS; // No rotation necessary, we are sampling the whole sphere.
+        }
+
+        const float DIFFERENTIAL_AREA = SOLID_ANGLE_SPHERE * SAMPLE_COUNT_INVERSE;
+        
+        float3 sampleIncomingRadiance = ProbeVolumeEvaluateIncomingRadiance(shIncomingIrradiance, sampleDirectionWS);
+
+        float3 basisAxisHitDirectionWS = mul(basisAxisHit.mean, probeVolumeLtw);
+        float3 basisAxisHitDirectionOS = basisAxisHit.mean;
+        basisAxisHit.mean = basisAxisHitDirectionWS;
+
+        incomingHitRadiance += sampleIncomingRadiance
+            * DIFFERENTIAL_AREA
+            * ComputeBasisAxisHitEvaluateFromDirection(basisAxisHit, sampleDirectionWS)
+            * ComputeBasisAxisHitEvaluateFromDirection(basisAxisNeighborHitWS, sampleDirectionWS);
+
+        basisAxisHit.mean = basisAxisHitDirectionOS; 
+    }
+
+    return incomingHitRadiance;
+}
+
+#define BASIS_FROM_SH_MODE_MONTE_CARLO_NAIVE_PROJECTION (0)
+#define BASIS_FROM_SH_MODE_SAMPLE_PEAKS (1)
+#define BASIS_FROM_SH_MODE_ZH_FIT (2)
+
+// See note in ProbePropagationInitialize.compute.
+// Sample Peaks turns out to be better than ZH Fit in practice for our combination of basis and SH order.
+#define BASIS_FROM_SH_MODE BASIS_FROM_SH_MODE_SAMPLE_PEAKS
+
 [numthreads(GROUP_SIZE, 1, 1)]
 void PropagateLight(uint3 id : SV_DispatchThreadID)
 {
@@ -286,8 +357,8 @@ void PropagateLight(uint3 id : SV_DispatchThreadID)
             if(hitIndex < (uint)_HitRadianceCacheAxisCount)
             {
                 // Hit
-                float weight = neighborAxisLookup.hitWeight * basisAxisNeighborHitIntegral;
-                incomingHitRadiance += DecodeRadiance(_HitRadianceCacheAxis[hitIndex]) * weight;
+                float hitWeight = neighborAxisLookup.hitWeight * basisAxisNeighborHitIntegral;
+                incomingHitRadiance += DecodeRadiance(_HitRadianceCacheAxis[hitIndex]) * hitWeight;
             }
             else
             {
@@ -304,8 +375,7 @@ void PropagateLight(uint3 id : SV_DispatchThreadID)
                         uint neighborProbeIndex = ProbeCoordinateToIndex(neighborProbeCoordinate);
                         float3 prevAxisRadiance = ReadPreviousPropagationAxis(neighborProbeIndex, axisIndex);
 
-                        float weight = neighborAxisLookup.propagationWeight;
-                        incomingMissRadiance += prevAxisRadiance * weight * InvalidScale(probeAxisValidity);
+                        incomingMissRadiance += prevAxisRadiance * neighborAxisLookup.propagationWeight * InvalidScale(probeAxisValidity);
                     #endif
                 }
 #if defined(SAMPLE_NEIGHBORS_DIRECTION_ONLY) || defined(SAMPLE_NEIGHBORS_POSITION_AND_DIRECTION)
@@ -320,17 +390,17 @@ void PropagateLight(uint3 id : SV_DispatchThreadID)
                         SHIncomingIrradiance shIncomingIrradiance = shIncomingIrradianceProbeCurrent;
                     #endif
 
-                    float weight = neighborAxisLookup.hitWeight * basisAxisNeighborHitIntegral;
-
-                    // Because convolution in SH space is cheap (zonal harmonic window component-wise multiplication), we are able to
-                    // more accurately capture the integral of incoming irradiance over the hit basis solid angle + weight.
-                    // It's important to note that this is an approximation - because our hit axis lobe is not a zonal harmonic - it is an SG or AmbientDice lobe,
-                    // but it can be approximately represented as an ZH.
-                    // This has been validated against a ground truth monte carlo sampling of the basis over the hemisphere and the results are almost completely indistinguishable.
-                    float3 basisAxisHitDirectionWS = mul(basisAxisHit.mean, probeVolumeLtw);
-                    ZHWindow basisAxisHitZHWindow = ComputeZHWindowFromBasisAxisHit(basisAxisHit);
-                    SHIncomingIrradianceConvolveZHWindowWithoutDeltaFunction(shIncomingIrradiance, basisAxisHitZHWindow);
-                    incomingHitRadiance += ProbeVolumeEvaluateIncomingRadiance(shIncomingIrradiance, basisAxisHitDirectionWS) * weight;
+                    float hitWeight = neighborAxisLookup.hitWeight * basisAxisNeighborHitIntegral;
+
+#if BASIS_FROM_SH_MODE == BASIS_FROM_SH_MODE_MONTE_CARLO_NAIVE_PROJECTION
+                    incomingHitRadiance += PropagateLightFromNeighborSHWithMonteCarloNaiveProjection(shIncomingIrradiance, neighborAxisLookup, probeVolumeLtw, hitWeight, basisAxisHit, basisAxisNeighborHit);
+
+#elif BASIS_FROM_SH_MODE == BASIS_FROM_SH_MODE_ZH_FIT
+                    incomingHitRadiance += PropagateLightFromNeighborSHWithZHFit(shIncomingIrradiance, neighborAxisLookup, probeVolumeLtw, hitWeight, basisAxisNeighborHit);
+
+#elif BASIS_FROM_SH_MODE == BASIS_FROM_SH_MODE_SAMPLE_PEAKS
+                    incomingHitRadiance += PropagateLightFromNeighborSHWithSamplePeaks(shIncomingIrradiance, neighborAxisLookup, probeVolumeLtw, hitWeight);
+#endif
                 }
 #endif
             }

com.unity.render-pipelines.high-definition/Runtime/Lighting/ProbeVolume/DynamicGI/ProbePropagationBasis.cs

@@ -73,11 +73,22 @@ public static void ComputeAmbientDiceSharpAmplitudeAndSharpnessFromAxisDirection
             componentNonZeroCount += Mathf.Abs(axisDirection.y) > 1e-3 ? 1 : 0;
             componentNonZeroCount += Mathf.Abs(axisDirection.z) > 1e-3 ? 1 : 0;
 
+#if true
+            // New coefficients tuned to to maximize unity across the sphere when integrating incoming radiance from monte carlo directions (i.e in a white furnace test).
+            float amplitudeEdge = 0.76f;
+            float amplitudeCorner = 0.52f;
+            float amplitudeOther = 0.73f;
+#else
+            float amplitudeEdge = 0.693f;
+            float amplitudeCorner = 0.3087f;
+            float amplitudeOther = 0.64575f;
+#endif
             amplitude = (componentNonZeroCount == 3)
-                ? 0.3087f // diagonal
+                ? amplitudeCorner // diagonal
                 : ((componentNonZeroCount == 2)
-                    ? 0.693f // edge
-                    : 0.64575f); // center
+                    ? amplitudeEdge // edge
+                    : amplitudeOther); // center
+
             sharpness = (componentNonZeroCount == 3)
                 ? 9f // diagonal
                 : 6f; // center

com.unity.render-pipelines.high-definition/Runtime/Lighting/ProbeVolume/DynamicGI/ProbePropagationBasis.hlsl

@@ -213,4 +213,24 @@ ZHWindow ComputeZHWindowFromBasisAxisHit(BasisAxisHit basisAxisHit)
     return zhWindow;
 }
 
+ZHWindow ComputeZHWindowFromBasisAxisMiss(BasisAxisMiss basisAxisMiss)
+{
+    ZHWindow zhWindow;
+    ZERO_INITIALIZE(ZHWindow, zhWindow);
+
+#if (defined(BASIS_PROPAGATION_OVERRIDE_NONE) && defined(BASIS_SPHERICAL_GAUSSIAN)) || defined(BASIS_PROPAGATION_OVERRIDE_SPHERICAL_GAUSSIAN)
+    zhWindow = ZHWindowComputeFromSphericalGaussian(basisAxisMiss.sharpness);
+#elif defined(BASIS_PROPAGATION_OVERRIDE_NONE) && defined(BASIS_SPHERICAL_GAUSSIAN_WINDOWED)
+    zhWindow = ZHWindowComputeFromSphericalGaussianCosineWindow(basisAxisMiss.sharpness);
+#elif defined(BASIS_PROPAGATION_OVERRIDE_NONE) && defined(BASIS_AMBIENT_DICE)
+    zhWindow = ZHWindowComputeFromAmbientDiceSharpness(basisAxisMiss.sharpness);
+#elif defined(BASIS_PROPAGATION_OVERRIDE_AMBIENT_DICE_WRAPPED_SOFTER) || defined(BASIS_PROPAGATION_OVERRIDE_AMBIENT_DICE_WRAPPED_SUPER_SOFT) || defined(BASIS_PROPAGATION_OVERRIDE_AMBIENT_DICE_WRAPPED_ULTRA_SOFT)
+    zhWindow = ZHWindowComputeFromAmbientDiceWrappedSharpness(basisAxisMiss.sharpness);
+#else
+#error "Undefined Probe Propagation Basis"
+#endif
+
+    return zhWindow;
+}
+
 #endif // endof PROBE_PROPAGATION_BASIS

com.unity.render-pipelines.high-definition/Runtime/Lighting/ProbeVolume/DynamicGI/ProbePropagationBasisInternal.hlsl

@@ -44,11 +44,23 @@ void ComputeAmbientDiceSharpAmplitudeAndSharpnessFromAxisDirectionBasis26Fit(out
     componentNonZeroCount += abs(axisDirection.y) > 1e-3 ? 1 : 0;
     componentNonZeroCount += abs(axisDirection.z) > 1e-3 ? 1 : 0;
 
+#if 1
+    // New coefficients tuned to to maximize unity across the sphere when integrating incoming radiance from monte carlo directions (i.e in a white furnace test).
+    float amplitudeEdge = 0.76;
+    float amplitudeDiagonal = 0.52;
+    float amplitudeCenter = 0.73;
+#else
+    float amplitudeEdge = 0.693;
+    float amplitudeDiagonal = 0.3087;
+    float amplitudeCenter = 0.64575;
+#endif
+
     amplitude = (componentNonZeroCount == 3)
-        ? 0.3087 // diagonal
+        ? amplitudeDiagonal // diagonal
         : ((componentNonZeroCount == 2)
-            ? 0.693 // edge
-            : 0.64575); // center
+            ? amplitudeEdge // edge
+            : amplitudeCenter); // center
+
     sharpness = (componentNonZeroCount == 3)
         ? 9.0 // diagonal
         : ((componentNonZeroCount == 2)
@@ -183,4 +195,16 @@ ZHWindow ZHWindowComputeFromAmbientDiceSharpness(float sharpness)
     return zhWindow;
 }
 
+ZHWindow ZHWindowComputeFromAmbientDiceWrappedSharpness(float sharpness)
+{
+    ZHWindow zhWindow;
+    
+    float3 coefficients = ComputeZonalHarmonicFromAmbientDiceWrappedSharpness(sharpness);
+    zhWindow.data[0] = coefficients.x;
+    zhWindow.data[1] = coefficients.y;
+    zhWindow.data[2] = coefficients.z;
+
+    return zhWindow;
+}
+
 #endif