IQ_Clouds

Author: IQ, website: https://www.shadertoy.com/view/XslGRr

Label: procedural, volumetric, lod

Render

The main function contains the conventional coordinate transformation and camera settings, so we directly enter the render function for analysis.

At the beginning, it is still the sky setting, including the color change of the sky and the setting of the sun.

// background sky     
float sun = clamp( dot(sundir,rd), 0.0, 1.0 );
vec3 col = vec3(0.6,0.71,0.75) - rd.y*0.2*vec3(1.0,0.5,1.0) + 0.15*0.5;
col += 0.2*vec3(1.0,.6,0.1)*pow( sun, 8.0 );

Then there is the core RayMarch function.

vec4 raymarch( in vec3 ro, in vec3 rd, in vec3 bgcol, in ivec2 px )
{
	vec4 sum = vec4(0.0);
	float t = 0.0;//0.05*texelFetch( iChannel0, px&255, 0 ).x;
    MARCH(40,map5);
    MARCH(40,map4);
    MARCH(30,map3);
    MARCH(30,map2);
    return clamp( sum, 0.0, 1.0 );
}

• For the MARCH macro definition, we need to pass in two parameters: the number of cycle steps and the type of FBM function (the number at the end represents the level of detail). Next, analyze the macro definition in detail

#define MARCH(STEPS,MAPLOD)
for(int i=0; i<STEPS; i++)
{
   vec3 pos = ro + t*rd;
   if( pos.y<-3.0 || pos.y>2.0 || sum.a>0.99 ) break;
   float den = MAPLOD( pos );
   if( den>0.01 )
   {
     float dif = clamp((den - MAPLOD(pos+0.3*sundir))/0.6, 0.0, 1.0 );
     vec3  lin = vec3(0.65,0.7,0.75)*1.4 + vec3(1.0,0.6,0.3)*dif;
     vec4  col = vec4( mix( vec3(1.0,0.95,0.8), vec3(0.25,0.3,0.35), den ), den );
     col.xyz *= lin;
     col.xyz = mix( col.xyz, bgcol, 1.0-exp(-0.003*t*t) );
     col.w *= 0.4;

     col.rgb *= col.a;
     sum += col*(1.0-sum.a);
   }
   t += max(0.05,0.02*t);
}

• First, Ray March's routine process: POS iteration and loop termination judgment. The accident here is - the ending condition is the Y coordinate (altitude) of POS.

  vec3 pos = ro + t*rd;
  if( pos.y<-3.0 || pos.y>2.0 || sum.a>0.99 ) break;
  

• Obtain the return value of the FBM function. If it is valid (i.e. greater than 0), stack and other related operations will be carried out. Firstly, the diffuse reflection coefficient dif is calculated. The core of this line of code is the subtraction of den. In terms of effect, not subtracting the latter will lead to too bright clouds. In terms of physics, FBM sampling will be carried out when the point moves to a certain position in the direction of the sun. If there are clouds in this position, it will obviously block the light of the sunlight hitting this point, so it needs to be subtracted to approximate the correct logic. This is a simple but ingenious simulation of physical laws

  float dif = clamp((den - MAPLOD(pos+0.3*sundir))/0.6, 0.0, 1.0 );
  

• Then there is the calculation of the parameter Lin. what's embarrassing is that I don't know the physical meaning of this parameter. At the same time, it has little impact on the effect, so I won't analyze it here.

  vec3  lin = vec3(0.65,0.7,0.75)*1.4 + vec3(1.0,0.6,0.3)*dif;
  

• For the setting of col, the FBM return value Den is used to make a Mix, and the w component is set to Den. The main effect is: the setting of cloud base color (gray and white).

  vec4  col = vec4( mix( vec3(1.0,0.95,0.8), vec3(0.25,0.3,0.35), den ), den );
  

• After Col and lin are multiplied, they are mixed according to the value of the current step distance and the background color. The general effect is: the smooth transition between distant clouds and the sky. Finally, the w component is multiplied by a 0.4 to scale: because the RGB of the later col needs to be multiplied by W, it needs to be scaled to avoid the final scene being too bright.

  col.xyz = mix( col.xyz, bgcol, 1.0-exp(-0.003*t*t) );
  col.w *= 0.4;
  col.rgb *= col.a;
  

