Planets (part 8)

Hybrid approach for procedural planets
Part VIII - Physically based clouds

The atmosphere is complex. Very complex. Now I know that's a bit like saying an atomic blast is a bit noisy, but it is very true. Whole arrays of super computers are used to attempt to model our planets atmosphere, and they still cannot tell you if it is going to rain tomorrow.
However we all can look at the sky and recognise features. There are many types of clouds, and we spend most of our life looking at them, so when we see a game with clouds in the sky we can tell if they look right or not.
Most people just use noise for clouds, Perlin is the favorite, but I wanted to see if we could do better. Which kicked off three weeks of intense coding. I am still not 100% happy with the results, but it is starting to look good and maybe someone out there can build on this. So here we go.

Building blocks

There are many physical factors that affect cloud formation, if you are good at math, get a copy of Cloud Dynamics by Houze and scare yourself.
We are going to take the most important ones and then cheat.
The basic building blocks we need to start with are...

Terrain
Heat
Friction
Pressure

These are pretty easy to generate so I will skip over them quite quickly.

Terrain
For this, just to make it more familiar for you, I have grabbed a heightmap of the Earth.

Heat
The heat map is a function of the terrain and the amount of star light that lands on the terrain.
We can do this very easily with a shader.
Note that water has a very high specific heat capacity, so I have adjusted the heat map in the water areas.
esun is a variable passed in from the planet definition and represents the amount of energy the star is throwing at the planet.
This texture is included for you to use in your own research. I ended up cheating in the simulation and I didn't need it, but I will explain it's potential uses later.

float4 PixelShaderFunction(float2 texCoord : TEXCOORD) : COLOR0
{
    float4 result=float4(0,0,0,1);
    float4 land = tex2D(HeightMapSampler,texCoord);
    if (land.r<=0.001)
    {
        land.r=0.5;
    }
    float  angle=3.141597*texCoord.y;
    float  intensity=sin(angle);
    result.r=(((1-land.r)*intensity*intensity)+(intensity/2))*esun;
    result.a=1;
    return result;
}

Friction
The terrain applies a frictional force to the atmosphere. This is a subtle effect, but if we want mountains to produce clouds we need to add it.

float4 PixelShaderFunction(float2 texCoord : TEXCOORD) : COLOR0
{
    float4 result=float4(0,0,0,1);
    float4 land = tex2D(HeightMapSampler,texCoord);

    float friction = lerp(cwater,cland,land.x);
    result.g=friction;
    result.a=1;
    return result;
}

Pressure
We define a new variable in the planet definition PlanetWeatherIntensity.
This is a representation of how violent the planets weather system is.
I am just using this as the number of weather systems to add to the planet.
We start by adding this number of random points
to a list.

public void Generate(GraphicsDevice graf, ContentManager cont, SpriteBatch sp)
{
    Random rand = new Random();
    HighPressureRegions.Clear();
    for (int i = 0; i < Game1.PlanetWeatherIntensity; i++)
    {
        HighPressureRegions.Add(new Vector2((float)rand.NextDouble(),(float)rand.NextDouble()));
    }
    graphics = graf;
    build();
}

To turn this into something we can use we need to add up all the pressure contributions from each point for each pixel.
This would require a lot of passes in a shader.
You would need to do a pass for each point, then recover the texture and use it in the next pass.
So I chose to do this in C#.

