Grid generation from scattered points?

I've made a program that controls a scanner, and the output file from this program contains the [X,Y] position of the data point and the depth Z. While the region that I'm scanning may be 200 mm by 500 mm in size, I only record the depth values at the most important points. The format of the output file (prior to mesh generation) might look like this (format is X Y Z):

884.3 -0.5 0.005
893.2 3.6 -0.066
898.4 -45.1 -0.965
865 -45.5 -0.838
846.6 -29.4 -1.158
838.8 -20.6 -0.279

After mesh generation, the file would be all filled in like this (in this example I told Surfer to use a 1mm grid spacing):

838.8 -45.5 -0.80678
839.793 -45.5 -0.813933
840.787 -45.5 -0.820559
841.78 -45.5 -0.826533
842.773 -45.5 -0.831739
…and so on …
893.433 3.6 -0.0660163
894.427 3.6 -0.0665912
895.42 3.6 -0.0682444
896.413 3.6 -0.0711039
897.407 3.6 -0.0751363
898.4 3.6 -0.0802335

Currently I'm using Surfer 6.1 (Golden Software) to generate a complete grid file (which can then be exported as a .txt), but I'd really like to move away from using this outside software and do everything within my scanner program. Basically, after I've finished scanning the points I'm interested in, I'd like to call a function that takes those points and returns the completed grid by interpolating all the points that I didn't scan. The method I use now in Surfer for interpolation is the default which is "inverse distance to a power" which seems to work fine, but most of the other methods seem to produce basically the same result.

I've scoured the internet including a number of different newsgroups looking for a simple C function that I could paste onto the end of my scanner program and call it within the main body of my code. Most of what I've found so far has been very complicated and didn't simply return a completed array. My C program runs in the console window, and really, my program so far is very very simple: if, then, else, while, select statements are pretty indicative of how uncomplicated it is.

I would assume that someone, somewhere, at some point has run into the same task and has already programmed a very simple little function that takes in an array of scattered points and returns a completed grid based on some interpolation method, and I just can't find it.

Can anyone help me?

did you attempt to change the good functions you found to make them return the array like you want?

Well this does not sound overly complicated. But you are going to have to select a 'sample size' for your linear interpolations between points. There are an infinite amount of points between point A and point B, so you must select a sample size.

The sample size would be based on how resolute you want the grid to be. The smaller it is the more detailed the grid will be. The larger it is, it will begin to look more like you just drew a straight line from point A to point B.

Here is a linear interpolation in pure assembly using the FPU opcodes.

Code:

---

_LI              proc
ARG  v1:QWORD,v2:QWORD,interpolant:QWORD

push  bp
mov    bp,sp

fld    [v2]
fsub  [v1]
fmul  [interpolant]
fadd  [v1]

pop  bp
ret

_LI            endp

---

Now you simply interpolate from x1,y1 to x2,y2 by your step size.

Code:

---

...
...
int newpoint_x=LI(x1,x2,interpolant_x);
int newpoint_y=LI(y1,y2,interpolant_y);

interpolant_x+=step_size_x;
interpolant_y+=step_size_y;
...
...

---

So if you want to step from x1,y1 to x2,y2 in .1 or 1/10th increments you simply set step_size_x and step_size_y to .1, run the loop, and draw a line from startx,starty to newpoint_x,newpoint_y -> then set startx and starty to newpoint_x and newpoint_y - wash, rinse, repeat.

Instead of drawing lines you could store them in an array. But you are going to have a larger array than you initially started out with. This is what interpolation is all about - you are finding an unknown value at x distance between two known values. Your new grid size will be this:

Code:

---

..
..
int sizex=x2-x1;
int sizey=y2-y1;
int num_points_x=(x2-x1)/step_size_x;
int num_points_y=(y2-y1)/step_size_y;

if ( _LINEAR_ARRAY_)
{
  //Use this value if you want to use a simple linear array
  int totalpoints=num_points_x*num_points_y;
  //Malloc array here
}
else
{
  //Use this for 2D array
  int grid_points[num_points_y][num_points_x];
  //Just an example - this variable will go out of scope so don't
  //do it like this
}

---

Hope this helps.
Given some more time I could code the entire thing for you - but not while just trying to answer this question.

PM me if you need more assistance and I'll be glad to help.

How do you use linear interpolation in 3 dimensions?

gg

Thanks for the tips, Bubba. In your examples, however, I didn't notice a reference to the depth (Z) values, which is what I need to fill in. Determining what the X and Y positions are is fairly straight forward.

As a real, complete (albeit very small) example, here is a series of 3 points sampled in a 5x5 grid (format is X Y Z)

0 0 1
2 2 0
4 4 1

The resulting grid, after interpolation using 'Inverse Distance to a Power' method, is as follows:

0 0 1
1 0 0.83871
2 0 0.545455
3 0 0.459364
4 0 0.5
0 1 0.83871
1 1 0.526316
2 1 0.216867
3 1 0.285714
4 1 0.459364
0 2 0.545455
1 2 0.216867
2 2 0
3 2 0.216867
4 2 0.545455
0 3 0.459364
1 3 0.285714
2 3 0.216867
3 3 0.526316
4 3 0.83871
0 4 0.5
1 4 0.459364
2 4 0.545455
3 4 0.83871
4 4 1

