🔒 Locked

Need help with some errors in my bilinear interpolation function

Started by polpoint Sep 28, 2009 at 4:49 PM 11 replies 2.3k views

Original Post

Sep 28, 2009 at 4:49 PM

Quote:
+---0,0---+---1,0---+ |? ? ? ? ?|? ? ? ? ?|... |? ? ? ? ?|? ? ? ? ?|... |? ? A ? ?|? ? B ? ?|... |? ? ? ? ?|? ? ? ? ?|... |? ? ? ? ?|? ? ? ? ?|... +---0,1---+---1,1---+ |? ? ? ? ?|? ? ? ? ?|... |? ? ? ? ?|? ? ? ? ?|... |? ? C ? ?|? ? D ? ?|... |? ? ? ? ?|? ? ? ? ?|... |? ? ? ? ?|? ? ? ? ?|... +---------+---------+ | ... | ... |...

I know the values of A, B, C & D and I want to interpolate between those values to fill in the rest of the ?'s Below is the source code as I currently have it, but I am doing something wrong somewhere and I can't figure it out. In the code below the NW quadrant == 0 NE quadrant == 1 SW quadrant == 2 SE quadrant == 3


	for(int pix = 0; pix < RMAP_X; pix++)
	{
		for(int piy = 0; piy < RMAP_Y; piy++)
		{
			quadrant = -1;
			if(pix < (halfRMAP_X) && piy <(halfRMAP_Y))
			{
				quadrant = 0;
			}
			else if(pix >= (halfRMAP_X) && piy < (halfRMAP_Y))
			{
				quadrant = 1;
			}
			else if(pix < (halfRMAP_X) && piy >= (halfRMAP_Y))
			{
				quadrant = 2;
			}
			else if(pix >= (halfRMAP_X) && piy >= (halfRMAP_Y))
			{
				quadrant = 3;
			}


			if(quadrant == 0)
			{
				//nw quadrant - 
				//pix < halfRMAP_X
				//piy < halfRMAP_Y

				float nwDist = 0;
				float nDist = 0;
				float wDist = 0;
				float centDist = 0;
				nDist = 1 / linearDist(pix,piy,pix + (halfRMAP_X - pix),piy - (halfRMAP_Y + piy));
				wDist = 1 / linearDist(pix,piy,pix - (halfRMAP_X + pix),piy + (halfRMAP_Y - piy));
				nwDist = 1 / linearDist(pix,piy, pix - (halfRMAP_X + pix), piy - (halfRMAP_Y - piy));
				centDist = 1 / linearDist(pix,piy, pix + (halfRMAP_X - pix), piy + (halfRMAP_Y - piy));
				
				float sum = nDist + wDist + nwDist + centDist;

				nDist = nDist / sum;
				wDist = wDist / sum;
				nwDist = nwDist / sum;
				centDist = centDist / sum;
				
				world.wMap[wX][wY].area.aMap[aX][aY].region.rMap[pix][piy].elevation = nw*nwDist + n*nDist + w*wDist + world.wMap[wX][wY].area.aMap[aX][aY].region.baseHt*centDist;

			}
			else if(quadrant == 1)
			{
				//ne quadrant
				//pix >= halfRMAP_X
				//piy < halfRMAP_Y

				float neDist = 0;
				float nDist = 0;
				float eDist = 0;
				float centDist = 0;
				nDist = 1 / linearDist(pix,piy, pix - (pix - halfRMAP_X), piy - (halfRMAP_Y + piy));
				eDist = 1 / linearDist(pix,piy, pix + (halfRMAP_X + (RMAP_X - pix)), piy + (halfRMAP_Y - piy));
				neDist = 1 / linearDist(pix,piy, pix + (halfRMAP_X + (RMAP_X - pix)), piy - (halfRMAP_Y + piy));
				centDist = 1 / linearDist(pix,piy, pix - (pix - halfRMAP_X), piy + (halfRMAP_Y - piy));

				float sum = nDist + eDist + neDist + centDist;

				nDist = nDist / sum;
				eDist = eDist / sum;
				neDist = neDist / sum;
				centDist = centDist / sum;

				world.wMap[wX][wY].area.aMap[aX][aY].region.rMap[pix][piy].elevation = ne*neDist + n*nDist + e*eDist + world.wMap[wX][wY].area.aMap[aX][aY].region.baseHt*centDist;

			}
			else if(quadrant == 2)
			{
				//sw quadrant
				//pix < halfRMAP_X
				//piy >= halfRMAP_Y

				float swDist = 0;
				float sDist = 0;
				float wDist = 0;
				float centDist = 0;

				wDist = 1 / linearDist(pix,piy,pix - (halfRMAP_X + pix), piy - (piy - halfRMAP_Y));
				sDist = 1 / linearDist(pix,piy, pix + (halfRMAP_X - pix), piy + (halfRMAP_Y + (RMAP_Y - piy)));
				swDist = 1 / linearDist(pix,piy, pix - (halfRMAP_X + pix), piy + (halfRMAP_Y + (RMAP_Y - piy)));
				centDist = 1 / linearDist(pix,piy, pix + (halfRMAP_X - pix), piy - (piy - halfRMAP_Y));

				float sum = sDist + wDist + swDist + centDist;

				sDist = sDist / sum;
				wDist = wDist / sum;
				swDist = swDist / sum;
				centDist = centDist / sum;


				world.wMap[wX][wY].area.aMap[aX][aY].region.rMap[pix][piy].elevation = sw*swDist + s*sDist + w*wDist + world.wMap[wX][wY].area.aMap[aX][aY].region.baseHt*centDist;

			}
			else if(quadrant == 3)
			{
				//se quadrant
				//pix >= halfRMAP_X
				//piy >= halfRMAP_Y

				float seDist = 0;
				float sDist = 0;
				float eDist = 0;
				float centDist = 0;

				eDist = 1 / linearDist(pix,piy, pix + (halfRMAP_X + (RMAP_X - pix)), piy - (piy - halfRMAP_Y));
				sDist = 1 / linearDist(pix,piy, pix - (pix - halfRMAP_X), piy + (halfRMAP_Y + (RMAP_Y - piy)));
				seDist = 1 / linearDist(pix,piy, pix + (halfRMAP_X + (RMAP_X - pix)), piy + (halfRMAP_Y + (RMAP_Y - piy)));
				centDist = 1 / linearDist(pix,piy, pix - (pix - halfRMAP_X), piy - (piy - halfRMAP_Y));

				float sum = sDist + eDist + seDist + centDist;

				sDist = sDist / sum;
				eDist = eDist / sum;
				seDist = seDist / sum;
				centDist = centDist / sum;



				world.wMap[wX][wY].area.aMap[aX][aY].region.rMap[pix][piy].elevation = se*seDist + s*sDist + e*eDist + world.wMap[wX][wY].area.aMap[aX][aY].region.baseHt*centDist;

			}
		}
	}

