Markyparky56/seamlessterrainblending.glsl

## seamlessterrainblending.glsl
// Check other page for Noise function
const float PI = 3.141592653589793238462643383279;
const float seed0 = 0.0;

// Rotations to hide some symmetry/bias, 4D rotations as the rotation must be applied to the
// points in 4D space to maintain seamlessness, rotations are arbitrary
// Matrices from here, only XZ and ZU are used here but any combination should work
const mat4 rot0 = mat4(
    cos(0.77), 0, -cos(0.77), 0,
    0, 1, 0, 0,
    sin(0.77), 0, cos(0.77), 0,
    0, 0, 0, 1);
const mat4 rot1 = mat4(
    1, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, cos(-0.23), -sin(-0.23),
    0, 0, sin(-0.23), cos(-0.23));

const int A = 0, // Billow
    	  B = 1, // Terrace
    	  C = 2, // Mountains/Ridged
    	  D = 3; // FBM

struct Biome
{
  float bias, scale;
};
Biome biomes[4] = Biome[4](
    	Biome(0.0, 0.5),
    	Biome(-0.25, 1.0),
    	Biome(0.0, 1.0),
    	Biome(-0.25, 0.75)
    );

struct BiomeSquare
{
  int biomes[4];
  float t0, t1;
  float m0, m1;
};
//const int numBiomeSquare = 1;
BiomeSquare biomeSquare[1] = BiomeSquare[1](
    BiomeSquare(int[4](A, B, C, D), 0.0, 1.0, 0.0, 1.0)
);

int findBiomes(float temp, float moisture)
{
    for(int i = 0; i < 1; i++)
    {
     	if((temp >= biomeSquare[i].t0 && temp <= biomeSquare[i].t1)
        && (moisture > biomeSquare[i].m0 && moisture <= biomeSquare[i].m1))
        {
        	return i;
        }
    }
}

float InterpQuinticFunc(float t) { return t*t*t*(t*(t * 6.0 - 15.0) + 10.0); }
float Lerp(float a, float b, float t) { return a + t * (b -a); }

float shiftRange(float v)
{
    return (v*0.5)+0.5;
}
vec3 shiftRange(vec3 v)
{
    return (v*0.5)+0.5;
}
float unshiftRange(float v)
{
    return (v-0.5)*2.0;
}

vec2 grad4to2(vec4 grad, float r, float theta, float phi)
{
    return vec2(
        	(-r * grad.x * sin(theta)) + (r * grad.y * cos(theta))
        , 	(-r * grad.z * sin(phi))   + (r * grad.w * cos(phi))
        );
}

// Have a non-gradient returning version of snoise!
float basicFbm(in vec2 x, float r, int octaves, float lac, float gain, float seed)
{
    float v = 0.0, amp = 1.0;
    vec4 gradDummy;

    float theta = (x.x) * 2.0 * PI;
    float phi = (x.y) * 2.0 * PI;

    for(int i = 0; i < octaves; i++)
    {
        vec4 p4 = vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
    	v += amp * snoise(p4, gradDummy, seed);
        amp *= gain;
        r *= lac;
    }

    return v;
}

vec4 basicFbm4(vec2 x, float r, int octaves, float lac, float gain)
{
    vec4 v = vec4(0.0); float amp = 1.0;
    vec4 gradDummy;

    float theta = (x.x) * 2.0 * PI;
    float phi = (x.y) * 2.0 * PI;

    for(int i = 0; i < octaves; i++)
    {
        vec4 p4 = vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
    	v += amp * vec4(  snoise(p4 + vec4( 0.1,  1.4, -5.2,  0.87), gradDummy, seed0)
                        , snoise(p4 + vec4( 2.4, -6.4,  3.4,  1.85), gradDummy, seed0)
                        , snoise(p4 + vec4( 2.7, -2.0,  5.9,  1.3 ), gradDummy, seed0)
                        , snoise(p4 + vec4(-3.3, -2.0,  1.7, -3.0 ), gradDummy, seed0));
        amp *= gain;
        r *= lac;
    }

    return v;
}

