Smoothing Depth blend in ray marching algorithm

Hi everyone! I’m trying to implement a raymarching algorithm in a post process materials.
I managed to get mostly what I wanted thanks to

Ryan Brucks - Creating a Volumetric Ray Marching

and

ShaderToy

Though I’m struggling on one of the last steps - Compositing final result with scene Depth.
I managed to get the ray killed at depth intersections by calculating max Steps relative to max ray distance (and also using a raytraced Axis Aligned Bouding Box for performance reasons) but i’m getting some banding artifacts I can’t get rid off.

Ryan Brucks mentionned this in his blog post but I didn’t managed to implement this correctly yet…

To fix those slicing artifacts requires just taking one additional step. We track how many steps would have fit up to the scene depth and then take one final step sized to fit the remainder. That assures we end up taking a final sample right at the depth location which smooths out those seams. In order to keep the main tracing loop as simple as possible, we do this outside of the main loop as an additional density/shadow pass.

Any help on this would be really welcomed, thanks!

#pragma once
#define noise_use_smoothstep

struct SDF_CommonLib
{
    float mod(float x, float y)
    {
        return x - y * floor(x / y); 
    }
    
    float2 mod(float2 x, float y)
    {
        return x - y * floor(x / y); 
    }
    
    float3 mod(float3 x, float y)
    {        
        return x - y * floor(x / y);         
    }
    
    float4 mod(float4 x, float y)
    {
        return x - y * floor(x / y); 
    }
    
    
    float hash(float x)
    {
    	return frac(sin(x*.0127863)*17143.321); //decent hash for noise generation
    }
    
    float hash(float2 x)
    {
    	return frac(cos(dot(x.xy,float2(2.31,53.21))*124.123)*412.0); 
    }
    
    float hashmix(float x0, float x1, float interp)
    {
    	x0 = hash(x0);
    	x1 = hash(x1);
    	#ifdef noise_use_smoothstep
    	    interp = smoothstep(0.0,1.0,interp);
    	#endif
    	return lerp(x0,x1,interp);
    }
    
    float hashmix(float2 p0, float2 p1, float2 interp)
    {
    	float v0 = hashmix(p0[0]+p0[1]*128.0,p1[0]+p0[1]*128.0,interp[0]);
    	float v1 = hashmix(p0[0]+p1[1]*128.0,p1[0]+p1[1]*128.0,interp[0]);
    	#ifdef noise_use_smoothstep
    	    interp = smoothstep(float2(0.0, 0.0),float2(1.0, 1.0),interp);
    	#endif
    	return lerp(v0,v1,interp[1]);
    }
    
    float hashmix(float3 p0, float3 p1, float3 interp)
    {
    	float v0 = hashmix(p0.xy+float2(p0.z*143.0,0.0),p1.xy+float2(p0.z*143.0,0.0),interp.xy);
    	float v1 = hashmix(p0.xy+float2(p1.z*143.0,0.0),p1.xy+float2(p1.z*143.0,0.0),interp.xy);
    	#ifdef noise_use_smoothstep
    	    interp = smoothstep(float3(0.0, 0.0, 0.0),float3(1.0, 1.0, 1.0),interp);
    	#endif
    	return lerp(v0,v1,interp[2]);
    }
    
    float hashmix(float4 p0, float4 p1, float4 interp)
    {
        float v0 = hashmix(p0.xyz + float3(p0.w * 17.0, 0.0, 0.0), p1.xyz + float3(p0.w * 17.0, 0.0, 0.0), interp.xyz);
        float v1 = hashmix(p0.xyz + float3(p1.w * 17.0, 0.0, 0.0), p1.xyz + float3(p1.w * 17.0, 0.0, 0.0), interp.xyz);
        #ifdef noise_use_smoothstep
            interp = smoothstep(float4(0.0, 0.0, 0.0, 0.0), float4(1.0, 1.0, 1.0, 1.0), interp);
    	#endif
    	return lerp(v0,v1,interp[3]);
    }
    
    float noise(float p) // 1D noise
    {
    	float pm = mod(p,1.0);
    	float pd = p-pm;
    	return hashmix(pd,pd+1.0,pm);
    }
    
    float noise(float2 p) // 2D noise
    {
    	float2 pm = mod(p,1.0);
    	float2 pd = p-pm;
    	return hashmix(pd,(pd+float2(1.0,1.0)), pm);
    }
    
    float noise(float3 p) // 3D noise
    {
        float3 pm = mod(p, 1.0);
        float3 pd = p - pm;
    	return hashmix(pd,(pd+float3(1.0,1.0,1.0)), pm);
    }
    