• Finally, col is superimposed on sum. The logic here is: col is scaled according to (1.0-sum). The rule is also simple: the weight of the cloud closest to the viewpoint must be greater. As far as the effect is concerned, if this scaling is not carried out, the final rendering effect will have serious artifacts - covered with topographic map?

  sum += col*(1.0-sum.a);
  

• Finally, update the step size t: here to avoid too large step size

  t += max(0.05,0.02*t);
  

• Continue to analyze the tool function map4. At the beginning, add the time attribute to p to make the clouds change constantly; Then there is the acyclic FBM process. The final return value is different, which is roughly equal to the return value − 0.5 + 1.75 ∗ f -0.5+1.75*f − 0.5 + 1.75 * f, and the result is limited to the unit interval. No debugging, visual inspection should be the result of adjusting parameters on the effect.

  float map4( in vec3 p )
  {
  	vec3 q = p - vec3(0.0,0.1,1.0)*iTime;
  	float f;
      f  = 0.50000*noise( q ); q = q*2.02;
      f += 0.25000*noise( q ); q = q*2.03;
      f += 0.12500*noise( q ); q = q*2.01;
      f += 0.06250*noise( q );
  	return clamp( 1.5 - p.y - 2.0 + 1.75*f, 0.0, 1.0 );
  }
  

Tool function - Noise

float noise( in vec3 x )
{
    vec3 p = floor(x);
    vec3 f = fract(x);
	f = f*f*(3.0-2.0*f);
    
#if 1
	vec2 uv = (p.xy+vec2(37.0,239.0)*p.z) + f.xy;
    vec2 rg = textureLod(iChannel0,(uv+0.5)/256.0,0.0).yx;
#else
    ivec3 q = ivec3(p);
	ivec2 uv = q.xy + ivec2(37,239)*q.z;

	vec2 rg = mix(mix(texelFetch(iChannel0,(uv           )&255,0),
				      texelFetch(iChannel0,(uv+ivec2(1,0))&255,0),f.x),
				  mix(texelFetch(iChannel0,(uv+ivec2(0,1))&255,0),
				      texelFetch(iChannel0,(uv+ivec2(1,1))&255,0),f.x),f.y).yx;
#endif    
	return -1.0+2.0*mix( rg.x, rg.y, f.z );
}

Tool function - FBM

float map5( in vec3 p )
{
	vec3 q = p - vec3(0.0,0.1,1.0)*iTime;
	float f;
    f  = 0.50000*noise( q ); q = q*2.02;
    f += 0.25000*noise( q ); q = q*2.03;
    f += 0.12500*noise( q ); q = q*2.01;
    f += 0.06250*noise( q ); q = q*2.02;
    f += 0.03125*noise( q );
	return clamp( 1.5 - p.y - 2.0 + 1.75*f, 0.0, 1.0 );
}
float map4( in vec3 p )
{
	vec3 q = p - vec3(0.0,0.1,1.0)*iTime;
	float f;
    f  = 0.50000*noise( q ); q = q*2.02;
    f += 0.25000*noise( q ); q = q*2.03;
    f += 0.12500*noise( q ); q = q*2.01;
    f += 0.06250*noise( q );
	return clamp( 1.5 - p.y - 2.0 + 1.75*f, 0.0, 1.0 );
}
float map3( in vec3 p )
{
	vec3 q = p - vec3(0.0,0.1,1.0)*iTime;
	float f;
    f  = 0.50000*noise( q ); q = q*2.02;
    f += 0.25000*noise( q ); q = q*2.03;
    f += 0.12500*noise( q );
	return clamp( 1.5 - p.y - 2.0 + 1.75*f, 0.0, 1.0 );
}
float map2( in vec3 p )
{
	vec3 q = p - vec3(0.0,0.1,1.0)*iTime;
	float f;
    f  = 0.50000*noise( q ); q = q*2.02;
    f += 0.25000*noise( q );;
	return clamp( 1.5 - p.y - 2.0 + 1.75*f, 0.0, 1.0 );
}