float terrain_fbm( in vec2 x, in vec4 pOffsets, in float r)
{
    float lac = 2.0, v = 0.0, amp = 1.0;
	vec2  dsum = vec2(0.0);

    float theta = (x.x) * 2.0 * PI;
    float phi = (x.y) * 2.0 * PI;

    for( int i=0; i<6; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);
        float n = snoise(p4, grad, seed0);
        dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
        v += amp*n/(1.0+dot(dsum,dsum));
		amp *= 0.65;
        r *= lac;
    }
	amp = 0.25;
    for( int i=0; i<5; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);
        float n = snoise(p4, grad, seed0)*0.15;
        dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
        v += amp*n;
		amp *= 0.25;
        r *= lac;
    }

	return v;
}

float terrain_ridged( in vec2 x, in vec4 pOffsets, in float r)
{
	float lac = 2.0, v = 0.0, amp = 1.0;
    vec2  dsum = vec2(0.0);

    float theta = (x.x) * 2.0 * PI;
    float phi = (x.y) * 2.0 * PI;

    for( int i=0; i<5; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);
        float n = (2.0 * pow(1.0 - abs(snoise(p4, grad, seed0)), 4.0)-1.0);
        dsum += grad4to2(grad, r, theta, phi)*(0.35/r);
        v += amp*n/(1.0+dot(dsum,dsum));
		amp *= 0.5;
        r *= lac;
    }
    r = 8.0; amp = 0.20;
    for( int i=0; i<5; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);
        float n = snoise(p4, grad, seed0);
        v += amp*n;
		amp *= 0.7;
        r *= lac;
    }

	return v/1.5;
}

float terrain_billow( in vec2 x, in vec4 pOffsets, in float r)
{
    //x *= rot;
	float lac = 1.8, v = -1.0, amp = 1.0;
	vec2  dsum = vec2(0.0);

    float theta = (x.x) * 2.0 * PI;
    float phi = (x.y) * 2.0 * PI;

    for( int i=0; i<4; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);

        float n = abs(snoise(p4, grad, seed0)) * 1.5;
        dsum += grad4to2(grad, r, theta, phi)*(0.125/r);
        v += amp*n/(1.0+dot(dsum,dsum));
		amp *= 0.6;
        r *= lac;

    }
    v = clamp(v, 0.0, v);
    r = 4.0; amp = 0.6;
    for( int i=0; i<6; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);

        float n = snoise(p4, grad, seed0);
        v += amp*n;
		amp *= 0.65;
        r *= lac;
    }

	return v;
}

float terrace_few( in vec2 x, in vec4 pOffsets, in float r)
{
	float v = 0.0, amp = 1.0, lac = 1.75;
    float bandWidth = 0.4;
	vec2  dsum = vec2(0.0);

    float theta = (x.x) * 2.0 * PI;
    float phi = (x.y) * 2.0 * PI;

    for( int i=0; i<5; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);

        float n = snoise(p4, grad, seed0);
        dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
        v += amp*n/(1.0+dot(dsum,dsum));
		amp *= 0.55;
        r *= lac;
    }
    if(v > -0.66 && v < 0.8)
    {
        float k = floor(v/bandWidth);
        float f = (v - k*bandWidth) / bandWidth;
        float s = min(2.0*f, 1.0);
        float blend = smoothstep(-0.66, -0.33, v) * (1.0 - smoothstep(0.66, 0.8, v));
        v = ((k + s) * bandWidth)*blend + (v * (1.0 - blend));
    }
    r = 6.0;
    for( int i=0; i<4; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);

        float n = snoise(p4, grad, seed0);
        v += amp*n;
		amp *= 0.5;
        r *= lac;
    }
	return v;
}

