CGAL users discussion list

[Adam Laža <>]

Hi,

no, I don"t have a grid of input points. I have scattered input points (as vector data) in a computational region, and I need to have an output raster grid with interpolated data for computational region (or convex hull of input points to be specific).

Adam

2015-05-22 10:54 GMT+02:00 Sebastien Loriot (GeometryFactory) <>:

IIUC, you have a grid of input points right?
Did you first checked that the value at the grid point are identical
because it's not clear to me on the pictures.

Sebastien.

On 05/21/2015 06:37 PM, Adam Laža wrote:

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

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss