Hi all,

Great work so far in this thread.  I've been implementing this technique over the past little while as I've had time, and thought I would dump what I've got so far in case it helps move the discussion along.  A few caveats:

1. This is my first big foray into volume rendering, so I'm all but positive that some of the math is off.
2. The atmospheric scattering below the clouds is completely fake and just trial and error values (you'll see in the shader what I'm doing).
3. Ambient light is also just trial and error values, nothing fancy.
4. My noise doesn't tile well - I get around this by using a mirrored sampler, but the downside is that the repeating pattern becomes clear pretty quickly if you looked at from the right distance.
5. No temporal reprojection (yet), so performance is about what you'd expect based on the articles.

I think items 2 and 3 can be somewhat addressed by looking into the way Uncharted did some of their scattering (mip fog - basically mip map the sky box and use that as a fog lookup based on distance).  FreneticPony mentioned this in Post 40.  I'm hoping some of you can help me with item 1 through code review.  Suggestions for item 4 are also welcome, of course.

Here's what's cooked up so far:

#include "../../ConstantBuffers/PerFrame.hlsli"
#include "../../Utils/DepthUtils.hlsli"
#include "../../Constants.hlsli"

Texture3D baseShapeLookup : register(t0);
Texture3D erosionLookup : register(t1);
Texture2D weatherLookup : register(t2);
Texture2D depthTexture : register(t3);

SamplerState texSampler : register(s0);

cbuffer cbVolumetricClouds : register(b0)
{
	float3 cb_lightDirection; // light direction world space
	float cb_groundRadius; // meters - for a ground/lower atmosphere only version, this could be smaller
			 
	float3 cb_sunColor; // color of sun light
	uint cb_baseShapeTextureBottomMipLevel;

	float cb_cloudVolumeStartHeight; // meters - height above ground level
	float cb_cloudVolumeHeight; // meters - height above ground level
	float cb_cloudSpeed;
	float cb_cloudTopOffset;

	float3 cb_windDirection;
	uint cb_erosionTextureBottomMipLevel;

	float3 cb_weatherTexMod; // scale(x), offset(y, z)
	float cb_windStrength;
};

static const float VOLUME_END_HEIGHT = cb_cloudVolumeStartHeight + cb_cloudVolumeHeight;
// planet center (world space)
static const float3 PLANET_CENTER = float3(0.0f, -cb_groundRadius - 100.0f, 0.0f); // TODO revisit - 100.0f offset is to match planet sky settings
// radius from the planet center to the bottom of the cloud volume
static const float PLANET_CENTER_TO_LOWER_CLOUD_RADIUS = cb_groundRadius + cb_cloudVolumeStartHeight;
// radius from the planet center to the top of the cloud volume
static const float PLANET_CENTER_TO_UPPER_CLOUD_RADIUS = cb_groundRadius + VOLUME_END_HEIGHT;
static const float CLOUD_SCALE = 1.0f / VOLUME_END_HEIGHT;
static const float3 WEATHER_TEX_MOD = float3(1.0f / (VOLUME_END_HEIGHT * cb_weatherTexMod.x), cb_weatherTexMod.y, cb_weatherTexMod.z); 
static const float2 WEATHER_TEX_MOVE_SPEED = float2(cb_windStrength * cb_windDirection.x, cb_windStrength * cb_windDirection.z); // this is modded by app run time
// samples based on shell thickness between inner and outer volume
static const uint2 SAMPLE_RANGE = uint2(64u, 128u);

static const float4 STRATUS_GRADIENT = float4(0.02f, 0.05f, 0.09f, 0.11f);
static const float4 STRATOCUMULUS_GRADIENT = float4(0.02f, 0.2f, 0.48f, 0.625f);
static const float4 CUMULUS_GRADIENT = float4(0.01f, 0.0625f, 0.78f, 1.0f); // these fractions would need to be altered if cumulonimbus are added to the same pass

/**
 * Perform a ray-sphere intersection test.
 * Returns the number of intersections in the direction of the ray (excludes intersections behind the ray origin), between 0 and 2.
 * In the case of more than one intersection, the nearest point will be returned in t1.
 * 
 * http://www.scratchapixel.com/lessons/3d-basic-rendering/minimal-ray-tracer-rendering-simple-shapes/ray-sphere-intersection
 */
uint intersectRaySphere(
	float3 rayOrigin,
	float3 rayDir, // must be normalized
	float3 sphereCenter,
	float  sphereRadius,
	out float2 t)
{
	float3 l = rayOrigin - sphereCenter;
	float a = 1.0f; // dot(rayDir, rayDir) where rayDir is normalized
	float b = 2.0f * dot(rayDir, l);
	float c = dot(l, l) - sphereRadius * sphereRadius;
	float discriminate = b * b - 4.0f * a * c;
	if(discriminate < 0.0f)
	{
		t.x = t.y = 0.0f;
		return 0u;
	}
	else if(abs(discriminate) - 0.00005f <= 0.0f)
	{
		t.x = t.y = -0.5f * b / a;
		return 1u;
	}
	else
	{
		float q = b > 0.0f ?
			-0.5f * (b + sqrt(discriminate)) :
			-0.5f * (b - sqrt(discriminate));
		float h1 = q / a;
		float h2 = c / q;
		t.x = min(h1, h2);
		t.y = max(h1, h2);
		if(t.x < 0.0f)
		{
			t.x = t.y;
			if(t.x < 0.0f)
			{
				return 0u;
			}
			return 1u;
		}
		return 2u;
	}
}

float remap(
	float value,
	float oldMin,
	float oldMax,
	float newMin,
	float newMax)
{
	return newMin + (value - oldMin) / (oldMax - oldMin) * (newMax - newMin);
}

float3 sampleWeather(float3 pos)
{
	return weatherLookup.SampleLevel(texSampler, pos.xz * WEATHER_TEX_MOD.x + WEATHER_TEX_MOD.yz + (WEATHER_TEX_MOVE_SPEED * cb_appRunTime), 0.0f).rgb;
}

float getCoverage(float3 weatherData)
{
	return weatherData.r;
}

float getPrecipitation(float3 weatherData)
{
	return weatherData.g;
}

float getCloudType(float3 weatherData)
{
	// weather b channel tells the cloud type 0.0 = stratus, 0.5 = stratocumulus, 1.0 = cumulus
	return weatherData.b;
}

float heightFraction(float3 pos)
{
	return saturate((distance(pos, PLANET_CENTER) - PLANET_CENTER_TO_LOWER_CLOUD_RADIUS) / cb_cloudVolumeHeight);
}