private void build()
{
    float[] map = new float[Game1.TextureHeight * Game1.TextureWidth];
    int pos = 0;
    float max = 0;
    for (int y = 0; y < Game1.TextureHeight; y++)
    {
        for (int x = 0; x < Game1.TextureWidth; x++)
        {
            float z = 0;
            Vector2 ppos = new Vector2((float)x / (float)Game1.TextureWidth, (float)y / (float)Game1.TextureHeight);
            for (int i = 0; i < Game1.PlanetWeatherIntensity; i++)
            {
                float dx = Math.Abs(ppos.X - HighPressureRegions[i].X);
                if (dx > 0.5)
                    dx = 1 - dx;
                float dy = Math.Abs(ppos.Y - HighPressureRegions[i].Y);
                if (dy > 0.5)
                    dy = 1 - dy;
                float d = (dx * dx) + (dy * dy);
                if (d < 0.0001f)
                    d = 0.0001f;
                z += 1.0f / d;
            }
            if (z > max) max = z;
            map[pos] = z;
            pos++;
        }
    }
    pos = 0;

    for (int y = 0; y < Game1.TextureHeight; y++)
    {
        for (int x = 0; x < Game1.TextureWidth; x++)
        {
            map[pos++] /= max;
        }
    }
    texture = new Texture2D(graphics, Game1.TextureWidth, Game1.TextureHeight, 1, TextureUsage.None, SurfaceFormat.Single);
    texture.SetData<float>(map);
}

Note that I have used a floating point format for this texture to give us some accuracy.
Now we have the basic building blocks, we need to calculate the velocity of the atmosphere.
The intial velocity is calculated based on the four maps above.
The primary force is due to differences in pressure, so we calculate the average pressure difference by looking at all the pixels around the current one

// atmospheric pressure at the current point
float4 here = tex2D(HeightMapSampler,texCoord);
float2 uv = float2(0,0);

for (int i = 0; i < SAMPLE_COUNT; i++)
{
    // for all the points around the current point,
    // calculate the pressure difference and use that as a force

    float2 samplePoint=texCoord + SampleOffsets[i];
    float4 localPressure=tex2D(HeightMapSampler,samplePoint);
    float pressureDifference=localPressure.x-here.x; // -1 to +1
    uv+=pressureDifference*SampleOffsets[i]; // -gridsize to +gridsize
}
uv*=200000;

The next factor we need to estimate is friction. We have the base value from the map above so we can add it fairly easily.

// estimate the friction force
float4 friction = tex2D(FrictionMapSampler,texCoord); // 0 to 1
float fl =length(uv);
float fx = (fl*fl*friction.y)+(fl*0.25*friction.y);
uv.x-=fx/128;
uv.y-=fx/128;

One of the major factors we need to add is the corriolis effect, this is what gives us our circular weather systems.
This is proportional to the speed relative to the axis of rotation of the planet, so we can cheat and just use the X velocity.
We also need to scale the result to fit in the range 0-1, so we might as well do that now.

// estimate the coriolis force
uv.y+=(uv.x-uv.y)*0.5;

uv/=2;
uv+=0.5;

The end result of this shader is the initial velocity of the atmosphere we can feed into the next stage of the simulation.

As you can see it has nice circular artifacts which will lead to weather systems in the finished map.

Now we start to get into the interesting bit. We treat the atmosphere as a great big ball of fluid. We use the velocity map to move this ball of fluid around. Fluid simulation is another complex subject, but thankfully there has been quite a lot of work published we can draw on to make this work. The first thing we have to do is advection. This is the process of moving the atmosphere around based purely on the velocities. We start with the initial velocity we calculated earlier and use this to advect the velocity, and a density map. Both results need to be stored in new textures. I will skip over the supporting C# code here, it's in the zip anyway so you can look at it at your leisure, but it basically just creates a render target, runs the shaders, and recovers the results. These are the pixel shaders for advection.

//===============================================================================================
//= This program performs a semi-lagrangian advection of a passive field by a moving
//= velocity field.
//=
//===============================================================================================
float4 AdvectVelocity(float2 Coords : TEXCOORD0) : COLOR0
{
    float2 pos = Coords;
    float2 vel = tex2D(VelocityMapSampler, Coords).xy;
    vel=2*(vel-0.5);
    pos -= timestep * rdx * vel;
    float4 xNew = dissipation * bilerp(VelocityMapSampler, pos); 
    xNew.w=1;
    return xNew;
}
float4 AdvectDensity(float2 Coords : TEXCOORD0) : COLOR0
{
    float2 vel = tex2D(VelocityMapSampler, Coords).xy;
    vel=2*(vel-0.5);

    float2 pos = Coords - timestep * rdx * vel;
    float4 xNew = dissipation * bilerp(DensityMapSampler, pos); 
    xNew.w=1;
    return xNew;
}