The chosen grid spacing was 1 for simplicity, but could have been anything. Also, the interpolation method could have been anything. But as you can see, the result is simply that the Z values (#'s in the 3rd column) have been interpolated.

I've attached a picture to this post which shows what the resulting surface looks like.

Then interpolate all three coords. I'm assuming your 'grid' is simply a type of heightmap. The z coordinate can be found by interpolating as well.

Code:

---

struct point3D
{
  double x;
  double y;
  double z;
};

inline double LI(double v1,double v2,double f1)
{
  return v1+f1*(v2-v1);
}

void Interpolate_3D(point3D Start, point3D End,point3D NewPoint,double xi,double yi,double zi)
{
  NewPoint.x=LI(Start.x,End.x,xi);
  NewPoint.y=LI(Start.y,End.y,yi);
  NewPoint.z=LI(Start.z,End.z,zi);
}

---

Simply call this inside of a loop that iterates through your known grid points. Your new point is stored in NewPoint.

How you do this also depends on your coordinate system. For instance I normally use x as left and right, y as up and down, and z as into the screen. So for a grid you iterate on x and z at the selected sample size (your step size) and change y to vary the heights. Since X and Z are simply being displaced, you really can just store height values - y is the only value we do not know. X and Z can be found by simply starting at the grid's origin and adding stepsize_x and stepsize_z to x and z respectively on each iteration.

So you could just save height values, which is what most games and simulators do. Saves a lot of space.

Thank you for your help. I managed to come up with something this afternoon which looks to be working fairly well. I used the "Inverse Distance to a Power" method as opposed to straight ahead linear interpolation.

The w[m] array is the "weighted" array which weights the contribution of each known point.

The code isn't really robust at all, but if you feed it the 3 line input file I showed a couple posts up, and tell it to use a grid increment of 1, it works perfectly. I had to use this new fangled (for me, anyways) dynamic memory allocation "new" and "delete" commands, which I'd never used before.

Code:

---

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <malloc.h>
#include <math.h>

using namespace std;

int main( void )
{
        ifstream in_stream;
        ofstream out_stream;
               
        char filename[50];
        char filename_out[50];
        int i=0, num_points, m;
        float X[10000], Y[10000], Z[10000];
        //float* grid;
        float X_min=0, Y_min=0, X_max=0, Y_max=0;
        float X_points, Y_points;
        float gridspace;
        float sum_w, w[3];
               
        out_stream.setf(ios::fixed);
        out_stream.precision(3);
       
        cout << "Assuming input file contains scattered data in the form:\n";
        cout << "X[a]        Y[b]        Z[a,a]\n";
        cout << "X[c]        Y[d]        Z[c,d]\n";
        cout << "X[e]        Y[f]        Z[e,f]...\n";
        cout << "\nOutput file is completed grid using 'Inverse Distance to a Power' to estimate Z values";
        cout << "\nEnter filename to convert: ";
        cin >> filename;
        cout << "\nEnter grid spacing: ";
        cin >> gridspace;

        strcat(filename_out, filename);
        strcat(filename_out, "_grid");
               
        out_stream.open(filename_out);
       
        in_stream.open(filename);
       
        // Read in 1st line
        in_stream >> X[0];
        in_stream >> Y[0];
        in_stream >> Z[0];

        // Initialize min and max values
        X_min = X[0];
        X_max = X[0];
        Y_min = Y[0];
        Y_max = Y[0];
       
        // Read scattered data file
        while (!in_stream.eof()) {
                i++;
                // Read in data
                in_stream >> X[i];
                in_stream >> Y[i];
                in_stream >> Z[i];
                // Set min and max values
                if (X[i] < X_min) X_min = X[i];
                if (X[i] > X_max) X_max = X[i];
                if (Y[i] < Y_min) Y_min = Y[i];
                if (Y[i] > Y_max) Y_max = Y[i]; }

        num_points = i;

        in_stream.close();

        X_points = (X_max - X_min) / gridspace + 1;
        Y_points = (Y_max - Y_min) / gridspace + 1;

        cout << "\nX_points: " << X_points;
        cout << "\nY_points: " << Y_points << "\n\n";

        // Dynamically allocate space for the grid
        float** grid;
        grid = new float*[X_points];
        for(int j = 0 ; j < X_points ; ++j)
                grid[j] = new float[Y_points];
       
        // Initialize the grid
        int x, y, n=0;
        for(x = 0 ; x < X_points ; ++x)
                for(y = 0 ; y < Y_points ; ++y) grid[x][y]=0;       
       
       
        // Assign values within the grid
        for(x = 0 ; x < X_points ; ++x) {
                for(y = 0 ; y < Y_points ; ++y) {
                        if ( (x==X[n]) & (y==Y[n]) ) { grid[x][y] = Z[n]; n++; } // If its a point we already have, don't recalculate it
                        //else grid[x][y] = 8;
                        else {
                                for(m = 0; m <= num_points; ++m) { // Look through all the scattered points
                                        w[m] = pow( pow( (x - X[m]), 2) + pow( (y - Y[m]), 2), -0.5 ); // Find linear distance to each point
                                        w[m] = pow( w[m], 2 );
                                        grid[x][y] += w[m] * Z[m];
                                        sum_w += w[m]; }
                                grid[x][y] /= sum_w; // Calculate the grid value
                                sum_w = 0;}

                        cout << x << "\t" << y << "\t" << grid[x][y] << "\n"; }}
       
        // Free up the memory for the grid                       
        for(j = 0 ; j < X_points ; ++j)
                delete[] grid[j];
        delete[] grid;

}

---