My results are below:

as you can see, things don't line up correctly and some other stuff isn't working out correctly. It looks like it is something with the se quadrant, because it seems from the shore in the 2nd picture that the sw & ne quadrants will line up, but I'm not sure. [Edited by - polpoint on September 28, 2009 7:08:48 PM]

polpoint

Sep 28, 2009 at 4:56 PM

The code seems to be interpolating correctly within each quadrant, but the 4 quadrants don't line up correctly. And there are those little low elevation squares (the little blue squares evenly spaced around) that mark the exact center of each square.

if I take into account all 8 of the surrounding elevation levels (as opposed to how I am doing it now and only taking into account the three closest) will that make everything smooth inbetween the 4 quadrants?

iMalc

Sep 29, 2009 at 1:37 AM

First thing you need to know is that your diagram and description of the interpolation is incorrect (apart from the missing [code][/code] tags). Bilinear interpolation only provides the values from between the samples. Also, one strip is usually left for calculation by the interpolation of the next "source pixel" in the row or column (to save duplication). Usually a diagram would represent the values takes the values from the "source pixel" corners. Here is your diagram modified to put X's in the values that are not calculated from an bilerp of A, B, C & D. Only the ?'s are calculated via interpolation of A, B, C & D.

+---0,0---+---1,0---+-|X X X X X|X X X X X|...|X X X X X|X X X X X|...|X X A ? ?|? ? B X X|...|X X ? ? ?|? ? X X X|...|X X ? ? ?|? ? X X X|...+---0,1---+---1,1---+-|X X ? ? ?|? ? X X X|...|X X ? ? ?|? ? X X X|...|X X C X X|X X D X X|...|X X X X X|X X X X X|...|X X X X X|X X X X X|...+---------+---------+-| ....... | ....... |

Now as you can see, the actual amount of samples taken from the one interpolation actually only covers an area the size of one logical "source pixel". Thus the more typical diagram would be:

+---0,0---+---1,0|A ? ? ? ?|B ...|? ? ? ? ?|X ...|? ? ? ? ?|X ...|? ? ? ? ?|X ...|? ? ? ? ?|X ...+---0,1---+---1,1|C X X X X|D ...| ....... | ....

Note now how you sample each "source pixel" from the top-left corner and use those four values to fill in an area corresponding to exactly one "source pixel". You may also notice that this also means you actually need one more row and one more column of source data than you have already. (Unless you scale your sampling across the entire row etc such that you sample each "source pixel" from a slightly different spot)

Sure it can be done without an extra sample by using extrapolation around the edges, but that's actually a lot more work and certainly not a good idea if you're already running into difficulties.

In light of this I think you unfortunately need to totally rewrite what you have. Having written plenty of bilerps myself, I don't even recognise this code as a bilerp of any kind. The whole thing with the quadrants has no place here IMHO.

Note that I use the terms "source pixel" and "destination pixel" very losely here. It only refers to the sampling frequencies of the source and desintation.

I'm trying to look up a useful source of the algorithm for you to start from...
Edit: I found this which might help.
Here's another link that might be useful.

[Edited by - iMalc on September 29, 2009 1:37:30 PM]

polpoint

Sep 30, 2009 at 12:22 AM

Thanks for the thorough description iMalc!

I don't think I explained the specifics of my problem very well or else I am not completely understanding your explanation.

I am not doing anything with graphics, I let openGL handle all of the resampling of my graphics. I am trying to smooth out the elevation of my heightmap after zooming in to a different "magnification" (but it seems that you understand that with the very loose usage of source and destination pixel)

I start out with this:
+-+-+-+
|5|1|4|
+-+-+-+
|9|5|2|
+-+-+-+
|3|8|7|
+-+-+-+

Then increase my magnification to this:
+---+---+---+
|555|111|444|
|555|111|444|
|555|111|444|
+---+---+---+
|999|555|222|
|999|555|222|
|999|555|222|
+---+---+---+
|333|888|777|
|333|888|777|
|333|888|777|
+---+---+---+

At this point I only have the center 5-value square loaded in memory for the new magnification and all of the squares at the previous magnification loaded in memory, so I basically know what the "base" surrounding values are supposed to be.

Then I want to smooth out the values so I get something close to this (the values in the center square are a close enough approximation for this example):
+---+---+---+
|555|111|444|
|555|111|444|
|555|111|444|
+---+---+---+
|999|635|222|
|999|754|222|
|999|566|222|
+---+---+---+
|333|888|777|
|333|888|777|
|333|888|777|
+---+---+---+

