CGAL users discussion list

[Adam Laža <>]

Hi all,

I use CGAL::linear_interpolation function in my GRASS GIS modul v.surf.nn that provides natural neighbor interpolation for 2D scattered data. Code works well, I write down the results into an output raster map. But the output differs a lot from another outputs.

For same data I use GRASS GIS modul v.surf.nnbathy (this module uses nn-c library) and ESRI ArcGIS tool Natural neighbor. Output of these tools are identical. But the output of v.surf.nn using CGAL is completely different.

Any idea what's going wrong?

Thanks in advance
Adam

ArcGIS output: http://oi60.tinypic.com/ou17c4.jpg

GRASS GIS (v.surf.nnbathy): http://oi60.tinypic.com/ml4jd5.jpg

GRASS GIS (v.surf.nn using CGAL): http://oi58.tinypic.com/208k9du.jpg

Code:

/* perform NN interpolation */

G_message(_("Computing..."));

/* triangulation */

Delaunay_triangulation T;

typedef CGAL::Data_access< std::map<Point, Coord_type, K::Less_xy_2 > >

Value_access;

T.insert(points.begin(), points.end());

//coordinate computation in grid

double coor_x, coor_y;

coor_x = window.west;

coor_y = window.north;

for (int rows=0 ; rows<window.rows ; rows++) {

  G_percent(rows, window.rows, 2);

  if (mask)

    Rast_get_c_row(maskfd, mask, rows);

    coor_x = window.west;

    for (int cols=0 ; cols<window.cols ; cols++) {

     /* don't interpolate outside of the mask */

     if (mask && mask[cols] == 0) {

       Rast_set_d_null_value(&dcell[cols], 1);

       continue;

     }

 

     K::Point_2 p(coor_x,coor_y);

     std::vector< std::pair< Point, Coord_type > > coords;

     Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p,std::back_inserter(coords)).second;

     Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(), norm, Value_access(function_values));

     G_debug(5, "x: %f y: %f -> res: %f (row=%d; col=%d)",

     coor_x, coor_y, res, rows, cols);

     coor_x += ewres;

     //std::cout << res << " ";

     dcell[cols] = (DCELL) res;

   }

   coor_y -= nsres;

   //std::cout << std::endl;

   Rast_put_d_row(fd, dcell);

 }

G_percent(1, 1, 1);