float4 mixGradients(
	float cloudType)
{
	float stratus = 1.0f - saturate(cloudType * 2.0f);
	float stratocumulus = 1.0f - abs(cloudType - 0.5f) * 2.0f;
	float cumulus = saturate(cloudType - 0.5f) * 2.0f;
	return STRATUS_GRADIENT * stratus + STRATOCUMULUS_GRADIENT * stratocumulus + CUMULUS_GRADIENT * cumulus;
}

float densityHeightGradient(
	float heightFrac,
	float cloudType)
{
	float4 cloudGradient = mixGradients(cloudType);
	return smoothstep(cloudGradient.x, cloudGradient.y, heightFrac) - smoothstep(cloudGradient.z, cloudGradient.w, heightFrac);
}

float sampleCloudDensity(
	float3 pos,
	float3 weatherData,
	float heightFrac,
	float lod)
{
	pos += heightFrac * cb_windDirection * cb_cloudTopOffset;
	pos += (cb_windDirection + float3(0.0f, -0.25f, 0.0f)) * cb_cloudSpeed * (cb_appRunTime/* * 25.0f*/); // the * 25.0f is just for testing to make the effect obvious
	pos *= CLOUD_SCALE;
	
	float4 lowFreqNoise = baseShapeLookup.SampleLevel(texSampler, pos, lerp(0.0f, cb_baseShapeTextureBottomMipLevel, lod));
	float lowFreqFBM =
		(lowFreqNoise.g * 0.625f) +
		(lowFreqNoise.b * 0.25f) +
		(lowFreqNoise.a * 0.125f);

	float baseCloud = remap(
		lowFreqNoise.r,
		-(1.0f - lowFreqFBM), 1.0f, // gets about the same results just using -lowFreqFBM
		0.0f, 1.0f);

	float densityGradient = densityHeightGradient(heightFrac, getCloudType(weatherData));
	baseCloud *= densityGradient;

	float cloudCoverage = getCoverage(weatherData);
	float baseCloudWithCoverage = remap(
		baseCloud,
		1.0f - cloudCoverage, 1.0f,
		0.0f, 1.0f);
	baseCloudWithCoverage *= cloudCoverage;

	//// TODO add curl noise
	//// pos += curlNoise.xy * (1.0f - heightFrac);

	float3 highFreqNoise = erosionLookup.SampleLevel(texSampler, pos * 0.1f, lerp(0.0f, cb_erosionTextureBottomMipLevel, lod)).rgb;
	float highFreqFBM =
		(highFreqNoise.r * 0.625f) +
		(highFreqNoise.g * 0.25f) +
		(highFreqNoise.b * 0.125f);

	float highFreqNoiseModifier = lerp(highFreqFBM, 1.0f - highFreqFBM, saturate(heightFrac * 10.0f));

	baseCloudWithCoverage = remap(
		baseCloudWithCoverage,
		highFreqNoiseModifier * 0.2f, 1.0f,
		0.0f, 1.0f);

	return saturate(baseCloudWithCoverage);
}

struct VertexOut
{
	float4 posH : SV_POSITION;
	float3 viewRay : VIEWRAY;
	float2 tex : TEXCOORD;
};

// random vectors on the unit sphere
static const float3 RANDOM_VECTORS[] =
{
	float3( 0.38051305f,  0.92453449f, -0.02111345f),
	float3(-0.50625799f, -0.03590792f, -0.86163418f),
	float3(-0.32509218f, -0.94557439f,  0.01428793f),
	float3( 0.09026238f, -0.27376545f,  0.95755165f),
	float3( 0.28128598f,  0.42443639f, -0.86065785f),
	float3(-0.16852403f,  0.14748697f,  0.97460106f)
};

static const uint LIGHT_RAY_ITERATIONS = 6u;
static const float RCP_LIGHT_RAY_ITERATIONS = 1.0f / float(LIGHT_RAY_ITERATIONS);

float beerLambert(float sampleDensity, float precipitation)
{
	return exp(-sampleDensity * precipitation);
}

float powder(float sampleDensity, float lightDotEye)
{
	float powd = 1.0f - exp(-sampleDensity * 2.0f);
	return lerp(
		1.0f,
		powd,
		saturate((-lightDotEye * 0.5f) + 0.5f) // [-1,1]->[0,1]
	);
}

float henyeyGreenstein(
	float lightDotEye,
	float g)
{
	float g2 = g * g;
	return ((1.0f - g2) / pow((1.0f + g2 - 2.0f * g * lightDotEye), 1.5f)) * 0.25f;
}

float lightEnergy(
	float lightDotEye,
	float densitySample,
	float originalDensity,
	float precipitation)
{
	return 2.0f *
		beerLambert(densitySample, precipitation) *
		powder(originalDensity, lightDotEye) * 
		lerp(henyeyGreenstein(lightDotEye, 0.8f), henyeyGreenstein(lightDotEye, -0.5f), 0.5f);
}

// TODO get from cb values - has to change as time of day changes
float3 ambientLight(float heightFrac)
{
	return lerp(
		float3(0.5f, 0.67f, 0.82f),
		float3(1.0f, 1.0f, 1.0f),
		heightFrac);
}

float sampleCloudDensityAlongCone(
	float3 startPos,
	float  stepSize,
	float  lightDotEye,
	float  originalDensity)
{
	float3 lightStep = stepSize * -cb_lightDirection;
	float3 pos = startPos;
	float coneRadius = 1.0f;
	float coneStep = RCP_LIGHT_RAY_ITERATIONS;
	float densityAlongCone = 0.0f;
	float lod = 0.0f;
	float lodStride = RCP_LIGHT_RAY_ITERATIONS;
	float3 weatherData = 0.0f;
	float rcpThickness = 1.0f / (stepSize * LIGHT_RAY_ITERATIONS);
	float density = 0.0f;

	for(uint i = 0u; i < LIGHT_RAY_ITERATIONS; ++i)
	{
		float3 conePos = pos + coneRadius * RANDOM_VECTORS[i] * float(i + 1u);
		float heightFrac = heightFraction(conePos);
		if(heightFrac <= 1.0f)
		{
			weatherData = sampleWeather(conePos);
			float cloudDensity = sampleCloudDensity(
				conePos,
				weatherData,
				heightFrac,
				lod);
			if(cloudDensity > 0.0f)
			{
				density += cloudDensity;
				float transmittance = 1.0f - (density * rcpThickness);
				densityAlongCone += (cloudDensity * transmittance);
			}
		}
		pos += lightStep;
		coneRadius += coneStep;
		lod += lodStride;
	}
	// take additional step at large distance away for shadowing from other clouds
	pos = pos + (lightStep * 8.0f);
	weatherData = sampleWeather(pos);
	float heightFrac = heightFraction(pos);
	if(heightFrac <= 1.0f)
	{
		float cloudDensity = sampleCloudDensity(
			pos,
			weatherData,
			heightFrac,
			0.8f);
		// no need to branch here since density variable is no longer used after this
		density += cloudDensity;
		float transmittance = 1.0f - saturate(density * rcpThickness);
		densityAlongCone += (cloudDensity * transmittance);
	}

	return saturate(lightEnergy(
		lightDotEye,
		densityAlongCone,
		originalDensity,
		lerp(1.0f, 2.0f, getPrecipitation(weatherData))));
}