So I am coming close to figuring it out, but I don't have access to the extra row & column and it seems like, from your example, I need a specific "start position," which would be the upper-left-most square (or whichever corner). I am using this to smooth out the elevation of a heightmap for a huge area that is divided up into smaller chunks.

So the code that I posted above (minus the stupid method I utilized to figure the coordinates for the surrounding center squares) works, but only within one quadrant and the surrounding square's bordering quadrants. I edited one of the pictures above to highlight this.

The arrows point to the center of the big squares (the values of that specific tile are messed up because I divided by 0 distance from the center) and the black square is outlining the se, sw, ne, and nw quadrant of 4 adjoining big squares.

As you can see everything is lining up within the black square, but it does not line up with that same big square's other 3 quadrants. I tried doing the same thing utilizing all 8 surrounding big square's values along with the central big square that I am modifying but that doesn't work either.

Ah! problems that I can not solve are frustrating!

iMalc

Sep 30, 2009 at 3:26 AM

Good, more diagrams, and you got the code tags - excellent![smile]
What I thought you might not fully understand is that bilinear interpolation of images is exactly the same process as it is for generating smoothed map data.
This is known as point sampling:

+---+---+---+|555|111|444||555|111|444||555|111|444|+---+---+---+|999|555|222||999|555|222||999|555|222|+---+---+---+|333|888|777||333|888|777||333|888|777|+---+---+---+

All well and good.

Here is the sample as you probably imagine it, without any of the interpolated pixels calculated:

+---+---+---+|???|???|???||?5?|?1?|?4?||???|???|???|+---+---+---+|???|???|???||?9?|?5?|?2?||???|???|???|+---+---+---+|???|???|???||?3?|?8?|?7?||???|???|???|+---+---+---+

Now, if you expect those grid points with numbers in them to be exactly as shown with the interpolation, how would you plan to calculate the numbers in say the position of the top-left question mark? It's not located between the 5 and the 1 horizontally, nor is it located between the 5 and the 9 vertically. You'd need to calculate it via extrapolation if this is how you imagined it. This is the difficult way, and it really doesn't work nicely for what you want.

On the other hand you can start out like:

+---+---+---+|5??|?1?|??4||???|???|???||???|???|???|+---+---+---+|???|???|???||9??|?5?|??2||???|???|???|+---+---+---+|???|???|???||???|???|???||3??|?8?|??7|+---+---+---+

What I've done here is stretch the source image out by almost one "source pixel". The numbers at the positions shown are the final values. Now you can clearly calculate the values for all of the question marks because they are all located between known values. For example the first three question marks on the top row would be 4, 3, and 2, correct?
Okay, if this is how you want to do it it's easy:
Iterate over every "destination pixel" (ahem grid tile or whatever) and calculate the exact fractional position that that position refers to on the source "image". when performing lookups. For the question mark in the sixth row and the second column it's source coordinates are 1.25 and 0.25. You use the whole numbers to look up the source pixels, and the whole numbers plus one. So that's source positions (1, 0) (1, 1) (2, 0) and (2, 1).
Now you perform the interpolation using the fractional parts. In this case it's 0.25 of the lower two samples and 0.75 of the upper two. Ditto for the horizontal interpolation.
("1"*0.75 + "4"*0.25)*0.75 + ("5"*0.75 + "2"*0.25)*0.25 = Bingo!

Or you can start with one extra row and column of data like this:

+---+---+---+|5??|1??|4??|8|???|???|???|X|???|???|???|X+---+---+---+-|9??|5??|2??|2|???|???|???|X|???|???|???|X+---+---+---+-|3??|8??|7??|6|???|???|???|X|???|???|???|X+---+---+---+-|0XX|3XX|2XX|5

Again you can calculate the value for every question mark square because it is actually between known values. This method puts each nown value in the same place and makes the interpolation potentially faster, but requires an extra row and column of data.
For the question mark in the sixth row and the second column it's source coordinates are now 1.667 and 0.333. The integer portions are the same here so you still read the values 1, 4, 5 and 2, but now you take slightly different portions of each to fill in that grid position.
Note that the easy way to get one extra row and column is to make it wrap around for the rightmost column such that it actually reads the value from the leftmost column again. I.e:

+---+---+---+|5??|1??|4??|5|???|???|???|X|???|???|???|X+---+---+---+-|9??|5??|2??|9|???|???|???|X|???|???|???|X+---+---+---+-|3??|8??|7??|3|???|???|???|X|???|???|???|X+---+---+---+-|5XX|1XX|4XX|5

This also makes the resulting map tileable, and it's the method usually used in texture bilinear filtering.
Hmm well that's all I have time for right now. I'll let you digest all that for now as you probably need to adjust your thinking somewhat. I suspect that you're either not sure which one of these methods your using, which has led the the problem, or you're trying the method that would require extrapolation.

Edit: Oh one last thing I just realised, it looks like you might be trying to interpolate the values in one destination pixel from the values of the source pixels above, below, left, and right of the desintation pixel. If so, that is definitely wrong and although it makes some sense to someone who is new to interpolation, it will not work and probably explains why you're getting what you're getting.

iMalc

Sep 30, 2009 at 1:14 PM

I've thought about your actual some more and I really can't tell what you ar actually doing, but it definitely isn't bilinear interpolation.

				nDist = 1 / linearDist(pix,piy,pix + (halfRMAP_X - pix),piy - (halfRMAP_Y + piy));				wDist = 1 / linearDist(pix,piy,pix - (halfRMAP_X + pix),piy + (halfRMAP_Y - piy));				nwDist = 1 / linearDist(pix,piy, pix - (halfRMAP_X + pix), piy - (halfRMAP_Y - piy));				centDist = 1 / linearDist(pix,piy, pix + (halfRMAP_X - pix), piy + (halfRMAP_Y - piy));//All of these can be simplified to just				nDist = 1 / linearDist(pix, piy, halfRMAP_X, -halfRMAP_Y));				wDist = 1 / linearDist(pix, piy, -halfRMAP_X, halfRMAP_Y);				nwDist = 1 / linearDist(pix, piy, -halfRMAP_X, -halfRMAP_Y);				centDist = 1 / linearDist(pix, piy, halfRMAP_X, halfRMAP_Y);

