Folks,

I am working on terrain rendering with elevation data by using SRTM. Does anyone have example program about bicubic/bilinear interpolation? I found some through google and vterrain website but they did not provide example C/C++ program.

Thanks,
Tim

	static float SampleVolumeLinear (const float *volume,
		const int *gridI, const float *gridF, const int *size)
	{
		const int &dimX = size[0];
		const int &dimY = size[1];
							
		int g00 = gridI[0] + gridI[1]*dimX + gridI[2]*(dimY*dimX);
		int g01 = dimX + g00;
		int g10 = dimY*dimX + g00;
		int g11 = dimY*dimX + dimX + g00;
							
		float gridG[3] = {
			1-gridF[0],
			1-gridF[1],
			1-gridF[2]};
		float sample = 0;
		sample += volume[g00] * gridG[0] * gridG[1] * gridG[2];
		sample += volume[g10] * gridG[0] * gridG[1] * gridF[2];
		sample += volume[g01] * gridG[0] * gridF[1] * gridG[2];
		sample += volume[g11] * gridG[0] * gridF[1] * gridF[2];
		sample += volume[g00+1] * gridF[0] * gridG[1] * gridG[2];
		sample += volume[g10+1] * gridF[0] * gridG[1] * gridF[2];
		sample += volume[g01+1] * gridF[0] * gridF[1] * gridG[2];
		sample += volume[g11+1] * gridF[0] * gridF[1] * gridF[2];
							
		return sample;
	}

	static float SampleVolumeCubic (const float *volume,
		const int *gridI, const float *gridF, const int *size)
	{
		auto Cubic = [] (float k0, float k1, float k2, float k3, float t)
		{
			float s = 1.0f - t;

			float g0 = k2 - k0;
			float g1 = k3 - k1;
			float gL = g0 * s + g1 * t;
			float x0 = k1 + (gL + g0) * t * 0.333333f;
			float x1 = k2 - (gL + g1) * s * 0.333333f;

			float tI = t * t * (3.0f - 2.0f*t); // cubic
			//float tI = t * t * t * (t * (t * 6 - 15) + 10); // quintic
			return x0 * (1.0f - tI) + x1 * tI;
		};

		const int &dimX = size[0];
		const int &dimY = size[1];

		float samplesX[16];			
		for (int z=0; z<4; z++)
		for (int y=0; y<4; y++)
		{
			int g = (gridI[0]-1) + 
					(gridI[1]-1+y) * dimX + 
					(gridI[2]-1+z) * (dimY*dimX);
				
			float s = Cubic(volume[g], volume[g+1], volume[g+2], volume[g+3], gridF[0]);
			samplesX[y + z*4] = s;
		}

		float samplesY[4];
		for (int i=0; i<4; i++)
		{
			int g = i*4;	
			float s = Cubic(samplesX[g], samplesX[g+1], samplesX[g+2], samplesX[g+3], gridF[1]);
			samplesY[i] = s;
		}

		return Cubic(samplesY[0], samplesY[1], samplesY[2], samplesY[3], gridF[2]);
	}

Above is both bilinear / bicubic filter lookup, but for volume data. Should be easy to port to 2D.

But not sure if what i did here really is what people mean when saying ‘bicubic’. There should be better resources. I know there are at least some shadertoy examples, but probably that's about optimization and obscures the initial idea.
What i did was to make tangents from neighboring voxel values, and interpolate them with smoothstep to get a curve instead a straight line:

To illustrate, get light blue / red lines from neighbor samples, transport them to sample values (dark blue and red), and interpolate both vectors for the 1D case. It's pretty simple.

Thanks. I now got it. I recently discovered another algorithm for different interpolation method. It called Lanczos algorithm.

I found some demonstration on YouTube about three different methods - bilinear, bicubic, and Lanzcos.

How about Lanzcos sample program?

Thanks,
Tim

I remember ImageMagick: https://imagemagick.org/index.php​

I think it has supports many filters, including Lancosz: https://legacy.imagemagick.org/Usage/filter/

Maybe you could use it to upscale some hightmaps and see which suits your needs best.

@sword7 Just a small warning and note (from years of experience with height data sets taken from satellites), keep in mind that those datasets contain errors (not just noise - it might be completely invalid data - like altitude 65535 meters).

There are some very good filtered data sets like http://www.viewfinderpanoramas.org/​ - but the resolution tends to be lower (3" … but has whole world), I think SRTM can be obtained up to 1" but not for whole world … also its precision is quite poor for most of Europe that is available at this resolution.

As for interpolation of data (and what we've used) - take samples from data, validate (i.e. min and max height - to throw off invalid samples) and push as a point set (at given x, y … height being value). Then we resampled it with a filter into 2D texture at discrete resolution.

The filter searched for N nearest points and used weighted average (weight being calculated based on distance of those points).

Why?

• You will often need area that crosses borders between tiles (assuming your data set is tiled - which most likely is)
• You may need area covering many tiles
• It is possible you will have part of your data available at high resolution of 1" and part at resolution of 3"
• Possibility of combining multiple data sets into one (thus possibility of fixing the errors)
• Both of these cases will most likely require resampling of the dataset into new height maps

Results may look like:

@joej Despite my inactivity on GameDev.net, I'm still not dead yet! Glad to see you're still active.

@Vilem Otte I now got it. Thanks. I got latest SRTM v3.0 global 1 arcsecond datasets (600 GB+) but it only cover 60 degrees from equator both hemispheres. I am now looking for elevation datasets above 60 degrees (both poles).

Well I now found ArcticDEM datasets by just searching on google. Thanks.