float4 traceClouds(
	float3 viewDirW,		// world space view direction
	float3 startPos,		// world space start position
	float3 endPos)			// world space end position
{
	float3 dir = endPos - startPos;
	float thickness = length(dir);
	float rcpThickness = 1.0f / thickness;
	uint sampleCount = lerp(SAMPLE_RANGE.x, SAMPLE_RANGE.y, saturate((thickness - cb_cloudVolumeHeight) / cb_cloudVolumeHeight));
	float stepSize = thickness / float(sampleCount);
	dir /= thickness;
	float3 posStep = stepSize * dir;

	float lightDotEye = -dot(cb_lightDirection, viewDirW);

	float3 pos = startPos;
	float3 weatherData = 0.0f;
	float4 result = 0.0f;
	float density = 0.0f;
	
	for(uint i = 0u; i < sampleCount; ++i)
	{
		float heightFrac = heightFraction(pos);
		weatherData = sampleWeather(pos);
		float cloudDensity = sampleCloudDensity(
			pos,
			weatherData,
			heightFrac,
			0.0f);

		if(cloudDensity > 0.0f)
		{
			density += cloudDensity;
			float transmittance = 1.0f - (density * rcpThickness);
			float lightDensity = sampleCloudDensityAlongCone(
				pos,
				stepSize,
				lightDotEye,
				cloudDensity);

			float3 ambientBadApprox = ambientLight(heightFrac) * min(1.0f, length(cb_sunColor.rgb * 0.0125f)) * transmittance;
			float4 source = float4((cb_sunColor.rgb * lightDensity) + ambientBadApprox/*+ ambientLight(heightFrac)*/, cloudDensity * transmittance); // TODO enable ambient when added to constant buffer
			source.rgb *= source.a;
			result = (1.0f - result.a) * source + result;
			if(result.a >= 1.0f) break;
		}

		pos += posStep;
	}

	// experimental fog - may not be needed if clouds are drawn before atmosphere - would have to draw sun by itself, then clouds, then atmosphere
	// fogAmt = 0 to disable
	float fogAmt = 1.0f - exp(-distance(startPos, cb_eyePositionW) * 0.00001f);
	float3 fogColor = float3(0.3f, 0.4f, 0.45f) * length(cb_sunColor.rgb * 0.125f) * 0.8f;
	float3 sunColor = normalize(cb_sunColor.rgb) * 4.0f * length(cb_sunColor.rgb * 0.125f);
	fogColor = lerp(fogColor, sunColor, pow(saturate(lightDotEye), 8.0f));
	return float4(clamp(lerp(result.rgb, fogColor, fogAmt), 0.0f, 1000.0f), saturate(result.a));
}

float4 main(VertexOut pIn) : SV_TARGET
{
	int3 loadIndices = int3(pIn.posH.xy, 0);
	float zwDepth = depthTexture.Load(loadIndices).r;
	bool depthPresent = zwDepth < 1.0f;
	float depth = linearizeDepth(zwDepth);
	float3 posV = pIn.viewRay * depth;
	float3 posW = mul(float4(posV, 1.0f), cb_inverseViewMatrix).xyz;
	float3 viewDirW = normalize(posW - cb_eyePositionW);

	// find nearest planet surface point
	float2 ph = 0.0f;
	uint planetHits = intersectRaySphere(
		cb_eyePositionW,
		viewDirW,
		PLANET_CENTER,
		cb_groundRadius,
		ph);

	// find nearest inner shell point
	float2 ih = 0.0f;
	uint innerShellHits = intersectRaySphere(
		cb_eyePositionW,
		viewDirW,
		PLANET_CENTER,
		PLANET_CENTER_TO_LOWER_CLOUD_RADIUS,
		ih);

	// find nearest outer shell point
	float2 oh = 0.0f;
	uint outerShellHits = intersectRaySphere(
		cb_eyePositionW,
		viewDirW,
		PLANET_CENTER,
		PLANET_CENTER_TO_UPPER_CLOUD_RADIUS,
		oh);

	// world space ray intersections
	float3 planetHitSpot = cb_eyePositionW + (viewDirW * ph.x);
	float3 innerShellHit = cb_eyePositionW + (viewDirW * ih.x);
	float3 outerShellHit = cb_eyePositionW + (viewDirW * oh.x);

	// eye radius from planet center
	float eyeRadius = distance(cb_eyePositionW, PLANET_CENTER);
	if(eyeRadius < PLANET_CENTER_TO_LOWER_CLOUD_RADIUS) // under inner shell
	{
		// exit if there's something in front of the start of the cloud volume
		if((depthPresent && (distance(posW, cb_eyePositionW) < distance(innerShellHit, cb_eyePositionW))) ||
			planetHits > 0u) // shell hits are guaranteed, but the ground may be occluding cloud layer
		{
			return float4(0.0f, 0.0f, 0.0f, 0.0f);
		}
		return traceClouds(
			viewDirW,
			innerShellHit,
			outerShellHit);
	}
	else if(eyeRadius > PLANET_CENTER_TO_UPPER_CLOUD_RADIUS) // over outer shell
	{
		// possibilities are
		// 1) enter outer shell, leave inner shell
		// 2) enter outer shell, leave outer shell
		float3 firstShellHit = outerShellHit;
		if(outerShellHits == 0u ||
			depthPresent && (distance(posW, cb_eyePositionW) < distance(firstShellHit, cb_eyePositionW)))
		{
			return float4(0.0f, 0.0f, 0.0f, 0.0f);
		}
		float3 secondShellHit = outerShellHits == 2u && innerShellHits == 0u ? cb_eyePositionW + (viewDirW * oh.y) : innerShellHit;
		return traceClouds(
			viewDirW,
			firstShellHit,
			secondShellHit);
	}
	else // between shells
	{
		/*
		 * From a practical viewpoint (properly scaled planet, atmosphere, etc.)
		 * only one shell will be hit.
		 * Start position is always eye position.
		 */
		float3 shellHit = innerShellHits > 0u ? innerShellHit : outerShellHit;
		// if there's something in the depth buffer that's closer, that's the end point
		if(depthPresent && (distance(posW, cb_eyePositionW) < distance(shellHit, cb_eyePositionW)))
		{
			shellHit = posW;
		}
		return traceClouds(
			viewDirW,
			cb_eyePositionW,
			shellHit);
	}
}