What does linearDist do?
I don't know where n, w, and nw come from. Are those samples form the corner of the map or something?

world.wMap[wX][wY].area.aMap[aX][aY].region.rMap[pix][piy].elevation = sw*swDist + s*sDist + w*wDist + world.wMap[wX][wY].area.aMap[aX][aY].region.baseHt*centDist;

This monstrosity screams to me "breaking encapsulation"!. Use a reference to the part that doesn't change during the loop.
Anyway, one thing remains clear here. Your code does not resemble bilinear interpolation at all. Bilinear interpolation would mean that every "destination pixel" is a weighted average of exactly four source pixels. There is normally no floating point used, no divisions, and nothing to do with quadrants.
If you could show me where the sample data comes from and where it is to go to, I can write the bilerp for you to try out.

polpoint

Sep 30, 2009 at 1:36 PM

I don't really understand -

so if I do something like your diagram to calculate only the central 5-value square:

+---+---+---+|5??|?1?|??4||???|???|???||???|???|???|+---+---+---+|???|???|???||9??|?5?|??2||???|???|???|+---+---+---+|???|???|???||???|???|???||3??|?8?|??7|+---+---+---+

Then if I move one square to the north and calculate the 1-value square like this:

+---+---+---+|7??|?2?|??3||???|???|???||???|???|???|+---+---+---+|???|???|???||5??|?1?|??4||???|???|???|+---+---+---+|???|???|???||???|???|???||9??|?5?|??2|+---+---+---+

will it line up with the calculated values for the 5-value square?

It doesn't seem like it would - it seems like it would be missing one row (the three ?'s that are right between the 5 and the 1)

polpoint

Sep 30, 2009 at 1:49 PM

Yeah, that code is pretty funked up -
the lowercase n, ne, e, se, s, sw, w & nw are the surrounding big square's base elevation values (or the A, B & D in the example where we try to calculate the values for the C square - the A would be n, the B would be ne and the D would be e)

	CORD N, NE, E, SE, S, SW, W, NW;	N.X = halfRMAP_X;	N.Y = -halfRMAP_Y;		NE.X = halfRMAP_X + RMAP_X;	NE.Y = -halfRMAP_Y;	E.X = halfRMAP_X + RMAP_X;	E.Y = halfRMAP_Y;	SE.X = halfRMAP_X + RMAP_X;	SE.Y = halfRMAP_Y + RMAP_Y;	S.X = halfRMAP_X;	S.Y = halfRMAP_Y + RMAP_Y;	SW.X = -halfRMAP_X;	SW.Y = halfRMAP_Y + RMAP_Y;	W.X = -halfRMAP_X;	W.Y = halfRMAP_Y;		NW.X = -halfRMAP_X;	NW.Y = -halfRMAP_Y;			//quadrants go:	//	0  1	//	2  3	for(int pix = 0; pix < RMAP_X; pix++)	{		for(int piy = 0; piy < RMAP_Y; piy++)		{			quadrant = -1;			if(pix < (halfRMAP_X) && piy <(halfRMAP_Y))			{				quadrant = 0;			}			else if(pix >= (halfRMAP_X) && piy < (halfRMAP_Y))			{				quadrant = 1;			}			else if(pix < (halfRMAP_X) && piy >= (halfRMAP_Y))			{				quadrant = 2;			}			else if(pix >= (halfRMAP_X) && piy >= (halfRMAP_Y))			{				quadrant = 3;			}			if(quadrant == 0)			{				float wDist = 1 / linearDist(pix,piy,W.X,W.Y);				float nwDist = 1 / linearDist(pix,piy,NW.X,NW.Y);				float nDist = 1 / linearDist(pix,piy,N.X,N.Y);				float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);				float sum = nDist + wDist + nwDist + centDist;				nDist = nDist / sum;				wDist = wDist / sum;				nwDist = nwDist / sum;				centDist = centDist / sum;				world.wMap[wX][wY].area.aMap[aX][aY].region.rMap[pix][piy].elevation = n*nDist + w*wDist + nw*nwDist + world.wMap[wX][wY].area.aMap[aX][aY].region.baseHt*centDist;							}			else if(quadrant == 1)			{				float nDist = 1 / linearDist(pix,piy,N.X,N.Y);				float neDist = 1 / linearDist(pix,piy,NE.X,NE.Y);				float eDist = 1 / linearDist(pix,piy,E.X,E.Y);				float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);				float sum = nDist + eDist + neDist + centDist;				nDist = nDist / sum;				eDist = eDist / sum;				neDist = neDist / sum;				centDist = centDist / sum;				world.wMap[wX][wY].area.aMap[aX][aY].region.rMap[pix][piy].elevation = n*nDist + e*eDist + ne*neDist + world.wMap[wX][wY].area.aMap[aX][aY].region.baseHt*centDist;							}			else if(quadrant == 2)			{				float sDist = 1 / linearDist(pix,piy,S.X,S.Y);				float swDist = 1 / linearDist(pix,piy,SW.X,SW.Y);				float wDist = 1 / linearDist(pix,piy,W.X,W.Y);					float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);				float sum = sDist + wDist + swDist + centDist;				sDist = sDist / sum;				wDist = wDist / sum;				swDist = swDist / sum;				centDist = centDist / sum;				world.wMap[wX][wY].area.aMap[aX][aY].region.rMap[pix][piy].elevation = s*sDist + w*wDist + sw*swDist + world.wMap[wX][wY].area.aMap[aX][aY].region.baseHt*centDist;							}			else if(quadrant == 3)			{				float eDist = linearDist(pix,piy,E.X,E.Y);				float seDist = linearDist(pix,piy,SE.X,SE.Y);				float sDist =  linearDist(pix,piy,S.X,S.Y);				float centDist = linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);								eDist = 1 / eDist;				seDist = 1 / seDist;				sDist = 1 / sDist;				if(centDist != 0)				{					centDist = 1 / centDist;				}				else				{					eDist = 0;					seDist = 0;					sDist = 0;					centDist = 1;				}													float sum = sDist + eDist + seDist + centDist;				sDist = sDist / sum;				eDist = eDist / sum;				seDist = seDist / sum;				centDist = centDist / sum; 				world.wMap[wX][wY].area.aMap[aX][aY].region.rMap[pix][piy].elevation = s*sDist + e*eDist + se*seDist + world.wMap[wX][wY].area.aMap[aX][aY].region.baseHt*centDist;							}                }         }

