Volumetric Clouds - Part 2 - Writing an Opacity Raymarcher

Quick update: A couple of days have passed since I wrote part one of the Volumetric Clouds and a few things have changed. We have decided to try and keep the clouds a bit more subtle to not distract as much from the buildings, so we will probably keep them low contrast and experiment a bit with the settings.

Clouds - Work in progress

Definitely needs tweaking

After figuring out the basics I had to dive into writing the Raymarcher, something I hadn't done before, in the Custom Node of Unreal's Material Editor. As mentioned before that meant learning HLSL from scratch. Fortunately it is very similar to C++ and the shader didn't require anything complex, so I was fine. Otherwise the time would have been way too short. Fortunately there are some great resources out there explaining the theory of raymarching. Here are some of the ones I used next to Guerilla Games Paper:

https://shaderbits.com/blog/creating-volumetric-ray-marcher (Ryan Brucks - Creating a Volumetric Ray Marcher)
https://computergraphics.stackexchange.com/questions/161/what-is-ray-marching-is-sphere-tracing-the-same-thing/163
https://www.youtube.com/watch?v=PGtv-dBi2wE (Art of Code - Raymarching for Dummies)
https://www.youtube.com/watch?v=Ff0jJyyiVyw (Art of Code - Raymarching simple Shapes)
https://www.youtube.com/watch?v=Cp5WWtMoeKg (Sebastian Lague - Coding Adventure - Raymarching)
http://www.diva-portal.org/smash/get/diva2:1223894/FULLTEXT01.pdf ( Fredrik Häggström - Real-Time Rendering of Volumetric Clouds)
http://jamie-wong.com/2016/07/15/ray-marching-signed-distance-functions/ (Jamie Wong - Raymarching and Signed Distance Functions)