Again, I'm sure some of the math and maybe even the general approach is off a bit, so please feel free to offer suggestions for improvement.  One thing I wish was better is that the clouds don't seem to taper well as they rise with my current implementation.  There may be some ideas from this thread I can use to help with that.  I also need to get curl noise in, but I'm more worried about getting the general solution working well first.

If there's anything I can help out with, let me know

Another interesting part of the system is the weather simulation.
I found that a good way to handle cloud formation/dissipation and movement is using Cellular Automata as Dobashi tried:

http://evasion.imag.fr/~Antoine.Bouthors/research/dea/sig00_cloud.pdf

A very brief summary of the paper is this:

https://graphics.stanford.edu/courses/cs348b-competition/cs348b-05/clouds/index.html

A working implementaton is here:

https://software.intel.com/en-us/articles/dynamic-volumetric-cloud-rendering-for-games-on-multi-core-platforms

All the simulation is a 3D grid, and it can be used to generate the weather map.

Accumulating cloudiness on the y-axis could create the cloudtype parameter, and if over a certain threshold the coverage too.

Probably another pass could see spots areas with taller clouds (maybe doing a downscaled version of the simulation map) to check precipitation values.

Just a train of thoughts!

@SourceDude22 have you finish it?

Hello again,

I'm going to be presenting a 40 Minute Talk at SIGGRAPH 2017 as part of the Advances course again. 

I am going to cover the noise generation and weather system as well as our authoring approach and some other updates. If it is allowed, I will share the Houdini tool that I made to generate our noise with the public. This way you can also experiment with other configurations.

Hope you can make it or review the slides when the come out. 

Andrew

PS this is awesome!!!

On 11/11/2016 at 8:19 PM, WFP said:

Hi all,

 

Great work so far in this thread.  I've been implementing this technique over the past little while as I've had time, and thought I would dump what I've got so far in case it helps move the discussion along.  A few caveats:

1. This is my first big foray into volume rendering, so I'm all but positive that some of the math is off.
2. The atmospheric scattering below the clouds is completely fake and just trial and error values (you'll see in the shader what I'm doing).
3. Ambient light is also just trial and error values, nothing fancy.
4. My noise doesn't tile well - I get around this by using a mirrored sampler, but the downside is that the repeating pattern becomes clear pretty quickly if you looked at from the right distance.
5. No temporal reprojection (yet), so performance is about what you'd expect based on the articles.

I think items 2 and 3 can be somewhat addressed by looking into the way Uncharted did some of their scattering (mip fog - basically mip map the sky box and use that as a fog lookup based on distance).  FreneticPony mentioned this in Post 40.  I'm hoping some of you can help me with item 1 through code review.  Suggestions for item 4 are also welcome, of course.

 

Here's what's cooked up so far:

#include "../../ConstantBuffers/PerFrame.hlsli"
#include "../../Utils/DepthUtils.hlsli"
#include "../../Constants.hlsli"

Texture3D baseShapeLookup : register(t0);
Texture3D erosionLookup : register(t1);
Texture2D weatherLookup : register(t2);
Texture2D depthTexture : register(t3);

SamplerState texSampler : register(s0);

cbuffer cbVolumetricClouds : register(b0)
{
	float3 cb_lightDirection; // light direction world space
	float cb_groundRadius; // meters - for a ground/lower atmosphere only version, this could be smaller
			 
	float3 cb_sunColor; // color of sun light
	uint cb_baseShapeTextureBottomMipLevel;

	float cb_cloudVolumeStartHeight; // meters - height above ground level
	float cb_cloudVolumeHeight; // meters - height above ground level
	float cb_cloudSpeed;
	float cb_cloudTopOffset;

	float3 cb_windDirection;
	uint cb_erosionTextureBottomMipLevel;

	float3 cb_weatherTexMod; // scale(x), offset(y, z)
	float cb_windStrength;
};

static const float VOLUME_END_HEIGHT = cb_cloudVolumeStartHeight + cb_cloudVolumeHeight;
// planet center (world space)
static const float3 PLANET_CENTER = float3(0.0f, -cb_groundRadius - 100.0f, 0.0f); // TODO revisit - 100.0f offset is to match planet sky settings
// radius from the planet center to the bottom of the cloud volume
static const float PLANET_CENTER_TO_LOWER_CLOUD_RADIUS = cb_groundRadius + cb_cloudVolumeStartHeight;
// radius from the planet center to the top of the cloud volume
static const float PLANET_CENTER_TO_UPPER_CLOUD_RADIUS = cb_groundRadius + VOLUME_END_HEIGHT;
static const float CLOUD_SCALE = 1.0f / VOLUME_END_HEIGHT;
static const float3 WEATHER_TEX_MOD = float3(1.0f / (VOLUME_END_HEIGHT * cb_weatherTexMod.x), cb_weatherTexMod.y, cb_weatherTexMod.z); 
static const float2 WEATHER_TEX_MOVE_SPEED = float2(cb_windStrength * cb_windDirection.x, cb_windStrength * cb_windDirection.z); // this is modded by app run time
// samples based on shell thickness between inner and outer volume
static const uint2 SAMPLE_RANGE = uint2(64u, 128u);

static const float4 STRATUS_GRADIENT = float4(0.02f, 0.05f, 0.09f, 0.11f);
static const float4 STRATOCUMULUS_GRADIENT = float4(0.02f, 0.2f, 0.48f, 0.625f);
static const float4 CUMULUS_GRADIENT = float4(0.01f, 0.0625f, 0.78f, 1.0f); // these fractions would need to be altered if cumulonimbus are added to the same pass

/**
 * Perform a ray-sphere intersection test.
 * Returns the number of intersections in the direction of the ray (excludes intersections behind the ray origin), between 0 and 2.
 * In the case of more than one intersection, the nearest point will be returned in t1.
 * 
 * http://www.scratchapixel.com/lessons/3d-basic-rendering/minimal-ray-tracer-rendering-simple-shapes/ray-sphere-intersection
 */
uint intersectRaySphere(
	float3 rayOrigin,
	float3 rayDir, // must be normalized
	float3 sphereCenter,
	float  sphereRadius,
	out float2 t)
{
	float3 l = rayOrigin - sphereCenter;
	float a = 1.0f; // dot(rayDir, rayDir) where rayDir is normalized
	float b = 2.0f * dot(rayDir, l);
	float c = dot(l, l) - sphereRadius * sphereRadius;
	float discriminate = b * b - 4.0f * a * c;
	if(discriminate < 0.0f)
	{
		t.x = t.y = 0.0f;
		return 0u;
	}
	else if(abs(discriminate) - 0.00005f <= 0.0f)
	{
		t.x = t.y = -0.5f * b / a;
		return 1u;
	}
	else
	{
		float q = b > 0.0f ?
			-0.5f * (b + sqrt(discriminate)) :
			-0.5f * (b - sqrt(discriminate));
		float h1 = q / a;
		float h2 = c / q;
		t.x = min(h1, h2);
		t.y = max(h1, h2);
		if(t.x < 0.0f)
		{
			t.x = t.y;
			if(t.x < 0.0f)
			{
				return 0u;
			}
			return 1u;
		}
		return 2u;
	}
}

