Skip to Content.
Sympa Menu

cgal-discuss - [cgal-discuss] CGAL::linear_interpolation

Subject: CGAL users discussion list

List archive

[cgal-discuss] CGAL::linear_interpolation


Chronological Thread 
  • From: Adam Laža <>
  • To:
  • Subject: [cgal-discuss] CGAL::linear_interpolation
  • Date: Thu, 21 May 2015 18:37:52 +0200

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);





Archive powered by MHonArc 2.6.18.

Top of Page