The function bilerp is here and the variable rdx is 1/texture size.

//===============================================================================================
//= These methods perform texture lookups at the four nearest neighbors of the position s and
//= bilinearly interpolate them.
//===============================================================================================
float4 bilerp(sampler tex, float2 s)
{
    float4 st;
    st.x = -rdx;
    st.y = -rdx;
    st.z = rdx;
    st.w = rdx;
    st.xy+=s;
    st.zw+=s;
 
    float2 t = float2(0.5,0.5);
   
    float4 tex11 = tex2D(tex, st.xy);
    float4 tex21 = tex2D(tex, st.zy);
    float4 tex12 = tex2D(tex, st.xw);
    float4 tex22 = tex2D(tex, st.zw);

    // bilinear interpolation
    return lerp(lerp(tex11, tex21, t.x), lerp(tex12, tex22, t.x), t.y);
}

When fluids move we get vortices. We want these in the clouds as it gives us extra detail. So the next stage is to find these.

//===============================================================================================
//= Vorticity
//=
//===============================================================================================
float4 vorticity(float2 Coords : TEXCOORD0) : COLOR0
{
  float4 uL, uR, uB, uT;
  neighbors(VelocityMapSampler, Coords);
 
  float vort = 10000*halfrdx * ((top.y - left.y) - (right.x - bottom.x));
  return float4(vort,0,0,1);
}

This will give you something like this.
Now we need to apply this to the velocity field.

//===============================================================================================
//= The second pass of vorticity confinement computes a vorticity confinement
//= force field and applies it to the velocity field to arrive at a new velocity field.
//===============================================================================================
float4 vortForce(float2 Coords : TEXCOORD0) : COLOR0
{
    float vL, vR, vB, vT, vC;
    h1neighbors(VorticityMapSampler, Coords, vL, vR, vB, vT);
 
    vC = tex2D(VorticityMapSampler, Coords).x;
 
    float2 force = halfrdx * float2(abs(vT) - abs(vB), abs(vR) - abs(vL));
 
    // safe normalize
    float magSqr = max(EPSILON, dot(force, force));
    force = force * sqrt(magSqr);
 
    force *= dxscale * vC * float2(1, -1);

    float4 uNew = tex2D(VelocityMapSampler, Coords);

    uNew.xy += timestep * force;
    return uNew;
}

Now we use the newly modified velocity map to calculate the divergence.
Divergence is the amount of stuff being moved by the velocity field.
This is the really important bit for me as I use the divergence directly as my clouds.
The idea I had was simple, if the planet started off with a uniform density field, then the change in density IS the divergence.
So I can cheat and just use this as the final density field.
I know this is a bit of a stretch, but it will do for me for the time being.

//===============================================================================================
//= This program computes the divergence of the specified vector field "velocity".
//= The divergence is defined as
//=     "grad dot v" = partial(v.x)/partial(x) + partial(v.y)/partial(y),
//=
//= and it represents the quantity of "stuff" flowing in and out of a parcel of fluid.
//= Incompressible fluids must be divergence-free.  In other words this quantity must be zero
//= everywhere.
//===============================================================================================
float4 divergence(float2 Coords : TEXCOORD0) : COLOR0
{
    float4 vL, vR, vB, vT;
    neighbors(VelocityMapSampler, Coords);
 
    float4 div =  (right.x - left.x + top.y - bottom.y);
    div*=10;
    div.w=1;
    return div;
}

The function neighbours is here. I had to use global variables in it, for some reason the shader compiler just wouldn't accept out float4 in the function definition.

//===============================================================================================
//= get neighbouring values
//=
//===============================================================================================
void neighbors(sampler tex, float2 s)
{
    left   = tex2D(tex, s - float2(rdx, 0));
    right  = tex2D(tex, s + float2(rdx, 0));
    bottom = tex2D(tex, s - float2(0, rdx));
    top    = tex2D(tex, s + float2(0, rdx));
}