float remap(
	float value,
	float oldMin,
	float oldMax,
	float newMin,
	float newMax)
{
	return newMin + (value - oldMin) / (oldMax - oldMin) * (newMax - newMin);
}

float3 sampleWeather(float3 pos)
{
	return weatherLookup.SampleLevel(texSampler, pos.xz * WEATHER_TEX_MOD.x + WEATHER_TEX_MOD.yz + (WEATHER_TEX_MOVE_SPEED * cb_appRunTime), 0.0f).rgb;
}

float getCoverage(float3 weatherData)
{
	return weatherData.r;
}

float getPrecipitation(float3 weatherData)
{
	return weatherData.g;
}

float getCloudType(float3 weatherData)
{
	// weather b channel tells the cloud type 0.0 = stratus, 0.5 = stratocumulus, 1.0 = cumulus
	return weatherData.b;
}

float heightFraction(float3 pos)
{
	return saturate((distance(pos, PLANET_CENTER) - PLANET_CENTER_TO_LOWER_CLOUD_RADIUS) / cb_cloudVolumeHeight);
}

float4 mixGradients(
	float cloudType)
{
	float stratus = 1.0f - saturate(cloudType * 2.0f);
	float stratocumulus = 1.0f - abs(cloudType - 0.5f) * 2.0f;
	float cumulus = saturate(cloudType - 0.5f) * 2.0f;
	return STRATUS_GRADIENT * stratus + STRATOCUMULUS_GRADIENT * stratocumulus + CUMULUS_GRADIENT * cumulus;
}

float densityHeightGradient(
	float heightFrac,
	float cloudType)
{
	float4 cloudGradient = mixGradients(cloudType);
	return smoothstep(cloudGradient.x, cloudGradient.y, heightFrac) - smoothstep(cloudGradient.z, cloudGradient.w, heightFrac);
}

float sampleCloudDensity(
	float3 pos,
	float3 weatherData,
	float heightFrac,
	float lod)
{
	pos += heightFrac * cb_windDirection * cb_cloudTopOffset;
	pos += (cb_windDirection + float3(0.0f, -0.25f, 0.0f)) * cb_cloudSpeed * (cb_appRunTime/* * 25.0f*/); // the * 25.0f is just for testing to make the effect obvious
	pos *= CLOUD_SCALE;
	
	float4 lowFreqNoise = baseShapeLookup.SampleLevel(texSampler, pos, lerp(0.0f, cb_baseShapeTextureBottomMipLevel, lod));
	float lowFreqFBM =
		(lowFreqNoise.g * 0.625f) +
		(lowFreqNoise.b * 0.25f) +
		(lowFreqNoise.a * 0.125f);

	float baseCloud = remap(
		lowFreqNoise.r,
		-(1.0f - lowFreqFBM), 1.0f, // gets about the same results just using -lowFreqFBM
		0.0f, 1.0f);

	float densityGradient = densityHeightGradient(heightFrac, getCloudType(weatherData));
	baseCloud *= densityGradient;

	float cloudCoverage = getCoverage(weatherData);
	float baseCloudWithCoverage = remap(
		baseCloud,
		1.0f - cloudCoverage, 1.0f,
		0.0f, 1.0f);
	baseCloudWithCoverage *= cloudCoverage;

	//// TODO add curl noise
	//// pos += curlNoise.xy * (1.0f - heightFrac);

	float3 highFreqNoise = erosionLookup.SampleLevel(texSampler, pos * 0.1f, lerp(0.0f, cb_erosionTextureBottomMipLevel, lod)).rgb;
	float highFreqFBM =
		(highFreqNoise.r * 0.625f) +
		(highFreqNoise.g * 0.25f) +
		(highFreqNoise.b * 0.125f);

	float highFreqNoiseModifier = lerp(highFreqFBM, 1.0f - highFreqFBM, saturate(heightFrac * 10.0f));

	baseCloudWithCoverage = remap(
		baseCloudWithCoverage,
		highFreqNoiseModifier * 0.2f, 1.0f,
		0.0f, 1.0f);

	return saturate(baseCloudWithCoverage);
}

struct VertexOut
{
	float4 posH : SV_POSITION;
	float3 viewRay : VIEWRAY;
	float2 tex : TEXCOORD;
};

// random vectors on the unit sphere
static const float3 RANDOM_VECTORS[] =
{
	float3( 0.38051305f,  0.92453449f, -0.02111345f),
	float3(-0.50625799f, -0.03590792f, -0.86163418f),
	float3(-0.32509218f, -0.94557439f,  0.01428793f),
	float3( 0.09026238f, -0.27376545f,  0.95755165f),
	float3( 0.28128598f,  0.42443639f, -0.86065785f),
	float3(-0.16852403f,  0.14748697f,  0.97460106f)
};

static const uint LIGHT_RAY_ITERATIONS = 6u;
static const float RCP_LIGHT_RAY_ITERATIONS = 1.0f / float(LIGHT_RAY_ITERATIONS);

float beerLambert(float sampleDensity, float precipitation)
{
	return exp(-sampleDensity * precipitation);
}

float powder(float sampleDensity, float lightDotEye)
{
	float powd = 1.0f - exp(-sampleDensity * 2.0f);
	return lerp(
		1.0f,
		powd,
		saturate((-lightDotEye * 0.5f) + 0.5f) // [-1,1]->[0,1]
	);
}

float henyeyGreenstein(
	float lightDotEye,
	float g)
{
	float g2 = g * g;
	return ((1.0f - g2) / pow((1.0f + g2 - 2.0f * g * lightDotEye), 1.5f)) * 0.25f;
}

float lightEnergy(
	float lightDotEye,
	float densitySample,
	float originalDensity,
	float precipitation)
{
	return 2.0f *
		beerLambert(densitySample, precipitation) *
		powder(originalDensity, lightDotEye) * 
		lerp(henyeyGreenstein(lightDotEye, 0.8f), henyeyGreenstein(lightDotEye, -0.5f), 0.5f);
}

// TODO get from cb values - has to change as time of day changes
float3 ambientLight(float heightFrac)
{
	return lerp(
		float3(0.5f, 0.67f, 0.82f),
		float3(1.0f, 1.0f, 1.0f),
		heightFrac);
}