float terrace_many( in vec2 x, in vec4 pOffsets, in float r)
{
	float v = 0.0, amp = 1.0, lac = 1.75;
    float bandWidth = 0.25;
	vec2  dsum = vec2(0.0);

    float theta = (x.x) * 2.0 * PI;
    float phi = (x.y) * 2.0 * PI;

    for( int i=0; i<5; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);

        float n = snoise(p4, grad, seed0);

        dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
        v += amp*n/(1.0+dot(dsum,dsum));
		amp *= 0.55;
        r *= lac;
    }
    if(v > -0.5 && v < 0.75)
    {
        float k = floor(v/bandWidth);
        float f = (v - k*bandWidth) / bandWidth;
        float s = min(2.0*f, 1.0);
        float blend = smoothstep(-0.5, -0.33, v) * (1.0 - smoothstep(0.66, 0.75, v));
        v = ((k + s) * bandWidth)*blend + (v * (1.0 - blend));
    }
    r = 6.0;
    for( int i=0; i<8; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);

        float n = snoise(p4, grad, seed0);

        dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
        v += amp*n;
		amp *= 0.5;
        r *= lac;
    }
	return v;
}

float terrace_many_warped( in vec2 x, in vec4 pOffsets, in float r)
{
	float v = 0.0, amp = 1.0, lac = 1.75;
    float bandWidth = 0.25;
	vec2  dsum = vec2(0.0);
    int w1o = 4;
    //int w2o = 2;
    //int w3o = 1;

    float theta = (x.x) * 2.0 * PI;
    float phi = (x.y) * 2.0 * PI;

    for( int i=0; i<2; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);

        vec4 warp4 = (p4 + basicFbm4(x, r, w1o, lac*0.5, 0.5));
        //vec4 warp4_2 = (warp4 + basicFbm4(x, r, w2o, lac, 0.5));
        //vec4 warp4_3 = (warp4_2 + basicFbm4(x, r, w3o, lac, 0.5));


        float n = snoise(warp4, grad, seed0);

        dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
        v += amp*n/(1.0+dot(dsum,dsum));
		amp *= 0.55;
        r *= lac;
    }
    if(v > -1.0 && v < 0.7)
    {
        float k = floor(v/bandWidth);
        float f = (v - k*bandWidth) / bandWidth;
        float s = min(2.0*f, 1.0);
        float blend = smoothstep(-1.0, -0.6, v) * (1.0 - smoothstep(0.6, 0.7, v));
        v = ((k + s) * bandWidth)*blend + (v * (1.0 - blend));
    }
    v = clamp(v, -1.0, 0.7);
    r = 10.0; amp = 0.05;
    for( int i=0; i<1; i++ )
    {
        vec4 p4 = pOffsets + vec4(
                (r * cos(theta))
            ,	(r * sin(theta))
            , 	(r * cos(phi))
            , 	(r * sin(phi))
        );
        p4 *= rot0;
        p4 *= rot1;
        vec4 grad = vec4(0.0);

        float n = snoise(p4, grad, seed0);

        dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
        v += amp*n;
		amp *= 0.5;
        r *= lac;
    }
	return v;
}

