CGAL users discussion list

[Adam Laža <>]

Hi Andreas,

here is a small data set:

I created a small region 10x10 with resolution 1m.
north:      10
south:      0
west:       0
east:       10
nsres:      1
ewres:      1
rows:       10
cols:       10
cells:      100

Then I generated 10 points:

x|y|z

6.1492718|1.38340775|91.12609400
5.13390161|2.73835363|12.07604600
8.07245632|0.4525116|89.91016900
4.09978228|4.81470035|91.63207700
5.64333425|5.35229498|5.53118800
4.17097588|2.71770808|3.44402800
1.18859736|6.30061561|36.33730100
3.60237755|4.66005842|16.61935800
9.63803784|0.10443555|99.67707600
2.71315208|0.58555986|58.36077900

Some examples of query points (represented as 1m x 1m cell):
4.5|3.5||30.6767022751186
4.5|1.5||27.4545361650472
6.5|1.5||62.9808584320274
5.5|3.5||27.2998580089716

What I expected:
4.5|3.5||10.5043057804237
4.5|1.5||62.9320720752062

6.5|1.5||84.2197781468109
5.5|3.5||15.7523200829975

Expected results are based on GRASS GIS module v.surf.nnbathy, that uses nn-c library[1]. These results(v.surf.nnbathy) are pretty similar to results from ArcGIS (these results should be correct).

Best regards,

Adam

[1]: https://code.google.com/p/nn-c/

2015-05-28 10:15 GMT+02:00 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

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