float sampleCloudDensityAlongCone(
	float3 startPos,
	float  stepSize,
	float  lightDotEye,
	float  originalDensity)
{
	float3 lightStep = stepSize * -cb_lightDirection;
	float3 pos = startPos;
	float coneRadius = 1.0f;
	float coneStep = RCP_LIGHT_RAY_ITERATIONS;
	float densityAlongCone = 0.0f;
	float lod = 0.0f;
	float lodStride = RCP_LIGHT_RAY_ITERATIONS;
	float3 weatherData = 0.0f;
	float rcpThickness = 1.0f / (stepSize * LIGHT_RAY_ITERATIONS);
	float density = 0.0f;

	for(uint i = 0u; i < LIGHT_RAY_ITERATIONS; ++i)
	{
		float3 conePos = pos + coneRadius * RANDOM_VECTORS * float(i + 1u);
		float heightFrac = heightFraction(conePos);
		if(heightFrac <= 1.0f)
		{
			weatherData = sampleWeather(conePos);
			float cloudDensity = sampleCloudDensity(
				conePos,
				weatherData,
				heightFrac,
				lod);
			if(cloudDensity > 0.0f)
			{
				density += cloudDensity;
				float transmittance = 1.0f - (density * rcpThickness);
				densityAlongCone += (cloudDensity * transmittance);
			}
		}
		pos += lightStep;
		coneRadius += coneStep;
		lod += lodStride;
	}
	// take additional step at large distance away for shadowing from other clouds
	pos = pos + (lightStep * 8.0f);
	weatherData = sampleWeather(pos);
	float heightFrac = heightFraction(pos);
	if(heightFrac <= 1.0f)
	{
		float cloudDensity = sampleCloudDensity(
			pos,
			weatherData,
			heightFrac,
			0.8f);
		// no need to branch here since density variable is no longer used after this
		density += cloudDensity;
		float transmittance = 1.0f - saturate(density * rcpThickness);
		densityAlongCone += (cloudDensity * transmittance);
	}

	return saturate(lightEnergy(
		lightDotEye,
		densityAlongCone,
		originalDensity,
		lerp(1.0f, 2.0f, getPrecipitation(weatherData))));
}

float4 traceClouds(
	float3 viewDirW,		// world space view direction
	float3 startPos,		// world space start position
	float3 endPos)			// world space end position
{
	float3 dir = endPos - startPos;
	float thickness = length(dir);
	float rcpThickness = 1.0f / thickness;
	uint sampleCount = lerp(SAMPLE_RANGE.x, SAMPLE_RANGE.y, saturate((thickness - cb_cloudVolumeHeight) / cb_cloudVolumeHeight));
	float stepSize = thickness / float(sampleCount);
	dir /= thickness;
	float3 posStep = stepSize * dir;

	float lightDotEye = -dot(cb_lightDirection, viewDirW);

	float3 pos = startPos;
	float3 weatherData = 0.0f;
	float4 result = 0.0f;
	float density = 0.0f;
	
	for(uint i = 0u; i < sampleCount; ++i)
	{
		float heightFrac = heightFraction(pos);
		weatherData = sampleWeather(pos);
		float cloudDensity = sampleCloudDensity(
			pos,
			weatherData,
			heightFrac,
			0.0f);

		if(cloudDensity > 0.0f)
		{
			density += cloudDensity;
			float transmittance = 1.0f - (density * rcpThickness);
			float lightDensity = sampleCloudDensityAlongCone(
				pos,
				stepSize,
				lightDotEye,
				cloudDensity);

			float3 ambientBadApprox = ambientLight(heightFrac) * min(1.0f, length(cb_sunColor.rgb * 0.0125f)) * transmittance;
			float4 source = float4((cb_sunColor.rgb * lightDensity) + ambientBadApprox/*+ ambientLight(heightFrac)*/, cloudDensity * transmittance); // TODO enable ambient when added to constant buffer
			source.rgb *= source.a;
			result = (1.0f - result.a) * source + result;
			if(result.a >= 1.0f) break;
		}

		pos += posStep;
	}

	// experimental fog - may not be needed if clouds are drawn before atmosphere - would have to draw sun by itself, then clouds, then atmosphere
	// fogAmt = 0 to disable
	float fogAmt = 1.0f - exp(-distance(startPos, cb_eyePositionW) * 0.00001f);
	float3 fogColor = float3(0.3f, 0.4f, 0.45f) * length(cb_sunColor.rgb * 0.125f) * 0.8f;
	float3 sunColor = normalize(cb_sunColor.rgb) * 4.0f * length(cb_sunColor.rgb * 0.125f);
	fogColor = lerp(fogColor, sunColor, pow(saturate(lightDotEye), 8.0f));
	return float4(clamp(lerp(result.rgb, fogColor, fogAmt), 0.0f, 1000.0f), saturate(result.a));
}

float4 main(VertexOut pIn) : SV_TARGET
{
	int3 loadIndices = int3(pIn.posH.xy, 0);
	float zwDepth = depthTexture.Load(loadIndices).r;
	bool depthPresent = zwDepth < 1.0f;
	float depth = linearizeDepth(zwDepth);
	float3 posV = pIn.viewRay * depth;
	float3 posW = mul(float4(posV, 1.0f), cb_inverseViewMatrix).xyz;
	float3 viewDirW = normalize(posW - cb_eyePositionW);

	// find nearest planet surface point
	float2 ph = 0.0f;
	uint planetHits = intersectRaySphere(
		cb_eyePositionW,
		viewDirW,
		PLANET_CENTER,
		cb_groundRadius,
		ph);

	// find nearest inner shell point
	float2 ih = 0.0f;
	uint innerShellHits = intersectRaySphere(
		cb_eyePositionW,
		viewDirW,
		PLANET_CENTER,
		PLANET_CENTER_TO_LOWER_CLOUD_RADIUS,
		ih);

	// find nearest outer shell point
	float2 oh = 0.0f;
	uint outerShellHits = intersectRaySphere(
		cb_eyePositionW,
		viewDirW,
		PLANET_CENTER,
		PLANET_CENTER_TO_UPPER_CLOUD_RADIUS,
		oh);

	// world space ray intersections
	float3 planetHitSpot = cb_eyePositionW + (viewDirW * ph.x);
	float3 innerShellHit = cb_eyePositionW + (viewDirW * ih.x);
	float3 outerShellHit = cb_eyePositionW + (viewDirW * oh.x);

	// eye radius from planet center
	float eyeRadius = distance(cb_eyePositionW, PLANET_CENTER);
	if(eyeRadius < PLANET_CENTER_TO_LOWER_CLOUD_RADIUS) // under inner shell
	{
		// exit if there's something in front of the start of the cloud volume
		if((depthPresent && (distance(posW, cb_eyePositionW) < distance(innerShellHit, cb_eyePositionW))) ||
			planetHits > 0u) // shell hits are guaranteed, but the ground may be occluding cloud layer
		{
			return float4(0.0f, 0.0f, 0.0f, 0.0f);
		}
		return traceClouds(
			viewDirW,
			innerShellHit,
			outerShellHit);
	}
	else if(eyeRadius > PLANET_CENTER_TO_UPPER_CLOUD_RADIUS) // over outer shell
	{
		// possibilities are
		// 1) enter outer shell, leave inner shell
		// 2) enter outer shell, leave outer shell
		float3 firstShellHit = outerShellHit;
		if(outerShellHits == 0u ||
			depthPresent && (distance(posW, cb_eyePositionW) < distance(firstShellHit, cb_eyePositionW)))
		{
			return float4(0.0f, 0.0f, 0.0f, 0.0f);
		}
		float3 secondShellHit = outerShellHits == 2u && innerShellHits == 0u ? cb_eyePositionW + (viewDirW * oh.y) : innerShellHit;
		return traceClouds(
			viewDirW,
			firstShellHit,
			secondShellHit);
	}
	else // between shells
	{
		/*
		 * From a practical viewpoint (properly scaled planet, atmosphere, etc.)
		 * only one shell will be hit.
		 * Start position is always eye position.
		 */
		float3 shellHit = innerShellHits > 0u ? innerShellHit : outerShellHit;
		// if there's something in the depth buffer that's closer, that's the end point
		if(depthPresent && (distance(posW, cb_eyePositionW) < distance(shellHit, cb_eyePositionW)))
		{
			shellHit = posW;
		}
		return traceClouds(
			viewDirW,
			cb_eyePositionW,
			shellHit);
	}
}