Nice?
You could just stop here and use this as your cloud map, but I find you have some artifacts in the display because the fluid hasn't settled down yet. So I do a bit more work. The next stage is to perform a Jacobi relaxation step on the divergence map and the pressure map. You need to do this in a loop, really you should do about 50 passes, but I find 10 works for me.

//===============================================================================================
//= This program performs a single Jacobi relaxation step for a poisson equation of the form
//=
//=              Laplacian(U) = b
//=
//= where U = (u, v) and Laplacian(U) is defined as
//=
//= grad(div x) = grad(grad dot x) = partial^2(u)/(partial(x))^2 + partial^2(v)/(partial(y))^2
//=
//= A solution of the equation can be found iteratively, by using this iteration:
//=
//=   U'(i,j) = (U(i-1,j) + U(i+1,j) + U(i,j-1) + U(i,j+1) + b) * 0.25
//=
//===============================================================================================
float4 jacobi(float2 Coords : TEXCOORD0) : COLOR0 
{
    float4 xNew=float4(0,0,0,0);
 
    float xL, xR, xB, xT;
    neighbors(DivergenceMapSampler, Coords);
 
    float bC = tex2D(PressureMapSampler, Coords).x;

    xNew.x = (left.x + right.x + bottom.x + top.x + alpha * bC) * rBeta;
    xNew.w=1;
    return xNew;
}

The final stage is to apply the changes in velocity caused by the changes in the pressure field back into the velocity field.
Yes another shader.

//===============================================================================================
//= This program implements the final step in the fluid simulation.  After the poisson solver
//= has iterated to find the pressure disturbance caused by the divergence of the velocity
//= field, the gradient of that pressure needs to be subtracted from this divergent velocity to
//= get a divergence-free velocity field:
//=
//= v-zero-divergence = v-divergent -  grad(p)
//=
//= The gradient(p) is defined:
//=     grad(p) = (partial(p)/partial(x), partial(p)/partial(y))
//=
//= The discrete form of this is:
//=     grad(p) = ((p(i+1,j) - p(i-1,j)) / 2dx, (p(i,j+1)-p(i,j-1)) / 2dy)
//=
//= where dx and dy are the dimensions of a grid cell.
//=
//= This program computes the gradient of the pressure and subtracts it from the velocity to
//= get a divergence free velocity.
//===============================================================================================
float4 gradient(float2 Coords : TEXCOORD0) : COLOR0 
{
    float pL, pR, pB, pT;
 
    h1neighbors(PressureMapSampler, Coords, pL, pR, pB, pT);

    float2 grad = float2(pR - pL, pT - pB) * rdx;

    float4 uNew = tex2D(VelocityMapSampler, Coords);
    uNew.xy -= grad;
    uNew.w=1;
    return uNew;
}

void h1neighbors(sampler tex, float2 s,out float left,out float right, out float bottom, out float top)
{
    left   = tex2D(tex, s - float2(rdx, 0)).x;
    right  = tex2D(tex, s + float2(rdx, 0)).x;
    bottom = tex2D(tex, s - float2(0, rdx)).x;
    top    = tex2D(tex, s + float2(0, rdx)).x;
}

Few.
To get the fluid simulation to settle into a fairly steady state, we need to loop over all those shaders a few times.
The more you loop the smoother the result. I use 50 iterations.
At the end of all that you have something like this.

A lot of work has gone into getting the simulation this close, but it's not right yet.
If you look at the pacific off the coast east of the Philipines you can see linear artifacts, I would like to get rid of them.
Secondly, although a lot of this work is based on the real physical systems, I have cheated by just using the divergence field.
What I would like to do at some stage is to use the height and heat maps to inject water vapour into the density map, then allow the system to run for a while and see what turns up.
Then I want to use the heat, pressure, and denisty maps to model rainfall.
But that can wait for another day.
New zip with all this rubbish integrated into the planet renderer. here
After all that I would like to move onto a nice easy subject, but sadly the next stage in our planetary creation quest is atmospheric scattering.

Next part of the tutorial

Hybrid approach for procedural planetsPart VIII - Physically based clouds

Building blocks

Hybrid approach for procedural planets
Part VIII - Physically based clouds