I was looking around the internet and couldn't find a perfect algorithm for this particular problem:

Our customer have a set of points and weight data along with each point as can be demonstrated by this image:

Of which, we have a GIS program that could generate a "heightmap" or a sort of terrain data from these points and their weight values but as we have near a thousand points of data and that these will change over time, we would like to create our own tools to auto-generate these heightmaps.

So far, I've tried to calculate the weight for each pixel from its distance to the closest data point with Sqrt((x1 - x2) ^ 2 + (y1 - y2) ^ 2) and applying weight and distance factor to the data point's color to produce the resulting gradient color for that particular pixel:

You can see that there are still problems with certain configuration of data points and the algorithm sometimes produce a rather polygonal image when there is a lot of data points. The ideal result should looks more like an ellipsis and less like a polygon.

Here is one example image from wikipedia article on gradient ascent which demonstrates the result I want:

The gradient ascent algorithm is not of my interest. What I'm interested in; is the algorithm to calculate the original function in that picture in the first place, provided data points with weights.

I've not taken any class in topological maths, but I can do some calculus. I think I may be missing something and am rather lost at what should I type in that Google search box.

How smooth do you want the terrain generated? Do you have negative as well as positive displacement. I assume the weight refers no the are a 'hill' directly occupies but instead it's maximum point.
–
ShuggyCoUkFeb 17 '09 at 16:08

@ShuggyCoUk 1. As smooth as possible without too much complications. 2. All values are positive. 3. Yes, it could interpreted as a maximum point of hills and still be correct.
–
chakritFeb 17 '09 at 16:12

9 Answers
9

The resulting function/spline/other mathematical construct can then be interrogated at the required resolution to supply the height map.

Your interpolation function

Sqrt((x1 - x2) ^ 2 + (y1 - y2) ^ 2)

Is similar to Inverse Distance Weighted methods except you are applying an arbitrary filter and discarding many of the other data points.

Most of these techniques rely on a reasonable number of samples and 'terrain-like' behaviour underpinning the values.

I suggest using the weight as the height sample and trying the simple Shepard's Method in the second link (do not filter any pixels to start with) by taking the proportion of a sample points contribution to the overall height value at an interpolation point you can blend the colours of the samples in those ratios to also colour the point. Use the intensity (roughly speaking the grayscale in simple RGB space) to display the height or add contour lines in black like the example image does.

This should be what I'm looking for... and I am already using ArcGIS.. thanks!
–
chakritFeb 17 '09 at 16:32

Also note that you can render at a lower resolution and do a simple (and quick) bilinear interpolation to produce your final height map, however your tool might do this itself if you supply a low res height map anyway.
–
ShuggyCoUkFeb 17 '09 at 17:00

This problem is not as easy as it looks on the surface. Your problem is that both sides of the border of two regions need to have the same height, which is to say, the height at a given pixel is determined by more than just one nearest neighbor.

If I understand it correctly, you need at least two algorithms (and a third piece of jargon).

You are probably going to want to use a kd-tree to help you find the nearest neighbor. Instead of taking O(n^2), this will bring it down to O(n log(n)) (the added benefit is that your Voronoi region generation phase will be fast enough in development to work on the height calculation phase).

Now that you have a 2-D map indexing each point to its nearest neighbor i, you need to walk across every x,y point on the map and calculate its height.

To do that for a given point x,y, first grab its nearest neighbor i and stick that into a list, then collect all contiguous regions on the Voronoi diagram. An easy way is to use flood fill to find all the points in the region, then look around the border and collect up the other identities.

Using this list of all the nearest neighbors, you now have a shot at interpolating correctly! (See other answers for interpolation schemes).

Just look at the picture on the Voronoi diagram article and it will be obvious that generating a map of nearest neighbors will make a picture like the one shown on that article.
–
Jared UpdikeFeb 17 '09 at 16:35

2

"This problem is not as easy as it looks on the surface" bravo :)
–
ShuggyCoUkFeb 17 '09 at 16:39

@ShuggyCoUk: thanks for sharing the Shepherd's method link. I've been puzzling over this for a few months and just knowing what to Google is a great help!
–
Jared UpdikeFeb 17 '09 at 16:44

You've asked for information on algorithms for 2-D interpolation of irregular data, which is quite a complex area. Since you say you have ArcGIS, I strongly advise you to interpolate automatically in ArcGIS using its built-in features for automatic calculations. I am sure that will be far easier than writing your own interpolation algorithm. I've done some automation of ArcGIS, it's fairly straightforward.

If you do write your own interpolation code - I advise you not to - the first thing is to choose the appropriate algorithm as there are several each with their own pluses and minuses. Here is some advice cribbed from the help for the excellent interpolation tool Surfer (which BTW can also be automated quite easily). There are more algorithms, these are just the ones I've tried.