Again, I'm sure some of the math and maybe even the general approach is off a bit, so please feel free to offer suggestions for improvement.  One thing I wish was better is that the clouds don't seem to taper well as they rise with my current implementation.  There may be some ideas from this thread I can use to help with that.  I also need to get curl noise in, but I'm more worried about getting the general solution working well first.

 

If there's anything I can help out with, let me know

 

[sharedmedia=gallery:albums:1149]

 

Edit:  Huge thank you to Andrew and the team at Guerrilla Games for their publications and willingness to help other developers work through implementing a solution.

Yaaaaaaaaaaayyyy!* Looking forward to it 

*Excited Kermit voice

Related, does anyone know how Reset does their variable distance volume marching for their precipitation effect/large scale god rays? EG: 

Near camera volume marching is fine, you set your z-slices and go. Clouds are fine if they're on a defined plane/sphere. But they also seem to use such to represent rain, a very cool effect. 

I ask because the most dramatic scattering, or god rays, often comes between clouds/distant mountains. And somehow that seems to be accomplished here. A logarithmic march? If everything (including the clouds) was in some sort of shadowmap like depth buffer you could get a summed area table, then do a gradient domain like march. EG skip space till sunlight, etc. Any other ideas?

Been working on this for a long time so I figured I’d share what I have, but I can’t find the right numbers to make the clouds look, well, good.

The highlights are overblown and the shadows are way too deep. The powder term is also making a mess here.

Here is the picture from the GPU pro book, for comparison.

Some more angles:

Silver lining effect, seems to reach too far into the interiors of the clouds, even with the absorption turned way up. Meanwhile the holdouts are too dark.

Aerial view, the powder term causes some very strange shadow casting. They actually make the clouds behind them lighter instead of darker. Clouds also look too dark and generally unrealistic.

This is the only view that looks okay, but the highlights and shadows still look bad.

Here is my source code (GLSL) if it helps. I’m brute forcing the lightmarching with 16 samples to debug the lighting model.

#version 330 core

layout (std140) uniform CameraDataBlock
{
    vec3  position;     // [ 0 ..< 12]
    float z;            // [12 ..< 16]
    vec3  direction;    // [16 ..< 28]
    float twice_size;   // [28 ..< 32]
    float half_h;       // [32 ..< 36]
    float half_k;       // [36 ..< 40]
    vec2  shift;        // [40 ..< 48]
    mat3  world;        // [48 ..< 96] size = 96
} camera;

uniform sampler3D low_freq_noise;

in Planet
{
    vec3  center;
    vec2  cloud_deck;
    float radius;
    vec3  sun;
} planet;

out vec4 color;

vec2 sphere_intersect(const vec3 ray, const vec3 camera_to_center, const float r)
{
    float midpoint = dot(camera_to_center, ray);
    float D = r*r - dot(camera_to_center, camera_to_center) + midpoint*midpoint;
    if (D <= 0)
    {
        return vec2(0, 0);
    }
    else
    {
        float half_depth = sqrt(D);
        return vec2(max(0, midpoint - half_depth), max(0, midpoint + half_depth));
    }
}

vec2 box_intersection(const vec3 p1, const vec3 p2, const vec3 start, const vec3 ray)
{
    vec3 ray_inverse = 1 / ray;
    vec3 v1 = (p1 - start) * ray_inverse;
    vec3 v2 = (p2 - start) * ray_inverse;
    vec3 n  = min(v1, v2);
    vec3 f  = max(v1, v2);
    vec2 slice = vec2(max(n.x, max(n.y, n.z)), min(f.x, min(f.y, f.z)));

    if (slice.x < slice.y)
    {
        return max(slice, 0);
    }
    else
    {
        return vec2(0, 0);
    }
}

float height_fraction(const vec3 position)
{
    return clamp(0.5 * (position.z + 1f), 0, 1);
}

float density_gradient_stratus(const float h)
{
    return max(smoothstep(0.00, 0.07, h) - smoothstep(0.07, 0.11, h), 0); // stratus, could be better
}

float density_gradient_cumulus(const float h)
{
    //return max(smoothstep(0.00, 0.22, h) - smoothstep(0.4, 0.62, h), 0); // cumulus
    return smoothstep(0.3, 0.35, h) - smoothstep(0.425, 0.7, h); // cumulus
}

float density_gradient_cumulonimbus(const float h)
{
    return smoothstep(0, 0.1, h) - smoothstep(0.7, 1, h); // cumulonimbus
}

float erode(const float x, const float e)
{
    return max(1 - (1 - x) / e, 0);
}

float cloud_density(const vec3 position, const float detail)
{
    vec4 low_freq_noise_sample = textureLod(low_freq_noise, 0.5 * (position + 1), detail);
    float low_freq_fbm = low_freq_noise_sample.g * 0.725 +
                         low_freq_noise_sample.b * 0.3  +
                         low_freq_noise_sample.a * 0.15;

    float base_cloud = erode(low_freq_noise_sample.r, density_gradient_cumulus(height_fraction(position)));
    base_cloud = erode(base_cloud, low_freq_fbm);
    base_cloud = clamp(6 * (base_cloud - 0.390625 - 0.1) + 0.5, 0, 1);

    return base_cloud;
}

