[GRASS-SVN] r65511 - grass-addons/grass7/vector/v.nnstat
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 22 01:03:48 PDT 2015
Author: evas
Date: 2015-06-22 01:03:48 -0700 (Mon, 22 Jun 2015)
New Revision: 65511
Added:
grass-addons/grass7/vector/v.nnstat/getval.c
grass-addons/grass7/vector/v.nnstat/mbb.c
grass-addons/grass7/vector/v.nnstat/mbr.c
grass-addons/grass7/vector/v.nnstat/nearest_neighbour.c
grass-addons/grass7/vector/v.nnstat/spatial_index.c
grass-addons/grass7/vector/v.nnstat/utils.c
Log:
v.nnstat: kd-tree (PCL library) replaced by R-tree (GRASS GIS)
Added: grass-addons/grass7/vector/v.nnstat/getval.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/getval.c (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/getval.c 2015-06-22 08:03:48 UTC (rev 65511)
@@ -0,0 +1,161 @@
+#include "local_proto.h"
+
+
+/* get array of values from attribute column (function based on part of v.buffer2 (Radim Blazek, Rosen Matev)) */
+double *get_col_values(struct Map_info *map, int field, const char *column)
+{
+ int i, nrec, ctype;
+ struct field_info *Fi;
+
+ dbCatValArray cvarr;
+ dbDriver *Driver;
+
+ double *values, *vals;
+
+ db_CatValArray_init(&cvarr); /* array of categories and values initialised */
+
+ Fi = Vect_get_field(map, field); /* info about call of DB */
+ if (Fi == NULL)
+ G_fatal_error(_("Database connection not defined for layer %d"),
+ field);
+ Driver = db_start_driver_open_database(Fi->driver, Fi->database); /* started connection to DB */
+ if (Driver == NULL)
+ G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+ Fi->database, Fi->driver);
+
+ /* Note do not check if the column exists in the table because it may be expression */
+
+ /* TODO: only select values we need instead of all in column */
+
+ /* Select pairs key/value to array, values are sorted by key (must be integer) */
+ nrec =
+ db_select_CatValArray(Driver, Fi->table, Fi->key, column, NULL,
+ &cvarr);
+ if (nrec < 0)
+ G_fatal_error(_("Unable to select data from table <%s>"), Fi->table);
+ G_message(_("Reading elevations from attribute table: %d records selected"), nrec);
+ ctype = cvarr.ctype;
+
+ db_close_database_shutdown_driver(Driver);
+
+ values = (double *) G_malloc(cvarr.n_values * sizeof(double));
+ vals = &values[0];
+
+ for (i = 0; i < cvarr.n_values; i++) {
+ if (ctype == DB_C_TYPE_INT) {
+ *vals = (double) cvarr.value[i].val.i;
+ }
+ if (ctype == DB_C_TYPE_DOUBLE) {
+ *vals = cvarr.value[i].val.d;
+ }
+ if (ctype == DB_C_TYPE_STRING) {
+ G_fatal_error(_("The column must be numeric..."));
+ }
+ vals++;
+ }
+
+ return values;
+}
+
+/* get coordinates of input points */
+void read_points(struct Map_info *map, int field, struct nna_par *xD,
+ const char *zcol, struct points *pnts)
+{
+ int ctrl, n, type, nskipped, pass;
+
+ struct line_pnts *Points; // structures to hold line *Points (map)
+ double *rx, *ry, *rz; // pointers to the coordinates
+ double *z_attr, *z; // pointers to attribute z values
+
+ Points = Vect_new_line_struct();
+ n = pnts->n = Vect_get_num_primitives(map, GV_POINT); /* topology required */
+ pnts->r = (double *) G_malloc(n*3 * sizeof(double));
+
+ rx = &pnts->r[0];
+ ry = &pnts->r[1];
+ rz = &pnts->r[2];
+
+ pnts->R_tree = create_spatial_index(xD); // create spatial index (R-tree)
+
+ /* Get 3rd coordinate of 2D points from attribute column -> 3D interpolation */
+ if (xD->v3 == FALSE && zcol != NULL) { // 2D input layer with z attribute column:
+ xD->zcol = (char *) G_malloc(strlen(zcol) * sizeof(char));
+ strcpy(xD->zcol, zcol);
+ z_attr = get_col_values(map, field, zcol); // read attribute z values
+ z = &z_attr[0];
+ }
+ else {
+ xD->zcol = NULL;
+ }
+
+ nskipped = ctrl = pass = 0;
+
+ while (TRUE) {
+ type = Vect_read_next_line(map, Points, NULL);
+ if (type == -1) {
+ G_fatal_error(_("Unable to read vector map"));
+ }
+ if (type == -2) {
+ break;
+ }
+
+ if (type != GV_POINT) {
+ nskipped++;
+ continue;
+ }
+
+ *rx = Points->x[0];
+ *ry = Points->y[0];
+
+ // 3D points or 2D points without attribute column -> 2D interpolation
+ if (xD->zcol == NULL) { // z attribute column not available:
+ if (xD->v3 == TRUE && xD->i3 == FALSE) { // 2D NNA:
+ *rz = 0.;
+ }
+ else { // 3D NNA:
+ *rz = Points->z[0];
+ }
+ }
+ else { // z attribute column available:
+ *rz = *z;
+ z++;
+ }
+ insert_rectangle(xD->i3, ctrl, pnts); // insert rectangle into spatial index
+
+ /* Find extends */
+ if (ctrl == 0) {
+ pnts->r_min = pnts->r_max = triple(*rx, *ry, *rz);
+ }
+ else {
+ pnts->r_min = triple(MIN(*pnts->r_min, *rx), MIN(*(pnts->r_min+1), *ry), MIN(*(pnts->r_min+2), *rz));
+ pnts->r_max = triple(MAX(*pnts->r_max, *rx), MAX(*(pnts->r_max+1), *ry), MAX(*(pnts->r_max+2), *rz));
+ }
+ if (ctrl < pnts->n - 1) {
+ rx += 3;
+ ry += 3;
+ rz += 3;
+ }
+ ctrl++;
+ }
+
+ if (nskipped > 0) {
+ G_warning(_("%d features skipped, only points are accepted"),
+ nskipped);
+ }
+
+ Vect_destroy_line_struct(Points);
+
+ if (ctrl <= 0) {
+ G_fatal_error(_("Unable to read coordinates of point layer"));
+ }
+
+ // Compute maximum distance of the point sample
+ double dx = pnts->r_max[0] - pnts->r_min[0];
+ double dy = pnts->r_max[1] - pnts->r_min[1];
+ double dz = pnts->r_max[2] - pnts->r_min[2];
+ pnts->max_dist = xD->i3 == TRUE ? sqrt(dx*dx + dy*dy + dz*dz) : sqrt(dx*dx + dy*dy);
+
+ G_message(_("Input coordinates have been read..."));
+
+ return;
+}
Added: grass-addons/grass7/vector/v.nnstat/mbb.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/mbb.c (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/mbb.c 2015-06-22 08:03:48 UTC (rev 65511)
@@ -0,0 +1,760 @@
+#include "local_proto.h"
+#include "macros.h"
+#include "mbb.h"
+
+/*******************************
+These functions are mostly taken from the module v.hull (Aime, A., Neteler, M., Ducke, B., Landa, M.)
+
+Minimum Bounding Block (MBB)
+ - obtain vertices of 3D convex hull,
+ - transform coordinates of vertices into coordinate system with axes parallel to hull's edges,
+ - find extents,
+ - compute volumes and find minimum of them.
+ *******************************
+ */
+
+/* Release all memory allocated for edges, faces and vertices */
+void freeMem(void)
+{
+ tEdge e; /* Primary index into edge list. */
+ tFace f; /* Primary pointer into face list. */
+ tVertex v;
+ tEdge te; /* Temporary edge pointer. */
+ tFace tf; /* Temporary face pointer. */
+ tVertex tv; /* Temporary vertex pointer. */
+
+ e = edges;
+ do {
+ te = e;
+ e = e->next;
+ DELETE(edges, te);
+ } while (e != edges);
+
+ f = faces;
+ do {
+ tf = f;
+ f = f->next;
+ DELETE(faces, tf);
+ } while (f != faces);
+
+ v = vertices;
+ do {
+ tv = v;
+ v = v->next;
+ DELETE(vertices, tv);
+ } while (v != vertices);
+
+ FREE(te);
+ FREE(tf);
+ FREE(tv);
+
+ DELETE(edges, e);
+ DELETE(faces, f);
+ DELETE(vertices, v);
+
+ FREE(edges);
+ FREE(faces);
+ FREE(vertices);
+
+ return;
+}
+
+
+/*-------------------------------------------------------------------*/
+int make3DHull(struct points *pnts, struct convex *hull)
+{
+ int error;
+
+ ReadVertices(pnts);
+
+ error = DoubleTriangle();
+ if (error < 0) {
+ G_fatal_error
+ ("All points of 3D input map are in the same plane.\n Cannot create a 3D hull.");
+ }
+
+ ConstructHull();
+
+ write_coord_faces(pnts, hull);
+
+ freeMem();
+
+ return (0);
+}
+
+/*---------------------------------------------------------------------
+ MakeNullVertex: Makes a vertex, nulls out fields.
+ ---------------------------------------------------------------------*/
+tVertex MakeNullVertex(void)
+{
+ tVertex v;
+
+ NEW(v, tsVertex); /* If memory not allocated => out of memory*/
+ v->duplicate = NULL;
+ v->onhull = !ONHULL;
+ v->mark = !PROCESSED;
+ ADD(vertices, v);
+
+ return v;
+}
+
+/*---------------------------------------------------------------------
+ ReadVertices: Reads in the vertices, and links them into a circular
+ list with MakeNullVertex. There is no need for the # of vertices to be
+ the first line: the function looks for EOF instead. Sets the global
+ variable vertices via the ADD macro.
+ ---------------------------------------------------------------------*/
+void ReadVertices(struct points *pnts)
+{
+ tVertex v;
+ int vnum = 0;
+ int i;
+
+ double *r;
+ r = &pnts->r[0];
+
+ G_message(_("Reading 3D vertices..."));
+ for (i = 0; i < pnts->n; i++) {
+ v = MakeNullVertex();
+ v->v[X] = *r;
+ v->v[Y] = *(r+1);
+ v->v[Z] = *(r+2);
+ v->vnum = vnum++;
+ r +=3;
+ G_percent(i, (pnts->n - 1), 1);
+ }
+ fflush(stdout);
+}
+
+/*---------------------------------------------------------------------
+ Outputs coordinates of 3D hull to matrix [n x 3]
+ ---------------------------------------------------------------------*/
+void write_coord_faces(struct points *pnts, struct convex *hull)
+{
+ int n = pnts->n;
+
+ /* Pointers to vertices, edges, faces. */
+ tVertex v;
+ tFace f;
+ int nv = 0, nf = 0;
+ double *hc, *hf;
+
+ hull->coord = (double *) malloc(n * 3 * sizeof(double));
+ hull->faces = (double *) malloc(3*n * 3 * sizeof(double));
+ hc = &hull->coord[0];
+ hf = &hull->faces[0];
+
+ v = vertices;
+ f = faces;
+
+ do {
+ /* Write vertex coordinates */
+ *hc = ((double)(v->v[X]));
+ *(hc+1) = ((double)(v->v[Y]));
+ *(hc+2) = ((double)(v->v[Z]));
+
+ v = v->next;
+ hc += 3;
+ nv += 3;
+
+ } while (v!=vertices);
+
+ do {
+ /* write one triangular face */
+ *hf = ((double)(f->vertex[0]->v[X]));
+ *(hf+1) = ((double)(f->vertex[0]->v[Y]));
+ *(hf+2) = ((double)(f->vertex[0]->v[Z]));
+
+ *(hf+3) = ((double)(f->vertex[1]->v[X]));
+ *(hf+4) = ((double)(f->vertex[1]->v[Y]));
+ *(hf+5) = ((double)(f->vertex[1]->v[Z]));
+
+ *(hf+6) = ((double)(f->vertex[2]->v[X]));
+ *(hf+7) = ((double)(f->vertex[2]->v[Y]));
+ *(hf+8) = ((double)(f->vertex[2]->v[Z]));
+
+ f = f->next;
+ hf += 9;
+ nf += 9;
+
+ } while (f!=faces);
+
+ /* reclaim uneeded memory */
+ hull->n = nv;
+ hull->n_faces = nf;
+ hull->coord = (double *) G_realloc(hull->coord, nv*3 * sizeof(double));
+ hull->faces = (double *) G_realloc(hull->faces, nf*3 * sizeof(double));
+
+ fflush(stdout);
+}
+
+
+/*---------------------------------------------------------------------
+ DoubleTriangle builds the initial double triangle. It first finds 3
+ noncollinear points and makes two faces out of them, in opposite order.
+ It then finds a fourth point that is not coplanar with that face. The
+ vertices are stored in the face structure in counterclockwise order so
+ that the volume between the face and the point is negative. Lastly, the
+ 3 newfaces to the fourth point are constructed and the data structures
+ are cleaned up.
+ ---------------------------------------------------------------------*/
+
+/* RETURN: 0 if OK */
+/* -1 if all points collinear */
+/* -2 if all points coplanar */
+
+int DoubleTriangle(void)
+{
+ tVertex v0, v1, v2, v3;
+ tFace f0, f1 = NULL;
+ long int vol;
+
+ /* Find 3 noncollinear points. */
+ v0 = vertices;
+ while (Collinear(v0, v0->next, v0->next->next)) {
+ if ((v0 = v0->next) == vertices) {
+ G_warning("DoubleTriangle: All points are collinear!\n");
+ return (-1);
+ }
+ }
+ v1 = v0->next;
+ v2 = v1->next;
+
+ /* Mark the vertices as processed. */
+ v0->mark = PROCESSED;
+ v1->mark = PROCESSED;
+ v2->mark = PROCESSED;
+
+ /* Create the two "twin" faces. */
+ f0 = MakeFace(v0, v1, v2, f1);
+ f1 = MakeFace(v2, v1, v0, f0);
+
+ /* Link adjacent face fields. */
+ f0->edge[0]->adjface[1] = f1;
+ f0->edge[1]->adjface[1] = f1;
+ f0->edge[2]->adjface[1] = f1;
+ f1->edge[0]->adjface[1] = f0;
+ f1->edge[1]->adjface[1] = f0;
+ f1->edge[2]->adjface[1] = f0;
+
+ /* Find a fourth, noncoplanar point to form tetrahedron. */
+ v3 = v2->next;
+ vol = VolumeSign(f0, v3);
+ while (!vol) {
+ if ((v3 = v3->next) == v0) {
+ G_warning("DoubleTriangle: All points are coplanar!\n");
+ return (-2);
+ }
+ vol = VolumeSign(f0, v3);
+ }
+
+ /* Insure that v3 will be the first added. */
+ vertices = v3;
+
+ return (0);
+}
+
+
+/*---------------------------------------------------------------------
+ ConstructHull adds the vertices to the hull one at a time. The hull
+ vertices are those in the list marked as onhull.
+ ---------------------------------------------------------------------*/
+int ConstructHull(void)
+{
+ tVertex v, vnext;
+ bool changed; /* T if addition changes hull; not used. */
+ int i;
+ int numVertices;
+
+
+ G_message(_("Constructing 3D hull..."));
+
+ v = vertices;
+ i = 0;
+ do {
+ vnext = v->next;
+ v = vnext;
+ i++;
+ } while (v != vertices);
+ numVertices = i;
+
+ v = vertices;
+ i = 0;
+ do {
+ vnext = v->next;
+ if (!v->mark) {
+ v->mark = PROCESSED;
+ changed = AddOne(v);
+ CleanUp();
+ }
+ v = vnext;
+ i++;
+
+ G_percent(i, numVertices, 1);
+
+ } while (v != vertices);
+
+ fflush(stdout);
+
+ return numVertices;
+
+}
+
+/*---------------------------------------------------------------------
+ AddOne is passed a vertex. It first determines all faces visible from
+ that point. If none are visible then the point is marked as not
+ onhull. Next is a loop over edges. If both faces adjacent to an edge
+ are visible, then the edge is marked for deletion. If just one of the
+ adjacent faces is visible then a new face is constructed.
+ ---------------------------------------------------------------------*/
+bool AddOne(tVertex p)
+{
+ tFace f;
+ tEdge e, temp;
+ long int vol;
+ bool vis = FALSE;
+
+
+ /* Mark faces visible from p. */
+ f = faces;
+ do {
+ vol = VolumeSign(f, p);
+
+ if (vol < 0) {
+ f->visible = VISIBLE;
+ vis = TRUE;
+ }
+ f = f->next;
+ } while (f != faces);
+
+ /* If no faces are visible from p, then p is inside the hull. */
+ if (!vis) {
+ p->onhull = !ONHULL;
+ return FALSE;
+ }
+
+ /* Mark edges in interior of visible region for deletion.
+ Erect a newface based on each border edge. */
+ e = edges;
+ do {
+ temp = e->next;
+ if (e->adjface[0]->visible && e->adjface[1]->visible)
+ /* e interior: mark for deletion. */
+ e->del = REMOVED;
+ else if (e->adjface[0]->visible || e->adjface[1]->visible)
+ /* e border: make a new face. */
+ e->newface = MakeConeFace(e, p);
+ e = temp;
+ } while (e != edges);
+ return TRUE;
+}
+
+/*---------------------------------------------------------------------
+ VolumeSign returns the sign of the volume of the tetrahedron determined by f
+ and p. VolumeSign is +1 iff p is on the negative side of f,
+ where the positive side is determined by the rh-rule. So the volume
+ is positive if the ccw normal to f points outside the tetrahedron.
+ The final fewer-multiplications form is due to Bob Williamson.
+ ---------------------------------------------------------------------*/
+int VolumeSign(tFace f, tVertex p)
+{
+ double vol;
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
+
+ ax = f->vertex[0]->v[X] - p->v[X];
+ ay = f->vertex[0]->v[Y] - p->v[Y];
+ az = f->vertex[0]->v[Z] - p->v[Z];
+ bx = f->vertex[1]->v[X] - p->v[X];
+ by = f->vertex[1]->v[Y] - p->v[Y];
+ bz = f->vertex[1]->v[Z] - p->v[Z];
+ cx = f->vertex[2]->v[X] - p->v[X];
+ cy = f->vertex[2]->v[Y] - p->v[Y];
+ cz = f->vertex[2]->v[Z] - p->v[Z];
+
+ vol = ax * (by * cz - bz * cy)
+ + ay * (bz * cx - bx * cz)
+ + az * (bx * cy - by * cx);
+
+ /* The volume should be an integer. */
+ if (vol > 0.0)
+ return 1;
+ else if (vol < -0.0)
+ return -1;
+ else
+ return 0;
+}
+
+
+/*---------------------------------------------------------------------
+ MakeConeFace makes a new face and two new edges between the
+ edge and the point that are passed to it. It returns a pointer to
+ the new face.
+ ---------------------------------------------------------------------*/
+tFace MakeConeFace(tEdge e, tVertex p)
+{
+ tEdge new_edge[2];
+ tFace new_face;
+ int i, j;
+
+ /* Make two new edges (if don't already exist). */
+ for (i = 0; i < 2; ++i)
+ /* If the edge exists, copy it into new_edge. */
+ if (!(new_edge[i] = e->endpts[i]->duplicate)) {
+ /* Otherwise (duplicate is NULL), MakeNullEdge. */
+ new_edge[i] = MakeNullEdge();
+ new_edge[i]->endpts[0] = e->endpts[i];
+ new_edge[i]->endpts[1] = p;
+ e->endpts[i]->duplicate = new_edge[i];
+ }
+
+ /* Make the new face. */
+ new_face = MakeNullFace();
+ new_face->edge[0] = e;
+ new_face->edge[1] = new_edge[0];
+ new_face->edge[2] = new_edge[1];
+ MakeCcw(new_face, e, p);
+
+ /* Set the adjacent face pointers. */
+ for (i = 0; i < 2; ++i)
+ for (j = 0; j < 2; ++j)
+ /* Only one NULL link should be set to new_face. */
+ if (!new_edge[i]->adjface[j]) {
+ new_edge[i]->adjface[j] = new_face;
+ break;
+ }
+
+ return new_face;
+}
+
+/*---------------------------------------------------------------------
+ MakeCcw puts the vertices in the face structure in counterclock wise
+ order. We want to store the vertices in the same
+ order as in the visible face. The third vertex is always p.
+ ---------------------------------------------------------------------*/
+void MakeCcw(tFace f, tEdge e, tVertex p)
+{
+ tFace fv; /* The visible face adjacent to e */
+ int i; /* Index of e->endpoint[0] in fv. */
+ tEdge s; /* Temporary, for swapping */
+
+ if (e->adjface[0]->visible)
+ fv = e->adjface[0];
+ else
+ fv = e->adjface[1];
+
+ /* Set vertex[0] & [1] of f to have the same orientation
+ as do the corresponding vertices of fv. */
+ for (i = 0; fv->vertex[i] != e->endpts[0]; ++i) ;
+ /* Orient f the same as fv. */
+ if (fv->vertex[(i + 1) % 3] != e->endpts[1]) {
+ f->vertex[0] = e->endpts[1];
+ f->vertex[1] = e->endpts[0];
+ }
+ else {
+ f->vertex[0] = e->endpts[0];
+ f->vertex[1] = e->endpts[1];
+ SWAP(s, f->edge[1], f->edge[2]);
+ }
+ /* This swap is tricky. e is edge[0]. edge[1] is based on endpt[0],
+ edge[2] on endpt[1]. So if e is oriented "forwards," we
+ need to move edge[1] to follow [0], because it precedes. */
+
+ f->vertex[2] = p;
+}
+
+/*---------------------------------------------------------------------
+ MakeNullEdge creates a new cell and initializes all pointers to NULL
+ and sets all flags to off. It returns a pointer to the empty cell.
+ ---------------------------------------------------------------------*/
+tEdge MakeNullEdge(void)
+{
+ tEdge e;
+
+ NEW(e, tsEdge);
+ e->adjface[0] = e->adjface[1] = e->newface = NULL;
+ e->endpts[0] = e->endpts[1] = NULL;
+ e->del = !REMOVED;
+ ADD(edges, e);
+ return e;
+}
+
+/*--------------------------------------------------------------------
+ MakeNullFace creates a new face structure and initializes all of its
+ flags to NULL and sets all the flags to off. It returns a pointer
+ to the empty cell.
+ ---------------------------------------------------------------------*/
+tFace MakeNullFace(void)
+{
+ tFace f;
+ int i;
+
+ NEW(f, tsFace);
+ for (i = 0; i < 3; ++i) {
+ f->edge[i] = NULL;
+ f->vertex[i] = NULL;
+ }
+ f->visible = !VISIBLE;
+ ADD(faces, f);
+ return f;
+}
+
+/*---------------------------------------------------------------------
+ MakeFace creates a new face structure from three vertices (in ccw
+ order). It returns a pointer to the face.
+ ---------------------------------------------------------------------*/
+tFace MakeFace(tVertex v0, tVertex v1, tVertex v2, tFace fold)
+{
+ tFace f;
+ tEdge e0, e1, e2;
+
+ /* Create edges of the initial triangle. */
+ if (!fold) {
+ e0 = MakeNullEdge();
+ e1 = MakeNullEdge();
+ e2 = MakeNullEdge();
+ }
+ else { /* Copy from fold, in reverse order. */
+ e0 = fold->edge[2];
+ e1 = fold->edge[1];
+ e2 = fold->edge[0];
+ }
+ e0->endpts[0] = v0;
+ e0->endpts[1] = v1;
+ e1->endpts[0] = v1;
+ e1->endpts[1] = v2;
+ e2->endpts[0] = v2;
+ e2->endpts[1] = v0;
+
+ /* Create face for triangle. */
+ f = MakeNullFace();
+ f->edge[0] = e0;
+ f->edge[1] = e1;
+ f->edge[2] = e2;
+ f->vertex[0] = v0;
+ f->vertex[1] = v1;
+ f->vertex[2] = v2;
+
+ /* Link edges to face. */
+ e0->adjface[0] = e1->adjface[0] = e2->adjface[0] = f;
+
+ return f;
+}
+
+/*---------------------------------------------------------------------
+ CleanUp goes through each data structure list and clears all
+ flags and NULLs out some pointers. The order of processing
+ (edges, faces, vertices) is important.
+ ---------------------------------------------------------------------*/
+void CleanUp(void)
+{
+ CleanEdges();
+ CleanFaces();
+ CleanVertices();
+
+ return;
+}
+
+/*---------------------------------------------------------------------
+ CleanEdges runs through the edge list and cleans up the structure.
+ If there is a newface then it will put that face in place of the
+ visible face and NULL out newface. It also deletes so marked edges.
+ ---------------------------------------------------------------------*/
+void CleanEdges(void)
+{
+ tEdge e; /* Primary index into edge list. */
+ tEdge t; /* Temporary edge pointer. */
+
+ /* Integrate the newface's into the data structure. */
+ /* Check every edge. */
+ e = edges;
+ do {
+ if (e->newface) {
+ if (e->adjface[0]->visible)
+ e->adjface[0] = e->newface;
+ else
+ e->adjface[1] = e->newface;
+ e->newface = NULL;
+ }
+ e = e->next;
+ } while (e != edges);
+
+ /* Delete any edges marked for deletion. */
+ while (edges && edges->del) {
+ e = edges;
+ DELETE(edges, e);
+ }
+ e = edges->next;
+ do {
+ if (e->del) {
+ t = e;
+ e = e->next;
+ DELETE(edges, t);
+ }
+ else
+ e = e->next;
+ } while (e != edges);
+
+ return;
+}
+
+/*---------------------------------------------------------------------
+ CleanFaces runs through the face list and deletes any face marked visible.
+ ---------------------------------------------------------------------*/
+void CleanFaces(void)
+{
+ tFace f; /* Primary pointer into face list. */
+ tFace t; /* Temporary pointer, for deleting. */
+
+
+ while (faces && faces->visible) {
+ f = faces;
+ DELETE(faces, f);
+ }
+ f = faces->next;
+ do {
+ if (f->visible) {
+ t = f;
+ f = f->next;
+ DELETE(faces, t);
+ }
+ else
+ f = f->next;
+ } while (f != faces);
+
+ return;
+}
+
+/*---------------------------------------------------------------------
+ CleanVertices runs through the vertex list and deletes the
+ vertices that are marked as processed but are not incident to any
+ undeleted edges.
+ ---------------------------------------------------------------------*/
+void CleanVertices(void)
+{
+ tEdge e;
+ tVertex v, t;
+
+ /* Mark all vertices incident to some undeleted edge as on the hull. */
+ e = edges;
+ do {
+ e->endpts[0]->onhull = e->endpts[1]->onhull = ONHULL;
+ e = e->next;
+ } while (e != edges);
+
+ /* Delete all vertices that have been processed but
+ are not on the hull. */
+ while (vertices && vertices->mark && !vertices->onhull) {
+ v = vertices;
+ DELETE(vertices, v);
+ }
+ v = vertices->next;
+ do {
+ if (v->mark && !v->onhull) {
+ t = v;
+ v = v->next;
+ DELETE(vertices, t)
+ }
+ else
+ v = v->next;
+ } while (v != vertices);
+
+ /* Reset flags. */
+ v = vertices;
+ do {
+ v->duplicate = NULL;
+ v->onhull = !ONHULL;
+ v = v->next;
+ } while (v != vertices);
+
+ return;
+}
+
+/*---------------------------------------------------------------------
+ Collinear checks to see if the three points given are collinear,
+ by checking to see if each element of the cross product is zero.
+ ---------------------------------------------------------------------*/
+bool Collinear(tVertex a, tVertex b, tVertex c)
+{
+ return
+ (c->v[Z] - a->v[Z]) * (b->v[Y] - a->v[Y]) -
+ (b->v[Z] - a->v[Z]) * (c->v[Y] - a->v[Y]) == 0
+ && (b->v[Z] - a->v[Z]) * (c->v[X] - a->v[X]) -
+ (b->v[X] - a->v[X]) * (c->v[Z] - a->v[Z]) == 0
+ && (b->v[X] - a->v[X]) * (c->v[Y] - a->v[Y]) -
+ (b->v[Y] - a->v[Y]) * (c->v[X] - a->v[X]) == 0;
+}
+
+void convexHull3d(struct points *pnts, struct convex *hull)
+{
+ int error;
+
+ error = make3DHull(pnts, hull); /* make 3D hull */
+ if (error < 0) {
+ G_fatal_error(_("Simple planar hulls not implemented yet"));
+ }
+
+ return;
+}
+
+ /* ----------------------------
+ * MBR area estimation
+ */
+double MBB(struct points *pnts)
+{
+ int i, k;
+ double usx, usy, cosusx, sinusx, cosusy, sinusy, V, V_min;
+ double *hc_k, *hf, *ht; // pointers to hull_trans and hull->faces
+ double *r_min, *r_max;
+ r_min = (double *) G_malloc(3 * sizeof(double));
+ r_max = (double *) G_malloc(3 * sizeof(double));
+
+ double *hull_trans; /* Coordinates of hull vertices transformed into coordinate system with axes parallel to hull's edges */
+ hull_trans = (double *) malloc(3 * sizeof(double));
+ ht = &hull_trans[0];
+
+ struct convex hull;
+ convexHull3d(pnts, &hull);
+ hc_k = &hull.coord[0];
+ hf = &hull.faces[0];
+
+ V_min = (pnts->r_max[0] - pnts->r_min[0]) * (pnts->r_max[1] - pnts->r_min[1]) * (pnts->r_max[2] - pnts->r_min[2]); /* Volume of extent */
+
+ for (i=0; i < hull.n_faces-2; i += 3) { /* n = number of vertices (n/3 = number of faces)*/
+ /* Bearings of hull edges */
+ usx = bearing(*(hull.faces+1), *(hull.faces+4), *(hull.faces+2), *(hull.faces+5));
+ usy = bearing(*hull.faces, *(hull.faces+6), *(hull.faces+2), *(hull.faces+8));
+
+ if (usx == -9999 || usy == -9999) /* Identical points */
+ continue;
+ cosusx = cos(usx);
+ sinusx = sin(usx);
+ cosusy = cos(usy);
+ sinusy = sin(usy);
+
+ hc_k = &hull.coord[0]; // original coords
+ for (k=0; k < hull.n; k++) {
+ /* Coordinate transformation */
+ *ht = *hc_k * cos(usy) + 0 - *(hc_k+2) * sin(usy);
+ *(ht+1) = *hc_k * sinusx * sinusy + *(hc_k+1) * cosusx + *(hc_k+2) * sinusx * cosusy;
+ *(ht+2) = *hc_k * cosusx * sinusy - *(hc_k+1) * sinusx + *(hc_k+2) * cosusx * cosusy;
+
+ /* Transformed extent */
+ switch (k) {
+ case 0: r_min = r_max = triple(*ht, *(ht+1), *(ht+2)); break;
+ default:
+ r_min = triple(MIN(*ht, *r_min), MIN(*(ht+1), *(r_min+1)), MIN(*(ht+2), *(r_min+2)));
+ r_max = triple(MAX(*ht, *r_max), MAX(*(ht+1), *(r_max+1)), MAX(*(ht+2), *(r_max+2)));
+ }
+ hc_k += 3;
+ } // end k
+
+ hf += 9;
+
+ V = (*r_max - *r_min) * (*(r_max+1) - *(r_min+1)) * (*(r_max+2) - *(r_min+2)); /* Area of transformed extent */
+ V_min = MIN(V, V_min);
+ } // end i
+
+ return V_min;
+}
+
+
Added: grass-addons/grass7/vector/v.nnstat/mbr.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/mbr.c (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/mbr.c 2015-06-22 08:03:48 UTC (rev 65511)
@@ -0,0 +1,205 @@
+#include "local_proto.h"
+#define PI M_PI
+
+/*******************************
+These functions are mostly taken from the module v.hull (Aime, A., Neteler, M., Ducke, B., Landa, M.)
+
+Minimum Bounding Rectangle (MBR)
+ - obtain vertices of convex hull,
+ - transform coordinates of vertices into coordinate system with axes parallel to hull's edges,
+ - find extents,
+ - compute areas and find minimum of them.
+ *******************************
+*/
+
+/* ----------------------------
+ * Obtain convex hull vertices: functions from hull.c (Andrea Aime)
+ * Theoretical background: Andrew's variant of Graham Scan
+ * Source: Mark Nelson, 2007: http://marknelson.us/2007/08/22/convex/
+ */
+
+/* Function to determine whether a point is above or below the line
+ * connecting the leftmost and the rightmost points */
+
+int rightTurn(double (*P)[3], int i, int j, int k)
+{
+ double a, b, c, d;
+
+ /* Coordinate translation: P[j] = 0.0 */
+ a = P[i][0] - P[j][0];
+ b = P[i][1] - P[j][1];
+ c = P[k][0] - P[j][0];
+ d = P[k][1] - P[j][1];
+ /* Determinant of P[i] and P[k] */
+ return a * d - b * c < 0; /* Returns |det| < 0 => P[k] is angled off in left direction */
+}
+
+/* Function to compare two points (for sorting points in convexHull) */
+int cmpPoints(const void *v1, const void *v2)
+{
+ double *p1, *p2;
+
+ p1 = (double *)v1;
+ p2 = (double *)v2;
+ if (p1[0] > p2[0]) {
+ return 1;
+ }
+ else if (p1[0] < p2[0]) {
+ return -1;
+ }
+ else {
+ return 0;
+ }
+}
+
+/* Function to obtain vertices of convex hull */
+int convexHull(struct points *pnts, struct convex *hull)
+{
+ int pointIdx, upPoints, loPoints;
+ int i, *upHull, *loHull;
+ int n = pnts->n;
+
+ double *ro, (*ri)[3], (*r)[3]; /* r = [xyz] */
+ r = (double (*)[3]) malloc(n * 3 * sizeof(double));
+ ri = &r[0];
+ ro = &pnts->r[0];
+
+ for (i=0; i<n; i++) {
+ (*ri)[0] = *ro;
+ (*ri)[1] = *(ro+1);
+ (*ri)[2] = *(ro+2);
+ ri++;
+ ro += 3;
+ }
+
+ /* sort points in ascending x order
+ * modified according to http://www.physicsforums.com/showthread.php?t=546209 */
+ qsort(&r[0][0], n, 3 * sizeof(double), cmpPoints);
+
+ hull->hull = (int *) G_malloc(n * 3 * sizeof(int));
+
+ /* compute upper hull */
+ upHull = hull->hull;
+ upHull[0] = 0;
+ upHull[1] = 1;
+ upPoints = 1; /* number of points in upper hull */
+ for (pointIdx = 2; pointIdx < n; pointIdx++) {
+ upPoints++;
+ upHull[upPoints] = pointIdx;
+ while (upPoints > 1 &&
+ !rightTurn(r, upHull[upPoints], upHull[upPoints - 1],
+ upHull[upPoints - 2])
+ ) {
+ upHull[upPoints - 1] = upHull[upPoints];
+ upPoints--;
+ }
+ }
+
+ /* compute lower hull, overwrite last point of upper hull */
+ loHull = &(upHull[upPoints]);
+ loHull[0] = n - 1;
+ loHull[1] = n - 2;
+ loPoints = 1; /* number of points in lower hull */
+ for (pointIdx = n - 3; pointIdx >= 0; pointIdx--) {
+ loPoints++;
+ loHull[loPoints] = pointIdx;
+ while (loPoints > 1 &&
+ !rightTurn(r, loHull[loPoints], loHull[loPoints - 1],
+ loHull[loPoints - 2])
+ ) {
+ loHull[loPoints - 1] = loHull[loPoints];
+ loPoints--;
+ }
+ }
+ hull->n = loPoints + upPoints;
+
+ G_debug(3, "numPoints:%d loPoints:%d upPoints:%d",
+ n, loPoints, upPoints);
+
+ /* reclaim uneeded memory */
+ hull->hull = (int *) G_realloc(hull->hull, (hull->n + 1) * sizeof(int));
+
+ /* Obtain coordinates of hull vertices */
+ hull->coord = (double *) G_malloc((hull->n + 1) * 3 * sizeof(double)); /* 1st = last pnt */
+
+ int *hh, *hh0;
+ double *hc, (*r0)[3];
+ hc = &hull->coord[0];
+ hh = &hull->hull[0];
+
+ for (i=0; i<=hull->n; i++) {
+ if (i < hull->n) {
+ *hc = (*(r+*hh))[0];
+ *(hc+1) = (*(r+*hh))[1];
+ *(hc+2) = (*(r+*hh))[2];
+ r++;
+ hc += 3;
+ hh++;
+ }
+ else { // coords of 1st equal to coords of last hull vertex
+ r0 = &r[0];
+ hh0 = &hull->hull[0];
+ *hc = (*(r0+*hh0))[0];
+ *(hc+1) = (*(r0+*hh0))[1];
+ *(hc+2) = (*(r0+*hh0))[2];
+ }
+ }
+ return hull->n;
+}
+
+/* ----------------------------
+ * MBR area estimation
+ * ----------------------------*/
+double MBR(struct points *pnts)
+{
+ int i, k;
+ double us, cosus, sinus, S, S_min;
+ double *hc, *hc_k; // original and to be transformed vertices of hull
+ double *r_min, *r_max;
+ r_min = (double *) G_malloc(3 * sizeof(double));
+ r_max = (double *) G_malloc(3 * sizeof(double));
+
+ double *hull_trans, *ht; /* Coordinates of hull vertices transformed into coordinate system with axes parallel to hull's edges */
+ hull_trans = (double *) G_malloc(3 * sizeof(double));
+ ht = &hull_trans[0];
+
+ struct convex hull;
+ hull.n = convexHull(pnts, &hull);
+ hc = &hull.coord[0];
+
+ S_min = (pnts->r_max[0] - pnts->r_min[0]) * (pnts->r_max[1] - pnts->r_min[1]); /* Area of extent */
+
+ for (i=0; i < hull.n; i++) {
+ /* Bearings of hull edges */
+ us = bearing(*hc, *(hc+3), *(hc+1), *(hc+4)); // x0, x1, y0, y1
+ if (us == -9999) // Identical points
+ continue;
+ cosus = cos(us);
+ sinus = sin(us);
+
+ hc_k = &hull.coord[0]; // original coords
+ for (k=0; k <= hull.n; k++) {
+ /* Coordinate transformation */
+ *ht = *hc_k * cosus + *(hc_k+1) * sinus;
+ *(ht+1) = -(*hc_k) * sinus + *(hc_k+1) * cosus;
+ *(ht+2) = *(hc_k+2);
+
+ /* Transformed extent */
+ switch (k) {
+ case 0:
+ r_min = r_max = triple(*ht, *(ht+1), *(ht+2)); break;
+ default:
+ r_min = triple(MIN(*ht, *r_min), MIN(*(ht+1), *(r_min+1)), MIN(*(ht+2), *(r_min+2)));
+ r_max = triple(MAX(*ht, *r_max), MAX(*(ht+1), *(r_max+1)), MAX(*(ht+2), *(r_max+2)));
+ }
+ hc_k += 3;
+ } // end k
+
+ hc += 3; // next point
+
+ S = (*r_max - *r_min) * (*(r_max+1) - *(r_min+1)); /* Area of transformed extent */
+ S_min = MIN(S, S_min);
+ } // end i
+
+ return S_min;
+}
Added: grass-addons/grass7/vector/v.nnstat/nearest_neighbour.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/nearest_neighbour.c (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/nearest_neighbour.c 2015-06-22 08:03:48 UTC (rev 65511)
@@ -0,0 +1,129 @@
+#include "local_proto.h"
+
+void nn_average_distance_real(struct nna_par *xD, struct points *pnts, struct nearest *nna)
+{
+ int i;
+ int i3 = xD->i3; // 2D or 3D NNA
+ int n = pnts->n; // # of points
+ double sum_r; // sum of the distances of the NNs
+
+ struct ilist *list;
+
+ G_message(_("Computing average distance between nearest neighbors..."));
+ sum_r = 0.;
+
+ for (i=0; i < n; i++) {
+ list = G_new_ilist(); // create list of overlapping rectangles
+ list = find_NNs(i3, i, pnts); // search NNs
+ sum_r += sum_NN(i3, i, list, pnts); // add NN distance to the sum
+ G_percent(i, n-1, 1); // progress bar
+ }
+ nna->rA = sum_r / n; // average NN distance
+
+ return;
+}
+
+void density(struct points *pnts, struct nna_par *xD, const char *Acol, struct nearest *nna)
+{
+ if (xD->i3 == TRUE) { // 3D NNA: Minimum Enclosing Block
+ nna->A = Acol != NULL ? atof(Acol) : MBB(pnts); // set volume specified by user or MBB volume
+ if (nna->A <= 0.) {
+ G_fatal_error(_("Volume must be greater than 0"));
+ }
+ }
+
+ if (xD->i3 == FALSE) { // 2D NNA: Minimum Enclosing Rectangle
+ nna->A = Acol != NULL ? atof(Acol) : MBR(pnts); // set area specified by user or MBR area
+ if (nna->A <= 0.) {
+ G_fatal_error(_("Area must be greater than 0"));
+ }
+ }
+
+ nna->rho = pnts->n / nna->A; // number of points per area/volume unit
+
+ return;
+}
+
+void nn_average_distance_expected(struct nna_par *xD, struct nearest *nna)
+{
+ if (xD->i3 == TRUE) { // 3D NNA:
+ nna->rE = 0.55396/pow(nna->rho,1./3.); // according to Chandrasekhar
+ }
+ if (xD->i3 == FALSE) { // 2D NNA:
+ nna->rE = 0.5/sqrt(nna->rho); // according to Clark & Evans
+ }
+}
+
+/* ---------------------------------
+ * Two-tailed Student's test of significance of the mean
+ * ---------------------------------*/
+
+void nn_results(struct points *pnts, struct nna_par *xD, struct nearest *nna)
+{
+ G_message(_("\n\n*** Nearest Neighbour Analysis results ***\n"));
+ if (xD->zcol != NULL) {
+ G_message(_("Input settings .. 3D layer: %d .. 3D NNA: %d .. zcolumn: %s"), xD->v3, xD->i3, xD->zcol);
+ }
+ else {
+ G_message(_("Input settings .. 3D layer: %d 3D NNA: %d"), xD->v3, xD->i3);
+ }
+ G_message(_("Number of points .......... %d"), pnts->n);
+ if (xD->i3 == TRUE) {
+ G_message(_("Volume .................... %f [units^3]"), nna->A);
+ }
+ else {
+ G_message(_("Area ...................... %f [units^2]"), nna->A);
+ }
+ G_message(_("Density of points ......... %f"), nna->rho);
+ G_message(_("Average distance between the nearest neighbours ........... %.3f [units]"), nna->rA);
+ G_message(_("Average expected distance between the nearest neighbours .. %.3f [units]"), nna->rE);
+ G_message(_("Ratio rA/rE ............... %f"), nna->R);
+ G_message(_("\n*** Results of two-tailed test of the mean ***"));
+ G_message(_("Null hypothesis: Point set is randomly distributed within the region."));
+ G_message(_("Standard variate of the normal curve> c = %f"), nna->c);
+
+ /* two-tailed test of the mean */
+ if (-1.96 < nna->c && nna->c < 1.96) {
+ G_message(_("Null hypothesis IS NOT REJECTED at the significance level alpha = 0.05"));
+ }
+ if (-2.58 < nna->c && nna->c <= -1.96) {
+ G_message(_("Null hypothesis IS NOT REJECTED at the significance level alpha = 0.01 => point set is clustered"));
+ G_message(_("Null hypothesis CAN BE REJECTED at the significance level alpha = 0.05 => point set is randomly distributed"));
+ }
+
+ if (1.96 <= nna->c && nna->c < 2.58) {
+ G_message(_("Null hypothesis IS NOT REJECTED at the significance level alpha = 0.01 => point set is dispersed"));
+ G_message(_("Null hypothesis CAN BE REJECTED at the significance level alpha = 0.05 => point set is randomly distributed"));
+ }
+
+ if (nna->c <= -2.58) {
+ G_message(_("Null hypothesis CAN BE REJECTED at the significance levels alpha = 0.05 and alpha = 0.01 => point set is clustered"));
+ }
+ if (nna->c >= 2.58) {
+ G_message(_("Null hypothesis CAN BE REJECTED at the significance levels alpha = 0.05 and alpha = 0.01 => point set is dispersed"));
+ }
+ G_message(" ");
+
+ return;
+}
+
+void nn_statistics(struct points *pnts, struct nna_par *xD, struct nearest *nna)
+{
+ double disp_rE, var_rE, sig_rE;
+
+ if (xD->i3 == TRUE) { // 3D NNA:
+ disp_rE = 0.347407742479038/pow(nna->rho,2./3.); // standard deviation of average distance between the NN in randomly distributed set of points; derived in (Stopkova, 2013)
+ var_rE = 0.0405357524727094/pow(nna->rho,2./3.); // variance of average distance between the NN in randomly distributed set of points; derived in (Stopkova, 2013)
+ }
+ if (xD->i3 == FALSE) { // 2D NNA:
+ disp_rE = 1./(nna->rho*PI); // standard deviation of average distance between the NN in randomly distributed set of points (Clark & Evans)
+ var_rE = (4.0 - PI)/(4.0 * PI * pnts->n * nna->rho); // variance of average distance between the NN in randomly distributed set of points (Clark & Evans)
+ }
+
+ sig_rE = sqrt(var_rE); // standard error of the mean distance
+ nna->c = (nna->rA - nna->rE) / sig_rE; // standard variate of the normal curve
+
+ nn_results(pnts, xD, nna); // write the result to the command line
+
+ return;
+}
Added: grass-addons/grass7/vector/v.nnstat/spatial_index.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/spatial_index.c (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/spatial_index.c 2015-06-22 08:03:48 UTC (rev 65511)
@@ -0,0 +1,146 @@
+#include "local_proto.h"
+
+struct RTree *create_spatial_index(struct nna_par *xD)
+{
+ struct RTree *R_tree; // spatial index
+ struct RTree_Rect *rect; // rectangle
+
+ if (xD->i3 == TRUE) { // 3D NNA:
+ R_tree = RTreeCreateTree(-1, 0, 3); // create 3D spatial index
+ }
+ if (xD->i3 == FALSE) { // 2D NNA:
+ R_tree = RTreeCreateTree(-1, 0, 2); // create 3D spatial index
+ }
+
+ return(R_tree);
+}
+
+void insert_rectangle(int i3, int i, struct points *pnts)
+{
+ struct RTree *R_tree; // spatial index
+ struct RTree_Rect *rect; // rectangle
+ double *r; // pointer to the input coordinates
+
+ r = &pnts->r[3 * i];
+ rect = RTreeAllocRect(pnts->R_tree); // allocate the rectangle
+
+ if (i3 == TRUE) { // 3D NNA:
+ RTreeSetRect3D(rect, pnts->R_tree, *r, *r, *(r+1), *(r+1), *(r+2), *(r+2)); // set 3D coordinates
+ }
+
+ if (i3 == FALSE) { // 2D NNA:
+ RTreeSetRect2D(rect, pnts->R_tree, *r, *r, *(r+1), *(r+1)); // set 2D coordinates
+ }
+
+ RTreeInsertRect(rect, i + 1, pnts->R_tree); // insert rectangle to the spatial index
+ RTreeFreeRect(rect); // free rectangle memory
+
+ return;
+}
+
+/* find neighbors in cubic (square) surroundings */
+struct ilist *spatial_search(int i3, int i, struct points *pnts, double max_dist)
+{
+ // local variables
+ struct RTree *R_tree = pnts->R_tree; // spatial index
+ double *r; // pointer to vector of the coordinates
+ r = &pnts->r[3 * i];
+
+ struct ilist *list;
+ struct RTree_Rect *search; // search rectangle
+
+ list = G_new_ilist(); // create new list
+ search = RTreeAllocRect(R_tree); // allocate new rectangle
+
+ if (i3 == TRUE) { // 3D NNA:
+ RTreeSetRect3D(search, R_tree,
+ *r - max_dist, *r + max_dist,
+ *(r+1) - max_dist, *(r+1) + max_dist,
+ *(r+2) - max_dist, *(r+2) + max_dist); // set up searching rectangle
+ }
+
+ if (i3 == FALSE) { // 2D NNA
+ RTreeSetRect2D(search, R_tree,
+ *r - max_dist, *r + max_dist,
+ *(r+1) - max_dist, *(r+1) + max_dist); // set up searching rectangle
+ }
+
+ RTreeSearch2(R_tree, search, list); // search the nearest rectangle
+ RTreeFreeRect(search);
+
+ return list;
+}
+
+/* find neighbors in spherical (circular) surroundings */
+struct ilist *find_NNs(int i3, int i, struct points *pnts)
+{
+ // local variables
+ double max_dist0 = 0.005 * pnts->max_dist; // initial search radius
+ struct ilist *list; // list of selected rectangles (points)
+
+ int n_vals = 0; // # of the nearest neighbours (NN)
+ int iter = 1; // iteration number (searching NN)
+ double max_dist = max_dist0; // iterated search radius
+
+ while (n_vals < 2) { // until some NN is found (2 points: identical and NN)
+ list = spatial_search(i3, i, pnts, max_dist);
+ n_vals = list->n_values; // set up control variable
+
+ if (n_vals < 2) {
+ iter += 1; // go to the next iteration
+ max_dist = iter * max_dist0; // enlarge search radius
+ }
+ }
+
+ // check spherical (circular) surrounding
+ list = spatial_search(i3, i, pnts, sqrt(2.) * max_dist);
+
+ return list;
+}
+
+/* sum of the distances between nearest neighbors */
+double sum_NN(int i3, int i, struct ilist *list, struct points *pnts)
+{
+ int j;
+ int n_vals = list->n_values; // # of selected points
+ int *value; // pointer to indices of selected points
+ value = &list->value[0];
+ *value -= 1; // decrease index of selected point by 1 (because rectangle IDs start by 1, not by 0 like point IDs)
+
+ double *r0 = &pnts->r[3 * i]; // pointer to search point coordinates
+ double *r = &pnts->r[3 * *value]; // pointer to the coordinates of selected point
+ double dx, dy, dz; // coordinate differences between search and selected point
+ double dist, *sqDist, *d; // vector of squared distances
+
+ sqDist = (double *) G_malloc(n_vals * sizeof(double));
+ d = &sqDist[0]; // pointer to vector of squared distances
+
+ for (j=0; j<n_vals; j++) { // for each rectangle:
+ dx = *r0 - *r;
+ dy = *(r0+1) - *(r+1);
+
+ if (i3 == TRUE) { // 3D NNA:
+ dz = *(r0+2) - *(r+2); // compute also 3rd difference
+ }
+
+ *d = dx*dx + dy*dy; // squared distance
+
+ if (i3 == TRUE) { // 3D NNA:
+ *d += dz*dz; // use also the difference in z
+ }
+
+ d++; // go to the next distance
+
+ value++; // go to the index of the next selected point
+ *value -= 1; // decrease the index by 1
+ r += 3 * (*value - *(value-1)); // go to the next selected point
+ }
+
+ qsort(&sqDist[0], n_vals, sizeof(double), cmpVals); // sort squared distances
+ dist = sqrt(sqDist[1]); // save the second one (1st NN is identical to the search point)
+
+ G_free(sqDist); // free memory of the vector of squared distances
+ G_free_ilist(list); // free list of selected points
+
+ return dist;
+}
Added: grass-addons/grass7/vector/v.nnstat/utils.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/utils.c (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/utils.c 2015-06-22 08:03:48 UTC (rev 65511)
@@ -0,0 +1,65 @@
+#include "local_proto.h"
+
+/* make triple of the coordinates */
+double *triple(double x, double y, double z)
+{
+ double *t;
+ t = (double *) G_malloc(3 * sizeof(double));
+
+ *t = x;
+ *(t+1) = y;
+ *(t+2) = z;
+
+ return t;
+}
+
+/* bearing */
+double bearing(double x0, double x1, double y0, double y1)
+{
+ double dy, dx, us;
+ dy = y1 - y0;
+ dx = x1 - x0;
+ if (dy == 0. && dx == 0.) {
+ return -9999;
+ }
+
+ us = atan(fabsf(dy/dx));
+ if (dx==0. && dy>0.) {
+ us = 0.5*PI;
+ }
+ if (dx<0. && dy>0.) {
+ us = PI - us;
+ }
+ if (dx<0. && dy==0.) {
+ us = PI;
+ }
+ if (dx<0. && dy<0.) {
+ us = PI + us;
+ }
+ if (dx==0. && dy<0.) {
+ us = 1.5*PI;
+ }
+ if (dx>0. && dy<0.) {
+ us = 2.0*PI - us;
+ }
+
+ return us;
+}
+
+/* According to two points comparison (sorting points in convexHull) */
+int cmpVals(const void *v1, const void *v2)
+{
+ double *p1, *p2;
+
+ p1 = (double *)v1;
+ p2 = (double *)v2;
+ if (p1[0] > p2[0]) {
+ return 1;
+ }
+ else if (p1[0] < p2[0]) {
+ return -1;
+ }
+ else {
+ return 0;
+ }
+}
More information about the grass-commit
mailing list