Thread: Grid generation from scattered points?

  1. #1
    Registered User
    Join Date
    Mar 2002

    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?

  2. #2
    Registered User Draco's Avatar
    Join Date
    Apr 2002
    did you attempt to change the good functions you found to make them return the array like you want?

  3. #3
    Super Moderator VirtualAce's Avatar
    Join Date
    Aug 2001
    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.

    _LI              proc
    ARG  v1:QWORD,v2:QWORD,interpolant:QWORD
    push   bp
    mov    bp,sp
    fld     [v2]
    fsub  [v1]
    fmul  [interpolant]
    fadd  [v1]
    pop   bp
    _LI             endp

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

    int newpoint_x=LI(x1,x2,interpolant_x);
    int newpoint_y=LI(y1,y2,interpolant_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:

    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
      //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.

  4. #4
    Registered User Codeplug's Avatar
    Join Date
    Mar 2003
    How do you use linear interpolation in 3 dimensions?


  5. #5
    Registered User
    Join Date
    Mar 2002
    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.
    Last edited by canada-paul; 12-12-2003 at 12:58 PM.

  6. #6
    Super Moderator VirtualAce's Avatar
    Join Date
    Aug 2001
    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.

    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)
    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.
    Last edited by VirtualAce; 12-12-2003 at 04:45 PM.

  7. #7
    Registered User
    Join Date
    Mar 2002
    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.

    #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];
    	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");
    	// 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()) {
    		// 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;
    	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;

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Help it won't compile!!!!!
    By esbo in forum C Programming
    Replies: 58
    Last Post: 01-04-2009, 03:22 PM
  2. Replies: 8
    Last Post: 11-03-2008, 09:48 PM
  3. Yahtzee C++ programme help
    By kenneth_888 in forum C++ Programming
    Replies: 13
    Last Post: 09-05-2007, 02:14 PM
  4. CProg Fantasy Football version pi
    By Govtcheez in forum A Brief History of
    Replies: 155
    Last Post: 12-26-2006, 04:30 PM
  5. More Windows Trouble
    By samGwilliam in forum Windows Programming
    Replies: 9
    Last Post: 04-12-2005, 10:02 AM