GLSL Vertex shader bilinear sampling heightmap

Question

I am creating a geomip-mapped terrain. So far I have it working fairly well. The terrain tessellation near the camera is very high and gets less so the further out the geometry is. The geometry of the terrain essentially follows the camera and samples a heightmap texture based on the position of the vertices. Because the geometry tessellation is very high, you can at times see each pixel in the texture when its sampled. It creates obvious pixel bumps. I figured I might be able to get around this by smoothing the sampling of the heightmap. However I seem to have a weird problem related to some bilinear sampling code. I am rendering the terrain by displacing each vertex according to a heightmap texture. To get the height of a vertex at a given UV coordinate I can use:

vec2 worldToMapSpace( vec2 worldPosition ) {
   return ( worldPosition / worldScale + 0.5 );
}

float getHeight( vec3 worldPosition )
{
        #ifdef USE_HEIGHTFIELD
        vec2 heightUv = worldToMapSpace(worldPosition.xz);
        vec2 tHeightSize = vec2( HEIGHTFIELD_SIZE_WIDTH, HEIGHTFIELD_SIZE_HEIGHT ); //both 512
        vec2 texel = vec2( 1.0 / tHeightSize );
        //float coarseHeight = texture2DBilinear( heightfield, heightUv,  texel, tHeightSize ).r;
        float coarseHeight = texture2D( heightfield, vUv ).r;
        return altitude * coarseHeight + heightOffset;
    #else
        return 0.0;
    #endif
}

Which produces this (notice how you can see each pixel):

Here is a wireframe:

I wanted to make the terrain sampling smoother. So I figured I could use some bilinear sampling instead of the standard texture2D function. So here is my bilinear sampling function:

vec4 texture2DBilinear( sampler2D textureSampler, vec2 uv, vec2 texelSize, vec2 textureSize )
{
    vec4 tl = texture2D(textureSampler, uv);
    vec4 tr = texture2D(textureSampler, uv + vec2( texelSize.x, 0.0 ));
    vec4 bl = texture2D(textureSampler, uv + vec2( 0.0, texelSize.y ));
    vec4 br = texture2D(textureSampler, uv + vec2( texelSize.x, texelSize.y ));
    vec2 f = fract( uv.xy * textureSize ); // get the decimal part
    vec4 tA = mix( tl, tr, f.x );
    vec4 tB = mix( bl, br, f.x );
    return mix( tA, tB, f.y ); 
}

The texelSize is calculated as 1 / heightmap size:

vec2 texel = vec2( 1.0 / tHeightSize );

and textureSize is the width and height of the heightmap. However, when I use this function I get this result:

float coarseHeight = texture2DBilinear( heightfield, heightUv,  texel, tHeightSize ).r;

That now seems worse :( Any ideas what I might be doing wrong? Or how I can get a smoother terrain sampling?

EDIT

Here is a vertical screenshot looking down at the terrain. You can see the layers work fine. Notice however that the outer layers that have less triangulation and look smoother while the ones with higher tessellation show each pixel. Im trying to find a way to smooth out the texture sampling.

Answer 1

I was able to find and implement a technique that uses catmulrom interpolation. Code is below.

// catmull works by specifying 4 control points p0, p1, p2, p3 and a weight. The function is used to calculate a point n between p1 and p2 based
// on the weight. The weight is normalized, so if it's a value of 0 then the return value will be p1 and if its 1 it will return p2. 
float catmullRom( float p0, float p1, float p2, float p3, float weight ) {
    float weight2 = weight * weight;
    return 0.5 * (
        p0 * weight * ( ( 2.0 - weight ) * weight - 1.0 ) +
        p1 * ( weight2 * ( 3.0 * weight - 5.0 ) + 2.0 ) +
        p2 * weight * ( ( 4.0 - 3.0 * weight ) * weight + 1.0 ) +
        p3 * ( weight - 1.0 ) * weight2 );
}

// Performs a horizontal catmulrom operation at a given V value.
float textureCubicU( sampler2D samp, vec2 uv00, float texel, float offsetV, float frac ) {
    return catmullRom(
        texture2DLod( samp, uv00 + vec2( -texel, offsetV ), 0.0 ).r,
        texture2DLod( samp, uv00 + vec2( 0.0, offsetV ), 0.0 ).r,
        texture2DLod( samp, uv00 + vec2( texel, offsetV ), 0.0 ).r,
        texture2DLod( samp, uv00 + vec2( texel * 2.0, offsetV ), 0.0 ).r,
    frac );
}

// Samples a texture using a bicubic sampling algorithm. This essentially queries neighbouring
// pixels to get an average value.
float textureBicubic( sampler2D samp, vec2 uv00, vec2 texel, vec2 frac ) {
    return catmullRom(
        textureCubicU( samp, uv00, texel.x, -texel.y, frac.x ),
        textureCubicU( samp, uv00, texel.x, 0.0, frac.x ),
        textureCubicU( samp, uv00, texel.x, texel.y, frac.x ),
        textureCubicU( samp, uv00, texel.x, texel.y * 2.0, frac.x ),
    frac.y );
}

    // Gets the  UV coordinates based on the world X Z position
    vec2 worldToMapSpace( vec2 worldPosition ) {
        return ( worldPosition / worldScale + 0.5 );
    }


// Gets the height at a location p (world space)
float getHeight( vec3 worldPosition )
{
    #ifdef USE_HEIGHTFIELD

        vec2 heightUv = worldToMapSpace(worldPosition.xz);
        vec2 tHeightSize = vec2( HEIGHTFIELD_WIDTH, HEIGHTFIELD_HEIGHT );

        // If we increase the smoothness factor, the terrain becomes a lot smoother.
        // This is because it has the effect of shrinking the texture size and increaing
        // the texel size. Which means when we do sampling the samples are from farther away - making
        // it smoother. However this means the terrain looks less like the original heightmap and so
        // terrain picking goes a bit off. 
        float smoothness = 1.1;
        tHeightSize /= smoothness;

        // The size of each texel
        vec2 texel = vec2( 1.0 / tHeightSize );

        // Find the top-left texel we need to sample.
        vec2 heightUv00 = ( floor( heightUv * tHeightSize ) ) / tHeightSize;

        // Determine the fraction across the 4-texel quad we need to compute.
        vec2 frac = vec2( heightUv - heightUv00 ) * tHeightSize;

        float coarseHeight = textureBicubic( heightfield, heightUv00, texel, frac );
        return altitude * coarseHeight + heightOffset;
    #else
        return 0.0;
    #endif
}