    float noise(float4 p) // 4D noise
    {
        float4 pm = mod(p, 1.0);
        float4 pd = p - pm;
        return hashmix(pd, (pd + float4(1.0, 1.0, 1.0, 1.0)), pm);
    }
    
    void RayTracedAABB(float3 rayOrigin, float3 rayDir, float3 boxMin, float3 boxMax, //Raytraced Axis aligned bouding box
                            inout float tNear, inout float tFar, 
                            inout float3 inPos, inout float3 outPos,
                            inout float hitTest, inout float insideTest)
    {    
        float3 tMin = (boxMin - rayOrigin) / rayDir;
        float3 tMax = (boxMax - rayOrigin) / rayDir;
        float3 t1 = min(tMin, tMax);
        float3 t2 = max(tMin, tMax);
        
        tNear = max(max(t1.x, t1.y), t1.z);
        tFar = min(min(t2.x, t2.y), t2.z);
        
        inPos = normalize(rayDir) * tNear + rayOrigin;
        outPos = normalize(rayDir) * tFar + rayOrigin;
        
        hitTest = step(tNear, tFar);
        
        float3 s = step(boxMin, rayOrigin) - step(boxMax, rayOrigin);        
        insideTest = s.x * s.y * s.z;
    }
};


struct FogRayMarcher
{  
    float3 _CameraPos;
    float3 _RayDir;
    float3 _Steps;
    float _VolumeRadius;
    float3 _VolumeCenter;
    float _SceneDepth;
    float _Density;
    float2 _UV;
    float _Time;
        
    SDF_CommonLib lib;
            
    float fogDensity(float3 p, float time, float3 sphereCenter, float sphereRadius) //density function for the cloud
    {
        float4 xp = float4(p * 0.4, time * 0.1 + lib.noise(p));
        
        float sphereMask = sphereRadius - distance(p, sphereCenter);
        float nv = pow(pow(lib.noise(xp), 2.0) * 2.1, 2.5) * sphereMask;
                        
        nv = max(0.0, nv); //negative density is illegal.
        nv = min(1.0, nv); //high density causes artifacts
        return nv;
    }
    
 
    float3 Main()
    {                           
        //Setting up our rays to propagate only inside an Axis-Aligned Bounding Box for performance reasons
        float tNear, tFar; 
        float3 inPos, outPos;
        float hitTest, insideTest;
        
        lib.RayTracedAABB(_CameraPos, _RayDir, -_VolumeRadius + _VolumeCenter, _VolumeRadius + _VolumeCenter, tNear, tFar, inPos, outPos, hitTest, insideTest);
           
        if(hitTest == 0) // Ray is not passing through bounding box so return
            return 0;
        
        float3 rayOrigin = lerp(inPos, _CameraPos, insideTest); // if camera is inside bounding box, ray origin becomes camera Pos.

     
        //Calculating maxSteps : Perf + depth test        
        float stepSize = (_VolumeRadius * 2) / (float) _Steps;
        float3 sceneWP = _SceneDepth * _RayDir + _CameraPos;
        
        float maxDist = min(distance(rayOrigin, outPos), distance(rayOrigin, sceneWP));

        int maxSteps = floor(maxDist / stepSize);
        maxSteps = min(maxSteps, _Steps); // clamping to avoid any miscalculation error
        
        
        
        float3 p = rayOrigin;
        float3 d = normalize(_RayDir);
                
        float dense_acc = 0.0;
        float density_coef = _Density * 0.1;
        //float density_coef = density * 0.35 + 0.23;
                       
        for (int i = 0; i < maxSteps; i++) //max 20 step raycast
        {         		
            float nv = fogDensity(p + d * lib.hash(_UV + float2(_Time, dense_acc)) * 0.25, _Time, _VolumeCenter, _VolumeRadius) * density_coef * stepSize;
		
            dense_acc += nv;
		
            if (dense_acc > 1.0)
            {
                dense_acc = 1.0; //break condition: following steps do not contribute 
                break; //to the color because it's occluded by the gas
            }
            
            p += d * stepSize;
        }
        
        p = maxDist * d + rayOrigin; // Trying to reproduce Ryan Brucks' last sample at maxDist for depth smoothing

        float nv = fogDensity(p + d * lib.hash(_UV + float2(_Time, dense_acc)) * 0.25, _Time, _VolumeCenter, _VolumeRadius) * density_coef * stepSize;
        
        dense_acc += nv;
                                
        float3 sky_color = float3(1, 1, 1);
        float3 dense_color = float3(1, 1, 1);
        float3 color = lerp(float3(0, 0, 0), dense_color, smoothstep(0.0, 1.0, dense_acc));
        color -= length(_UV * 0.2);
        color += lib.hash(color.xy + _UV) * 0.01;
            
        return color;
    } 
};