vec3 terrainBlend(in vec2 uv)
{
    vec4 pOffsets = vec4(
        	 123.456
        ,	-432.912
        ,	-198.023
        ,	 543.298
        );
    float terrainHeight = (basicFbm(uv, 1.0, 5, 2.2, 0.55, 24.0)); // Just basic fbm for terrain height guide

    float terrainTemp;
    float my = mod(uv.y, 1.0);
    //terrainTemp = (my < 0.5) ? (my - 0.25)*4.0 : (my - 0.75)*-4.0; // Gradient centred on 0.5

    terrainTemp = (my < 0.5) ? (my*2.0) : (2.0-(my*2.0)); // Gradient centred on 0.5
	terrainTemp *= shiftRange(basicFbm(uv, 2.0, 6, 1.75, 0.5, 24.0));

    terrainTemp = unshiftRange(terrainTemp);
    terrainTemp -= (0.3  * terrainHeight)*smoothstep(-0.15, 0.10, terrainHeight)*(1.0 - smoothstep(0.10, 0.3, terrainHeight));
    terrainTemp -= (0.4  * terrainHeight)*smoothstep( 0.10, 0.25, terrainHeight)*(1.0 - smoothstep(0.45, 0.65, terrainHeight));
    terrainTemp -= (0.6  * terrainHeight)*smoothstep( 0.45, 0.65, terrainHeight)*(1.0 - smoothstep(0.7, 0.85, terrainHeight));
    terrainTemp -= (0.7  * terrainHeight)*smoothstep( 0.7, 0.85, terrainHeight);
    terrainTemp = shiftRange(terrainTemp);
    terrainTemp = terrainTemp*2.25; // Strengthen gradient around centre
    terrainTemp = clamp(terrainTemp, 0.0, 1.0);

    float terrainMoisture = shiftRange(basicFbm(uv, 2.0, 4, 2.0, 0.5, 12.0));
    terrainMoisture += 1.0 - smoothstep(-1.0, -0.1, terrainHeight);
    terrainMoisture = clamp(terrainMoisture, 0.0, 1.0);

    // Calculate biomes & respective weights
    //int biome = findBiomes(terrainTemp, terrainMoisture);
    float tBlend, mBlend;
    //tBlend = InterpQuinticFunc(terrainTemp - biomeSquare[biome].t0);
    //mBlend = InterpQuinticFunc(terrainMoisture - biomeSquare[biome].t1);
    tBlend = InterpQuinticFunc(terrainTemp);
    mBlend = InterpQuinticFunc(terrainMoisture);

    // Sample four nearest biomes
    float values[4];
    //for(int i = 0; i < 4; i++)
    //{
    //   switch(biomeSquare[biome].biomes[i])
    //   {
    //   case A: values[i] = terrain_billow(uv, pOffsets);
    //   case B: values[i] = terrace_many(uv, pOffsets);
    //   case C: values[i] = terrain_ridged(uv, pOffsets);
    //   case D: values[i] = terrain_fbm(uv, pOffsets);
    //   default: values[i] = 0.0;
    //   }
    //}

    values[0] = terrain_billow(uv, pOffsets, 3.0) * 0.5;
    values[1] = (terrace_few(uv, pOffsets, 2.5) - 0.25);
    values[2] = terrain_ridged(uv, pOffsets, 3.0);
    values[3] = (terrain_fbm(uv, pOffsets, 2.0) - 0.25) * 0.75;

    // Calculate new height map
    float heightmap = Lerp( Lerp(values[0], values[1], tBlend)
                           ,Lerp(values[2], values[3], tBlend), mBlend);

    //return vec3(terrainTemp, 0.0, 0.0);
    //return vec3(terrainTemp, terrainMoisture, heightmap);
    //return vec3(tBlend, mBlend, heightmap);
    //return vec3(terrainHeight, terrainHeight, terrainHeight);
    return vec3(heightmap, heightmap, heightmap);


    //return terrainHeight;
    //return terrainTemp;
    //return terrainMoisture;
    //return heightmap;
    //return 0.0;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (fragCoord)/iResolution.y;
    // Set to 1 to enable cyclic zoom in to show increased detail
#if 0
    uv.x += sin(iTime*0.35)*0.5;
    uv.y += cos(iTime*0.35)*0.5;
    uv = uv*((sin(iTime*0.25)*0.5 + 0.5) + 0.125);
#else
    uv *= 1.125;
    //uv *= 0.5;
    uv.y -= 0.06;
  	uv.x -= iTime*0.1;
#endif
    vec4 pOffsets = vec4(
        	 123.456
        ,	-432.912
        ,	-198.023
        ,	 543.298
        ); //pOffsets += iTime * 0.1;

	//float v = basicFbm(uv, 1.0, 7, 1.55, 0.6);
    //float v = terrainBlend(uv);
    //v = shiftRange(v);
    vec3 v= terrainBlend(uv);
    v = shiftRange(v);
    //v.z = 0.0;

    // Output to screen
    if((mod(uv.x, 1.0025) <= 1.0 || mod(uv.x, 1.0025) >= 1.0025)
    && (mod(uv.y, 1.0025) <= 1.0 || mod(uv.y, 1.0025) >= 1.0025)) fragColor = vec4(v, 1.0);
    else fragColor = vec4(1.0, 0.0, 1.0, 1.0);

}
	// Check other page for Noise function
	const float PI = 3.141592653589793238462643383279;
	const float seed0 = 0.0;

	// Rotations to hide some symmetry/bias, 4D rotations as the rotation must be applied to the
	// points in 4D space to maintain seamlessness, rotations are arbitrary
	// Matrices from here, only XZ and ZU are used here but any combination should work
	const mat4 rot0 = mat4(
	cos(0.77), 0, -cos(0.77), 0,
	0, 1, 0, 0,
	sin(0.77), 0, cos(0.77), 0,
	0, 0, 0, 1);
	const mat4 rot1 = mat4(
	1, 0, 0, 0,
	0, 1, 0, 0,
	0, 0, cos(-0.23), -sin(-0.23),
	0, 0, sin(-0.23), cos(-0.23));

	const int A = 0, // Billow
	B = 1, // Terrace
	C = 2, // Mountains/Ridged
	D = 3; // FBM

	struct Biome
	{
	float bias, scale;
	};
	Biome biomes[4] = Biome[4](
	Biome(0.0, 0.5),
	Biome(-0.25, 1.0),
	Biome(0.0, 1.0),
	Biome(-0.25, 0.75)
	);

	struct BiomeSquare
	{
	int biomes[4];
	float t0, t1;
	float m0, m1;
	};
	//const int numBiomeSquare = 1;
	BiomeSquare biomeSquare[1] = BiomeSquare[1](
	BiomeSquare(int[4](A, B, C, D), 0.0, 1.0, 0.0, 1.0)
	);

	int findBiomes(float temp, float moisture)
	{
	for(int i = 0; i < 1; i++)
	{
	if((temp >= biomeSquare[i].t0 && temp <= biomeSquare[i].t1)
	&& (moisture > biomeSquare[i].m0 && moisture <= biomeSquare[i].m1))
	{
	return i;
	}
	}
	}

	float InterpQuinticFunc(float t) { return ttt(t(t * 6.0 - 15.0) + 10.0); }
	float Lerp(float a, float b, float t) { return a + t * (b -a); }

	float shiftRange(float v)
	{
	return (v*0.5)+0.5;
	}
	vec3 shiftRange(vec3 v)
	{
	return (v*0.5)+0.5;
	}
	float unshiftRange(float v)
	{
	return (v-0.5)*2.0;
	}

	vec2 grad4to2(vec4 grad, float r, float theta, float phi)
	{
	return vec2(
	(-r * grad.x * sin(theta)) + (r * grad.y * cos(theta))
	, (-r * grad.z * sin(phi)) + (r * grad.w * cos(phi))
	);
	}

	// Have a non-gradient returning version of snoise!
	float basicFbm(in vec2 x, float r, int octaves, float lac, float gain, float seed)
	{
	float v = 0.0, amp = 1.0;
	vec4 gradDummy;

	float theta = (x.x) * 2.0 * PI;
	float phi = (x.y) * 2.0 * PI;

	for(int i = 0; i < octaves; i++)
	{
	vec4 p4 = vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	v += amp * snoise(p4, gradDummy, seed);
	amp *= gain;
	r *= lac;
	}

	return v;
	}

	vec4 basicFbm4(vec2 x, float r, int octaves, float lac, float gain)
	{
	vec4 v = vec4(0.0); float amp = 1.0;
	vec4 gradDummy;

	float theta = (x.x) * 2.0 * PI;
	float phi = (x.y) * 2.0 * PI;

	for(int i = 0; i < octaves; i++)
	{
	vec4 p4 = vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	v += amp * vec4( snoise(p4 + vec4( 0.1, 1.4, -5.2, 0.87), gradDummy, seed0)
	, snoise(p4 + vec4( 2.4, -6.4, 3.4, 1.85), gradDummy, seed0)
	, snoise(p4 + vec4( 2.7, -2.0, 5.9, 1.3 ), gradDummy, seed0)
	, snoise(p4 + vec4(-3.3, -2.0, 1.7, -3.0 ), gradDummy, seed0));
	amp *= gain;
	r *= lac;
	}

	return v;
	}

	float terrain_fbm( in vec2 x, in vec4 pOffsets, in float r)
	{
	float lac = 2.0, v = 0.0, amp = 1.0;
	vec2 dsum = vec2(0.0);

	float theta = (x.x) * 2.0 * PI;
	float phi = (x.y) * 2.0 * PI;

	for( int i=0; i<6; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);
	float n = snoise(p4, grad, seed0);
	dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
	v += amp*n/(1.0+dot(dsum,dsum));
	amp *= 0.65;
	r *= lac;
	}
	amp = 0.25;
	for( int i=0; i<5; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);
	float n = snoise(p4, grad, seed0)*0.15;
	dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
	v += amp*n;
	amp *= 0.25;
	r *= lac;
	}

	return v;
	}

	float terrain_ridged( in vec2 x, in vec4 pOffsets, in float r)
	{
	float lac = 2.0, v = 0.0, amp = 1.0;
	vec2 dsum = vec2(0.0);

	float theta = (x.x) * 2.0 * PI;
	float phi = (x.y) * 2.0 * PI;

	for( int i=0; i<5; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);
	float n = (2.0 * pow(1.0 - abs(snoise(p4, grad, seed0)), 4.0)-1.0);
	dsum += grad4to2(grad, r, theta, phi)*(0.35/r);
	v += amp*n/(1.0+dot(dsum,dsum));
	amp *= 0.5;
	r *= lac;
	}
	r = 8.0; amp = 0.20;
	for( int i=0; i<5; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);
	float n = snoise(p4, grad, seed0);
	v += amp*n;
	amp *= 0.7;
	r *= lac;
	}

	return v/1.5;
	}

	float terrain_billow( in vec2 x, in vec4 pOffsets, in float r)
	{
	//x *= rot;
	float lac = 1.8, v = -1.0, amp = 1.0;
	vec2 dsum = vec2(0.0);

	float theta = (x.x) * 2.0 * PI;
	float phi = (x.y) * 2.0 * PI;

	for( int i=0; i<4; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);

	float n = abs(snoise(p4, grad, seed0)) * 1.5;
	dsum += grad4to2(grad, r, theta, phi)*(0.125/r);
	v += amp*n/(1.0+dot(dsum,dsum));
	amp *= 0.6;
	r *= lac;

	}
	v = clamp(v, 0.0, v);
	r = 4.0; amp = 0.6;
	for( int i=0; i<6; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);

	float n = snoise(p4, grad, seed0);
	v += amp*n;
	amp *= 0.65;
	r *= lac;
	}

	return v;
	}

	float terrace_few( in vec2 x, in vec4 pOffsets, in float r)
	{
	float v = 0.0, amp = 1.0, lac = 1.75;
	float bandWidth = 0.4;
	vec2 dsum = vec2(0.0);

	float theta = (x.x) * 2.0 * PI;
	float phi = (x.y) * 2.0 * PI;

	for( int i=0; i<5; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);

	float n = snoise(p4, grad, seed0);
	dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
	v += amp*n/(1.0+dot(dsum,dsum));
	amp *= 0.55;
	r *= lac;
	}
	if(v > -0.66 && v < 0.8)
	{
	float k = floor(v/bandWidth);
	float f = (v - k*bandWidth) / bandWidth;
	float s = min(2.0*f, 1.0);
	float blend = smoothstep(-0.66, -0.33, v) * (1.0 - smoothstep(0.66, 0.8, v));
	v = ((k + s) * bandWidth)blend + (v (1.0 - blend));
	}
	r = 6.0;
	for( int i=0; i<4; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);

	float n = snoise(p4, grad, seed0);
	v += amp*n;
	amp *= 0.5;
	r *= lac;
	}
	return v;
	}

	float terrace_many( in vec2 x, in vec4 pOffsets, in float r)
	{
	float v = 0.0, amp = 1.0, lac = 1.75;
	float bandWidth = 0.25;
	vec2 dsum = vec2(0.0);

	float theta = (x.x) * 2.0 * PI;
	float phi = (x.y) * 2.0 * PI;

	for( int i=0; i<5; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);

	float n = snoise(p4, grad, seed0);

	dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
	v += amp*n/(1.0+dot(dsum,dsum));
	amp *= 0.55;
	r *= lac;
	}
	if(v > -0.5 && v < 0.75)
	{
	float k = floor(v/bandWidth);
	float f = (v - k*bandWidth) / bandWidth;
	float s = min(2.0*f, 1.0);
	float blend = smoothstep(-0.5, -0.33, v) * (1.0 - smoothstep(0.66, 0.75, v));
	v = ((k + s) * bandWidth)blend + (v (1.0 - blend));
	}
	r = 6.0;
	for( int i=0; i<8; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);

	float n = snoise(p4, grad, seed0);

	dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
	v += amp*n;
	amp *= 0.5;
	r *= lac;
	}
	return v;
	}

	float terrace_many_warped( in vec2 x, in vec4 pOffsets, in float r)
	{
	float v = 0.0, amp = 1.0, lac = 1.75;
	float bandWidth = 0.25;
	vec2 dsum = vec2(0.0);
	int w1o = 4;
	//int w2o = 2;
	//int w3o = 1;

	float theta = (x.x) * 2.0 * PI;
	float phi = (x.y) * 2.0 * PI;

	for( int i=0; i<2; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);

	vec4 warp4 = (p4 + basicFbm4(x, r, w1o, lac*0.5, 0.5));
	//vec4 warp4_2 = (warp4 + basicFbm4(x, r, w2o, lac, 0.5));
	//vec4 warp4_3 = (warp4_2 + basicFbm4(x, r, w3o, lac, 0.5));


	float n = snoise(warp4, grad, seed0);

	dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
	v += amp*n/(1.0+dot(dsum,dsum));
	amp *= 0.55;
	r *= lac;
	}
	if(v > -1.0 && v < 0.7)
	{
	float k = floor(v/bandWidth);
	float f = (v - k*bandWidth) / bandWidth;
	float s = min(2.0*f, 1.0);
	float blend = smoothstep(-1.0, -0.6, v) * (1.0 - smoothstep(0.6, 0.7, v));
	v = ((k + s) * bandWidth)blend + (v (1.0 - blend));
	}
	v = clamp(v, -1.0, 0.7);
	r = 10.0; amp = 0.05;
	for( int i=0; i<1; i++ )
	{
	vec4 p4 = pOffsets + vec4(
	(r * cos(theta))
	, (r * sin(theta))
	, (r * cos(phi))
	, (r * sin(phi))
	);
	p4 *= rot0;
	p4 *= rot1;
	vec4 grad = vec4(0.0);

	float n = snoise(p4, grad, seed0);

	dsum += grad4to2(grad, r, theta, phi)*(0.5/r);
	v += amp*n;
	amp *= 0.5;
	r *= lac;
	}
	return v;
	}

	vec3 terrainBlend(in vec2 uv)
	{
	vec4 pOffsets = vec4(
	123.456
	, -432.912
	, -198.023
	, 543.298
	);
	float terrainHeight = (basicFbm(uv, 1.0, 5, 2.2, 0.55, 24.0)); // Just basic fbm for terrain height guide

	float terrainTemp;
	float my = mod(uv.y, 1.0);
	//terrainTemp = (my < 0.5) ? (my - 0.25)4.0 : (my - 0.75)-4.0; // Gradient centred on 0.5

	terrainTemp = (my < 0.5) ? (my2.0) : (2.0-(my2.0)); // Gradient centred on 0.5
	terrainTemp *= shiftRange(basicFbm(uv, 2.0, 6, 1.75, 0.5, 24.0));

	terrainTemp = unshiftRange(terrainTemp);
	terrainTemp -= (0.3 * terrainHeight)smoothstep(-0.15, 0.10, terrainHeight)(1.0 - smoothstep(0.10, 0.3, terrainHeight));
	terrainTemp -= (0.4 * terrainHeight)smoothstep( 0.10, 0.25, terrainHeight)(1.0 - smoothstep(0.45, 0.65, terrainHeight));
	terrainTemp -= (0.6 * terrainHeight)smoothstep( 0.45, 0.65, terrainHeight)(1.0 - smoothstep(0.7, 0.85, terrainHeight));
	terrainTemp -= (0.7 * terrainHeight)*smoothstep( 0.7, 0.85, terrainHeight);
	terrainTemp = shiftRange(terrainTemp);
	terrainTemp = terrainTemp*2.25; // Strengthen gradient around centre
	terrainTemp = clamp(terrainTemp, 0.0, 1.0);

	float terrainMoisture = shiftRange(basicFbm(uv, 2.0, 4, 2.0, 0.5, 12.0));
	terrainMoisture += 1.0 - smoothstep(-1.0, -0.1, terrainHeight);
	terrainMoisture = clamp(terrainMoisture, 0.0, 1.0);

	// Calculate biomes & respective weights
	//int biome = findBiomes(terrainTemp, terrainMoisture);
	float tBlend, mBlend;
	//tBlend = InterpQuinticFunc(terrainTemp - biomeSquare[biome].t0);
	//mBlend = InterpQuinticFunc(terrainMoisture - biomeSquare[biome].t1);
	tBlend = InterpQuinticFunc(terrainTemp);
	mBlend = InterpQuinticFunc(terrainMoisture);

	// Sample four nearest biomes
	float values[4];
	//for(int i = 0; i < 4; i++)
	//{
	// switch(biomeSquare[biome].biomes[i])
	// {
	// case A: values[i] = terrain_billow(uv, pOffsets);
	// case B: values[i] = terrace_many(uv, pOffsets);
	// case C: values[i] = terrain_ridged(uv, pOffsets);
	// case D: values[i] = terrain_fbm(uv, pOffsets);
	// default: values[i] = 0.0;
	// }
	//}

	values[0] = terrain_billow(uv, pOffsets, 3.0) * 0.5;
	values[1] = (terrace_few(uv, pOffsets, 2.5) - 0.25);
	values[2] = terrain_ridged(uv, pOffsets, 3.0);
	values[3] = (terrain_fbm(uv, pOffsets, 2.0) - 0.25) * 0.75;

	// Calculate new height map
	float heightmap = Lerp( Lerp(values[0], values[1], tBlend)
	,Lerp(values[2], values[3], tBlend), mBlend);

	//return vec3(terrainTemp, 0.0, 0.0);
	//return vec3(terrainTemp, terrainMoisture, heightmap);
	//return vec3(tBlend, mBlend, heightmap);
	//return vec3(terrainHeight, terrainHeight, terrainHeight);
	return vec3(heightmap, heightmap, heightmap);


	//return terrainHeight;
	//return terrainTemp;
	//return terrainMoisture;
	//return heightmap;
	//return 0.0;
	}

	void mainImage( out vec4 fragColor, in vec2 fragCoord )
	{
	vec2 uv = (fragCoord)/iResolution.y;
	// Set to 1 to enable cyclic zoom in to show increased detail
	#if 0
	uv.x += sin(iTime0.35)0.5;
	uv.y += cos(iTime0.35)0.5;
	uv = uv((sin(iTime0.25)*0.5 + 0.5) + 0.125);
	#else
	uv *= 1.125;
	//uv *= 0.5;
	uv.y -= 0.06;
	uv.x -= iTime*0.1;
	#endif
	vec4 pOffsets = vec4(
	123.456
	, -432.912
	, -198.023
	, 543.298
	); //pOffsets += iTime * 0.1;

	//float v = basicFbm(uv, 1.0, 7, 1.55, 0.6);
	//float v = terrainBlend(uv);
	//v = shiftRange(v);
	vec3 v= terrainBlend(uv);
	v = shiftRange(v);
	//v.z = 0.0;

	// Output to screen
	if((mod(uv.x, 1.0025) <= 1.0 \|\| mod(uv.x, 1.0025) >= 1.0025)
	&& (mod(uv.y, 1.0025) <= 1.0 \|\| mod(uv.y, 1.0025) >= 1.0025)) fragColor = vec4(v, 1.0);
	else fragColor = vec4(1.0, 0.0, 1.0, 1.0);

	}