[SCM] saga branch, dfsg-clean, updated. 2ad85c3fb4774bc7bec3fc660a9732103295e416
Johan Van de Wauw
johan.vandewauw at gmail.com
Sat Mar 26 00:27:38 UTC 2011
The following commit has been merged in the dfsg-clean branch:
commit 2ad85c3fb4774bc7bec3fc660a9732103295e416
Author: Johan Van de Wauw <johan.vandewauw at gmail.com>
Date: Fri Mar 25 22:35:23 2011 +0100
Reinclude natural neighbour routines
- originally deleted because they relied on triangle.c, but depend now on
qhull
diff --git a/src/modules/grid/grid_gridding/Makefile.am b/src/modules/grid/grid_gridding/Makefile.am
index 6a69dfa..ce33855 100644
--- a/src/modules/grid/grid_gridding/Makefile.am
+++ b/src/modules/grid/grid_gridding/Makefile.am
@@ -12,7 +12,7 @@ UC_DEFS = -D_SAGA_UNICODE
endif
DEF_SAGA = -D_SAGA_LINUX -D_TYPEDEF_BYTE -D_TYPEDEF_WORD
CXX_INCS = -I$(top_srcdir)/src/saga_core
-AM_CXXFLAGS = -fPIC $(CXX_INCS) $(DEP_DEFS) $(DEF_SAGA) $(UC_DEFS) $(DBGFLAGS)
+AM_CXXFLAGS = -fPIC $(CXX_INCS) $(DEP_DEFS) $(DEF_SAGA) $(UC_DEFS) $(DBGFLAGS) -lqhull
AM_LDFLAGS = -fPIC -shared -avoid-version
pkglib_LTLIBRARIES = libgrid_gridding.la
libgrid_gridding_la_SOURCES =\
diff --git a/src/modules/grid/grid_gridding/Makefile.in b/src/modules/grid/grid_gridding/Makefile.in
index 036283f..0479bc0 100644
--- a/src/modules/grid/grid_gridding/Makefile.in
+++ b/src/modules/grid/grid_gridding/Makefile.in
@@ -230,7 +230,7 @@ top_srcdir = @top_srcdir@
@SAGA_UNICODE_TRUE at UC_DEFS = -D_SAGA_UNICODE
DEF_SAGA = -D_SAGA_LINUX -D_TYPEDEF_BYTE -D_TYPEDEF_WORD
CXX_INCS = -I$(top_srcdir)/src/saga_core
-AM_CXXFLAGS = -fPIC $(CXX_INCS) $(DEP_DEFS) $(DEF_SAGA) $(UC_DEFS) $(DBGFLAGS)
+AM_CXXFLAGS = -fPIC $(CXX_INCS) $(DEP_DEFS) $(DEF_SAGA) $(UC_DEFS) $(DBGFLAGS) -lqhull
AM_LDFLAGS = -fPIC -shared -avoid-version
pkglib_LTLIBRARIES = libgrid_gridding.la
libgrid_gridding_la_SOURCES = \
diff --git a/src/modules/grid/grid_gridding/nn/delaunay.c b/src/modules/grid/grid_gridding/nn/delaunay.c
new file mode 100644
index 0000000..f511a63
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/delaunay.c
@@ -0,0 +1,731 @@
+/******************************************************************************
+ *
+ * File: delaunay.c
+ *
+ * Created: 04/08/2000
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: Delaunay triangulation - a wrapper to triangulate()
+ *
+ * Description: None
+ *
+ * Revisions: 10/06/2003 PS: delaunay_build(); delaunay_destroy();
+ * struct delaunay: from now on, only shallow copy of the
+ * input data is contained in struct delaunay. This saves
+ * memory and is consistent with libcsa.
+ *
+ * Modified: Joao Cardoso, 4/2/2003
+ * Adapted for use with Qhull instead of "triangle".
+ *
+ *****************************************************************************/
+
+#define USE_QHULL
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <string.h>
+#include <limits.h>
+#include <float.h>
+#ifdef USE_QHULL
+#include <qhull/qhull_a.h>
+#else
+#include "triangle.h"
+#endif
+#include "istack.h"
+#include "nan.h"
+#include "delaunay.h"
+
+int circle_build(circle* c, point* p0, point* p1, point* p2);
+int circle_contains(circle* c, point* p);
+
+#ifdef USE_QHULL
+static int cw(delaunay *d, triangle *t);
+#endif
+
+#ifndef USE_QHULL
+static void tio_init(struct triangulateio* tio)
+{
+ tio->pointlist = NULL;
+ tio->pointattributelist = NULL;
+ tio->pointmarkerlist = NULL;
+ tio->numberofpoints = 0;
+ tio->numberofpointattributes = 0;
+ tio->trianglelist = NULL;
+ tio->triangleattributelist = NULL;
+ tio->trianglearealist = NULL;
+ tio->neighborlist = NULL;
+ tio->numberoftriangles = 0;
+ tio->numberofcorners = 0;
+ tio->numberoftriangleattributes = 0;
+ tio->segmentlist = 0;
+ tio->segmentmarkerlist = NULL;
+ tio->numberofsegments = 0;
+ tio->holelist = NULL;
+ tio->numberofholes = 0;
+ tio->regionlist = NULL;
+ tio->numberofregions = 0;
+ tio->edgelist = NULL;
+ tio->edgemarkerlist = NULL;
+ tio->normlist = NULL;
+ tio->numberofedges = 0;
+}
+
+static void tio_destroy(struct triangulateio* tio)
+{
+ if (tio->pointlist != NULL)
+ free(tio->pointlist);
+ if (tio->pointattributelist != NULL)
+ free(tio->pointattributelist);
+ if (tio->pointmarkerlist != NULL)
+ free(tio->pointmarkerlist);
+ if (tio->trianglelist != NULL)
+ free(tio->trianglelist);
+ if (tio->triangleattributelist != NULL)
+ free(tio->triangleattributelist);
+ if (tio->trianglearealist != NULL)
+ free(tio->trianglearealist);
+ if (tio->neighborlist != NULL)
+ free(tio->neighborlist);
+ if (tio->segmentlist != NULL)
+ free(tio->segmentlist);
+ if (tio->segmentmarkerlist != NULL)
+ free(tio->segmentmarkerlist);
+ if (tio->holelist != NULL)
+ free(tio->holelist);
+ if (tio->regionlist != NULL)
+ free(tio->regionlist);
+ if (tio->edgelist != NULL)
+ free(tio->edgelist);
+ if (tio->edgemarkerlist != NULL)
+ free(tio->edgemarkerlist);
+ if (tio->normlist != NULL)
+ free(tio->normlist);
+}
+
+static delaunay* delaunay_create()
+{
+ delaunay* d = (delaunay *)malloc(sizeof(delaunay));
+
+ d->npoints = 0;
+ d->points = NULL;
+ d->xmin = DBL_MAX;
+ d->xmax = -DBL_MAX;
+ d->ymin = DBL_MAX;
+ d->ymax = -DBL_MAX;
+ d->ntriangles = 0;
+ d->triangles = NULL;
+ d->circles = NULL;
+ d->neighbours = NULL;
+ d->n_point_triangles = NULL;
+ d->point_triangles = NULL;
+ d->nedges = 0;
+ d->edges = NULL;
+ d->flags = NULL;
+ d->first_id = -1;
+ d->t_in = NULL;
+ d->t_out = NULL;
+
+ return d;
+}
+
+static void tio2delaunay(struct triangulateio* tio_out, delaunay* d)
+{
+ int i, j;
+
+ /*
+ * I assume that all input points appear in tio_out in the same order as
+ * they were written to tio_in. I have seen no exceptions so far, even
+ * if duplicate points were presented. Just in case, let us make a couple
+ * of checks.
+ */
+ assert(tio_out->numberofpoints == d->npoints);
+ assert(tio_out->pointlist[2 * d->npoints - 2] == d->points[d->npoints - 1].x && tio_out->pointlist[2 * d->npoints - 1] == d->points[d->npoints - 1].y);
+
+ for (i = 0, j = 0; i < d->npoints; ++i) {
+ point* p = &d->points[i];
+
+ if (p->x < d->xmin)
+ d->xmin = p->x;
+ if (p->x > d->xmax)
+ d->xmax = p->x;
+ if (p->y < d->ymin)
+ d->ymin = p->y;
+ if (p->y > d->ymax)
+ d->ymax = p->y;
+ }
+ if (nn_verbose) {
+ fprintf(stderr, "input:\n");
+ for (i = 0, j = 0; i < d->npoints; ++i) {
+ point* p = &d->points[i];
+
+ fprintf(stderr, " %d: %15.7g %15.7g %15.7g\n", i, p->x, p->y, p->z);
+ }
+ }
+
+ d->ntriangles = tio_out->numberoftriangles;
+ if (d->ntriangles > 0) {
+ d->triangles = (triangle *)malloc(d->ntriangles * sizeof(triangle));
+ d->neighbours = (triangle_neighbours *)malloc(d->ntriangles * sizeof(triangle_neighbours));
+ d->circles = (circle *)malloc(d->ntriangles * sizeof(circle));
+ d->n_point_triangles = (int *)calloc(d->npoints, sizeof(int));
+ d->point_triangles = (int **)malloc(d->npoints * sizeof(int*));
+ d->flags = (int *)calloc(d->ntriangles, sizeof(int));
+ }
+
+ if (nn_verbose)
+ fprintf(stderr, "triangles:\n");
+ for (i = 0; i < d->ntriangles; ++i) {
+ int offset = i * 3;
+ triangle* t = &d->triangles[i];
+ triangle_neighbours* n = &d->neighbours[i];
+ circle* c = &d->circles[i];
+
+ t->vids[0] = tio_out->trianglelist[offset];
+ t->vids[1] = tio_out->trianglelist[offset + 1];
+ t->vids[2] = tio_out->trianglelist[offset + 2];
+
+ n->tids[0] = tio_out->neighborlist[offset];
+ n->tids[1] = tio_out->neighborlist[offset + 1];
+ n->tids[2] = tio_out->neighborlist[offset + 2];
+
+ circle_build(c, &d->points[t->vids[0]], &d->points[t->vids[1]], &d->points[t->vids[2]]);
+
+ if (nn_verbose)
+ fprintf(stderr, " %d: (%d,%d,%d)\n", i, t->vids[0], t->vids[1], t->vids[2]);
+ }
+
+ for (i = 0; i < d->ntriangles; ++i) {
+ triangle* t = &d->triangles[i];
+
+ for (j = 0; j < 3; ++j)
+ d->n_point_triangles[t->vids[j]]++;
+ }
+ if (d->ntriangles > 0) {
+ for (i = 0; i < d->npoints; ++i) {
+ if (d->n_point_triangles[i] > 0)
+ d->point_triangles[i] = (int *)malloc(d->n_point_triangles[i] * sizeof(int));
+ else
+ d->point_triangles[i] = NULL;
+ d->n_point_triangles[i] = 0;
+ }
+ }
+ for (i = 0; i < d->ntriangles; ++i) {
+ triangle* t = &d->triangles[i];
+
+ for (j = 0; j < 3; ++j) {
+ int vid = t->vids[j];
+
+ d->point_triangles[vid][d->n_point_triangles[vid]] = i;
+ d->n_point_triangles[vid]++;
+ }
+ }
+
+ if (tio_out->edgelist != NULL) {
+ d->nedges = tio_out->numberofedges;
+ d->edges = (int *)malloc(d->nedges * 2 * sizeof(int));
+ memcpy(d->edges, tio_out->edgelist, d->nedges * 2 * sizeof(int));
+ }
+}
+#endif
+
+/* Builds Delaunay triangulation of the given array of points.
+ *
+ * @param np Number of points
+ * @param points Array of points [np] (input)
+ * @param ns Number of forced segments
+ * @param segments Array of (forced) segment endpoint indices [2*ns]
+ * @param nh Number of holes
+ * @param holes Array of hole (x,y) coordinates [2*nh]
+ * @return Delaunay triangulation structure with triangulation results
+ */
+delaunay* delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[])
+#ifndef USE_QHULL
+{
+ delaunay* d = delaunay_create();
+ struct triangulateio tio_in;
+ struct triangulateio tio_out;
+ char cmd[64] = "eznC";
+ int i, j;
+
+ assert(sizeof(REAL) == sizeof(double));
+
+ tio_init(&tio_in);
+
+ if (np == 0) {
+ free(d);
+ return NULL;
+ }
+
+ tio_in.pointlist = (double *)malloc(np * 2 * sizeof(double));
+ tio_in.numberofpoints = np;
+ for (i = 0, j = 0; i < np; ++i) {
+ tio_in.pointlist[j++] = points[i].x;
+ tio_in.pointlist[j++] = points[i].y;
+ }
+
+ if (ns > 0) {
+ tio_in.segmentlist = (int *)malloc(ns * 2 * sizeof(int));
+ tio_in.numberofsegments = ns;
+ memcpy(tio_in.segmentlist, segments, ns * 2 * sizeof(int));
+ }
+
+ if (nh > 0) {
+ tio_in.holelist = (double *)malloc(nh * 2 * sizeof(double));
+ tio_in.numberofholes = nh;
+ memcpy(tio_in.holelist, holes, nh * 2 * sizeof(double));
+ }
+
+ tio_init(&tio_out);
+
+ if (!nn_verbose)
+ strcat(cmd, "Q");
+ else if (nn_verbose > 1)
+ strcat(cmd, "VV");
+ if (ns != 0)
+ strcat(cmd, "p");
+
+ if (nn_verbose)
+ fflush(stderr);
+
+ /*
+ * climax
+ */
+ triangulate(cmd, &tio_in, &tio_out, NULL);
+
+ if (nn_verbose)
+ fflush(stderr);
+
+ d->npoints = np;
+ d->points = points;
+
+ tio2delaunay(&tio_out, d);
+
+ tio_destroy(&tio_in);
+ tio_destroy(&tio_out);
+
+ return d;
+}
+#else /* USE_QHULL */
+{
+ delaunay* d = (delaunay *)malloc(sizeof(delaunay));
+
+ coordT *qpoints; /* array of coordinates for each point */
+ boolT ismalloc = False; /* True if qhull should free points */
+ char flags[64] = "qhull d Qbb Qt"; /* option flags for qhull */
+ facetT *facet,*neighbor,**neighborp; /* variables to walk through facets */
+ vertexT *vertex, **vertexp; /* variables to walk through vertex */
+
+ int curlong, totlong; /* memory remaining after qh_memfreeshort */
+ FILE *outfile = stdout;
+ FILE *errfile = stderr; /* error messages from qhull code */
+
+ int i, j;
+ int exitcode;
+ int dim, ntriangles;
+ int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars;
+
+ dim = 2;
+
+ assert(sizeof(realT) == sizeof(double)); /* Qhull was compiled with doubles? */
+
+ if (np == 0 || ns > 0 || nh > 0) {
+ fprintf(stderr, "segments=%d holes=%d\n, aborting Qhull implementation, use 'triangle' instead.\n", ns, nh);
+ free(d);
+ return NULL;
+ }
+
+ qpoints = (coordT *) malloc(np * (dim+1) * sizeof(coordT));
+
+ for (i=0; i<np; i++) {
+ qpoints[i*dim] = points[i].x;
+ qpoints[i*dim+1] = points[i].y;
+ }
+
+ if (!nn_verbose)
+ outfile = NULL;
+ if (nn_verbose)
+ strcat(flags, " s");
+ if (nn_verbose > 1)
+ strcat(flags, " Ts");
+
+ if (nn_verbose)
+ fflush(stderr);
+
+ /*
+ * climax
+ */
+
+ exitcode = qh_new_qhull (dim, np, qpoints, ismalloc,
+ flags, outfile, errfile);
+
+ if(!exitcode) {
+
+ if (nn_verbose)
+ fflush(stderr);
+
+ d->xmin = DBL_MAX;
+ d->xmax = -DBL_MAX;
+ d->ymin = DBL_MAX;
+ d->ymax = -DBL_MAX;
+
+ d->npoints = np;
+ d->points = malloc(np * sizeof(point));
+ for (i = 0; i < np; ++i) {
+ point* p = &d->points[i];
+
+ p->x = points[i].x;
+ p->y = points[i].y;
+ p->z = points[i].z;
+
+ if (p->x < d->xmin)
+ d->xmin = p->x;
+ if (p->x > d->xmax)
+ d->xmax = p->x;
+ if (p->y < d->ymin)
+ d->ymin = p->y;
+ if (p->y > d->ymax)
+ d->ymax = p->y;
+ }
+
+ if (nn_verbose) {
+ fprintf(stderr, "input:\n");
+ for (i = 0; i < np; ++i) {
+ point* p = &d->points[i];
+
+ fprintf(stderr, " %d: %15.7g %15.7g %15.7g\n",
+ i, p->x, p->y, p->z);
+ }
+ }
+
+ qh_findgood_all (qh facet_list);
+ qh_countfacets (qh facet_list, NULL, !qh_ALL, &numfacets,
+ &numsimplicial, &totneighbors, &numridges,
+ &numcoplanars, &numtricoplanars);
+
+ ntriangles = 0;
+ FORALLfacets {
+ if (!facet->upperdelaunay && facet->simplicial)
+ ntriangles++;
+ }
+
+ d->ntriangles = ntriangles;
+ d->triangles = malloc(d->ntriangles * sizeof(triangle));
+ d->neighbours = malloc(d->ntriangles * sizeof(triangle_neighbours));
+ d->circles = malloc(d->ntriangles * sizeof(circle));
+
+ if (nn_verbose)
+ fprintf(stderr, "triangles:\tneighbors:\n");
+
+ i = 0;
+ FORALLfacets {
+ if (!facet->upperdelaunay && facet->simplicial) {
+ triangle* t = &d->triangles[i];
+ triangle_neighbours* n = &d->neighbours[i];
+ circle* c = &d->circles[i];
+
+ j = 0;
+ FOREACHvertex_(facet->vertices)
+ t->vids[j++] = qh_pointid(vertex->point);
+
+ j = 0;
+ FOREACHneighbor_(facet)
+ n->tids[j++] = neighbor->visitid ? neighbor->visitid - 1 : - 1;
+
+ /* Put triangle vertices in counterclockwise order, as
+ * 'triangle' do.
+ * The same needs to be done with the neighbors.
+ *
+ * The following works, i.e., it seems that Qhull maintains a
+ * relationship between the vertices and the neighbors
+ * triangles, but that is not said anywhere, so if this stop
+ * working in a future Qhull release, you know what you have
+ * to do, reorder the neighbors.
+ */
+
+ if(cw(d, t)) {
+ int tmp = t->vids[1];
+ t->vids[1] = t->vids[2];
+ t->vids[2] = tmp;
+
+ tmp = n->tids[1];
+ n->tids[1] = n->tids[2];
+ n->tids[2] = tmp;
+ }
+
+ circle_build(c, &d->points[t->vids[0]], &d->points[t->vids[1]],
+ &d->points[t->vids[2]]);
+
+ if (nn_verbose)
+ fprintf(stderr, " %d: (%d,%d,%d)\t(%d,%d,%d)\n",
+ i, t->vids[0], t->vids[1], t->vids[2], n->tids[0],
+ n->tids[1], n->tids[2]);
+
+ i++;
+ }
+ }
+
+ d->flags = calloc(d->ntriangles, sizeof(int));
+
+ d->n_point_triangles = calloc(d->npoints, sizeof(int));
+ for (i = 0; i < d->ntriangles; ++i) {
+ triangle* t = &d->triangles[i];
+
+ for (j = 0; j < 3; ++j)
+ d->n_point_triangles[t->vids[j]]++;
+ }
+ d->point_triangles = malloc(d->npoints * sizeof(int*));
+ for (i = 0; i < d->npoints; ++i) {
+ if (d->n_point_triangles[i] > 0)
+ d->point_triangles[i] = malloc(d->n_point_triangles[i] * sizeof(int));
+ else
+ d->point_triangles[i] = NULL;
+ d->n_point_triangles[i] = 0;
+ }
+ for (i = 0; i < d->ntriangles; ++i) {
+ triangle* t = &d->triangles[i];
+
+ for (j = 0; j < 3; ++j) {
+ int vid = t->vids[j];
+
+ d->point_triangles[vid][d->n_point_triangles[vid]] = i;
+ d->n_point_triangles[vid]++;
+ }
+ }
+
+ d->nedges = 0;
+ d->edges = NULL;
+
+ d->t_in = NULL;
+ d->t_out = NULL;
+ d->first_id = -1;
+
+ } else {
+ free(d);
+ d = NULL;
+ }
+
+ free(qpoints);
+ qh_freeqhull(!qh_ALL); /* free long memory */
+ qh_memfreeshort (&curlong, &totlong); /* free short memory and memory allocator */
+ if (curlong || totlong)
+ fprintf (errfile,
+ "qhull: did not free %d bytes of long memory (%d pieces)\n",
+ totlong, curlong);
+
+ return d;
+}
+
+ /* returns 1 if a,b,c are clockwise ordered */
+static int cw(delaunay *d, triangle *t)
+{
+ point* pa = &d->points[t->vids[0]];
+ point* pb = &d->points[t->vids[1]];
+ point* pc = &d->points[t->vids[2]];
+
+ return ((pb->x - pa->x)*(pc->y - pa->y) <
+ (pc->x - pa->x)*(pb->y - pa->y));
+}
+
+#endif
+
+/* Releases memory engaged in the Delaunay triangulation structure.
+ *
+ * @param d Structure to be destroyed
+ */
+void delaunay_destroy(delaunay* d)
+{
+ if (d == NULL)
+ return;
+
+ if (d->point_triangles != NULL) {
+ int i;
+
+ for (i = 0; i < d->npoints; ++i)
+ if (d->point_triangles[i] != NULL)
+ free(d->point_triangles[i]);
+ free(d->point_triangles);
+ }
+ if (d->nedges > 0)
+ free(d->edges);
+#ifdef USE_QHULL
+ /* This is a shallow copy if we're not using qhull so we don't
+ * need to free it */
+ if (d->points != NULL)
+ free(d->points);
+#endif
+ if (d->n_point_triangles != NULL)
+ free(d->n_point_triangles);
+ if (d->flags != NULL)
+ free(d->flags);
+ if (d->circles != NULL)
+ free(d->circles);
+ if (d->neighbours != NULL)
+ free(d->neighbours);
+ if (d->triangles != NULL)
+ free(d->triangles);
+ if (d->t_in != NULL)
+ istack_destroy(d->t_in);
+ if (d->t_out != NULL)
+ istack_destroy(d->t_out);
+ free(d);
+}
+
+/* Returns whether the point p is on the right side of the vector (p0, p1).
+ */
+static int on_right_side(point* p, point* p0, point* p1)
+{
+ return (p1->x - p->x) * (p0->y - p->y) > (p0->x - p->x) * (p1->y - p->y);
+}
+
+/* Finds triangle specified point belongs to (if any).
+ *
+ * @param d Delaunay triangulation
+ * @param p Point to be mapped
+ * @param seed Triangle index to start with
+ * @return Triangle id if successful, -1 otherwhile
+ */
+int delaunay_xytoi(delaunay* d, point* p, int id)
+{
+ triangle* t;
+ int i;
+
+ if (p->x < d->xmin || p->x > d->xmax || p->y < d->ymin || p->y > d->ymax)
+ return -1;
+
+ if (id < 0 || id > d->ntriangles)
+ id = 0;
+ t = &d->triangles[id];
+ do {
+ for (i = 0; i < 3; ++i) {
+ int i1 = (i + 1) % 3;
+
+ if (on_right_side(p, &d->points[t->vids[i]], &d->points[t->vids[i1]])) {
+ id = d->neighbours[id].tids[(i + 2) % 3];
+ if (id < 0)
+ return id;
+ t = &d->triangles[id];
+ break;
+ }
+ }
+ } while (i < 3);
+
+ return id;
+}
+
+/* Finds all tricircles specified point belongs to.
+ *
+ * @param d Delaunay triangulation
+ * @param p Point to be mapped
+ * @param n Pointer to the number of tricircles within `d' containing `p'
+ * (output)
+ * @param out Pointer to an array of indices of the corresponding triangles
+ * [n] (output)
+ *
+ * There is a standard search procedure involving search through triangle
+ * neighbours (not through vertex neighbours). It must be a bit faster due to
+ * the smaller number of triangle neighbours (3 per triangle) but can fail
+ * for a point outside convex hall.
+ *
+ * We may wish to modify this procedure in future: first check if the point
+ * is inside the convex hall, and depending on that use one of the two
+ * search algorithms. It not 100% clear though whether this will lead to a
+ * substantial speed gains because of the check on convex hall involved.
+ */
+void delaunay_circles_find(delaunay* d, point* p, int* n, int** out)
+{
+ int i;
+
+ if (d->t_in == NULL) {
+ d->t_in = istack_create();
+ d->t_out = istack_create();
+ }
+
+ /*
+ * It is important to have a reasonable seed here. If the last search
+ * was successful -- start with the last found tricircle, otherwhile (i)
+ * try to find a triangle containing (x,y); if fails then (ii) check
+ * tricircles from the last search; if fails then (iii) make linear
+ * search through all tricircles
+ */
+ if (d->first_id < 0 || !circle_contains(&d->circles[d->first_id], p)) {
+ /*
+ * if any triangle contains (x,y) -- start with this triangle
+ */
+ d->first_id = delaunay_xytoi(d, p, d->first_id);
+
+ /*
+ * if no triangle contains (x,y), there still is a chance that it is
+ * inside some of circumcircles
+ */
+ if (d->first_id < 0) {
+ int nn = d->t_out->n;
+ int tid = -1;
+
+ /*
+ * first check results of the last search
+ */
+ for (i = 0; i < nn; ++i) {
+ tid = d->t_out->v[i];
+ if (circle_contains(&d->circles[tid], p))
+ break;
+ }
+ /*
+ * if unsuccessful, search through all circles
+ */
+ if (tid < 0 || tid == nn) {
+ double nt = d->ntriangles;
+
+ for (tid = 0; tid < nt; ++tid) {
+ if (circle_contains(&d->circles[tid], p))
+ break;
+ }
+ if (tid == nt) {
+ istack_reset(d->t_out);
+ *n = 0;
+ *out = NULL;
+ return; /* failed */
+ }
+ }
+ d->first_id = tid;
+ }
+ }
+
+ istack_reset(d->t_in);
+ istack_reset(d->t_out);
+
+ istack_push(d->t_in, d->first_id);
+ d->flags[d->first_id] = 1;
+
+ /*
+ * main cycle
+ */
+ while (d->t_in->n > 0) {
+ int tid = istack_pop(d->t_in);
+ triangle* t = &d->triangles[tid];
+
+ if (circle_contains(&d->circles[tid], p)) {
+ istack_push(d->t_out, tid);
+ for (i = 0; i < 3; ++i) {
+ int vid = t->vids[i];
+ int nt = d->n_point_triangles[vid];
+ int j;
+
+ for (j = 0; j < nt; ++j) {
+ int ntid = d->point_triangles[vid][j];
+
+ if (d->flags[ntid] == 0) {
+ istack_push(d->t_in, ntid);
+ d->flags[ntid] = 1;
+ }
+ }
+ }
+ }
+ }
+
+ *n = d->t_out->n;
+ *out = d->t_out->v;
+}
diff --git a/src/modules/grid/grid_gridding/nn/delaunay.h b/src/modules/grid/grid_gridding/nn/delaunay.h
new file mode 100644
index 0000000..27c161a
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/delaunay.h
@@ -0,0 +1,76 @@
+/******************************************************************************
+ *
+ * File: delaunay.h
+ *
+ * Created: 04/08/2000
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: Header for delaunay triangulation wrapper
+ *
+ * Description: None
+ *
+ * Revisions: None
+ *
+ *****************************************************************************/
+
+#if !defined(_DELAUNAY_H)
+#define _DELAUNAY_H
+
+#include "nn.h"
+
+typedef struct {
+ int vids[3];
+} triangle;
+
+typedef struct {
+ int tids[3];
+} triangle_neighbours;
+
+typedef struct {
+ double x;
+ double y;
+ double r;
+} circle;
+
+#if !defined(_ISTACK_H)
+struct istack;
+typedef struct istack istack;
+#endif
+
+struct delaunay {
+ int npoints;
+ point* points;
+ double xmin;
+ double xmax;
+ double ymin;
+ double ymax;
+
+ int ntriangles;
+ triangle* triangles;
+ circle* circles;
+ triangle_neighbours* neighbours; /* for delaunay_xytoi() */
+
+ int* n_point_triangles; /* n_point_triangles[i] is number of
+ * triangles i-th point belongs to */
+ int** point_triangles; /* point_triangles[i][j] is index of j-th
+ * triangle i-th point belongs to */
+
+ int nedges;
+ int* edges; /* n-th edge is formed by points[edges[n*2]]
+ * and points[edges[n*2+1]] */
+
+ /*
+ * Work data for delaunay_circles_find(). Placed here for efficiency
+ * reasons. Should be moved to the procedure if parallelizable code
+ * needed.
+ */
+ int* flags;
+ int first_id; /* last search result, used in start up of a
+ * new search */
+ istack* t_in;
+ istack* t_out;
+};
+
+#endif
diff --git a/src/modules/grid/grid_gridding/nn/hash.c b/src/modules/grid/grid_gridding/nn/hash.c
new file mode 100644
index 0000000..bc1aa09
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/hash.c
@@ -0,0 +1,682 @@
+/******************************************************************************
+ *
+ * File: hash.c
+ *
+ * Purpose: Hash table implementation
+ *
+ * Author: Jerry Coffin
+ *
+ * Description: Public domain code by Jerry Coffin, with improvements by
+ * HenkJan Wolthuis.
+ * Date last modified: 05-Jul-1997
+ *
+ * Revisions: 18-09-2002 -- modified by Pavel Sakov
+ *
+ *****************************************************************************/
+
+#include <string.h>
+#include <stdlib.h>
+#include <assert.h>
+#include "hash.h"
+
+#define INT_PER_DOUBLE 2
+
+/** A hash table consists of an array of these buckets.
+ */
+typedef struct ht_bucket {
+ void* key;
+ void* data;
+ int id; /* unique id -- just in case */
+ struct ht_bucket* next;
+} ht_bucket;
+
+/** Hash table structure.
+ * Note that more nodes than `size' can be inserted in the table,
+ * but performance degrades as this happens.
+ */
+struct hashtable {
+ int size; /* table size */
+ int n; /* current number of entries */
+ int naccum; /* number of inserted entries */
+ int nhash; /* number of used table elements */
+ ht_keycp cp;
+ ht_keyeq eq;
+ ht_key2hash hash;
+ ht_bucket** table;
+};
+
+/* Creates a hashtable of specified size.
+ */
+hashtable* ht_create(int size, ht_keycp cp, ht_keyeq eq, ht_key2hash hash)
+{
+ hashtable* table = (hashtable *)malloc(sizeof(hashtable));
+ ht_bucket** bucket;
+ int i;
+
+ assert(sizeof(double) == INT_PER_DOUBLE * sizeof(int));
+ /*
+ * (used in d1hash() and d2hash())
+ */
+
+ if (table == NULL)
+ return NULL;
+
+ if (size <= 0) {
+ free(table);
+ return NULL;
+ }
+
+ table->size = size;
+ table->table = (ht_bucket **)malloc(sizeof(ht_bucket*) * size);
+ bucket = table->table;
+
+ if (bucket == NULL) {
+ free(table);
+ return NULL;
+ }
+
+ for (i = 0; i < size; ++i)
+ bucket[i] = NULL;
+ table->n = 0;
+ table->naccum = 0;
+ table->nhash = 0;
+ table->eq = eq;
+ table->cp = cp;
+ table->hash = hash;
+
+ return table;
+}
+
+/* Destroys a hash table.
+ * (Take care of deallocating data by ht_process() prior to destroying the
+ * table if necessary.)
+ *
+ * @param table Hash table to be destroyed
+ */
+void ht_destroy(hashtable* table)
+{
+ int i;
+
+ if (table == NULL)
+ return;
+
+ for (i = 0; i < table->size; ++i) {
+ ht_bucket* bucket;
+
+ for (bucket = (table->table)[i]; bucket != NULL;) {
+ ht_bucket* prev = bucket;
+
+ free(bucket->key);
+ bucket = bucket->next;
+ free(prev);
+ }
+ }
+
+ free(table->table);
+ free(table);
+}
+
+/* Inserts a new entry into the hash table.
+ *
+ * @param table The hash table
+ * @param key Ponter to entry's key
+ * @param data Pointer to associated data
+ * @return Pointer to the old data associated with the key, NULL if the key
+ * wasn't in the table previously
+ */
+void* ht_insert(hashtable* table, void* key, void* data)
+{
+ unsigned int val = table->hash(key) % table->size;
+ ht_bucket* bucket;
+
+ /*
+ * NULL means this bucket hasn't been used yet. We'll simply allocate
+ * space for our new bucket and put our data there, with the table
+ * pointing at it.
+ */
+ if ((table->table)[val] == NULL) {
+ bucket = (ht_bucket *)malloc(sizeof(ht_bucket));
+ if (bucket == NULL)
+ return NULL;
+
+ bucket->key = table->cp(key);
+ bucket->next = NULL;
+ bucket->data = data;
+ bucket->id = table->naccum;
+
+ (table->table)[val] = bucket;
+ table->n++;
+ table->naccum++;
+ table->nhash++;
+
+ return bucket->data;
+ }
+
+ /*
+ * This spot in the table is already in use. See if the current string
+ * has already been inserted, and if so, return corresponding data.
+ */
+ for (bucket = (table->table)[val]; bucket != NULL; bucket = bucket->next)
+ if (table->eq(key, bucket->key) == 1) {
+ void* old_data = bucket->data;
+
+ bucket->data = data;
+ bucket->id = table->naccum;
+ table->naccum++;
+
+ return old_data;
+ }
+
+ /*
+ * This key must not be in the table yet. We'll add it to the head of
+ * the list at this spot in the hash table. Speed would be slightly
+ * improved if the list was kept sorted instead. In this case, this
+ * code would be moved into the loop above, and the insertion would take
+ * place as soon as it was determined that the present key in the list
+ * was larger than this one.
+ */
+ bucket = (ht_bucket*) malloc(sizeof(ht_bucket));
+ if (bucket == NULL)
+ return 0;
+ bucket->key = table->cp(key);
+ bucket->data = data;
+ bucket->next = (table->table)[val];
+ bucket->id = table->naccum;
+
+ (table->table)[val] = bucket;
+ table->n++;
+ table->naccum++;
+
+ return data;
+}
+
+/* Returns a pointer to the data associated with a key. If the key has
+ * not been inserted in the table, returns NULL.
+ *
+ * @param table The hash table
+ * @param key The key
+ * @return The associated data or NULL
+ */
+void* ht_find(hashtable* table, void* key)
+{
+ unsigned int val = table->hash(key) % table->size;
+ ht_bucket* bucket;
+
+ if ((table->table)[val] == NULL)
+ return NULL;
+
+ for (bucket = (table->table)[val]; bucket != NULL; bucket = bucket->next)
+ if (table->eq(key, bucket->key) == 1)
+ return bucket->data;
+
+ return NULL;
+}
+
+/* Deletes an entry from the table. Returns a pointer to the data that
+ * was associated with the key so that the calling code can dispose it
+ * properly.
+ *
+ * @param table The hash table
+ * @param key The key
+ * @return The associated data or NULL
+ */
+void* ht_delete(hashtable* table, void* key)
+{
+ unsigned int val = table->hash(key) % table->size;
+ ht_bucket* prev;
+ ht_bucket* bucket;
+ void* data;
+
+ if ((table->table)[val] == NULL)
+ return NULL;
+
+ /*
+ * Traverse the list, keeping track of the previous node in the list.
+ * When we find the node to delete, we set the previous node's next
+ * pointer to point to the node after ourself instead. We then delete
+ * the key from the present node, and return a pointer to the data it
+ * contains.
+ */
+ for (prev = NULL, bucket = (table->table)[val]; bucket != NULL; prev = bucket, bucket = bucket->next) {
+ if (table->eq(key, bucket->key) == 1) {
+ data = bucket->data;
+ if (prev != NULL)
+ prev->next = bucket->next;
+ else {
+ /*
+ * If 'prev' still equals NULL, it means that we need to
+ * delete the first node in the list. This simply consists
+ * of putting our own 'next' pointer in the array holding
+ * the head of the list. We then dispose of the current
+ * node as above.
+ */
+ (table->table)[val] = bucket->next;
+ table->nhash--;
+ }
+ free(bucket->key);
+ free(bucket);
+ table->n--;
+
+ return data;
+ }
+ }
+
+ /*
+ * If we get here, it means we didn't find the item in the table. Signal
+ * this by returning NULL.
+ */
+ return NULL;
+}
+
+/* For each entry, calls a specified function with corresponding data as a
+ * parameter.
+ *
+ * @param table The hash table
+ * @param func The action function
+ */
+void ht_process(hashtable* table, void (*func) (void*))
+{
+ int i;
+
+ for (i = 0; i < table->size; ++i)
+ if ((table->table)[i] != NULL) {
+ ht_bucket* bucket;
+
+ for (bucket = (table->table)[i]; bucket != NULL; bucket = bucket->next)
+ func(bucket->data);
+ }
+}
+
+/*
+ * functions for for string keys
+ */
+
+static unsigned int strhash(void* key)
+{
+ char* str = (char*) key;
+ unsigned int hashvalue = 0;
+
+ while (*str != 0) {
+ hashvalue ^= *(unsigned int*) str;
+ hashvalue <<= 1;
+ str++;
+ }
+
+ return hashvalue;
+}
+
+static void* strcp(void* key)
+{
+ return strdup((const char *)key);
+}
+
+static int streq(void* key1, void* key2)
+{
+ return !strcmp((const char *)key1, (const char *)key2);
+}
+
+/* functions for for double keys */
+
+static unsigned int d1hash(void* key)
+{
+ unsigned int* v = (unsigned int*) key;
+
+#if INT_PER_DOUBLE == 2
+ return v[0] + v[1];
+#else
+#error not implemented
+#endif
+}
+
+static void* d1cp(void* key)
+{
+ double* newkey = (double *)malloc(sizeof(double));
+
+ *newkey = *(double*) key;
+
+ return newkey;
+}
+
+int d1eq(void* key1, void* key2)
+{
+ return *(double*) key1 == *(double*) key2;
+}
+
+/*
+ * functions for for double[2] keys
+ */
+
+#include "math.h"
+
+static unsigned int d2hash(void* key)
+{
+ unsigned int* v = (unsigned int*) key;
+
+#if INT_PER_DOUBLE == 2
+ /*
+ * PS: here multiplications suppose to make (a,b) and (b,a) generate
+ * different hash values
+ */
+ return v[0] + v[1] + v[2] * 3 + v[3] * 7;
+#else
+#error not implemented
+#endif
+}
+
+static void* d2cp(void* key)
+{
+ double* newkey = (double *)malloc(sizeof(double) * 2);
+
+ newkey[0] = ((double*) key)[0];
+ newkey[1] = ((double*) key)[1];
+
+ return newkey;
+}
+
+static int d2eq(void* key1, void* key2)
+{
+ return (((double*) key1)[0] == ((double*) key2)[0]) && (((double*) key1)[1] == ((double*) key2)[1]);
+}
+
+hashtable* ht_create_d1(int size)
+{
+ return ht_create(size, d1cp, d1eq, d1hash);
+}
+
+hashtable* ht_create_d2(int size)
+{
+ return ht_create(size, d2cp, d2eq, d2hash);
+}
+
+hashtable* ht_create_str(int size)
+{
+ return ht_create(size, strcp, streq, strhash);
+}
+
+#ifdef HT_TEST
+
+#include <stdio.h>
+#include <limits.h>
+
+#define BUFSIZE 1024
+
+static void print_double(void* data)
+{
+ printf(" \"%d\"", (int)* (double*) data);
+}
+
+static void print_string(void* data)
+{
+ printf(" \"%s\"", (char*) data);
+}
+
+int main()
+{
+ double points[] = {
+ 922803.7855, 7372394.688, 0,
+ 922849.2037, 7372307.027, 1,
+ 922894.657, 7372219.306, 2,
+ 922940.1475, 7372131.528, 3,
+ 922985.6777, 7372043.692, 4,
+ 923031.2501, 7371955.802, 5,
+ 923076.8669, 7371867.857, 6,
+ 923122.5307, 7371779.861, 7,
+ 923168.2439, 7371691.816, 8,
+ 923214.0091, 7371603.722, 9,
+ 923259.8288, 7371515.583, 10,
+ 922891.3958, 7372440.117, 11,
+ 922936.873, 7372352.489, 12,
+ 922982.3839, 7372264.804, 13,
+ 923027.9308, 7372177.064, 14,
+ 923073.5159, 7372089.268, 15,
+ 923119.1415, 7372001.42, 16,
+ 923164.8099, 7371913.521, 17,
+ 923210.5233, 7371825.572, 18,
+ 923256.2841, 7371737.575, 19,
+ 923302.0946, 7371649.534, 20,
+ 923347.9572, 7371561.45, 21,
+ 922978.9747, 7372485.605, 22,
+ 923024.5085, 7372398.009, 23,
+ 923070.0748, 7372310.358, 24,
+ 923115.6759, 7372222.654, 25,
+ 923161.3136, 7372134.897, 26,
+ 923206.9903, 7372047.09, 27,
+ 923252.7079, 7371959.233, 28,
+ 923298.4686, 7371871.33, 29,
+ 923344.2745, 7371783.381, 30,
+ 923390.1279, 7371695.389, 31,
+ 923436.0309, 7371607.357, 32,
+ 923066.5232, 7372531.148, 33,
+ 923112.1115, 7372443.583, 34,
+ 923157.7311, 7372355.966, 35,
+ 923203.3842, 7372268.296, 36,
+ 923249.0725, 7372180.577, 37,
+ 923294.7981, 7372092.808, 38,
+ 923340.5628, 7372004.993, 39,
+ 923386.3686, 7371917.132, 40,
+ 923432.2176, 7371829.229, 41,
+ 923478.1116, 7371741.284, 42,
+ 923524.0527, 7371653.302, 43,
+ 923154.0423, 7372576.746, 44,
+ 923199.6831, 7372489.211, 45,
+ 923245.3541, 7372401.625, 46,
+ 923291.0572, 7372313.989, 47,
+ 923336.7941, 7372226.305, 48,
+ 923382.5667, 7372138.574, 49,
+ 923428.3766, 7372050.798, 50,
+ 923474.2256, 7371962.978, 51,
+ 923520.1155, 7371875.118, 52,
+ 923566.0481, 7371787.218, 53,
+ 923612.0252, 7371699.282, 54,
+ 923241.533, 7372622.396, 55,
+ 923287.2244, 7372534.889, 56,
+ 923332.9449, 7372447.334, 57,
+ 923378.6963, 7372359.731, 58,
+ 923424.4801, 7372272.081, 59,
+ 923470.2979, 7372184.385, 60,
+ 923516.1513, 7372096.646, 61,
+ 923562.0418, 7372008.866, 62,
+ 923607.9709, 7371921.046, 63,
+ 923653.9402, 7371833.188, 64,
+ 923699.9514, 7371745.296, 65,
+ 923328.9962, 7372668.095, 66,
+ 923374.7365, 7372580.617, 67,
+ 923420.5049, 7372493.091, 68,
+ 923466.303, 7372405.519, 69,
+ 923512.1321, 7372317.901, 70,
+ 923557.9936, 7372230.24, 71,
+ 923603.8889, 7372142.536, 72,
+ 923649.8192, 7372054.793, 73,
+ 923695.786, 7371967.011, 74,
+ 923741.7905, 7371879.193, 75,
+ 923787.8341, 7371791.342, 76,
+ 923416.4327, 7372713.844, 77,
+ 923462.2204, 7372626.393, 78,
+ 923508.0353, 7372538.895, 79,
+ 923553.8787, 7372451.353, 80,
+ 923599.7517, 7372363.766, 81,
+ 923645.6555, 7372276.137, 82,
+ 923691.5914, 7372188.467, 83,
+ 923737.5603, 7372100.757, 84,
+ 923783.5634, 7372013.011, 85,
+ 923829.6017, 7371925.231, 86,
+ 923875.6763, 7371837.419, 87,
+ 923503.8433, 7372759.64, 88,
+ 923549.6771, 7372672.214, 89,
+ 923595.5372, 7372584.744, 90,
+ 923641.4246, 7372497.23, 91,
+ 923687.3404, 7372409.673, 92,
+ 923733.2855, 7372322.074, 93,
+ 923779.2608, 7372234.436, 94,
+ 923825.2672, 7372146.759, 95,
+ 923871.3056, 7372059.047, 96,
+ 923917.3766, 7371971.301, 97,
+ 923963.4812, 7371883.524, 98,
+ 923591.2288, 7372805.481, 99,
+ 923637.1076, 7372718.081, 100,
+ 923683.0118, 7372630.638, 101,
+ 923728.9423, 7372543.151, 102,
+ 923774.8998, 7372455.622, 103,
+ 923820.8852, 7372368.052, 104,
+ 923866.8991, 7372280.443, 105,
+ 923912.9422, 7372192.797, 106,
+ 923959.015, 7372105.116, 107,
+ 924005.118, 7372017.402, 108,
+ 924051.2518, 7371929.657, 109,
+ 923678.5898, 7372851.367, 110,
+ 923724.5126, 7372763.992, 111,
+ 923770.46, 7372676.574, 112,
+ 923816.4328, 7372589.113, 113,
+ 923862.4314, 7372501.611, 114,
+ 923908.4564, 7372414.069, 115,
+ 923954.5083, 7372326.488, 116,
+ 924000.5875, 7372238.87, 117,
+ 924046.6941, 7372151.218, 118,
+ 924092.8286, 7372063.533, 119,
+ 924138.9911, 7371975.818, 120
+ };
+
+ int size = sizeof(points) / sizeof(double) / 3;
+ hashtable* ht;
+ int i;
+
+ /*
+ * double[2] key
+ */
+
+ printf("\n1. Testing a table with key of double[2] type\n\n");
+
+ printf(" creating a table...");
+ ht = ht_create_d2(size);
+ printf("done\n");
+
+ printf(" inserting %d values from a file...", size);
+ for (i = 0; i < size; ++i)
+ ht_insert(ht, &points[i * 3], &points[i * 3 + 2]);
+ printf("done\n");
+
+ printf(" stats:\n");
+ printf(" %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash);
+ printf(" %f entries per hash value in use\n", (double) ht->n / ht->nhash);
+
+ printf(" finding and printing each 10th data:\n");
+ for (i = 0; i < size; i += 10) {
+ double* point = &points[i * 3];
+ double* data = ht_find(ht, point);
+
+ if (data != NULL)
+ printf(" i = %d; data = \"%d\"\n", i, (int)* data);
+ else
+ printf(" i = %d; data = <none>\n", i);
+ }
+
+ printf(" removing every 3rd element...");
+ for (i = 0; i < size; i += 3) {
+ double* point = &points[i * 3];
+ ht_delete(ht, point);
+ }
+ printf("done\n");
+
+ printf(" stats:\n");
+ printf(" %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash);
+ printf(" %f entries per hash value in use\n", (double) ht->n / ht->nhash);
+
+ printf(" finding and printing each 10th data:\n");
+ for (i = 0; i < size; i += 10) {
+ double* point = &points[i * 3];
+ double* data = ht_find(ht, point);
+
+ if (data != NULL)
+ printf(" i = %d; data = \"%d\"\n", i, (int)* data);
+ else
+ printf(" i = %d; data = <none>\n", i);
+ }
+
+ printf(" printing all data by calling ht_process():\n ");
+ ht_process(ht, print_double);
+
+ printf("\n destroying the hash table...");
+ ht_destroy(ht);
+ printf("done\n");
+
+ /*
+ * char* key
+ */
+
+ printf("\n2. Testing a table with key of char* type\n\n");
+
+ printf(" creating a table...");
+ ht = ht_create_str(size);
+ printf("done\n");
+
+ printf(" inserting %d elements with deep copy of each data string...", size);
+ for (i = 0; i < size; ++i) {
+ char key[BUFSIZE];
+ char str[BUFSIZE];
+ char* data;
+
+ sprintf(key, "%d-th key", i);
+ sprintf(str, "%d-th data", i);
+ data = strdup(str);
+ ht_insert(ht, key, data);
+ }
+ printf("done\n");
+
+ printf(" stats:\n");
+ printf(" %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash);
+ printf(" %f entries per hash value in use\n", (double) ht->n / ht->nhash);
+
+ printf(" finding and printing each 10th data:\n");
+ for (i = 0; i < size; i += 10) {
+ char key[BUFSIZE];
+ char* data;
+
+ sprintf(key, "%d-th key", i);
+ data = ht_find(ht, key);
+ if (data != NULL)
+ printf(" i = %d; data = \"%s\"\n", i, data);
+ else
+ printf(" i = %d; data = <none>\n", i);
+ }
+
+ printf(" removing every 3rd element...");
+ for (i = 0; i < size; i += 3) {
+ char key[BUFSIZE];
+
+ sprintf(key, "%d-th key", i);
+ free(ht_delete(ht, key));
+ }
+ printf("done\n");
+
+ printf(" stats:\n");
+ printf(" %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash);
+ printf(" %f entries per hash value in use\n", (double) ht->n / ht->nhash);
+
+ printf(" finding and printing each 10th data:\n");
+ for (i = 0; i < size; i += 10) {
+ char key[BUFSIZE];
+ char* data;
+
+ sprintf(key, "%d-th key", i);
+ data = ht_find(ht, key);
+ if (data != NULL)
+ printf(" i = %d; data = \"%s\"\n", i, data);
+ else
+ printf(" i = %d; data = <none>\n", i);
+ }
+
+ printf(" printing all data by calling ht_process():\n ");
+ ht_process(ht, print_string);
+
+ printf("\n freeing the remaining data by calling ht_process()...");
+ ht_process(ht, free);
+ printf("done\n");
+
+ printf(" destroying the hash table...");
+ ht_destroy(ht);
+ printf("done\n");
+
+ return 0;
+}
+
+#endif /* HT_TEST */
diff --git a/src/modules/grid/grid_gridding/nn/hash.h b/src/modules/grid/grid_gridding/nn/hash.h
new file mode 100644
index 0000000..3aa0c13
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/hash.h
@@ -0,0 +1,95 @@
+/******************************************************************************
+ *
+ * File: hash.h
+ *
+ * Purpose: Hash table header
+ *
+ * Author: Jerry Coffin
+ *
+ * Description: Public domain code by Jerry Coffin, with improvements by
+ * HenkJan Wolthuis.
+ * Date last modified: 05-Jul-1997
+ *
+ * Revisions: 18-09-2002 -- modified by Pavel Sakov
+ *
+ *****************************************************************************/
+
+#ifndef _HASH_H
+#define _HASH_H
+
+struct hashtable;
+typedef struct hashtable hashtable;
+
+/** Copies a key. The key must be able to be deallocated by free().
+ */
+typedef void* (*ht_keycp) (void*);
+
+/** Returns 1 if two keys are equal, 0 otherwise.
+ */
+typedef int (*ht_keyeq) (void*, void*);
+
+/** Converts key to an unsigned integer (not necessarily unique).
+ */
+typedef unsigned int (*ht_key2hash) (void*);
+
+/** Creates a hash table of specified size.
+ *
+ * @param size Size of hash table for output points
+ * @param cp Key copy function
+ * @param eq Key equality check function
+ * @param hash Hash value calculation function
+ */
+hashtable* ht_create(int size, ht_keycp cp, ht_keyeq eq, ht_key2hash hash);
+
+/** Create a hash table of specified size and key type.
+ */
+hashtable* ht_create_d1(int size); /* double[1] */
+hashtable* ht_create_d2(int size); /* double[2] */
+hashtable* ht_create_str(int size); /* char* */
+
+/** Destroys a hash table.
+ * (Take care of deallocating data by ht_process() prior to destroying the
+ * table if necessary.)
+ *
+ * @param table Hash table to be destroyed
+ */
+void ht_destroy(hashtable* table);
+
+/** Inserts a new entry into the hash table.
+ *
+ * @param table The hash table
+ * @param key Ponter to entry's key
+ * @param data Pointer to associated data
+ * @return Pointer to the old data associated with the key, NULL if the key
+ * wasn't in the table previously
+ */
+void* ht_insert(hashtable* table, void* key, void* data);
+
+/** Returns a pointer to the data associated with a key. If the key has
+ * not been inserted in the table, returns NULL.
+ *
+ * @param table The hash table
+ * @param key The key
+ * @return The associated data or NULL
+ */
+void* ht_find(hashtable* table, void* key);
+
+/** Deletes an entry from the table. Returns a pointer to the data that
+ * was associated with the key so that the calling code can dispose it
+ * properly.
+ *
+ * @param table The hash table
+ * @param key The key
+ * @return The associated data or NULL
+ */
+void* ht_delete(hashtable* table, void* key);
+
+/** For each entry, calls a specified function with corresponding data as a
+ * parameter.
+ *
+ * @param table The hash table
+ * @param func The action function
+ */
+void ht_process(hashtable* table, void (*func) (void*));
+
+#endif /* _HASH_H */
diff --git a/src/modules/grid/grid_gridding/nn/istack.c b/src/modules/grid/grid_gridding/nn/istack.c
new file mode 100644
index 0000000..1ced5d3
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/istack.c
@@ -0,0 +1,78 @@
+/******************************************************************************
+ *
+ * File: istack.c
+ *
+ * Created: 06/06/2001
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: Handling stack of integers
+ *
+ * Description: None
+ *
+ * Revisions: None
+ *
+ *****************************************************************************/
+
+#define STACK_NSTART 50
+#define STACK_NINC 50
+
+#include <stdlib.h>
+#include <string.h>
+#include "istack.h"
+
+static void istack_init(istack* s)
+{
+ s->n = 0;
+ s->nallocated = STACK_NSTART;
+ s->v = (int *)malloc(STACK_NSTART * sizeof(int));
+}
+
+istack* istack_create()
+{
+ istack* s = (istack *)malloc(sizeof(istack));
+
+ istack_init(s);
+ return s;
+}
+
+void istack_reset(istack* s)
+{
+ s->n = 0;
+}
+
+int istack_contains(istack* s, int v)
+{
+ int i;
+
+ for (i = 0; i < s->n; ++i)
+ if (s->v[i] == v)
+ return 1;
+ return 0;
+}
+
+void istack_push(istack* s, int v)
+{
+ if (s->n == s->nallocated) {
+ s->v = (int *)realloc(s->v, (s->nallocated + STACK_NINC) * sizeof(int));
+ s->nallocated += STACK_NINC;
+ }
+
+ s->v[s->n] = v;
+ s->n++;
+}
+
+int istack_pop(istack* s)
+{
+ s->n--;
+ return s->v[s->n];
+}
+
+void istack_destroy(istack* s)
+{
+ if (s != NULL) {
+ free(s->v);
+ free(s);
+ }
+}
diff --git a/src/modules/grid/grid_gridding/nn/istack.h b/src/modules/grid/grid_gridding/nn/istack.h
new file mode 100644
index 0000000..ff2c5f3
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/istack.h
@@ -0,0 +1,34 @@
+/******************************************************************************
+ *
+ * File: istack.h
+ *
+ * Created: 06/06/2001
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: Header for handling stack of integers.
+ *
+ * Description: None
+ *
+ * Revisions: None
+ *
+ *****************************************************************************/
+
+#if !defined(_ISTACK_H)
+#define _ISTACK_H
+
+typedef struct {
+ int n;
+ int nallocated;
+ int* v;
+} istack;
+
+int istack_contains(istack* s, int v);
+istack* istack_create();
+void istack_destroy(istack* s);
+void istack_push(istack* s, int v);
+int istack_pop(istack* s);
+void istack_reset(istack* s);
+
+#endif
diff --git a/src/modules/grid/grid_gridding/nn/lpi.c b/src/modules/grid/grid_gridding/nn/lpi.c
new file mode 100644
index 0000000..37801ea
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/lpi.c
@@ -0,0 +1,156 @@
+/******************************************************************************
+ *
+ * File: linear.c
+ *
+ * Created: 04/08/2000
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: 2D linear interpolation
+ *
+ * Description: `lpi' -- "Linear Point Interpolator" -- is
+ * a structure for conducting linear interpolation on a given
+ * data on a "point-to-point" basis. It interpolates linearly
+ * within each triangle resulted from the Delaunay
+ * triangluation of input data. `lpi' is much
+ * faster than all Natural Neighbours interpolators in `nn'
+ * library.
+ *
+ * Revisions: None
+ *
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include "nan.h"
+#include "delaunay.h"
+
+typedef struct {
+ double w[3];
+} lweights;
+
+struct lpi {
+ delaunay* d;
+ lweights* weights;
+};
+
+int delaunay_xytoi(delaunay* d, point* p, int seed);
+
+/* Builds linear interpolator.
+ *
+ * @param d Delaunay triangulation
+ * @return Linear interpolator
+ */
+lpi* lpi_build(delaunay* d)
+{
+ int i;
+ lpi* l = (lpi *)malloc(sizeof(lpi));
+
+ l->d = d;
+ l->weights = (lweights *)malloc(d->ntriangles * sizeof(lweights));
+
+ for (i = 0; i < d->ntriangles; ++i) {
+ triangle* t = &d->triangles[i];
+ lweights* lw = &l->weights[i];
+ double x0 = d->points[t->vids[0]].x;
+ double y0 = d->points[t->vids[0]].y;
+ double z0 = d->points[t->vids[0]].z;
+ double x1 = d->points[t->vids[1]].x;
+ double y1 = d->points[t->vids[1]].y;
+ double z1 = d->points[t->vids[1]].z;
+ double x2 = d->points[t->vids[2]].x;
+ double y2 = d->points[t->vids[2]].y;
+ double z2 = d->points[t->vids[2]].z;
+ double x02 = x0 - x2;
+ double y02 = y0 - y2;
+ double z02 = z0 - z2;
+ double x12 = x1 - x2;
+ double y12 = y1 - y2;
+ double z12 = z1 - z2;
+
+ if (y12 != 0.0) {
+ double y0212 = y02 / y12;
+
+ lw->w[0] = (z02 - z12 * y0212) / (x02 - x12 * y0212);
+ lw->w[1] = (z12 - lw->w[0] * x12) / y12;
+ lw->w[2] = (z2 - lw->w[0] * x2 - lw->w[1] * y2);
+ } else {
+ double x0212 = x02 / x12;
+
+ lw->w[1] = (z02 - z12 * x0212) / (y02 - y12 * x0212);
+ lw->w[0] = (z12 - lw->w[1] * y12) / x12;
+ lw->w[2] = (z2 - lw->w[0] * x2 - lw->w[1] * y2);
+ }
+ }
+
+ return l;
+}
+
+/* Destroys linear interpolator.
+ *
+ * @param l Structure to be destroyed
+ */
+void lpi_destroy(lpi* l)
+{
+ free(l->weights);
+ free(l);
+}
+
+/* Finds linearly interpolated value in a point.
+ *
+ * @param l Linear interpolation
+ * @param p Point to be interpolated (p->x, p->y -- input; p->z -- output)
+ */
+void lpi_interpolate_point(lpi* l, point* p)
+{
+ delaunay* d = l->d;
+ int tid = delaunay_xytoi(d, p, d->first_id);
+
+ if (tid >= 0) {
+ lweights* lw = &l->weights[tid];
+
+ d->first_id = tid;
+ p->z = p->x * lw->w[0] + p->y * lw->w[1] + lw->w[2];
+ } else
+ p->z = NaN;
+}
+
+/* Linearly interpolates data from one array of points for another array of
+ * points.
+ *
+ * @param nin Number of input points
+ * @param pin Array of input points [pin]
+ * @param nout Number of ouput points
+ * @param pout Array of output points [nout]
+ */
+void lpi_interpolate_points(int nin, point pin[], int nout, point pout[])
+{
+ delaunay* d = delaunay_build(nin, pin, 0, NULL, 0, NULL);
+ lpi* l = lpi_build(d);
+ int seed = 0;
+ int i;
+
+ if (nn_verbose) {
+ fprintf(stderr, "xytoi:\n");
+ for (i = 0; i < nout; ++i) {
+ point* p = &pout[i];
+
+ fprintf(stderr, "(%.7g,%.7g) -> %d\n", p->x, p->y, delaunay_xytoi(d, p, seed));
+ }
+ }
+
+ for (i = 0; i < nout; ++i)
+ lpi_interpolate_point(l, &pout[i]);
+
+ if (nn_verbose) {
+ fprintf(stderr, "output:\n");
+ for (i = 0; i < nout; ++i) {
+ point* p = &pout[i];;
+ fprintf(stderr, " %d:%15.7g %15.7g %15.7g\n", i, p->x, p->y, p->z);
+ }
+ }
+
+ lpi_destroy(l);
+ delaunay_destroy(d);
+}
diff --git a/src/modules/grid/grid_gridding/nn/nan.h b/src/modules/grid/grid_gridding/nn/nan.h
new file mode 100644
index 0000000..7292fa1
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/nan.h
@@ -0,0 +1,40 @@
+/******************************************************************************
+ *
+ * File: nan.h
+ *
+ * Created: 18/10/2001
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: NaN definition
+ *
+ * Description: Should cover machines with 64 bit doubles or other machines
+ * with GCC
+ *
+ * Revisions: None
+ *
+ *****************************************************************************/
+
+#if !defined(_NAN_H)
+#define _NAN_H
+
+
+#if defined(__GNUC__)
+ static const double NaN = 0.0 / 0.0;
+
+#elif defined(BIG_ENDIAN) || defined(_BIG_ENDIAN)
+ static const long long lNaN = 0x7fffffffffffffff;
+ #define NaN (*(double*)&lNaN)
+
+#elif defined(_SAGA_VC)
+ static const __int64 lNaN = 0xfff8000000000000;
+ #define NaN (*(double*)&lNaN)
+
+#else
+ static const long long lNaN = 0xfff8000000000000;
+ #define NaN (*(double*)&lNaN)
+#endif
+
+
+#endif
diff --git a/src/modules/grid/grid_gridding/nn/nn.h b/src/modules/grid/grid_gridding/nn/nn.h
new file mode 100644
index 0000000..1560b4d
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/nn.h
@@ -0,0 +1,339 @@
+/******************************************************************************
+ *
+ * File: nn.h
+ *
+ * Created: 04/08/2000
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: Header file for nn library
+ *
+ * Description: None
+ *
+ * Revisions: None
+ *
+ *****************************************************************************/
+
+#if !defined(_NN_H)
+#define _NN_H
+
+//---------------------------------------------------------
+#ifdef __cplusplus
+extern "C" {
+#endif
+//---------------------------------------------------------
+
+typedef enum { SIBSON, NON_SIBSONIAN } NN_RULE;
+
+#if !defined(_POINT_STRUCT)
+#define _POINT_STRUCT
+typedef struct {
+ double x;
+ double y;
+ double z;
+} point;
+#endif
+
+/** Smoothes the input point array by averaging the input x,y and z values
+ ** for each cell within virtual rectangular nx by ny grid. The corners of the
+ ** grid are created from min and max values of the input array. It also frees
+ ** the original array and returns results and new dimension via original
+ ** data and size pointers.
+ *
+ * @param pn Pointer to number of points (input/output)
+ * @param ppoints Pointer to array of points (input/output) [*pn]
+ * @param nx Number of x nodes in decimation
+ * @param ny Number of y nodes in decimation
+ */
+void points_thin(int* n, point** points, int nx, int ny);
+
+/** Generates rectangular grid nx by ny using min and max x and y values from
+ ** the input point array. Allocates space for the output point array, be sure
+ ** to free it when necessary!
+ *
+ * @param n Number of points
+ * @param points Array of points [n]
+ * @param nx Number of x nodes
+ * @param ny Number of y nodes
+ * @param nout Pointer to number of output points
+ * @param pout Ppointer to array of output points [*nout]
+ */
+void points_generate1(int n, point points[], int nx, int ny, double zoom, int* nout, point** pout);
+
+/** Generates rectangular grid nx by ny using specified min and max x and y
+ ** values. Allocates space for the output point array, be sure to free it
+ ** when necessary!
+ *
+ * @param xmin Min x value
+ * @param xmax Max x value
+ * @param ymin Min y value
+ * @param ymax Max y value
+ * @param nx Number of x nodes
+ * @param ny Number of y nodes
+ * @param zoom Zoom coefficient
+ * @param nout Pointer to number of output points
+ * @param pout Pointer to array of output points [*nout]
+ */
+void points_generate2(double xmin, double xmax, double ymin, double ymax, int nx, int ny, int* nout, point** pout);
+
+/** Reads array of points from a columnar file.
+ *
+ * @param fname File name (can be "stdin" dor stndard input)
+ * @param dim Number of dimensions (must be 2 or 3)
+ * @param n Pointer to number of points (output)
+ * @param points Pointer to array of points [*n] (output)
+ */
+void points_read(char* fname, int dim, int* n, point** points);
+
+/** Scales Y coordinate so that the resulting set fits into square:
+ ** xmax - xmin = ymax - ymin
+ *
+ * @param n Number of points
+ * @param points The points to scale
+ * @return Y axis compression coefficient
+ */
+double points_scaletosquare(int n, point* points);
+
+/** Compresses Y domain by a given multiple.
+ *
+ * @param n Number of points
+ * @param points The points to scale
+ * @param Y axis compression coefficient as returned by points_scaletosquare()
+ */
+void points_scale(int n, point* points, double k);
+
+/** Structure to perform the Delaunay triangulation of a given array of points.
+ *
+ * Contains a deep copy of the input array of points.
+ * Contains triangles, circles and edges resulted from the triangulation.
+ * Contains neighbour triangles for each triangle.
+ * Contains point to triangle map.
+ */
+struct delaunay;
+typedef struct delaunay delaunay;
+
+/** Builds Delaunay triangulation of the given array of points.
+ *
+ * @param np Number of points
+ * @param points Array of points [np] (input)
+ * @param ns Number of forced segments
+ * @param segments Array of (forced) segment endpoint indices [2*ns]
+ * @param nh Number of holes
+ * @param holes Array of hole (x,y) coordinates [2*nh]
+ * @return Delaunay triangulation with triangulation results
+ */
+delaunay* delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[]);
+
+/** Destroys Delaunay triangulation.
+ *
+ * @param d Structure to be destroyed
+ */
+void delaunay_destroy(delaunay* d);
+
+/** `lpi' -- "linear point interpolator" is a structure for
+ * conducting linear interpolation on a given data on a "point-to-point" basis.
+ * It interpolates linearly within each triangle resulted from the Delaunay
+ * triangluation of input data. `lpi' is much faster than all
+ * Natural Neighbours interpolators below.
+ */
+struct lpi;
+typedef struct lpi lpi;
+
+/** Builds linear interpolator.
+ *
+ * @param d Delaunay triangulation
+ * @return Linear interpolator
+ */
+lpi* lpi_build(delaunay* d);
+
+/** Destroys linear interpolator.
+ *
+ * @param l Structure to be destroyed
+ */
+void lpi_destroy(lpi* l);
+
+/** Finds linearly interpolated value in a point.
+ *
+ * @param l Linear point interpolator
+ * @param p Point to be interpolated (p->x, p->y -- input; p->z -- output)
+ */
+void lpi_interpolate_point(lpi* l, point* p);
+
+/* Linearly interpolates data from one array of points for another array of
+ * points.
+ *
+ * @param nin Number of input points
+ * @param pin Array of input points [pin]
+ * @param nout Number of ouput points
+ * @param pout Array of output points [nout]
+ */
+void lpi_interpolate_points(int nin, point pin[], int nout, point pout[]);
+
+/** `nnpi' -- "Natural Neighbours point interpolator" is a
+ * structure for conducting Natural Neighbours interpolation on a given data on
+ * a "point-to-point" basis. Because it involves weight calculation for each
+ * next output point, it is not particularly suitable for consequitive
+ * interpolations on the same set of observation points -- use
+ * `nnhpi' or `nnai' in these cases.
+ */
+struct nnpi;
+typedef struct nnpi nnpi;
+
+/** Creates Natural Neighbours point interpolator.
+ *
+ * @param d Delaunay triangulation
+ * @return Natural Neighbours interpolation
+ */
+nnpi* nnpi_create(delaunay* d);
+
+/** Destroys Natural Neighbours point interpolation.
+ *
+ * @param nn Structure to be destroyed
+ */
+void nnpi_destroy(nnpi* nn);
+
+/** Finds Natural Neighbours-interpolated value in a point.
+ *
+ * @param nn NN point interpolator
+ * @param p Point to be interpolated (p->x, p->y -- input; p->z -- output)
+ */
+void nnpi_interpolate_point(nnpi* nn, point* p);
+
+/** Natural Neighbours-interpolates data in one array of points for another
+ ** array of points.
+ *
+ * @param nin Number of input points
+ * @param pin Array of input points [pin]
+ * @param wmin Minimal allowed weight
+ * @param nout Number of output points
+ * @param pout Array of output points [nout]
+ */
+void nnpi_interpolate_points(int nin, point pin[], double wmin, int nout, point pout[]);
+
+/** Sets minimal allowed weight for Natural Neighbours interpolation.
+ * @param nn Natural Neighbours point interpolator
+ * @param wmin Minimal allowed weight
+ */
+void nnpi_setwmin(nnpi* nn, double wmin);
+
+/** `nnhpi' is a structure for conducting consequitive
+ * Natural Neighbours interpolations on a given spatial data set in a random
+ * sequence of points from a set of finite size, taking advantage of repeated
+ * interpolations in the same point. It allows to modify Z
+ * coordinate of data between interpolations.
+ */
+struct nnhpi;
+typedef struct nnhpi nnhpi;
+
+/** Creates Natural Neighbours hashing point interpolator.
+ *
+ * @param d Delaunay triangulation
+ * @param size Hash table size (should be of order of number of output points)
+ * @return Natural Neighbours interpolation
+ */
+nnhpi* nnhpi_create(delaunay* d, int size);
+
+/** Destroys Natural Neighbours hashing point interpolation.
+ *
+ * @param nn Structure to be destroyed
+ */
+void nnhpi_destroy(nnhpi* nn);
+
+/** Finds Natural Neighbours-interpolated value in a point.
+ *
+ * @param nnhpi NN hashing point interpolator
+ * @param p Point to be interpolated (p->x, p->y -- input; p->z -- output)
+ */
+void nnhpi_interpolate(nnhpi* nn, point* p);
+
+/** Modifies interpolated data.
+ * Finds point* pd in the underlying Delaunay triangulation such that
+ * pd->x = p->x and pd->y = p->y, and copies p->z to pd->z. Exits with error
+ * if the point is not found.
+ *
+ * @param nn Natural Neighbours hashing point interpolator
+ * @param p New data
+ */
+void nnhpi_modify_data(nnhpi* nn, point* p);
+
+/** Sets minimal allowed weight for Natural Neighbours interpolation.
+ * @param nn Natural Neighbours point hashing interpolator
+ * @param wmin Minimal allowed weight
+ */
+void nnhpi_setwmin(nnhpi* nn, double wmin);
+
+/* `nnai' is a tructure for conducting consequitive Natural
+ * Neighbours interpolations on a given spatial data set in a given array of
+ * points. It allows to modify Z coordinate of data between interpolations.
+ * `nnai' is the fastest of the three Natural Neighbours
+ * interpolators here.
+ */
+struct nnai;
+typedef struct nnai nnai;
+
+/** Builds Natural Neighbours array interpolator. This includes calculation of
+ * weights used in nnai_interpolate().
+ *
+ * @param d Delaunay triangulation
+ * @return Natural Neighbours interpolation
+ */
+nnai* nnai_build(delaunay* d, int n, double* x, double* y);
+
+/** Destroys Natural Neighbours array interpolator.
+ *
+ * @param nn Structure to be destroyed
+ */
+void nnai_destroy(nnai* nn);
+
+/** Conducts NN interpolation in a fixed array of output points using
+ * data specified for a fixed array of input points. Uses pre-calculated
+ * weights.
+ *
+ * @param nn NN array interpolator
+ * @param zin input data [nn->d->npoints]
+ * @param zout output data [nn->n]. Must be pre-allocated!
+ */
+void nnai_interpolate(nnai* nn, double* zin, double* zout);
+
+/** Sets minimal allowed weight for Natural Neighbours interpolation.
+ * @param nn Natural Neighbours array interpolator
+ * @param wmin Minimal allowed weight
+ */
+void nnai_setwmin(nnai* nn, double wmin);
+
+/* Sets the verbosity level within nn package.
+ * 0 (default) - silent
+ * 1 - verbose
+ * 2 - very verbose
+ */
+extern int nn_verbose;
+
+/* Switches between weight calculation methods.
+ * SIBSON -- classic Sibson method
+ * NON_SIBSONIAN -- simpler and (I think) more robust method
+ */
+extern NN_RULE nn_rule;
+
+/* Contains version string for the nn package.
+ */
+extern char* nn_version;
+
+/* Limits verbose information to a particular vertex (used mainly for
+ * debugging purposes).
+ */
+extern int nn_test_vertice;
+
+//---------------------------------------------------------
+#include <float.h>
+
+#include "delaunay.h"
+#include "nan.h"
+
+#ifdef __cplusplus
+} // extern "C" {
+#endif
+//---------------------------------------------------------
+
+
+#endif /* _NN_H */
diff --git a/src/modules/grid/grid_gridding/nn/nnai.c b/src/modules/grid/grid_gridding/nn/nnai.c
new file mode 100644
index 0000000..1cc0a49
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/nnai.c
@@ -0,0 +1,420 @@
+/******************************************************************************
+ *
+ * File: nnai.c
+ *
+ * Created: 15/11/2002
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: Code for:
+ * -- Natural Neighbours Array Interpolator
+ *
+ * Description: `nnai' is a tructure for conducting
+ * consequitive Natural Neighbours interpolations on a given
+ * spatial data set in a given array of points. It allows to
+ * modify Z coordinate of data in between interpolations.
+ * `nnai' is the fastest of the three Natural
+ * Neighbours interpolators in `nn' library.
+ *
+ * Revisions: None
+ *
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include "nn.h"
+#include "delaunay.h"
+#include "nan.h"
+
+typedef struct {
+ int nvertices;
+ int* vertices; /* vertex indices [nvertices] */
+ double* weights; /* vertex weights [nvertices] */
+} nn_weights;
+
+struct nnai {
+ delaunay* d;
+ double wmin;
+ double n; /* number of output points */
+ double* x; /* [n] */
+ double* y; /* [n] */
+ nn_weights* weights;
+};
+
+void nn_quit(char* format, ...);
+void nnpi_calculate_weights(nnpi* nn);
+int nnpi_get_nvertices(nnpi* nn);
+int* nnpi_get_vertices(nnpi* nn);
+double* nnpi_get_weights(nnpi* nn);
+void nnpi_normalize_weights(nnpi* nn);
+void nnpi_reset(nnpi* nn);
+void nnpi_set_point(nnpi* nn, point* p);
+
+/* Builds Natural Neighbours array interpolator. This includes calculation of
+ * weights used in nnai_interpolate().
+ *
+ * @param d Delaunay triangulation
+ * @return Natural Neighbours interpolation
+ */
+nnai* nnai_build(delaunay* d, int n, double* x, double* y)
+{
+ nnai* nn = (nnai *)malloc(sizeof(nnai));
+ nnpi* nnpi = nnpi_create(d);
+ int* vertices;
+ double* weights;
+ int i;
+
+ if (n <= 0)
+ nn_quit("nnai_create(): n = %d\n", n);
+
+ nn->d = d;
+ nn->n = n;
+ nn->x = (double *)malloc(n * sizeof(double));
+ memcpy(nn->x, x, n * sizeof(double));
+ nn->y = (double *)malloc(n * sizeof(double));
+ memcpy(nn->y, y, n * sizeof(double));
+ nn->weights = (nn_weights *)malloc(n * sizeof(nn_weights));
+
+ for (i = 0; i < n; ++i) {
+ nn_weights* w = &nn->weights[i];
+ point p;
+
+ p.x = x[i];
+ p.y = y[i];
+
+ nnpi_reset(nnpi);
+ nnpi_set_point(nnpi, &p);
+ nnpi_calculate_weights(nnpi);
+ nnpi_normalize_weights(nnpi);
+
+ vertices = nnpi_get_vertices(nnpi);
+ weights = nnpi_get_weights(nnpi);
+
+ w->nvertices = nnpi_get_nvertices(nnpi);
+ w->vertices = (int *)malloc(w->nvertices * sizeof(int));
+ memcpy(w->vertices, vertices, w->nvertices * sizeof(int));
+ w->weights = (double *)malloc(w->nvertices * sizeof(double));
+ memcpy(w->weights, weights, w->nvertices * sizeof(double));
+ }
+
+ nnpi_destroy(nnpi);
+
+ return nn;
+}
+
+/* Destroys Natural Neighbours array interpolator.
+ *
+ * @param nn Structure to be destroyed
+ */
+void nnai_destroy(nnai* nn)
+{
+ int i;
+
+ for (i = 0; i < nn->n; ++i) {
+ nn_weights* w = &nn->weights[i];
+
+ free(w->vertices);
+ free(w->weights);
+ }
+
+ free(nn->x);
+ free(nn->y);
+ free(nn->weights);
+ free(nn);
+}
+
+/* Conducts NN interpolation in a fixed array of output points using
+ * data specified for a fixed array of input points. Uses pre-calculated
+ * weights.
+ *
+ * @param nn NN array interpolator
+ * @param zin input data [nn->d->npoints]
+ * @param zout output data [nn->n]. Must be pre-allocated!
+ */
+void nnai_interpolate(nnai* nn, double* zin, double* zout)
+{
+ int i;
+
+ for (i = 0; i < nn->n; ++i) {
+ nn_weights* w = &nn->weights[i];
+ double z = 0.0;
+ int j;
+
+ for (j = 0; j < w->nvertices; ++j) {
+ double weight = w->weights[j];
+
+ if (weight < nn->wmin) {
+ z = NaN;
+ break;
+ }
+ z += weight * zin[w->vertices[j]];
+ }
+
+ zout[i] = z;
+ }
+}
+
+/** Sets minimal allowed weight for Natural Neighbours interpolation.
+ * @param nn Natural Neighbours array interpolator
+ * @param wmin Minimal allowed weight
+ */
+void nnai_setwmin(nnai* nn, double wmin)
+{
+ nn->wmin = wmin;
+}
+
+/* The rest of this file contains a number of test programs.
+ */
+#if defined(NNAI_TEST)
+
+#include <sys/time.h>
+
+#define NPOINTSIN 10000
+#define NMIN 10
+#define NX 101
+#define NXMIN 1
+
+#define SQ(x) ((x) * (x))
+
+static double franke(double x, double y)
+{
+ x *= 9.0;
+ y *= 9.0;
+ return 0.75 * exp((-SQ(x - 2.0) - SQ(y - 2.0)) / 4.0)
+ + 0.75 * exp(-SQ(x - 2.0) / 49.0 - (y - 2.0) / 10.0)
+ + 0.5 * exp((-SQ(x - 7.0) - SQ(y - 3.0)) / 4.0)
+ - 0.2 * exp(-SQ(x - 4.0) - SQ(y - 7.0));
+}
+
+/* *INDENT-OFF* */
+static void usage()
+{
+ printf(
+"Usage: nn_test [-v|-V] [-n <nin> <nxout>]\n"
+"Options:\n"
+" -a -- use non-Sibsonian interpolation rule\n"
+" -n <nin> <nout>:\n"
+" <nin> -- number of input points (default = 10000)\n"
+" <nout> -- number of output points per side (default = 64)\n"
+" -v -- verbose\n"
+" -V -- very verbose\n"
+);
+}
+/* *INDENT-ON* */
+
+int main(int argc, char* argv[])
+{
+ int nin = NPOINTSIN;
+ int nx = NX;
+ int nout = 0;
+ point* pin = NULL;
+ delaunay* d = NULL;
+ point* pout = NULL;
+ nnai* nn = NULL;
+ double* zin = NULL;
+ double* xout = NULL;
+ double* yout = NULL;
+ double* zout = NULL;
+ int cpi = -1; /* control point index */
+ struct timeval tv0, tv1, tv2;
+ struct timezone tz;
+ int i;
+
+ i = 1;
+ while (i < argc) {
+ switch (argv[i][1]) {
+ case 'a':
+ i++;
+ nn_rule = NON_SIBSONIAN;
+ break;
+ case 'n':
+ i++;
+ if (i >= argc)
+ nn_quit("no number of data points found after -i\n");
+ nin = atoi(argv[i]);
+ i++;
+ if (i >= argc)
+ nn_quit("no number of ouput points per side found after -i\n");
+ nx = atoi(argv[i]);
+ i++;
+ break;
+ case 'v':
+ i++;
+ nn_verbose = 1;
+ break;
+ case 'V':
+ i++;
+ nn_verbose = 2;
+ break;
+ default:
+ usage();
+ break;
+ }
+ }
+
+ if (nin < NMIN)
+ nin = NMIN;
+ if (nx < NXMIN)
+ nx = NXMIN;
+
+ printf("\nTest of Natural Neighbours array interpolator:\n\n");
+ printf(" %d data points\n", nin);
+ printf(" %d output points\n", nx * nx);
+
+ /*
+ * generate data
+ */
+ printf(" generating data:\n");
+ fflush(stdout);
+ pin = (point *)malloc(nin * sizeof(point));
+ zin = (double *)malloc(nin * sizeof(double));
+ for (i = 0; i < nin; ++i) {
+ point* p = &pin[i];
+
+ p->x = (double) random() / RAND_MAX;
+ p->y = (double) random() / RAND_MAX;
+ p->z = franke(p->x, p->y);
+ zin[i] = p->z;
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ /*
+ * triangulate
+ */
+ printf(" triangulating:\n");
+ fflush(stdout);
+ d = delaunay_build(nin, pin, 0, NULL, 0, NULL);
+
+ /*
+ * generate output points
+ */
+ points_generate2(-0.1, 1.1, -0.1, 1.1, nx, nx, &nout, &pout);
+ xout = (double *)malloc(nout * sizeof(double));
+ yout = (double *)malloc(nout * sizeof(double));
+ zout = (double *)malloc(nout * sizeof(double));
+ for (i = 0; i < nout; ++i) {
+ point* p = &pout[i];
+
+ xout[i] = p->x;
+ yout[i] = p->y;
+ zout[i] = NaN;
+ }
+ cpi = (nx / 2) * (nx + 1);
+
+ gettimeofday(&tv0, &tz);
+
+ /*
+ * create interpolator
+ */
+ printf(" creating interpolator:\n");
+ fflush(stdout);
+ nn = nnai_build(d, nout, xout, yout);
+
+ fflush(stdout);
+ gettimeofday(&tv1, &tz);
+ {
+ long dt = 1000000 * (tv1.tv_sec - tv0.tv_sec) + tv1.tv_usec - tv0.tv_usec;
+
+ printf(" interpolator creation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout);
+ }
+
+ /*
+ * interpolate
+ */
+ printf(" interpolating:\n");
+ fflush(stdout);
+ nnai_interpolate(nn, zin, zout);
+ if (nn_verbose)
+ for (i = 0; i < nout; ++i)
+ printf(" (%f, %f, %f)\n", xout[i], yout[i], zout[i]);
+
+ fflush(stdout);
+ gettimeofday(&tv2, &tz);
+ {
+ long dt = 1000000.0 * (tv2.tv_sec - tv1.tv_sec) + tv2.tv_usec - tv1.tv_usec;
+
+ printf(" interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout);
+ }
+
+ if (!nn_verbose)
+ printf(" control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], franke(xout[cpi], yout[cpi]));
+
+ printf(" interpolating one more time:\n");
+ fflush(stdout);
+ nnai_interpolate(nn, zin, zout);
+ if (nn_verbose)
+ for (i = 0; i < nout; ++i)
+ printf(" (%f, %f, %f)\n", xout[i], yout[i], zout[i]);
+
+ fflush(stdout);
+ gettimeofday(&tv0, &tz);
+ {
+ long dt = 1000000.0 * (tv0.tv_sec - tv2.tv_sec) + tv0.tv_usec - tv2.tv_usec;
+
+ printf(" interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout);
+ }
+
+ if (!nn_verbose)
+ printf(" control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], franke(xout[cpi], yout[cpi]));
+
+ printf(" entering new data:\n");
+ fflush(stdout);
+ for (i = 0; i < nin; ++i) {
+ point* p = &pin[i];
+
+ p->z = p->x * p->x - p->y * p->y;
+ zin[i] = p->z;
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ printf(" interpolating:\n");
+ fflush(stdout);
+ nnai_interpolate(nn, zin, zout);
+ if (nn_verbose)
+ for (i = 0; i < nout; ++i)
+ printf(" (%f, %f, %f)\n", xout[i], yout[i], zout[i]);
+
+ if (!nn_verbose)
+ printf(" control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], xout[cpi] * xout[cpi] - yout[cpi] * yout[cpi]);
+
+ printf(" restoring data:\n");
+ fflush(stdout);
+ for (i = 0; i < nin; ++i) {
+ point* p = &pin[i];
+
+ p->z = franke(p->x, p->y);
+ zin[i] = p->z;
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ printf(" interpolating:\n");
+ fflush(stdout);
+ nnai_interpolate(nn, zin, zout);
+ if (nn_verbose)
+ for (i = 0; i < nout; ++i)
+ printf(" (%f, %f, %f)\n", xout[i], yout[i], zout[i]);
+
+ if (!nn_verbose)
+ printf(" control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], franke(xout[cpi], yout[cpi]));
+
+ printf("\n");
+
+ nnai_destroy(nn);
+ free(zin);
+ free(xout);
+ free(yout);
+ free(zout);
+ free(pout);
+ delaunay_destroy(d);
+ free(pin);
+
+ return 0;
+}
+
+#endif
diff --git a/src/modules/grid/grid_gridding/nn/nncommon.c b/src/modules/grid/grid_gridding/nn/nncommon.c
new file mode 100644
index 0000000..7d4580e
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/nncommon.c
@@ -0,0 +1,498 @@
+/******************************************************************************
+ *
+ * File: nncommon.c
+ *
+ * Created: 04/08/2000
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: Common stuff for NN interpolation library
+ *
+ * Description: None
+ *
+ * Revisions: 15/11/2002 PS: Changed name from "utils.c"
+ * 28/02/2003 PS: Modified points_read() to do the job without
+ * rewinding the file. This allows to read from stdin when
+ * necessary.
+ * 09/04/2003 PS: Modified points_read() to read from a
+ * file specified by name, not by handle.
+ *
+ *****************************************************************************/
+#ifdef _SAGA_MSW
+ #define isnan _isnan
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include <assert.h>
+#include <math.h>
+#include <limits.h>
+#include <float.h>
+#include <string.h>
+#include <errno.h>
+#include "nan.h"
+#include "delaunay.h"
+
+#define BUFSIZE 1024
+
+int nn_verbose = 0;
+int nn_test_vertice = -1;
+NN_RULE nn_rule = SIBSON;
+
+#include "version.h"
+
+void nn_quit(char* format, ...)
+{
+ va_list args;
+
+ fflush(stdout); /* just in case, to have the exit message
+ * last */
+
+ fprintf(stderr, "error: nn: ");
+ va_start(args, format);
+ vfprintf(stderr, format, args);
+ va_end(args);
+
+ exit(1);
+}
+
+int circle_build(circle* c, point* p1, point* p2, point* p3)
+{
+ double x1sq = p1->x * p1->x;
+ double x2sq = p2->x * p2->x;
+ double x3sq = p3->x * p3->x;
+ double y1sq = p1->y * p1->y;
+ double y2sq = p2->y * p2->y;
+ double y3sq = p3->y * p3->y;
+ double t1 = x3sq - x2sq + y3sq - y2sq;
+ double t2 = x1sq - x3sq + y1sq - y3sq;
+ double t3 = x2sq - x1sq + y2sq - y1sq;
+ double D = (p1->x * (p2->y - p3->y) + p2->x * (p3->y - p1->y) + p3->x * (p1->y - p2->y)) * 2.0;
+
+ if (D == 0.0)
+ return 0;
+
+ c->x = (p1->y * t1 + p2->y * t2 + p3->y * t3) / D;
+ c->y = -(p1->x * t1 + p2->x * t2 + p3->x * t3) / D;
+ c->r = hypot(c->x - p1->x, c->y - p1->y);
+
+ return 1;
+}
+
+/* This procedure has taken it final shape after a number of tries. The problem
+ * was to have the calculated and stored radii being the same if (x,y) is
+ * exactly on the circle border (i.e. not to use FCPU extended precision in
+ * the radius calculation). This may have little effect in practice but was
+ * important in some tests when both input and output data were placed
+ * in rectangular grid nodes.
+ */
+int circle_contains(circle* c, point* p)
+{
+ return hypot(c->x - p->x, c->y - p->y) <= c->r;
+}
+
+/* Smoothes the input point array by averaging the input x,y and z values
+ * for each cell within virtual rectangular nx by ny grid. The corners of the
+ * grid are created from min and max values of the input array. It also frees
+ * the original array and returns results and new dimension via original
+ * data and size pointers.
+ *
+ * @param pn Pointer to number of points (input/output)
+ * @param ppoints Pointer to array of points (input/output) [*pn]
+ * @param nx Number of x nodes in decimation
+ * @param ny Number of y nodes in decimation
+ */
+void points_thin(int* pn, point** ppoints, int nx, int ny)
+{
+ int n = *pn;
+ point* points = *ppoints;
+ double xmin = DBL_MAX;
+ double xmax = -DBL_MAX;
+ double ymin = DBL_MAX;
+ double ymax = -DBL_MAX;
+ int nxy = nx * ny;
+ double* sumx = (double *)calloc(nxy, sizeof(double));
+ double* sumy = (double *)calloc(nxy, sizeof(double));
+ double* sumz = (double *)calloc(nxy, sizeof(double));
+ int* count = (int *)calloc(nxy, sizeof(int));
+ double stepx = 0.0;
+ double stepy = 0.0;
+ int nnew = 0;
+ point* pointsnew = NULL;
+ int i, j, ii;
+
+ if (nn_verbose)
+ fprintf(stderr, "thinned: %d points -> ", *pn);
+
+ if (nx < 1 || ny < 1) {
+ free(points);
+ *ppoints = NULL;
+ *pn = 0;
+ if (nn_verbose)
+ fprintf(stderr, "0 points");
+ return;
+ }
+
+ for (ii = 0; ii < n; ++ii) {
+ point* p = &points[ii];
+
+ if (p->x < xmin)
+ xmin = p->x;
+ if (p->x > xmax)
+ xmax = p->x;
+ if (p->y < ymin)
+ ymin = p->y;
+ if (p->y > ymax)
+ ymax = p->y;
+ }
+
+ stepx = (nx > 1) ? (xmax - xmin) / nx : 0.0;
+ stepy = (ny > 1) ? (ymax - ymin) / ny : 0.0;
+
+ for (ii = 0; ii < n; ++ii) {
+ point* p = &points[ii];
+ int index;
+
+ /*
+ * Following is the portion of the code which really depends on the
+ * floating point particulars. Do not be surprised if different
+ * compilers/options give different results here.
+ */
+ i = (nx == 1) ? 0 : (int)((p->x - xmin) / stepx);
+ j = (ny == 1) ? 0 : (int)((p->y - ymin) / stepy);
+
+ if (i == nx)
+ i--;
+ if (j == ny)
+ j--;
+ index = i + j * nx;
+ sumx[index] += p->x;
+ sumy[index] += p->y;
+ sumz[index] += p->z;
+ count[index]++;
+ }
+
+ for (j = 0; j < ny; ++j) {
+ for (i = 0; i < nx; ++i) {
+ int index = i + j * nx;
+
+ if (count[index] > 0)
+ nnew++;
+ }
+ }
+
+ pointsnew = (point *)malloc(nnew * sizeof(point));
+
+ ii = 0;
+ for (j = 0; j < ny; ++j) {
+ for (i = 0; i < nx; ++i) {
+ int index = i + j * nx;
+ int nn = count[index];
+
+ if (nn > 0) {
+ point* p = &pointsnew[ii];
+
+ p->x = sumx[index] / nn;
+ p->y = sumy[index] / nn;
+ p->z = sumz[index] / nn;
+ ii++;
+ }
+ }
+ }
+
+ if (nn_verbose)
+ fprintf(stderr, "%d points\n", nnew);
+
+ free(sumx);
+ free(sumy);
+ free(sumz);
+ free(count);
+
+ free(points);
+ *ppoints = pointsnew;
+ *pn = nnew;
+}
+
+/* Generates rectangular grid nx by ny using min and max x and y values from
+ * the input point array. Allocates space for the output point array, be sure
+ * to free it when necessary!
+ *
+ * @param n Number of points
+ * @param points Array of points [n]
+ * @param nx Number of x nodes
+ * @param ny Number of y nodes
+ * @param zoom Zoom coefficient
+ * @param nout Pointer to number of output points
+ * @param pout Pointer to array of output points [*nout]
+ */
+void points_generate1(int nin, point pin[], int nx, int ny, double zoom, int* nout, point** pout)
+{
+ double xmin = DBL_MAX;
+ double xmax = -DBL_MAX;
+ double ymin = DBL_MAX;
+ double ymax = -DBL_MAX;
+ double stepx, stepy;
+ double x0, xx, yy;
+ int i, j, ii;
+
+ if (nx < 1 || ny < 1) {
+ *pout = NULL;
+ *nout = 0;
+ return;
+ }
+
+ for (ii = 0; ii < nin; ++ii) {
+ point* p = &pin[ii];
+
+ if (p->x < xmin)
+ xmin = p->x;
+ if (p->x > xmax)
+ xmax = p->x;
+ if (p->y < ymin)
+ ymin = p->y;
+ if (p->y > ymax)
+ ymax = p->y;
+ }
+
+ if (isnan(zoom) || zoom <= 0.0)
+ zoom = 1.0;
+
+ if (zoom != 1.0) {
+ double xdiff2 = (xmax - xmin) / 2.0;
+ double ydiff2 = (ymax - ymin) / 2.0;
+ double xav = (xmax + xmin) / 2.0;
+ double yav = (ymax + ymin) / 2.0;
+
+ xmin = xav - xdiff2 * zoom;
+ xmax = xav + xdiff2 * zoom;
+ ymin = yav - ydiff2 * zoom;
+ ymax = yav + ydiff2 * zoom;
+ }
+
+ *nout = nx * ny;
+ *pout = (point *)malloc(*nout * sizeof(point));
+
+ stepx = (nx > 1) ? (xmax - xmin) / (nx - 1) : 0.0;
+ stepy = (ny > 1) ? (ymax - ymin) / (ny - 1) : 0.0;
+ x0 = (nx > 1) ? xmin : (xmin + xmax) / 2.0;
+ yy = (ny > 1) ? ymin : (ymin + ymax) / 2.0;
+
+ ii = 0;
+ for (j = 0; j < ny; ++j) {
+ xx = x0;
+ for (i = 0; i < nx; ++i) {
+ point* p = &(*pout)[ii];
+
+ p->x = xx;
+ p->y = yy;
+ xx += stepx;
+ ii++;
+ }
+ yy += stepy;
+ }
+}
+
+/* Generates rectangular grid nx by ny using specified min and max x and y
+ * values. Allocates space for the output point array, be sure to free it
+ * when necessary!
+ *
+ * @param xmin Min x value
+ * @param xmax Max x value
+ * @param ymin Min y value
+ * @param ymax Max y value
+ * @param nx Number of x nodes
+ * @param ny Number of y nodes
+ * @param nout Pointer to number of output points
+ * @param pout Pointer to array of output points [*nout]
+ */
+void points_generate2(double xmin, double xmax, double ymin, double ymax, int nx, int ny, int* nout, point** pout)
+{
+ double stepx, stepy;
+ double x0, xx, yy;
+ int i, j, ii;
+
+ if (nx < 1 || ny < 1) {
+ *pout = NULL;
+ *nout = 0;
+ return;
+ }
+
+ *nout = nx * ny;
+ *pout = (point *)malloc(*nout * sizeof(point));
+
+ stepx = (nx > 1) ? (xmax - xmin) / (nx - 1) : 0.0;
+ stepy = (ny > 1) ? (ymax - ymin) / (ny - 1) : 0.0;
+ x0 = (nx > 1) ? xmin : (xmin + xmax) / 2.0;
+ yy = (ny > 1) ? ymin : (ymin + ymax) / 2.0;
+
+ ii = 0;
+ for (j = 0; j < ny; ++j) {
+ xx = x0;
+ for (i = 0; i < nx; ++i) {
+ point* p = &(*pout)[ii];
+
+ p->x = xx;
+ p->y = yy;
+ xx += stepx;
+ ii++;
+ }
+ yy += stepy;
+ }
+}
+
+static int str2double(char* token, double* value)
+{
+ char* end = NULL;
+
+ if (token == NULL) {
+ *value = NaN;
+ return 0;
+ }
+
+ *value = strtod(token, &end);
+
+ if (end == token) {
+ *value = NaN;
+ return 0;
+ }
+
+ return 1;
+}
+
+#define NALLOCATED_START 1024
+
+/* Reads array of points from a columnar file.
+ *
+ * @param fname File name (can be "stdin" for standard input)
+ * @param dim Number of dimensions (must be 2 or 3)
+ * @param n Pointer to number of points (output)
+ * @param points Pointer to array of points [*n] (output) (to be freed)
+ */
+void points_read(char* fname, int dim, int* n, point** points)
+{
+ FILE* f = NULL;
+ int nallocated = NALLOCATED_START;
+ char buf[BUFSIZE];
+ char seps[] = " ,;\t";
+ char* token;
+
+ if (dim < 2 || dim > 3) {
+ *n = 0;
+ *points = NULL;
+ return;
+ }
+
+ if (fname == NULL)
+ f = stdin;
+ else {
+ if (strcmp(fname, "stdin") == 0 || strcmp(fname, "-") == 0)
+ f = stdin;
+ else {
+ f = fopen(fname, "r");
+ if (f == NULL)
+ nn_quit("%s: %s\n", fname, strerror(errno));
+ }
+ }
+
+ *points = (point *)malloc(nallocated * sizeof(point));
+ *n = 0;
+ while (fgets(buf, BUFSIZE, f) != NULL) {
+ point* p;
+
+ if (*n == nallocated) {
+ nallocated *= 2;
+ *points = (point *)realloc(*points, nallocated * sizeof(point));
+ }
+
+ p = &(*points)[*n];
+
+ if (buf[0] == '#')
+ continue;
+ if ((token = strtok(buf, seps)) == NULL)
+ continue;
+ if (!str2double(token, &p->x))
+ continue;
+ if ((token = strtok(NULL, seps)) == NULL)
+ continue;
+ if (!str2double(token, &p->y))
+ continue;
+ if (dim == 2)
+ p->z = NaN;
+ else {
+ if ((token = strtok(NULL, seps)) == NULL)
+ continue;
+ if (!str2double(token, &p->z))
+ continue;
+ }
+ (*n)++;
+ }
+
+ if (*n == 0) {
+ free(*points);
+ *points = NULL;
+ } else
+ *points = (point *)realloc(*points, *n * sizeof(point));
+
+ if (f != stdin)
+ if (fclose(f) != 0)
+ nn_quit("%s: %s\n", fname, strerror(errno));
+}
+
+/** Scales Y coordinate so that the resulting set fits into square:
+ ** xmax - xmin = ymax - ymin
+ *
+ * @param n Number of points
+ * @param points The points to scale
+ * @return Y axis compression coefficient
+ */
+double points_scaletosquare(int n, point* points)
+{
+ double xmin, ymin, xmax, ymax;
+ double k;
+ int i;
+
+ if (n <= 0)
+ return NaN;
+
+ xmin = xmax = points[0].x;
+ ymin = ymax = points[0].y;
+
+ for (i = 1; i < n; ++i) {
+ point* p = &points[i];
+
+ if (p->x < xmin)
+ xmin = p->x;
+ else if (p->x > xmax)
+ xmax = p->x;
+ if (p->y < ymin)
+ ymin = p->y;
+ else if (p->y > ymax)
+ ymax = p->y;
+ }
+
+ if (xmin == xmax || ymin == ymax)
+ return NaN;
+ else
+ k = (ymax - ymin) / (xmax - xmin);
+
+ for (i = 0; i < n; ++i)
+ points[i].y /= k;
+
+ return k;
+}
+
+/** Compresses Y domain by a given multiple.
+ *
+ * @param n Number of points
+ * @param points The points to scale
+ * @param Y axis compression coefficient as returned by points_scaletosquare()
+ */
+void points_scale(int n, point* points, double k)
+{
+ int i;
+
+ for (i = 0; i < n; ++i)
+ points[i].y /= k;
+}
diff --git a/src/modules/grid/grid_gridding/nn/nnpi.c b/src/modules/grid/grid_gridding/nn/nnpi.c
new file mode 100644
index 0000000..81cdeda
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/nnpi.c
@@ -0,0 +1,857 @@
+ /*****************************************************************************
+ *
+ * File: nnpi.c
+ *
+ * Created: 15/11/2002
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: Code for:
+ * -- Natural Neighbours Point Interpolator
+ * -- Natural Neighbours Point Hashing Interpolator
+ *
+ * Description: `nnpi' -- "Natural Neighbours Point
+ * Interpolator" -- is a structure for conducting Natural
+ * Neighbours interpolation on a given data on a
+ * "point-to-point" basis. Because it involves weight
+ * calculation for each next output point, it is not
+ * particularly suitable for consequitive interpolations on
+ * the same set of observation points -- use
+ * `nnhpi' or `nnai'
+ * in these cases.
+ *
+ * `nnhpi' is a structure for
+ * conducting consequitive Natural Neighbours interpolations
+ * on a given spatial data set in a random sequence of points
+ * from a set of finite size, taking advantage of repeated
+ * interpolations in the same point. It allows to modify Z
+ * coordinate of data in between interpolations.
+ *
+ *
+ * Revisions: 01/04/2003 PS: modified nnpi_triangle_process(): for
+ * Sibson interpolation, if circle_build fails(), now a
+ * local copy of a point is moved slightly rather than the
+ * data point itself. The later approach have found leading
+ * to inconsistencies of the new point position with the
+ * earlier built triangulation.
+ *
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <limits.h>
+#include <float.h>
+#include <string.h>
+#include <assert.h>
+#include <math.h>
+#include "nn.h"
+#include "delaunay.h"
+#include "nan.h"
+#include "hash.h"
+
+struct nnpi {
+ delaunay* d;
+ point* p;
+ double wmin;
+ /*
+ * work variables
+ */
+ int nvertices;
+ int nallocated;
+ int* vertices; /* vertex indices */
+ double* weights;
+ int n; /* number of points processed */
+};
+
+int circle_build(circle* c, point* p0, point* p1, point* p2);
+int circle_contains(circle* c, point* p);
+void delaunay_circles_find(delaunay* d, point* p, int* n, int** out);
+int delaunay_xytoi(delaunay* d, point* p, int seed);
+void nn_quit(char* format, ...);
+
+#define NSTART 10
+#define NINC 10
+#define EPS_SHIFT 1.0e-9
+#define N_SEARCH_TURNON 20
+#define BIGNUMBER 1.0e+100
+
+#ifndef min
+#define min(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#ifndef max
+#define max(x,y) ((x) > (y) ? (x) : (y))
+#endif
+
+/* Creates Natural Neighbours point interpolator.
+ *
+ * @param d Delaunay triangulation
+ * @return Natural Neighbours interpolation
+ */
+nnpi* nnpi_create(delaunay* d)
+{
+ nnpi* nn = (nnpi *)malloc(sizeof(nnpi));
+
+ nn->d = d;
+ nn->wmin = -DBL_MAX;
+ nn->vertices = (int *)calloc(NSTART, sizeof(int));
+ nn->weights = (double *)calloc(NSTART, sizeof(double));
+ nn->nvertices = 0;
+ nn->nallocated = NSTART;
+ nn->p = NULL;
+ nn->n = 0;
+
+ return nn;
+}
+
+/* Destroys Natural Neighbours point interpolator.
+ *
+ * @param nn Structure to be destroyed
+ */
+void nnpi_destroy(nnpi* nn)
+{
+ free(nn->weights);
+ free(nn->vertices);
+ free(nn);
+}
+
+void nnpi_reset(nnpi* nn)
+{
+ nn->nvertices = 0;
+ nn->p = NULL;
+ memset(nn->d->flags, 0, nn->d->ntriangles * sizeof(int));
+}
+
+static void nnpi_add_weight(nnpi* nn, int vertex, double w)
+{
+ int i;
+
+ /*
+ * find whether the vertex is already in the list
+ */
+ for (i = 0; i < nn->nvertices; ++i)
+ if (nn->vertices[i] == vertex)
+ break;
+
+ if (i == nn->nvertices) { /* not in the list */
+ /*
+ * get more memory if necessary
+ */
+ if (nn->nvertices == nn->nallocated) {
+ nn->vertices = (int *)realloc(nn->vertices, (nn->nallocated + NINC) * sizeof(int));
+ nn->weights = (double *)realloc(nn->weights, (nn->nallocated + NINC) * sizeof(double));
+ nn->nallocated += NINC;
+ }
+
+ /*
+ * add the vertex to the list
+ */
+ nn->vertices[i] = vertex;
+ nn->weights[i] = w;
+ nn->nvertices++;
+ } else { /* in the list */
+
+ if (nn_rule == SIBSON)
+ nn->weights[i] += w;
+ else if (w > nn->weights[i])
+ nn->weights[i] = w;
+ }
+}
+
+static double triangle_scale_get(delaunay* d, triangle* t)
+{
+ double x0 = d->points[t->vids[0]].x;
+ double x1 = d->points[t->vids[1]].x;
+ double x2 = d->points[t->vids[2]].x;
+ double y0 = d->points[t->vids[0]].y;
+ double y1 = d->points[t->vids[1]].y;
+ double y2 = d->points[t->vids[2]].y;
+ double xmin = min(min(x0, x1), x2);
+ double xmax = max(max(x0, x1), x2);
+ double ymin = min(min(y0, y1), y2);
+ double ymax = max(max(y0, y1), y2);
+
+ return xmax - xmin + ymax - ymin;
+}
+
+/* This is a central procedure for the Natural Neighbours interpolation. It
+ * uses the Watson's algorithm for the required areas calculation and implies
+ * that the vertices of the delaunay triangulation are listed in uniform
+ * (clockwise or counterclockwise) order.
+ */
+static void nnpi_triangle_process(nnpi* nn, point* p, int i)
+{
+ delaunay* d = nn->d;
+ triangle* t = &d->triangles[i];
+ circle* c = &d->circles[i];
+ circle cs[3];
+ int j;
+
+ assert(circle_contains(c, p));
+
+ if (nn_rule == SIBSON) {
+ point pp;
+
+ pp.x = p->x;
+ pp.y = p->y;
+ /*
+ * Sibson interpolation by using Watson's algorithm
+ */
+ do {
+ for (j = 0; j < 3; ++j) {
+ int j1 = (j + 1) % 3;
+ int j2 = (j + 2) % 3;
+ int v1 = t->vids[j1];
+ int v2 = t->vids[j2];
+
+ if (!circle_build(&cs[j], &d->points[v1], &d->points[v2], &pp)) {
+ double scale = triangle_scale_get(d, t);
+
+ if (d->points[v1].y == d->points[v2].y)
+ pp.y += EPS_SHIFT * scale;
+ else
+ pp.x += EPS_SHIFT * scale;
+ break;
+ }
+ }
+ } while (j != 3);
+
+ for (j = 0; j < 3; ++j) {
+ int j1 = (j + 1) % 3;
+ int j2 = (j + 2) % 3;
+ double det = ((cs[j1].x - c->x) * (cs[j2].y - c->y) - (cs[j2].x - c->x) * (cs[j1].y - c->y));
+
+ nnpi_add_weight(nn, t->vids[j], det);
+ }
+ } else if (nn_rule == NON_SIBSONIAN) {
+ double d1 = c->r - hypot(p->x - c->x, p->y - c->y);
+
+ for (i = 0; i < 3; ++i) {
+ int vid = t->vids[i];
+ point* pp = &d->points[vid];
+ double d2 = hypot(p->x - pp->x, p->y - pp->y);
+
+ if (d2 == 0.0)
+ nnpi_add_weight(nn, vid, BIGNUMBER);
+ else
+ nnpi_add_weight(nn, vid, d1 / d2);
+ }
+ } else
+ nn_quit("unknown rule\n");
+}
+
+void nnpi_calculate_weights(nnpi* nn)
+{
+ point* p = nn->p;
+ int n = nn->d->ntriangles;
+ int i;
+
+ if (n > N_SEARCH_TURNON) {
+ int* tids;
+
+ delaunay_circles_find(nn->d, p, &n, &tids);
+ for (i = 0; i < n; ++i)
+ nnpi_triangle_process(nn, p, tids[i]);
+ } else
+ for (i = 0; i < n; ++i)
+ if (circle_contains(&nn->d->circles[i], p))
+ nnpi_triangle_process(nn, p, i);
+}
+
+void nnpi_normalize_weights(nnpi* nn)
+{
+ int n = nn->nvertices;
+ double sum = 0.0;
+ int i;
+
+ for (i = 0; i < n; ++i)
+ sum += nn->weights[i];
+
+ for (i = 0; i < n; ++i)
+ nn->weights[i] /= sum;
+}
+
+/* Finds Natural Neighbours-interpolated value for a point.
+ *
+ * @param nn NN interpolation
+ * @param p Point to be interpolated (p->x, p->y -- input; p->z -- output)
+ */
+void nnpi_interpolate_point(nnpi* nn, point* p)
+{
+ delaunay* d = nn->d;
+ int i;
+
+ nnpi_reset(nn);
+ nn->p = p;
+ nnpi_calculate_weights(nn);
+ nnpi_normalize_weights(nn);
+
+ if (nn_verbose) {
+ if (nn_test_vertice == -1) {
+ if (nn->n == 0)
+ fprintf(stderr, "weights:\n");
+ fprintf(stderr, " %d: {", nn->n);
+ for (i = 0; i < nn->nvertices; ++i) {
+ fprintf(stderr, "(%d,%.5g)", nn->vertices[i], nn->weights[i]);
+ if (i < nn->nvertices - 1)
+ fprintf(stderr, ", ");
+ }
+ fprintf(stderr, "}\n");
+ } else {
+ double w = 0.0;
+
+ if (nn->n == 0)
+ fprintf(stderr, "weights for vertex %d:\n", nn_test_vertice);
+ for (i = 0; i < nn->nvertices; ++i) {
+ if (nn->vertices[i] == nn_test_vertice) {
+ w = nn->weights[i];
+ break;
+ }
+ }
+ fprintf(stderr, "%15.7g %15.7g %15.7g\n", p->x, p->y, w);
+ }
+ }
+
+ nn->n++;
+
+ if (nn->nvertices == 0) {
+ p->z = NaN;
+ return;
+ }
+
+ p->z = 0.0;
+ for (i = 0; i < nn->nvertices; ++i) {
+ double weight = nn->weights[i];
+
+ if (weight < nn->wmin) {
+ p->z = NaN;
+ return;
+ }
+ p->z += d->points[nn->vertices[i]].z * weight;
+ }
+}
+
+/* Performs Natural Neighbours interpolation for an array of points.
+ *
+ * @param nin Number of input points
+ * @param pin Array of input points [pin]
+ * @param wmin Minimal allowed weight
+ * @param nout Number of output points
+ * @param pout Array of output points [nout]
+ */
+void nnpi_interpolate_points(int nin, point pin[], double wmin, int nout, point pout[])
+{
+ delaunay* d = delaunay_build(nin, pin, 0, NULL, 0, NULL);
+ nnpi* nn = nnpi_create(d);
+ int seed = 0;
+ int i;
+
+ nn->wmin = wmin;
+
+ if (nn_verbose) {
+ fprintf(stderr, "xytoi:\n");
+ for (i = 0; i < nout; ++i) {
+ point* p = &pout[i];
+
+ fprintf(stderr, "(%.7g,%.7g) -> %d\n", p->x, p->y, delaunay_xytoi(d, p, seed));
+ }
+ }
+
+ for (i = 0; i < nout; ++i)
+ nnpi_interpolate_point(nn, &pout[i]);
+
+ if (nn_verbose) {
+ fprintf(stderr, "output:\n");
+ for (i = 0; i < nout; ++i) {
+ point* p = &pout[i];
+
+ fprintf(stderr, " %d:%15.7g %15.7g %15.7g\n", i, p->x, p->y, p->z);
+ }
+ }
+
+ nnpi_destroy(nn);
+ delaunay_destroy(d);
+}
+
+/* Sets minimal allowed weight for Natural Neighbours interpolation.
+ * @param nn Natural Neighbours point interpolator
+ * @param wmin Minimal allowed weight
+ */
+void nnpi_setwmin(nnpi* nn, double wmin)
+{
+ nn->wmin = wmin;
+}
+
+/* Sets point to interpolate in.
+ * @param nn Natural Neighbours point interpolator
+ * @param p Point to interpolate in
+ */
+void nnpi_set_point(nnpi* nn, point* p)
+{
+ nn->p = p;
+}
+
+/* Gets number of data points involved in current interpolation.
+ * @return Number of data points involved in current interpolation
+ */
+int nnpi_get_nvertices(nnpi* nn)
+{
+ return nn->nvertices;
+}
+
+/* Gets indices of data points involved in current interpolation.
+ * @return indices of data points involved in current interpolation
+ */
+int* nnpi_get_vertices(nnpi* nn)
+{
+ return nn->vertices;
+}
+
+/* Gets weights of data points involved in current interpolation.
+ * @return weights of data points involved in current interpolation
+ */
+double* nnpi_get_weights(nnpi* nn)
+{
+ return nn->weights;
+}
+
+/*
+ * nnhpi
+ */
+
+struct nnhpi {
+ struct nnpi* nnpi;
+ hashtable* ht_data;
+ hashtable* ht_weights;
+ int n; /* number of points processed */
+};
+
+typedef struct {
+ int nvertices;
+ int* vertices; /* vertex indices [nvertices] */
+ double* weights; /* vertex weights [nvertices] */
+} nn_weights;
+
+/* Creates Natural Neighbours hashing point interpolator.
+ *
+ * @param d Delaunay triangulation
+ * @param size Hash table size (should be of order of number of output points)
+ * @return Natural Neighbours interpolation
+ */
+nnhpi* nnhpi_create(delaunay* d, int size)
+{
+ nnhpi* nn = (nnhpi *)malloc(sizeof(nnhpi));
+ int i;
+
+ nn->nnpi = nnpi_create(d);
+
+ nn->ht_data = ht_create_d2(d->npoints);
+ nn->ht_weights = ht_create_d2(size);
+ nn->n = 0;
+
+ for (i = 0; i < d->npoints; ++i)
+ ht_insert(nn->ht_data, &d->points[i], &d->points[i]);
+
+ return nn;
+}
+
+static void free_nn_weights(void* data)
+{
+ nn_weights* weights = (nn_weights*) data;
+
+ free(weights->vertices);
+ free(weights->weights);
+ free(weights);
+}
+
+/* Destroys Natural Neighbours hashing point interpolation.
+ *
+ * @param nn Structure to be destroyed
+ */
+void nnhpi_destroy(nnhpi* nn)
+{
+ ht_destroy(nn->ht_data);
+ ht_process(nn->ht_weights, free_nn_weights);
+ ht_destroy(nn->ht_weights);
+ nnpi_destroy(nn->nnpi);
+}
+
+/* Finds Natural Neighbours-interpolated value in a point.
+ *
+ * @param nnhpi NN point hashing interpolator
+ * @param p Point to be interpolated (p->x, p->y -- input; p->z -- output)
+ */
+void nnhpi_interpolate(nnhpi* nnhpi, point* p)
+{
+ nnpi* nnpi = nnhpi->nnpi;
+ delaunay* d = nnpi->d;
+ hashtable* ht_weights = nnhpi->ht_weights;
+ nn_weights* weights;
+ int i;
+
+ if (ht_find(ht_weights, p) != NULL) {
+ weights = (nn_weights *)ht_find(ht_weights, p);
+ if (nn_verbose)
+ fprintf(stderr, " <hashtable>\n");
+ } else {
+ nnpi_reset(nnpi);
+ nnpi->p = p;
+ nnpi_calculate_weights(nnpi);
+ nnpi_normalize_weights(nnpi);
+
+ weights = (nn_weights *)malloc(sizeof(nn_weights));
+ weights->vertices = (int *)malloc(sizeof(int) * nnpi->nvertices);
+ weights->weights = (double *)malloc(sizeof(double) * nnpi->nvertices);
+
+ weights->nvertices = nnpi->nvertices;
+
+ for (i = 0; i < nnpi->nvertices; ++i) {
+ weights->vertices[i] = nnpi->vertices[i];
+ weights->weights[i] = nnpi->weights[i];
+ }
+
+ ht_insert(ht_weights, p, weights);
+
+ if (nn_verbose) {
+ if (nn_test_vertice == -1) {
+ if (nnpi->n == 0)
+ fprintf(stderr, "weights:\n");
+ fprintf(stderr, " %d: {", nnpi->n);
+
+ for (i = 0; i < nnpi->nvertices; ++i) {
+ fprintf(stderr, "(%d,%.5g)", nnpi->vertices[i], nnpi->weights[i]);
+
+ if (i < nnpi->nvertices - 1)
+ fprintf(stderr, ", ");
+ }
+ fprintf(stderr, "}\n");
+ } else {
+ double w = 0.0;
+
+ if (nnpi->n == 0)
+ fprintf(stderr, "weights for vertex %d:\n", nn_test_vertice);
+ for (i = 0; i < nnpi->nvertices; ++i) {
+ if (nnpi->vertices[i] == nn_test_vertice) {
+ w = nnpi->weights[i];
+
+ break;
+ }
+ }
+ fprintf(stderr, "%15.7g %15.7g %15.7g\n", p->x, p->y, w);
+ }
+ }
+
+ nnpi->n++;
+ }
+
+ nnhpi->n++;
+
+ if (weights->nvertices == 0) {
+ p->z = NaN;
+ return;
+ }
+
+ p->z = 0.0;
+ for (i = 0; i < weights->nvertices; ++i) {
+ if (weights->weights[i] < nnpi->wmin) {
+ p->z = NaN;
+ return;
+ }
+ p->z += d->points[weights->vertices[i]].z * weights->weights[i];
+ }
+}
+
+/* Modifies interpolated data.
+ * Finds point* pd in the underlying Delaunay triangulation such that
+ * pd->x = p->x and pd->y = p->y, and copies p->z to pd->z. Exits with error
+ * if the point is not found.
+ *
+ * @param nnhpi Natural Neighbours hashing point interpolator
+ * @param p New data
+ */
+void nnhpi_modify_data(nnhpi* nnhpi, point* p)
+{
+ point* orig = (point *)ht_find(nnhpi->ht_data, p);
+
+ assert(orig != NULL);
+ orig->z = p->z;
+}
+
+/* Sets minimal allowed weight for Natural Neighbours interpolation.
+ * @param nn Natural Neighbours point hashing interpolator
+ * @param wmin Minimal allowed weight
+ */
+void nnhpi_setwmin(nnhpi* nn, double wmin)
+{
+ nn->nnpi->wmin = wmin;
+}
+
+#if defined(NNPHI_TEST)
+
+#include <sys/time.h>
+
+#define NPOINTSIN 10000
+#define NMIN 10
+#define NX 101
+#define NXMIN 1
+
+#define SQ(x) ((x) * (x))
+
+static double franke(double x, double y)
+{
+ x *= 9.0;
+ y *= 9.0;
+ return 0.75 * exp((-SQ(x - 2.0) - SQ(y - 2.0)) / 4.0)
+ + 0.75 * exp(-SQ(x - 2.0) / 49.0 - (y - 2.0) / 10.0)
+ + 0.5 * exp((-SQ(x - 7.0) - SQ(y - 3.0)) / 4.0)
+ - 0.2 * exp(-SQ(x - 4.0) - SQ(y - 7.0));
+}
+
+static void usage()
+{
+ printf("Usage: nnhpi_test [-a] [-n <nin> <nxout>] [-v|-V]\n");
+ printf("Options:\n");
+ printf(" -a -- use non-Sibsonian interpolation rule\n");
+ printf(" -n <nin> <nout>:\n");
+ printf(" <nin> -- number of input points (default = 10000)\n");
+ printf(" <nout> -- number of output points per side (default = 64)\n");
+ printf(" -v -- verbose\n");
+ printf(" -V -- very verbose\n");
+
+ exit(0);
+}
+
+int main(int argc, char* argv[])
+{
+ int nin = NPOINTSIN;
+ int nx = NX;
+ int nout = 0;
+ point* pin = NULL;
+ delaunay* d = NULL;
+ point* pout = NULL;
+ nnhpi* nn = NULL;
+ int cpi = -1; /* control point index */
+ struct timeval tv0, tv1;
+ struct timezone tz;
+ int i;
+
+ i = 1;
+ while (i < argc) {
+ switch (argv[i][1]) {
+ case 'a':
+ i++;
+ nn_rule = NON_SIBSONIAN;
+ break;
+ case 'n':
+ i++;
+ if (i >= argc)
+ nn_quit("no number of data points found after -n\n");
+ nin = atoi(argv[i]);
+ i++;
+ if (i >= argc)
+ nn_quit("no number of ouput points per side found after -i\n");
+ nx = atoi(argv[i]);
+ i++;
+ break;
+ case 'v':
+ i++;
+ nn_verbose = 1;
+ break;
+ case 'V':
+ i++;
+ nn_verbose = 2;
+ break;
+ default:
+ usage();
+ break;
+ }
+ }
+
+ if (nin < NMIN)
+ nin = NMIN;
+ if (nx < NXMIN)
+ nx = NXMIN;
+
+ printf("\nTest of Natural Neighbours hashing point interpolator:\n\n");
+ printf(" %d data points\n", nin);
+ printf(" %d output points\n", nx * nx);
+
+ /*
+ * generate data
+ */
+ printf(" generating data:\n");
+ fflush(stdout);
+ pin = (point *)malloc(nin * sizeof(point));
+ for (i = 0; i < nin; ++i) {
+ point* p = &pin[i];
+
+ p->x = (double) random() / RAND_MAX;
+ p->y = (double) random() / RAND_MAX;
+ p->z = franke(p->x, p->y);
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ /*
+ * triangulate
+ */
+ printf(" triangulating:\n");
+ fflush(stdout);
+ d = delaunay_build(nin, pin, 0, NULL, 0, NULL);
+
+ /*
+ * generate output points
+ */
+ points_generate2(-0.1, 1.1, -0.1, 1.1, nx, nx, &nout, &pout);
+ cpi = (nx / 2) * (nx + 1);
+
+ gettimeofday(&tv0, &tz);
+
+ /*
+ * create interpolator
+ */
+ printf(" creating interpolator:\n");
+ fflush(stdout);
+ nn = nnhpi_create(d, nout);
+
+ fflush(stdout);
+ gettimeofday(&tv1, &tz);
+ {
+ long dt = 1000000 * (tv1.tv_sec - tv0.tv_sec) + tv1.tv_usec - tv0.tv_usec;
+
+ printf(" interpolator creation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout);
+ }
+
+ /*
+ * interpolate
+ */
+ printf(" interpolating:\n");
+ fflush(stdout);
+ gettimeofday(&tv1, &tz);
+ for (i = 0; i < nout; ++i) {
+ point* p = &pout[i];
+
+ nnhpi_interpolate(nn, p);
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ fflush(stdout);
+ gettimeofday(&tv0, &tz);
+ {
+ long dt = 1000000.0 * (tv0.tv_sec - tv1.tv_sec) + tv0.tv_usec - tv1.tv_usec;
+
+ printf(" interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout);
+ }
+
+ if (!nn_verbose)
+ printf(" control point: (%f, %f, %f) (expected z = %f)\n", pout[cpi].x, pout[cpi].y, pout[cpi].z, franke(pout[cpi].x, pout[cpi].y));
+
+ printf(" interpolating one more time:\n");
+ fflush(stdout);
+ gettimeofday(&tv0, &tz);
+ for (i = 0; i < nout; ++i) {
+ point* p = &pout[i];
+
+ nnhpi_interpolate(nn, p);
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ fflush(stdout);
+ gettimeofday(&tv1, &tz);
+ {
+ long dt = 1000000.0 * (tv1.tv_sec - tv0.tv_sec) + tv1.tv_usec - tv0.tv_usec;
+
+ printf(" interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout);
+ }
+
+ if (!nn_verbose)
+ printf(" control point: (%f, %f, %f) (expected z = %f)\n", pout[cpi].x, pout[cpi].y, pout[cpi].z, franke(pout[cpi].x, pout[cpi].y));
+
+ printf(" entering new data:\n");
+ fflush(stdout);
+ for (i = 0; i < nin; ++i) {
+ point* p = &pin[i];
+
+ p->z = p->x * p->x - p->y * p->y;
+ nnhpi_modify_data(nn, p);
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ printf(" interpolating:\n");
+ fflush(stdout);
+ gettimeofday(&tv1, &tz);
+ for (i = 0; i < nout; ++i) {
+ point* p = &pout[i];
+
+ nnhpi_interpolate(nn, p);
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ fflush(stdout);
+ gettimeofday(&tv0, &tz);
+ {
+ long dt = 1000000.0 * (tv0.tv_sec - tv1.tv_sec) + tv0.tv_usec - tv1.tv_usec;
+
+ printf(" interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout);
+ }
+
+ if (!nn_verbose)
+ printf(" control point: (%f, %f, %f) (expected z = %f)\n", pout[cpi].x, pout[cpi].y, pout[cpi].z, pout[cpi].x * pout[cpi].x - pout[cpi].y * pout[cpi].y);
+
+ printf(" restoring data:\n");
+ fflush(stdout);
+ for (i = 0; i < nin; ++i) {
+ point* p = &pin[i];
+
+ p->z = franke(p->x, p->y);
+ nnhpi_modify_data(nn, p);
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ printf(" interpolating:\n");
+ fflush(stdout);
+ gettimeofday(&tv0, &tz);
+ for (i = 0; i < nout; ++i) {
+ point* p = &pout[i];
+
+ nnhpi_interpolate(nn, p);
+ if (nn_verbose)
+ printf(" (%f, %f, %f)\n", p->x, p->y, p->z);
+ }
+
+ fflush(stdout);
+ gettimeofday(&tv1, &tz);
+ {
+ long dt = 1000000.0 * (tv1.tv_sec - tv0.tv_sec) + tv1.tv_usec - tv0.tv_usec;
+
+ printf(" interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout);
+ }
+
+ if (!nn_verbose)
+ printf(" control point: (%f, %f, %f) (expected z = %f)\n", pout[cpi].x, pout[cpi].y, pout[cpi].z, franke(pout[cpi].x, pout[cpi].y));
+
+ printf(" hashtable stats:\n");
+ fflush(stdout);
+ {
+ hashtable* ht = nn->ht_data;
+
+ printf(" input points: %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash);
+ ht = nn->ht_weights;
+ printf(" weights: %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash);
+ }
+ printf("\n");
+
+ nnhpi_destroy(nn);
+ free(pout);
+ delaunay_destroy(d);
+ free(pin);
+
+ return 0;
+}
+
+#endif
diff --git a/src/modules/grid/grid_gridding/nn/version.h b/src/modules/grid/grid_gridding/nn/version.h
new file mode 100644
index 0000000..2205ddc
--- /dev/null
+++ b/src/modules/grid/grid_gridding/nn/version.h
@@ -0,0 +1,19 @@
+/******************************************************************************
+ *
+ * File: version.h
+ *
+ * Created: 18/10/2001
+ *
+ * Author: Pavel Sakov
+ * CSIRO Marine Research
+ *
+ * Purpose: Contains version string
+ *
+ *****************************************************************************/
+
+#if !defined(_VERSION_H)
+#define _VERSION_H
+
+char* nn_version = "1.38";
+
+#endif
--
Saga GIS
More information about the Pkg-grass-devel
mailing list