float linearDist(float x1,float y1,float x2,float y2) {  float dx = x2 - x1;  float dy = y2 - y1;  dx = dx * dx + dy * dy;  if(dx != 0)  {	return sqrt(dx);  }  else  {	return 0;  }}

polpoint

Sep 30, 2009 at 5:30 PM

I think I am just missing some minor step in my function somewhere -

I need 'e' to only take into account 5 & 2, because at that point, both 4 & 7 and 1 & 8 should have 0% influence on the value

I only have odd sized maps at this magnification so I think it might produce the results I want if I calculate the 8 cardinal directions first and then interpolate inbetween those values to fill in the remainder. I'll have to write that up when I get the chance & test it out.

+---+---+---+|???|???|???||?5?|?1?|?4?||???|???|???|+---+---+---+|???|abc|???||?9?|d5e|?2?||???|fgh|???|+---+---+---+|???|???|???||?3?|?8?|?7?||???|???|???|+---+---+---+

iMalc

Sep 30, 2009 at 11:47 PM

Quote:
Original post by polpoint
I don't really understand -

so if I do something like your diagram to calculate only the central 5-value square:
+---+---+---+|5??|?1?|??4||???|???|???||???|???|???|+---+---+---+|???|???|???||9??|?5?|??2||???|???|???|+---+---+---+|???|???|???||???|???|???||3??|?8?|??7|+---+---+---+
Then if I move one square to the north and calculate the 1-value square like this:
+---+---+---+|7??|?2?|??3||???|???|???||???|???|???|+---+---+---+|???|???|???||5??|?1?|??4||???|???|???|+---+---+---+|???|???|???||???|???|???||9??|?5?|??2|+---+---+---+
will it line up with the calculated values for the 5-value square?

It doesn't seem like it would - it seems like it would be missing one row (the three ?'s that are right between the 5 and the 1)

Ah, that's because for that example I was assuming that the entire source board was 3x3, i.e. you don't scroll that part at all. I wasn't aware that you were not calculating the entire interpolated board at once, beforehand. (excuse the double-negation there)
From the extra code posted it seems that you're not actually interpolating the source data at all, as I suspected. Your code is only sampling from one spot and is creating all the point values in the "big square" based solely on that.
I think I have enough to have a shot at implementing actual bilinear interpolation now so I'll give it a go...

iMalc

Oct 1, 2009 at 12:14 AM

