Improve the quality of trilinear filtering of density volume textures

com.unity.render-pipelines.core/CoreRP/ShaderLibrary/Common.hlsl

     return ComputeTextureLOD(uv);
 }
 
-float ComputeTextureLOD(float3 Px, float3 Py, float3 Pz)
+float ComputeTextureLOD(float3 duvw_dx, float3 duvw_dy, float3 duvw_dz, float scale)
-    float d = max(dot(Px, Px), max(dot(Py, Py), dot(Pz, Pz)));
-    return max(0.0, 0.5 * log2(d));
+    float d = Max3(dot(duvw_dx, duvw_dx), dot(duvw_dy, duvw_dy), dot(duvw_dz, duvw_dz));
+    return max(0.0, 0.5 * log2(d * (scale * scale)));
 }
 
 

com.unity.render-pipelines.high-definition/CHANGELOG.md

 - SSS and Transmission code have been refactored to be able to share it between various material. Guidelines are in SubsurfaceScattering.hlsl
 - Change code in area light with LTC for Lit shader. Magnitude is now take from FGD texture instead of a separate texture
 - Improve camera relative rendering: We now apply camera translation on the model matrix, so before the TransformObjectToWorld(). Note: unity_WorldToObject and unity_ObjectToWorld must never be used directly.
-- Rename positionWS to positionRWS (Camera relative world position) at a lot of places (mainly in interpolator and FragInputs). In case of custom shader user will be required to update their code. 
+- Rename positionWS to positionRWS (Camera relative world position) at a lot of places (mainly in interpolator and FragInputs). In case of custom shader user will be required to update their code.
+- Improve the quality of trilinear filtering of density volume textures.
 
 ### Fixed
 - Fix contact shadows applied on transmission

com.unity.render-pipelines.high-definition/HDRP/Lighting/Volumetrics/VolumeVoxelization.compute

 //--------------------------------------------------------------------------------------------------
 // Implementation
 //--------------------------------------------------------------------------------------------------
-float SampleVolumeMask(DensityVolumeData volumeData, float3 voxelCenterUV, float3 VxUV, float3 VyUV, float3 VzUV)
+float SampleVolumeMask(DensityVolumeData volumeData, float3 voxelCenterUV, float3 duvw_dx, float3 duvw_dy, float3 duvw_dz)
-    float clampBorder = 0.5f * _VolumeMaskDimensions.y;
 
     //scale and bias the UVs and then take fractional part, will be in [0,1] range
     voxelCenterUV = frac(voxelCenterUV * volumeData.textureTiling + volumeData.textureScroll);
 
     // TODO: bugfix.
     // Note that this clamping to edge doesn't quite work.
-    // First of all, the distance to the edge should depend on the LoD. In the current implementation, it does not.
+    // First of all, the distance to the edge should depend on the LoD.
-    float lod = ComputeTextureLOD(VxUV * _VolumeMaskDimensions.z, VyUV * _VolumeMaskDimensions.z, VzUV * _VolumeMaskDimensions.z);
-    voxelCenterUV.z = clamp(voxelCenterUV.z, offset + clampBorder, offset + _VolumeMaskDimensions.x - clampBorder);    
+    // For now, we choose the second option.
+    float lod         = ComputeTextureLOD(duvw_dx, duvw_dy, duvw_dz, _VolumeMaskDimensions.z);
+    float lodScale    = exp2(ceil(lod));
+    float clampBorder = 0.5f * lodScale * _VolumeMaskDimensions.y; // 0.5 * (lodScale / textureSize)
+    voxelCenterUV.z   = clamp(voxelCenterUV.z, offset + clampBorder, offset + _VolumeMaskDimensions.x - clampBorder);
     
     // TODO: expose the LoD bias parameter.
     // TODO: can we use anisotropic filtering to improve the quality? SampleGrad() doesn't support ddz().
         float  z             = z0 + halfDZ;
         float3 voxelCenterWS = rayOriginWS + z * rayUnDirWS; // Works due to the length of of the dir
 