Here is the full code I am going to break down, it is still subject to changes and improvements, and the Material setup in Unreal.

	//CONSTANTS
	static const float PI = 3.14159265f;
	//VARIABLES
	float sceneDepth = CalcSceneDepth(ScreenAlignedPosition(GetScreenPosition(Parameters)));
	float time = View.GameTime;
	//Camera
	float3 cameraVec = - normalize(Parameters.CameraVector);
	float3 cameraPos = ResolvedView.WorldCameraOrigin;
	//RayMarcher Settings
	int raySteps = 30;
	int lightSteps = 4;
	lightVec = -normalize(lightVec);
	float3 worldPosition = Parameters.AbsoluteWorldPosition;
	float3 sunPos = -lightVec * 10000000;
	float3 customLightVec = normalize(sunPos - worldPosition);
	//BlueNoise
	float2 blueNoiseScale = 5;
	float blueNoiseMult = .2;
	float2 viewportUVs = GetViewportUV(Parameters) * blueNoiseScale;
	float bNoise = (ProcessMaterialColorTextureLookup(Texture2DSample(noiseTex,GetMaterialSharedSampler(noiseTexSampler,View.MaterialTextureBilinearWrapedSampler), viewportUVs)).x - .5) * 2;
	//Distortion
	float2 distortionUV0 = (worldPosition * 0.001 * distortionScale).xy + time * 0.1 * distortionSpeed * -windDirection.xy;
	float2 distortionUV1 = (worldPosition * 0.001 * distortionScale).xy + time * 0.1 * distortionSpeed * windDirection.xy;
	float2 distortionUV2 = (worldPosition * 0.001 * distortionScale).xy + time * 0.1 * distortionSpeed * float2(-windDirection.y, windDirection.x);
	float2 distortionUV3 = (worldPosition * 0.001 * distortionScale).xy + time * 0.1 * distortionSpeed * float2(windDirection.y, -windDirection.x);
	distortionStrength *= 0.01;
	float2 distortion0 = (ProcessMaterialColorTextureLookup(Texture2DSample(DistortionMap,GetMaterialSharedSampler(noiseTexSampler,View.MaterialTextureBilinearWrapedSampler), distortionUV0)).rg - .25) * 2 * distortionStrength;
	float2 distortion1 = (ProcessMaterialColorTextureLookup(Texture2DSample(DistortionMap,GetMaterialSharedSampler(noiseTexSampler,View.MaterialTextureBilinearWrapedSampler), distortionUV1)).rg - .25) * 2 * distortionStrength;
	float2 distortion2 = (ProcessMaterialColorTextureLookup(Texture2DSample(DistortionMap,GetMaterialSharedSampler(noiseTexSampler,View.MaterialTextureBilinearWrapedSampler), distortionUV2)).rg - .25) * 2 * distortionStrength;
	float2 distortion3 = (ProcessMaterialColorTextureLookup(Texture2DSample(DistortionMap,GetMaterialSharedSampler(noiseTexSampler,View.MaterialTextureBilinearWrapedSampler), distortionUV3)).rg - .25) * 2 * distortionStrength;
	float3 distortion = float3((distortion0 + distortion1 + distortion2 + distortion3)/4, 0);
	//Wind
	windSpeed *= .01;
	//FUNCTIONS
	struct function {
	//Sample Density
	float4 sampleDensity(Texture3D Tex, float3 position, float3 cloudScale, float3 cloudOffset, float4 channel_scales){
	float3 uvw0 = position * cloudScale * channel_scales.r * .0001 + cloudOffset + distortion + windDirection * windSpeed.r * time;
	float3 uvw1 = position * cloudScale * channel_scales.g * .0001 + cloudOffset + distortion + windDirection * windSpeed.g * time;
	float3 uvw2 = position * cloudScale * channel_scales.b * .0001 + cloudOffset + distortion + windDirection * windSpeed.b * time;
	float3 uvw3 = position * cloudScale * channel_scales.a * .0001 + cloudOffset + distortion + windDirection * windSpeed.a * time;
	float sample0 = ProcessMaterialColorTextureLookup(Texture3DSample(Tex,GetMaterialSharedSampler(ShapeTexSampler,Material.Wrap_WorldGroupSettings),uvw0)).r;
	float sample1 = ProcessMaterialColorTextureLookup(Texture3DSample(Tex,GetMaterialSharedSampler(ShapeTexSampler,Material.Wrap_WorldGroupSettings),uvw1)).g;
	float sample2 = ProcessMaterialColorTextureLookup(Texture3DSample(Tex,GetMaterialSharedSampler(ShapeTexSampler,Material.Wrap_WorldGroupSettings),uvw2)).b;
	float sample3 = ProcessMaterialColorTextureLookup(Texture3DSample(Tex,GetMaterialSharedSampler(ShapeTexSampler,Material.Wrap_WorldGroupSettings),uvw3)).a;
	float4 sample = float4(sample0, sample1, sample2, sample3);
	float4 density = max(sample - densityThreshold, 0) * densityMultiplier;
	return density;
	}
	float combineDensity(float4 sample, float4 weights){
	weights = max(0, weights);
	float combinedWeights = weights.r + weights.g + weights.b + weights.a;
	float combinedDensities = sample.r * weights.r + sample.g * weights.g + sample.b * weights.b + sample.a * weights.a;
	combinedDensities /= combinedWeights;
	return combinedDensities;
	}
	//Raymarching Box
	float4 boxRayInfo(float3 rayOrigin, float3 rayVec){
	//Bounding Box
	float3 boundMin = TransformLocalPositionToWorld(Parameters, MaterialFloat4(Primitive.LocalObjectBoundsMin,1.00000000).xyz).xyz;
	float3 boundMax = TransformLocalPositionToWorld(Parameters, MaterialFloat4(Primitive.LocalObjectBoundsMax,1.00000000).xyz).xyz;
	rayVec = normalize(rayVec);
	float3 dist1 = (boundMin - rayOrigin)/rayVec;
	float3 dist2 = (boundMax - rayOrigin)/rayVec;
	float3 closeDist = min(dist1, dist2);
	float3 farDist = min(max(dist1, dist2), sceneDepth); //don't draw parts of volume, that's occluded by something
	//Check if ray intersects box
	float distA = max(max(closeDist.x, closeDist.y), closeDist.z); //If distA < 0, camera is in Volume
	float distB = min( min(farDist.x, farDist.y), farDist.z); //If distB <= distA, no hit
	float distToBox = max(0, distA);
	float distInsideBox = max(0, distB - distToBox);
	float3 entrypos = rayOrigin + (distToBox * rayVec);
	return float4(entrypos, distInsideBox);
	}
	//Henyey - Greenstein Phase Function
	float HenyeyGreenstein(float g, float angle){
	float hg = (1 - pow(g,2)) / pow((1 + pow(g,2) - 2*g * degrees(cos(angle))), 1.5 ) / (4 * PI);
	return hg;
	}
	//Extra scatter based on sun position
	float ISExtra(float angle, float csi, float cse){
	float inScatter = csi * pow(saturate(angle), cse);
	return inScatter;
	}
	//Combined HenyeyGreenstein and ISExtra
	float scatter(float ivo, float ins, float outs, float3 cameraVec, float3 lightVec){
	float angle = dot(cameraVec, lightVec)/(length(cameraVec)*length(lightVec));
	float inScatter = max(HenyeyGreenstein(ins, angle), ISExtra(angle, .3, 1));
	float outScatter = HenyeyGreenstein(outs, angle);
	float ios = lerp(inScatter, outScatter, ivo);
	return ios;
	}
	};
	function f;
	//MAIN
	float distInsideBox = f.boxRayInfo(cameraPos, cameraVec).w;
	float3 entrypos = f.boxRayInfo(cameraPos, cameraVec).xyz;
	entrypos += bNoise * cameraVec * blueNoiseMult;
	float distLight = f.boxRayInfo(entrypos, customLightVec).w;
	float stepSize = distInsideBox/raySteps;
	float3 rayPos = entrypos;
	float3 lpos = entrypos;
	float lightStepSize = distLight/raySteps;
	float accumDens = 0;
	float accumShadow = 0;
	customLightVec = normalize(customLightVec);
	//RayMarch
	for(int i = 0; i < raySteps; i++) {
	float4 sample = saturate(f.sampleDensity(ShapeTex, rayPos, cloud_Scale, cloud_Offset, shapeScales));
	float4 lsample = saturate(f.sampleDensity(ShapeTex, lpos, cloud_Scale, cloud_Offset, shapeScales));
	float combSample = f.combineDensity(sample, shapeWeights);
	float combLSample = f.combineDensity(lsample, shapeWeights);
	accumDens += combSample * stepSize / 1000;
	accumShadow += combLSample * lightStepSize / 1000 ;
	rayPos += cameraVec * stepSize;
	lpos += customLightVec * stepSize;
	}
	float scatter = f.scatter(IVO, inScatter, 0, normalize(cameraVec), normalize(customLightVec));
	return float4(accumDens, accumShadow, (scatter * accumShadow * accumDens)/100, exp(-accumShadow));

