CGAL users discussion list

[Andreas Fabri <>]

Hi Adam,

can you provide us with a small data set, that demonstrates
what you say.  Ideally something like < 10 input points and
a single query point, plus what you expect as result.

Best regards,

Andreas

On 21/05/2015 18:37, 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);

--
Andreas Fabri, PhD
Chief Officer, GeometryFactory
Editor, The CGAL Project

phone: +33.492.954.912    skype: andreas.fabri