-        // Dimensions of the voxel as we step along the ray.
-        float3 voxelRightSize = z * voxelAxisRight;
-        float3 voxelUpSize = z * voxelAxisUp;
-        float3 voxelDepthSize = halfDZ * voxelAxisForward;
-
-
         // TODO: define a function ComputeGlobalFogCoefficients(float3 voxelCenterWS),
         // which allows procedural definition of extinction and scattering.
         float3 voxelScattering = _GlobalScattering;
             {
     #endif // USE_CLUSTERED_LIGHTLIST
 
-                OrientedBBox obb = _VolumeBounds[volumeIndex];
+                const OrientedBBox obb = _VolumeBounds[volumeIndex];
-                float3x3 obbFrame   = float3x3(obb.right, obb.up, cross(obb.up, obb.right));
-                float3   obbExtents = float3(obb.extentX, obb.extentY, obb.extentZ);
+                const float3x3 obbFrame   = float3x3(obb.right, obb.up, cross(obb.up, obb.right));
+                const float3   obbExtents = float3(obb.extentX, obb.extentY, obb.extentZ);
-                float3 voxelCenterBS = mul(voxelCenterWS - obb.center, transpose(obbFrame));
-                float3 voxelCenterUV = (voxelCenterBS / obbExtents);
+                const float3 voxelCenterBS = mul(voxelCenterWS - obb.center, transpose(obbFrame));
+                const float3 voxelCenterUV = (voxelCenterBS / obbExtents);
+
+                const float3 voxelAxisRightBS   = mul(voxelAxisRight,   transpose(obbFrame));
+                const float3 voxelAxisUpBS      = mul(voxelAxisUp,      transpose(obbFrame));
+                const float3 voxelAxisForwardBS = mul(voxelAxisForward, transpose(obbFrame));
 
             #if SOFT_VOXELIZATION
                 // We need to determine which is the face closest to 'voxelCenterBS'.
                 // We have determined the normal of the closest face.
                 // We now have to construct the diagonal of the voxel with the longest extent along this normal.
                 float3 minDiagPointBS, maxDiagPointBS;
-
-                float3 voxelAxisRightBS   = mul(voxelAxisRight,   transpose(obbFrame));
-                float3 voxelAxisUpBS      = mul(voxelAxisUp,      transpose(obbFrame));
-                float3 voxelAxisForwardBS = mul(voxelAxisForward, transpose(obbFrame));
 
                 // Start at the center of the voxel.
                 minDiagPointBS = maxDiagPointBS = voxelCenterBS;
 
                 if (overlapFraction > 0)
                 {
-                    float scatteringAndExtinctionMask = 1.0f;
+                    float densityMask = 1.0f;
-                        float3 voxelRightSizeBS = mul(voxelRightSize, transpose(obbFrame));
-                        float3 voxelRightSizeUV = (voxelRightSizeBS / obbExtents);
-                        
-                        float3 voxelUpSizeBS = mul(voxelUpSize, transpose(obbFrame));
-                        float3 voxelUpSizeUV = (voxelUpSizeBS / obbExtents);
-                        
-                        float3 voxelDepthSizeBS = mul(voxelDepthSize, transpose(obbFrame));
-                        float3 voxelDepthSizeUV = (voxelDepthSizeBS / obbExtents);
+                        // We divide extents (half-sizes) by extents here, obtaining full-sized gradients.
+                        float3 voxelGradRightUV   = z      * voxelAxisRightBS   / obbExtents;
+                        float3 voxelGradUpUV      = z      * voxelAxisUpBS      / obbExtents;
+                        float3 voxelGradForwardUV = halfDZ * voxelAxisForwardBS / obbExtents;
-                        scatteringAndExtinctionMask = SampleVolumeMask(_VolumeData[volumeIndex], voxelCenterUV * 0.5 + 0.5, voxelRightSizeUV, voxelUpSizeUV, voxelDepthSizeUV);
+                        densityMask = SampleVolumeMask(_VolumeData[volumeIndex], voxelCenterUV * 0.5 + 0.5, voxelGradRightUV, voxelGradUpUV, voxelGradForwardUV);
-                    voxelScattering += overlapFraction * _VolumeData[volumeIndex].scattering * scatteringAndExtinctionMask;
-                    voxelExtinction += overlapFraction * _VolumeData[volumeIndex].extinction * scatteringAndExtinctionMask;
+                    voxelScattering += overlapFraction * _VolumeData[volumeIndex].scattering * densityMask;
+                    voxelExtinction += overlapFraction * _VolumeData[volumeIndex].extinction * densityMask;
                 }
 
     #ifndef USE_CLUSTERED_LIGHTLIST
     float3 upDirWS     = mul(-float3(upCoord,     1), (float3x3)_VBufferCoordToViewDirWS);
 
     // Compute the axes of the voxel. These are not normalized, but rather computed to scale with Z.
+    // This coordinate system is generally not orthogonal.
     float3 voxelAxisForward = centerDirWS;
     float3 voxelAxisUp      = 0.5 * (upDirWS - centerDirWS);
     float3 voxelAxisRight   = 0.5 * (centerDirWS - leftDirWS);