Alright, here we go:

	const float XSCALE = (1.0f/RMAP_X);	const float YSCALE = (1.0f/RMAP_Y);	someTypeHere &sourceArea = world.wMap[wX][wY].area;	someOtherTypeHere &destRegion = sourceArea.aMap[aX][aY].region;	for(int pix = 0; pix < RMAP_X; pix++)	{		for(int piy = 0; piy < RMAP_Y; piy++)		{			int xval = pix + halfRMAP_X, yval = piy + halfRMAP_Y;			int xwhole = aX + xval/RMAP_X, ywhole = aY + yval/RMAP_Y;			float xfrac = (xval%RMAP_X) * XSCALE, yfrac = (yval%RMAP_Y) * YSCALE;			// wrapping			if (xwhole < 1) xwhole += SOME_AREA_XMAX_CONSTANT_HERE-1;			else if (xwhole >= SOME_AREA_XMAX_CONSTANT_HERE) xwhole -= SOME_AREA_XMAX_CONSTANT_HERE-1;			if (ywhole < 1) xwhole += SOME_AREA_YMAX_CONSTANT_HERE-1;			else if (ywhole >= SOME_AREA_YMAX_CONSTANT_HERE) ywhole -= SOME_AREA_YMAX_CONSTANT_HERE-1;			// sampling			float sample1 = sourceArea.aMap[xwhole-1][ywhole-1].region.baseHt;			float sample2 = sourceArea.aMap[xwhole  ][ywhole-1].region.baseHt;			float sample3 = sourceArea.aMap[xwhole-1][ywhole  ].region.baseHt;			float sample4 = sourceArea.aMap[xwhole  ][ywhole  ].region.baseHt;			// interpolating			float xinterp1 = sample1 * xfrac +  sample2 * (1.0f - xfrac);			float xinterp2 = sample3 * xfrac +  sample4 * (1.0f - xfrac);			float yinterp = xinterp1 * yfrac + xinterp2 * (1.0f - yfrac);			destRegion.elevation = yinterp;		}	}

This code wraps the edges of the area map.
What's missing:
You need to fill in someTypeHere and someOtherTypeHere. I cant work out these types from your code.
SOME_AREA_XMAX_CONSTANT_HERE should be replaced with whatever constant you have for the maximum value aX can go up to, similar for SOME_AREA_YMAX_CONSTANT_HERE.
I'm assuming C++, but if it's C then change the references to pointers.
Note that this code is very much written for clarity, not for performance. To optimse it you could start with hoisting some calculations out of the innermost loop. However it's proably already faster than what you had since it has no sqrt.

Do let me know how you get on trying this code out, and any questions that arise from it. It it gives correct looking but not perfectly visually appealing enough output then we could try going for a bicubic interpolation next (which uses 16 samples rather than 4).

polpoint

Oct 2, 2009 at 12:05 AM

I've got to thank you for all your help, iMalc!
This is what I wound up doing (I broke up the function into a couple of different code quotes because it is kind of long):

void smoothifyRegion(int wX, int wY, int aX, int aY){	std::stringstream ddp;	//this figures the elevation from the surrounding areamaps and smoothifies everything	int quadrant = -1;	ACORD cd;	cd.wX = wX;	cd.wY = wY;	cd.aX = aX;	cd.aY = aY;	//n	cd.aY--;        //if aY < 0 then we need to modify wY until the coords are valid!	cd = fixCoordinates(cd);        //n is the elevation of the areamap to the north	int n = world.wMap[cd.wX][cd.wY].area.aMap[cd.aX][cd.aY].region.baseHt;	cd.wX = wX;	cd.wY = wY;	cd.aX = aX;	cd.aY = aY;	//ne	cd.aY--;	cd.aX++;	cd = fixCoordinates(cd);	int ne = world.wMap[cd.wX][cd.wY].area.aMap[cd.aX][cd.aY].region.baseHt;	cd.wX = wX;	cd.wY = wY;	cd.aX = aX;	cd.aY = aY;		//e	cd.aX++;	cd = fixCoordinates(cd);	int e = world.wMap[cd.wX][cd.wY].area.aMap[cd.aX][cd.aY].region.baseHt;	cd.wX = wX;	cd.wY = wY;	cd.aX = aX;	cd.aY = aY;		//se	cd.aX++;	cd.aY++;	cd = fixCoordinates(cd);	int se = world.wMap[cd.wX][cd.wY].area.aMap[cd.aX][cd.aY].region.baseHt;	cd.wX = wX;	cd.wY = wY;	cd.aX = aX;	cd.aY = aY;		//s	cd.aY++;	cd = fixCoordinates(cd);	int s = world.wMap[cd.wX][cd.wY].area.aMap[cd.aX][cd.aY].region.baseHt;	cd.wX = wX;	cd.wY = wY;	cd.aX = aX;	cd.aY = aY;		//sw	cd.aY++;	cd.aX--;	cd = fixCoordinates(cd);	int sw = world.wMap[cd.wX][cd.wY].area.aMap[cd.aX][cd.aY].region.baseHt;	cd.wX = wX;	cd.wY = wY;	cd.aX = aX;	cd.aY = aY;		//w	cd.aX--;	cd = fixCoordinates(cd);	int w = world.wMap[cd.wX][cd.wY].area.aMap[cd.aX][cd.aY].region.baseHt;	cd.wX = wX;	cd.wY = wY;	cd.aX = aX;	cd.aY = aY;		//nw	cd.aX--;	cd.aY--;	cd = fixCoordinates(cd);	int nw = world.wMap[cd.wX][cd.wY].area.aMap[cd.aX][cd.aY].region.baseHt;

	CORD N, NE, E, SE, S, SW, W, NW;        //capital directions represent the coordinates of the center of the surrounding areamaps	N.X = halfRMAP_X;	N.Y = -halfRMAP_Y;		NE.X = halfRMAP_X + RMAP_X;	NE.Y = -halfRMAP_Y;	E.X = halfRMAP_X + RMAP_X;	E.Y = halfRMAP_Y;	SE.X = halfRMAP_X + RMAP_X;	SE.Y = halfRMAP_Y + RMAP_Y;	S.X = halfRMAP_X;	S.Y = halfRMAP_Y + RMAP_Y;	SW.X = -halfRMAP_X;	SW.Y = halfRMAP_Y + RMAP_Y;	W.X = -halfRMAP_X;	W.Y = halfRMAP_Y;		NW.X = -halfRMAP_X;	NW.Y = -halfRMAP_Y;	        //haaa less typing!	AREAMAP &sourceArea = world.wMap[wX][wY].area;	REGIONMAP &destRegion = sourceArea.aMap[aX][aY].region;	//do the center of the square first! 	destRegion.rMap[halfRMAP_X][halfRMAP_Y].elevation = destRegion.baseHt;    	//ok - do the 8 cardinal directions next!		int pix;	int piy;	//n	pix = halfRMAP_X;	for(piy = halfRMAP_Y-1; piy >= 0; piy--)	{	 		float nDist = 1 / linearDist(pix,piy,N.X,N.Y);		float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);		float sum = nDist + centDist;		nDist = nDist / sum;		centDist = centDist / sum;		destRegion.rMap[pix][piy].elevation = n*nDist + destRegion.baseHt*centDist;					}	//ne	for(pix = halfRMAP_X+1, piy = halfRMAP_Y-1; pix < RMAP_X && piy >= 0; pix++, piy--)	{				float neDist = 1 / linearDist(pix,piy,NE.X,NE.Y);				float nDist = 1 / linearDist(pix,piy,N.X,N.Y);				float eDist = 1 / linearDist(pix,piy,E.X,E.Y);				float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);				float sum = neDist + nDist + eDist + centDist;				neDist = neDist / sum;				nDist = nDist / sum;				eDist = eDist / sum;				centDist = centDist / sum;				destRegion.rMap[pix][piy].elevation = ne*neDist + n*nDist + e*eDist + destRegion.baseHt*centDist;							}	//e	piy = halfRMAP_Y;	for(pix = halfRMAP_X+1; pix < RMAP_X; pix++)	{	 		float eDist = 1 / linearDist(pix,piy,E.X,E.Y);		float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);		float sum = eDist + centDist;		eDist = eDist / sum;		centDist = centDist / sum;		destRegion.rMap[pix][piy].elevation = e*eDist + destRegion.baseHt*centDist;					}	//se	for(pix = halfRMAP_X+1, piy = halfRMAP_Y+1; pix < RMAP_X && piy < RMAP_Y; pix++, piy++)	{				float seDist = 1 / linearDist(pix,piy,SE.X,SE.Y);				float eDist = 1 / linearDist(pix,piy,E.X,E.Y);				float sDist = 1 / linearDist(pix,piy,S.X,S.Y);				float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);				float sum = seDist + centDist + sDist + eDist;				seDist = seDist / sum;				eDist = eDist / sum;				sDist = sDist / sum;				centDist = centDist / sum;				destRegion.rMap[pix][piy].elevation = se*seDist + s*sDist + e*eDist + destRegion.baseHt*centDist;							}	//s	pix = halfRMAP_X;	for(piy = halfRMAP_Y+1; piy < RMAP_Y; piy++)	{	 		float sDist = 1 / linearDist(pix,piy,S.X,S.Y);		float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);		float sum = sDist + centDist;		sDist = sDist / sum;		centDist = centDist / sum;		destRegion.rMap[pix][piy].elevation = s*sDist + destRegion.baseHt*centDist;					}	//sw	for(pix = halfRMAP_X-1, piy = halfRMAP_Y+1; pix >= 0 && piy < RMAP_Y; pix--, piy++)	{				float swDist = 1 / linearDist(pix,piy,SW.X,SW.Y);				float sDist = 1 / linearDist(pix,piy,S.X,S.Y);				float wDist = 1 / linearDist(pix,piy,W.X,W.Y);				float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);				float sum = swDist + sDist + wDist + centDist;				swDist = swDist / sum;				sDist = sDist / sum;				wDist = wDist / sum;				centDist = centDist / sum;				destRegion.rMap[pix][piy].elevation = sw*swDist + s*sDist + w*wDist + destRegion.baseHt*centDist;							}	//w	piy = halfRMAP_Y;	for(pix = halfRMAP_X-1; pix >= 0; pix--)	{	 		float wDist = 1 / linearDist(pix,piy,W.X,W.Y);		float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);		float sum = wDist + centDist;		wDist = wDist / sum;		centDist = centDist / sum;		destRegion.rMap[pix][piy].elevation = w*wDist + destRegion.baseHt*centDist;					}	//nw	for(pix = halfRMAP_X-1, piy = halfRMAP_Y-1; pix >= 0 && piy >= 0; pix--, piy--)	{				float nwDist = 1 / linearDist(pix,piy,NW.X,NW.Y);				float wDist = 1 / linearDist(pix,piy,W.X,W.Y);				float nDist = 1 / linearDist(pix,piy,N.X,N.Y);				float centDist = 1 / linearDist(pix,piy,halfRMAP_X,halfRMAP_Y);				float sum = nwDist + wDist + nDist + centDist;				nwDist = nwDist / sum;				wDist = wDist / sum;				nDist = nDist / sum;				centDist = centDist / sum;				destRegion.rMap[pix][piy].elevation = nw*nwDist + w*wDist + n*nDist + destRegion.baseHt*centDist;							}

	//octodrants go	//\7  |	0 /	// \  |	 /	//  \ |	/	//6  \|/  1	//----+----	//5  /|\  2	//  / | \	// /  |  \	/// 4 |3  \	//0!	//starts at halfRMAP_X+2, halfRMAP_Y-2	//	int pipY;	int pipX;	int popp;	int counter;	//0 	//starts at halfRMAP_X+1, halfRMAP_Y-2	pipX = halfRMAP_X + 1;	counter = 1;	for(pipY = halfRMAP_Y - 2; pipY >= 0; pipY--, counter++)	{		for(int popp = 0; popp < counter; popp++)		{			float nDist = 1 / linearDist(pipX + popp, pipY, halfRMAP_X, pipY);			float neDist = 1 / linearDist(pipX + popp, pipY, halfRMAP_X + counter + 1, pipY);			 float sum = nDist + neDist;			 nDist = nDist / sum;			 neDist = neDist / sum;			destRegion.rMap[pipX + popp][pipY].elevation = destRegion.rMap[halfRMAP_X][pipY].elevation * nDist + destRegion.rMap[halfRMAP_X + counter + 1][pipY].elevation * neDist;		}	}	//1	//starts at halfRMAP_X + 2, halfRMAP_Y -1	pipY = halfRMAP_Y-1;	counter = 1;	for(pipX = halfRMAP_X + 2; pipX < RMAP_X; pipX++, counter++)	{		for(int popp = 0; popp < counter; popp++)		{		 float eDist = 1 / linearDist(pipX, pipY - popp, pipX, halfRMAP_Y);		 float neDist = 1 / linearDist(pipX, pipY - popp, pipX, halfRMAP_Y - counter - 1);		 float sum = eDist + neDist;		 eDist = eDist / sum;		 neDist = neDist / sum;		 destRegion.rMap[pipX][pipY - popp].elevation = destRegion.rMap[pipX][halfRMAP_Y].elevation * eDist + destRegion.rMap[pipX][halfRMAP_Y - counter - 1].elevation * neDist;		}	}	//2	//starts at halfRMAP_X + 2, halfRMAP_Y+1	pipY = halfRMAP_Y+1;	counter = 1;	for(pipX = halfRMAP_X + 2; pipX < RMAP_X; pipX++, counter++)	{		for(int popp = 0; popp < counter; popp++)		{			float eDist = 1 / linearDist(pipX, pipY + popp, pipX, halfRMAP_Y);			float seDist = 1 / linearDist(pipX, pipY + popp, pipX, halfRMAP_Y + counter + 1);			float sum = eDist + seDist;			eDist = eDist / sum;			seDist = seDist / sum;			destRegion.rMap[pipX][pipY + popp].elevation = destRegion.rMap[pipX][halfRMAP_Y].elevation * eDist + destRegion.rMap[pipX][halfRMAP_Y + counter + 1].elevation * seDist;		}			}	 	//3!	//starts at halfRMAP_X+1, halfRMAP_Y + 2  	pipX = halfRMAP_X+1;	counter = 1;	for(pipY = halfRMAP_Y + 2; pipY < RMAP_Y; pipY++, counter++)	{		for(int popp = 0; popp < counter; popp++)		{			float sDist = 1 / linearDist(pipX + popp, pipY, halfRMAP_X, pipY);			float seDist = 1 / linearDist(pipX + popp, pipY, halfRMAP_X + counter + 1, pipY);			float sum = sDist + seDist;			sDist = sDist / sum;			seDist = seDist / sum;			destRegion.rMap[pipX + popp][pipY].elevation = destRegion.rMap[halfRMAP_X][pipY].elevation * sDist + destRegion.rMap[halfRMAP_X + counter + 1][pipY].elevation * seDist;		}			}	//4!	//starts at halfRMAP_X-1, halfRMAP_Y + 2 	pipX = halfRMAP_X-1;	counter = 1;	for(pipY = halfRMAP_Y + 2; pipY < RMAP_Y; pipY++, counter++)	{		for(int popp = 0; popp < counter; popp++)		{			float sDist = 1 / linearDist(pipX - popp, pipY, halfRMAP_X, pipY);			float swDist = 1 / linearDist(pipX - popp, pipY, halfRMAP_X - counter - 1, pipY);			float sum = sDist + swDist;			sDist = sDist / sum;			swDist = swDist / sum;			destRegion.rMap[pipX - popp][pipY].elevation = destRegion.rMap[halfRMAP_X][pipY].elevation * sDist + destRegion.rMap[halfRMAP_X - counter - 1][pipY].elevation * swDist;		}			}	//5	//starts at halfRMAP_X-2, halfRMAP_Y-1	pipY = halfRMAP_Y + 1;	counter = 1;	for(pipX = halfRMAP_X - 2; pipX >= 0; pipX--, counter++)	{		for(int popp = 0; popp < counter; popp++)		{			float wDist = 1 / linearDist(pipX, pipY + popp, pipX, halfRMAP_Y);			float swDist = 1 / linearDist(pipX, pipY + popp, pipX, halfRMAP_Y + counter + 1);			 float sum = wDist + swDist;			 wDist = wDist / sum;			 swDist = swDist / sum;			destRegion.rMap[pipX][pipY + popp].elevation = destRegion.rMap[pipX][halfRMAP_Y].elevation * wDist + destRegion.rMap[pipX][halfRMAP_Y + counter + 1].elevation * swDist;		}	}	//6	//starts at halfRMAP_X-2, halfRMAP_Y-1	pipY = halfRMAP_Y - 1;	counter = 1;	for(pipX = halfRMAP_X - 2; pipX >= 0; pipX--, counter++)	{		for(int popp = 0; popp < counter; popp++)		{			float wDist = 1 / linearDist(pipX, pipY - popp, pipX, halfRMAP_Y);			float nwDist = 1 / linearDist(pipX, pipY - popp, pipX, halfRMAP_Y - counter - 1);			 float sum = wDist + nwDist;			 wDist = wDist / sum;			 nwDist = nwDist / sum;			destRegion.rMap[pipX][pipY - popp].elevation = destRegion.rMap[pipX][halfRMAP_Y].elevation * wDist + destRegion.rMap[pipX][halfRMAP_Y - counter - 1].elevation * nwDist;		}	}	//7	//starts at halfRMAP_X-1, halfRMAP_Y-2	pipX = halfRMAP_X - 1;	counter = 1;	for(pipY = halfRMAP_Y - 2; pipY >= 0; pipY--, counter++)	{		for(int popp = 0; popp < counter; popp++)		{			float nDist = 1 / linearDist(pipX - popp, pipY, halfRMAP_X, pipY);			float nwDist = 1 / linearDist(pipX - popp, pipY, halfRMAP_X - counter - 1, pipY);			 float sum = nDist + nwDist;			 nDist = nDist / sum;			 nwDist = nwDist / sum;			destRegion.rMap[pipX - popp][pipY].elevation = destRegion.rMap[halfRMAP_X][pipY].elevation * nDist + destRegion.rMap[halfRMAP_X - counter - 1][pipY].elevation * nwDist;		}	}	world.wMap[wX][wY].area.aMap[aX][aY].region.reflectREle(wX,wY,aX,aY,3);}

and I went from this:

to this:

Which is good enough for my needs, even though it produces little strips every once in a while - since I am going to be modifying this map slightly, the extra processing will hide whatever little artifacts show up.

iMalc

Oct 2, 2009 at 2:07 AM

Ok well that's good to hear that you've got something you're fairly happy with. It may not be what I was thinking you'd end up with, but you've learnt some tricks anyway so you've certainly gotten somewhere.

That's a lot of code though. What I had posted was actually to replace the whole lot. Admitedly the method you're using probably does have the advantages of looking less ... linear. That's why I would have gone for bicubic myself.

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.