view raw RM_Clouds.hlsl hosted with ❤ by GitHub

M_Clouds Overview

 

1 - Writing an Opacity Raymarcher

To sample the opacity of the volume textures a vector/ray between the camera and volume is used to define the marching direction and then the texture gets sampled whilst moving along that ray in increments. At the end the the samples get combined to get the accumulated density as seen from the camera's view.

Raymarching Visualisation

The code to do this looks something like this:

	float3 cameraVec = -normalize(Parameters.CameraVector);
	float3 cameraPos = ResolvedView.WorldCameraOrigin;
	int raySteps = 30;
	float marchingDist = 10;
	float stepSize = marchingDist/raySteps;
	float totalDensity = 0;
	for(int i = 0; i < raySteps, i++){
	float3 pos = cameraPos + cameraVec * stepSize * i;
	float sample = Texture3DSample(Tex,GetMaterialSharedSampler(ShapeTexSampler,Material.Wrap_WorldGroupSettings),pos)).r;
	totalDensity += sample / stepSize;
	}
	return totalDensity;

view raw RM_Basic_Algorithm.hlsl hosted with ❤ by GitHub

However to actually make sure, that the Raymarcher samples the points within the box containing the volume, we need to figure out  where it needs to start sampling and how far it needs to move within the box. A lot of the resources I found make use of the Ray-Box intersection Algorithm. To keep it simple I stuck with the Bounding Box being aligned to the world axis.

The following resources were really useful for understanding how the algorithm works:
https://www.youtube.com/watch?v=4h-jlOBsndU (SketchpunkLab -WebGL 2.0:046:Ray intersects Bounding Box (AABB))
https://www.scratchapixel.com/lessons/3d-basic-rendering/minimal-ray-tracer-rendering-simple-shapes/ray-box-intersection (Scratchapixel - Ray-Box Intersection)
https://developer.arm.com/docs/100140/0302/advanced-graphics-techniques/implementing-reflections-with-a-local-cubemap/ray-box-intersection-algorithm (ARM Developer - Ray-Box Intersection Algorithm)
https://medium.com/@bromanz/another-view-on-the-classic-ray-aabb-intersection-algorithm-for-bvh-traversal-41125138b525 (Roman Wiche - Another View on the Classic Ray-AABB Intersection Algorithm for BVH Traversal)

These all so a far better job at explaining how it works, than I could ever do, but I will try to give a brief overview. I have also annotated the code to hopefully make a few things a bit clearer.

So essentially the bounding box is defined by a minimum and maximum (X, Y, Z) position, which become the origin of corresponding new axis, as shown in the image below.

Bounding Box