Kriging is one of the more flexible methods and is useful for gridding almost any type of data set. With most data sets, Kriging with the default linear variogram is quite effective. In general, we would most often recommend this method. Kriging is the default gridding method because it generates a good map for most data sets. For larger data sets, Kriging can be rather slow. Kriging can extrapolate grid values beyond your data's Z range.

Inverse distance weighting is fast but has the tendency to generate "bull's-eye" patterns of concentric contours around the data points. Inverse Distance to a Power does not extrapolate Z values beyond the range of data. A simple inverse distance weighting algorithm is easy to implement but will be slowish.

Triangulation with linear interpolation is fast. When you use small data sets, Triangulation with Linear Interpolation generates distinct triangular faces between data points. Triangulation with Linear Interpolation does not extrapolate Z values beyond the range of data.

Shephard's Method is similar to Inverse Distance to a Power but does not tend to generate "bull's eye" patterns, especially when a smoothing factor is used. Shepard's Method can extrapolate values beyond your data's Z range.

To implement the algorithms: you can try Googling, or follow the links in some of the other answers. There are some open-source GIS packages which include interpolation, so maybe you can extract the algorithms from them if you like potholing through C++. Or this book by David Watson is apparently considered a classic, although it is a tricky read and the sample code is spaghetti Basic!! But, from what I hear, it's the best available. If anyone else on Stack Overflow knows better, please do correct me as I can't believe it either.

Kriging is one of the heavyweight methods for doing this, particularly within the field of GIS. It has several nice mathematical properties - the downside is that it can be slow depending on your variogram.

If you want something simpler, there are many interpolation routines which handle this quite well. If you can get ahold of a copy of Numerical Recipes, Chapter 3 is dedicated to explaining many variants for interpolation, and includes code examples and descriptions of their functional properties.

From what I remember Numerical Recipes in C 2nd Edition only includes a routine for 2D linear interpolation. A bit limiting, you might also want to consider Kriging or inverse distance weighting or one of the other methods suggested in other answers.
–
MarkJFeb 21 '09 at 20:15

Well, it depends on what you want the outcome to be. Ideally you'd take into account every datapoint in the universe for every pixel, but that is processing intensive, and may not be what you really want. However, it might be what you need (magnetic fields, for instance)
–
Adam DavisFeb 17 '09 at 16:17

Beware of iterating over all n points for each of the m pixels in your output map. This is O(n*m), and for a 1000x1000 image and 1000 data points that is a billion operations. This will not scale. Use Shepherd's algorithm to tile the plane or something similar.
–
Jared UpdikeFeb 17 '09 at 16:36

I implemented something like this in Winamp AVS a while ago. It uses a "metaballs" type approach of calculating the inverse squared distance (to avoid the sqrt for speed) from each data point, capping it (e.g. to 1.0), and taking a sum of those distances for each point on the 2D grid. This will give a smoothly varying colour/height map.

If you want to look at the code, its in the "Glowy" preset from my J10 AVS pack.

EDIT: Just looking at it I added some other jazz to make it look prettier, the part that is most important is:

Which takes the sum for the 6 points. Everything else done to the red, green and blue output values is to make it look prettier. 6 points isn't much but bear in mind I was trying to make this run in real-time on a 320x200 grid on a 400MHz machine when it was new (which it does at ~20fps). :)

Replace the red =, green = and blue = ... lines with red = d; etc... to see what I mean. All the prettiness goes away and you are left with a greyscale image of smoothly varying blobs around the data points.

Another edit: I forgot to say "s" is the shared weight for all of the points, changing it for each one gives weights individual to each point, e.g. d1 = 2/(...) and d2 = 1/(...) would give d1 twice as much height at its centre as d2. You might also want to cap the expression at the bottom with something like d1 = 2/max(..., 1.0) to smooth off the tops of the points so that they don't peak at infinity in the middle. :)

Sorry for the messiness of the answer... I thought posting the code example would be good enough but on inspection my code is confusing and difficult to read. :(

Your objects are cones which stick out of the ground. They are all parabolas and the weight tells how far they stick out of the ground. Alternatively, make them all the same height and adjust the "flatness" of the parabola accordingly, so a big weight makes the cone very wide while a low weight makes it sharp. Maybe even both to a certain degree.

I suggest that you implement this and see how it looks.

Next, you need to hang a cloth or rubber sheet over the result. The cloth will stretch by a certain amount and it will generally hang down due to gravity. The cones keep it up.

As long as you are close to the center of a cone, the Z coordinate is just the position on the surface of the cone. As you leave the cone center, gravity starts to pull down and the influence of other cones grows.