float exponential_integral(const float z)
{
	return 0.5772156649015328606065 + log( 1e-4 + abs(z) ) + z * (1.0 + z * (0.25 + z * ( (1.0/18.0) + z * ( (1.0/96.0) + z * (1.0/600.0) ) ) ) ); // For x!=0
}

vec3 ambient_color(float h, float extinction_coeff)
{
    // For now make ambient color white (less instructions)
    float ambient_term = -extinction_coeff * (1.0 - h);
    vec3 isotropic_scattering_top = vec3(0.8, 0.9, 1.0) * max(0.0, exp(ambient_term) - ambient_term * exponential_integral(ambient_term));

    ambient_term = -extinction_coeff * h;
    vec3 isotropic_scattering_bottom = vec3(0.6, 0.8, 1.0) * max(0.0, exp(ambient_term) - ambient_term * exponential_integral(ambient_term));

    // Modulate the ambient term by the altitude
    isotropic_scattering_top *= h * 0.5;

    return 7*(isotropic_scattering_top + isotropic_scattering_bottom);
}

float henyey_greenstein(const float phase, const float g)
{
    float g2 = g * g;
    return 1f / (3.14159 * 4f) * (1f - g2) / pow(1f + g2 - 2f * g * phase, 1.5f);
}

float lightmarch(vec3 start, float absorption, vec3 light)
{
    vec2 slice = box_intersection(vec3(-1, -1, -0.4), vec3(1, 1, 0.4), start, -light);

    float ds = (slice.y - slice.x) / 16f;
    float s  = slice.x + 0.5f * ds;

    //float T         = 1f;
    float sigma_ds  = absorption * ds;
    float optical_depth = 0;
    for (uint i = 0u; i < 16f; ++i)
    {
        vec3 position = start - s * light;
        optical_depth += sigma_ds * cloud_density(position, 0);

        /* early exit
        if (T < 1e-8)
        {
            break;
        }
        */
        s += ds;
    }

    return optical_depth;
}

vec4 raymarch(vec3 ray, float absorption, vec2 slice)
{
    float ds = (slice.y - slice.x) / 128f;
    float s  = slice.x + 0.5f * ds;

    vec4  color     = vec4(0);
    float T         = 1f;
    float sigma_ds  = -3*absorption * ds;
    float phase     = -dot(planet.sun, ray);

    for (uint i = 0u; i < 128f; ++i)
    {
        vec3 position = camera.position + s * ray;
        float density = cloud_density(position, 0);
        float Ti = exp(sigma_ds * density);
        T *= Ti;

        if (T < 1e-8)
        {
            break;
        }

        float lightmarch_depth = lightmarch(position, absorption, planet.sun);
        float powder_factor = mix(2f*(1 - exp(-2f*lightmarch_depth)), 1f, 0.5f * (phase + 1));
        float phase_factor = mix(henyey_greenstein(phase, 0.06f), henyey_greenstein(phase, 0.9f), 0.05f);
        vec3 light_color = powder_factor * exp(-lightmarch_depth) * 75f*vec3(5f, 4f, 3f) * mix(1f/(4 * 3.14159f), phase_factor, exp(-lightmarch_depth)) +
                            ambient_color(height_fraction(position), sigma_ds * density);

        color.rgb += T * light_color * vec3(1.0f) * density * ds;
        color.a   += (1.0f - Ti) * (1.0f - color.a);
        s += ds;
    }

    return color;
}



void main()
{
    vec3 ray   = camera.world *
                 normalize(vec3(camera.twice_size * ((gl_FragCoord.x - camera.shift.x) - camera.half_h),
                                camera.twice_size * ((gl_FragCoord.y - camera.shift.y) - camera.half_k),
                               -camera.z));

    vec2 slice = box_intersection(vec3(-1, -1, -0.4), vec3(1, 1, 0.4), camera.position, ray);
    vec4 cloud_color;
    if (slice.x == slice.y)
    {
        cloud_color = vec4(0, 0, 0, 0);
    }
    else
    {
        cloud_color = raymarch(ray, 7, slice);
    }


    vec2 horizon = sphere_intersect(ray, planet.center - camera.position, planet.cloud_deck.y + 1000000);
    float optical_depth = clamp(0.00000025 * (horizon.y - horizon.x), 0, 1);
    vec3 sky = vec3(mix(optical_depth*optical_depth*optical_depth*optical_depth, 1, 0.1),
                    mix(pow(optical_depth, 1.7), 1, 0.12),
                    mix(optical_depth, 1, 0.1));

    //float fade = 1.0 / (1 + 0.0000001*slice.x*slice.x);
    sky = cloud_color.rgb + sky * (1 - cloud_color.a);
    color = vec4(sky, 1);
}

PS to the person having problems with making the clouds taper, you have to feed the altitude gradient into a remap function, not just multiply it into the cloud density. That makes the clouds shrink as you go up rather than just get lighter.

Andrew, looking forward to your SIGGRAPH 2017 Talk.

thanks for your share.

ps. the gallery pic can't show correctly now, is it because of gamedev revision? have chance to recover it?

Yeah, I think something with the site updates may have messed up the gallery link. You can find them by going to my profile and clicking the Albums tab, or use this link: 

My results look a bit like yours, taylorswift:

If I drop the absorption rate to get rid of the dark areas then the clouds get blown out. But I'm not modeling scattering yet, so I think this makes sense. The only light is that which comes directly from the sun but isn't absorbed on the way. I'll see what happens when I add ambient lighting and the henyey greenstein term but right now I'm trying to make sure what I already have is correct-ish. But it looks like you already have both. What happens if you remove them?

I'm not very happy with my cloud shapes either. I don't get very clean edges even at really high sample counts. It ends up looking pretty fuzzy. Your first image looks really good in that respect. Maybe I need some kind of exponential ramp on the density.

I also don't like how the clouds appear near the horizon but in that case I might just need some fog.

Really looking forward to Andrew's talk in a couple weeks.

2 minutes ago, wild_pointer said:

My results look a bit like yours, taylorswift:

If I drop the absorption rate to get rid of the dark areas then the clouds get blown out. But I'm not modeling scattering yet, so I think this makes sense. The only light is that which comes directly from the sun but isn't absorbed on the way. I'll see what happens when I add ambient lighting and the henyey greenstein term but right now I'm trying to make sure what I already have is correct-ish. But it looks like you already have both. What happens if you remove them?

I'm not very happy with my cloud shapes either. I don't get very clean edges even at really high sample counts. It ends up looking pretty fuzzy. Your first image looks really good in that respect. Maybe I need some kind of exponential ramp on the density.

I also don't like how the clouds appear near the horizon but in that case I might just need some fog.

Really looking forward to Andrew's talk in a couple weeks.

Removing henyey greenstein takes away the bright fringes but brightens the clouds when viewed from the sun since it’s the source of the anisotropic shading in the clouds. Ambient lighting doesn’t do much except brighten the whole scene.