The ray is defined by the following formula: f(t) = origin + Direction * t. The Intersection algorithm then checks where that ray intersects with the bounding axis and if the intersection is actually on the surface of the bounding box.
I further check, if part of the bounding box is occluded by something else by comparing Scene Depth to the further away intersection. Otherwise the volume is drawn on top of all other objects, without regarding depth. Here is a comparison between the Raymarcher with and without Depth Check.

Depth Check OFF vs Depth Check ON

Here is the HLSL Ray-Box Intersection Code:

	//Axis Aligned Bounding Box AABB Intersection
	float4 boxRayInfo(float3 rayOrigin, float3 rayVec){
	//Bounding Box
	float3 boundMin = TransformLocalPositionToWorld(Parameters, MaterialFloat4(Primitive.LocalObjectBoundsMin,1.00000000).xyz).xyz;
	float3 boundMax = TransformLocalPositionToWorld(Parameters, MaterialFloat4(Primitive.LocalObjectBoundsMax,1.00000000).xyz).xyz;
	rayVec = normalize(rayVec);
	//This determines where the ray intersects the bounds
	float3 tMin = (boundMin - rayOrigin)/rayVec;
	float3 tMax = (boundMax - rayOrigin)/rayVec;
	tMin = min(tMin, tMax);
	tMax = min(max(tMin, tMax), sceneDepth); //don't draw parts of volume, that's occluded by something
	//Check if ray intersects box
	float distA = max(max(tMin.x, tMin.y), tMin.z); //If distA < 0, camera is in Volume
	float distB = min( min(tMax.x, tMax.y), tMax.z); //If distB <= distA, no hit
	float distToBox = max(0, distA);
	float distInsideBox = max(0, distB - distToBox);
	float3 entrypos = rayOrigin + (distToBox * rayVec);
	return float4(entrypos, distInsideBox);
	}

view raw Ray_Box_Intersection.hlsl hosted with ❤ by GitHub

The Algorithm returns both the point to start the sampling from, as well as the distance to travel within the box. I can now use the later to calculate the step size. I have visualised that step length throughout the box, which gets used by the Raymarcher, being shorter towards the thinner sections from the camera's view point.

I can now update the Raymarching code to this:

	float3 cameraVec = -normalize(Parameters.CameraVector);
	float3 cameraPos = ResolvedView.WorldCameraOrigin;
	struct function {
	//Axis Aligned Bounding Box AABB Intersection
	float4 boxRayInfo(float3 rayOrigin, float3 rayVec){
	//Bounding Box
	float3 boundMin = TransformLocalPositionToWorld(Parameters, MaterialFloat4(Primitive.LocalObjectBoundsMin,1.00000000).xyz).xyz;
	float3 boundMax = TransformLocalPositionToWorld(Parameters, MaterialFloat4(Primitive.LocalObjectBoundsMax,1.00000000).xyz).xyz;
	rayVec = normalize(rayVec);
	//This determines where the ray intersects the bounds
	float3 tMin = (boundMin - rayOrigin)/rayVec;
	float3 tMax = (boundMax - rayOrigin)/rayVec;
	tMin = min(tMin, tMax);
	tMax = min(max(tMin, tMax), sceneDepth); //don't draw parts of volume, that's occluded by something
	//Check if ray intersects box
	float distA = max(max(tMin.x, tMin.y), tMin.z); //If distA < 0, camera is in Volume
	float distB = min( min(tMax.x, tMax.y), tMax.z); //If distB <= distA, no hit
	float distToBox = max(0, distA);
	float distInsideBox = max(0, distB - distToBox);
	float3 entrypos = rayOrigin + (distToBox * rayVec);
	return float4(entrypos, distInsideBox);
	}
	};
	function f;
	//MAIN
	int raySteps = 64;
	float distInsideBox = f.boxRayInfo(cameraPos, cameraVec).w;
	float3 entrypos = f.boxRayInfo(cameraPos, cameraVec).xyz;
	float stepSize = distInsideBox/raySteps;
	float3 rayPos = entrypos;
	float totalDensity = 0;
	for(int i = 0; i < raySteps; i++){
	float sample = Texture3DSample(Tex,GetMaterialSharedSampler(TexSampler,Material.Wrap_WorldGroupSettings),rayPos).r;
	totalDensity += sample / stepSize;
	rayPos += cameraVec * stepSize;
	}
	return float4(totalDensity, 0, 0, 0);

view raw RM_AABB.hlsl hosted with ❤ by GitHub

Unreal's Custom Node doesn't allow having functions, so the functions need to be put into a struct, that way they can be called.