[GRASS-SVN] r62314 - in grass-addons/grass7/vector: . v.nnstat
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Oct 21 06:01:49 PDT 2014
Author: evas
Date: 2014-10-21 06:01:49 -0700 (Tue, 21 Oct 2014)
New Revision: 62314
Added:
grass-addons/grass7/vector/v.nnstat/
grass-addons/grass7/vector/v.nnstat/Makefile
grass-addons/grass7/vector/v.nnstat/getval.cpp
grass-addons/grass7/vector/v.nnstat/hull.h
grass-addons/grass7/vector/v.nnstat/local_proto.h
grass-addons/grass7/vector/v.nnstat/macros.h
grass-addons/grass7/vector/v.nnstat/main.cpp
grass-addons/grass7/vector/v.nnstat/mbb.cpp
grass-addons/grass7/vector/v.nnstat/mbb.h
grass-addons/grass7/vector/v.nnstat/mbr.cpp
grass-addons/grass7/vector/v.nnstat/nearest_neighbour.cpp
grass-addons/grass7/vector/v.nnstat/utils.cpp
grass-addons/grass7/vector/v.nnstat/v.nnstat.html
Log:
new module v.nnstat
Added: grass-addons/grass7/vector/v.nnstat/Makefile
===================================================================
--- grass-addons/grass7/vector/v.nnstat/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/Makefile 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,16 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.nnstat
+
+LIBES = $(VECTORLIB) $(DBMILIB) $(GISLIB) $(GMATHLIB) $(IOSTREAMLIB) $(MATHLIB) -lpcl_apps -lpcl_common -lpcl_features -lpcl_filters -lpcl_io -lpcl_kdtree -lpcl_keypoints -lpcl_octree -lpcl_registration -lpcl_sample_consensus -lpcl_search -lpcl_segmentation -lpcl_surface -lpcl_visualization
+DEPENDENCIES= $(VECTORDEP) $(DBMIDEP) $(GISDEP) $(IOSTREAMDEP)
+EXTRA_INC = $(VECT_INC) -I/usr/include/pcl-1.7 -I/usr/include/eigen3 $(PCL_INCLUDE_DIRS)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+LINK = $(CXX)
+
+ifneq ($(strip $(CXX)),)
+default: cmd
+endif
Added: grass-addons/grass7/vector/v.nnstat/getval.cpp
===================================================================
--- grass-addons/grass7/vector/v.nnstat/getval.cpp (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/getval.cpp 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,134 @@
+#include "local_proto.h"
+
+extern "C" {
+ /* get array of values from attribute column (function based on part of v.buffer2 (Radim Blazek, Rosen Matev)) */
+ void *get_col_values(struct Map_info *map, int field, const char *column, union dat *ifs)
+ {
+ int i, nrec, ctype;
+ struct field_info *Fi;
+
+ dbCatValArray cvarr;
+ dbDriver *Driver;
+
+ int *iv;
+ double *fv;
+ string *sv;
+
+ 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(_("%d records selected from table"), nrec);
+ ctype = cvarr.ctype;
+
+ db_close_database_shutdown_driver(Driver);
+
+ switch (ctype) {
+ case DB_C_TYPE_INT: ifs->i_dat = (int *) G_malloc(cvarr.n_values * sizeof(int)); iv = &ifs->i_dat[0]; break;
+ case DB_C_TYPE_DOUBLE: ifs->f_dat = (double *) G_malloc(cvarr.n_values * sizeof(double)); fv = &ifs->f_dat[0]; break;
+ case DB_C_TYPE_STRING: ifs->s_dat = (char **) G_malloc(cvarr.n_values * sizeof(char)); sv = &ifs->s_dat[0];
+ }
+
+ for (i = 0; i < cvarr.n_values; i++) {
+ switch (ctype) {
+ case DB_C_TYPE_INT: *iv = cvarr.value[i].val.i; iv++; break;
+ case DB_C_TYPE_DOUBLE: *fv = cvarr.value[i].val.d; fv++; break;
+ case DB_C_TYPE_STRING:
+ G_message(_("cvarr = %s"), db_get_string(cvarr.value[i].val.s));
+ *sv = db_get_string(cvarr.value[i].val.s);
+ G_message(_("ifs = %s"), *sv);
+ sv++;
+ }
+ }
+ }
+
+/* get coordinates of input points */
+ void read_points(struct Map_info *map, int field, struct nna_par *xD,
+ const char *zcol, union dat *ifs, struct points *point)
+ {
+ int ctrl, n, type, nskipped, pass;
+
+ struct line_pnts *Points; /* structures to hold line *Points (map) */
+ double *rx, *ry, *rz, *z_attr;
+
+ Points = Vect_new_line_struct();
+ n = point->n = Vect_get_num_primitives(map, GV_POINT); /* topology required */
+ point->r = (double *) G_malloc(n*3 * sizeof(double));
+
+ rx = &point->r[0];
+ ry = &point->r[1];
+ rz = &point->r[2];
+
+ /* Get 3rd coordinate of 2D points from attribute column -> 3D interpolation */
+ /*if (xD.v3 == FALSE && zcol != NULL)
+ z_attr = get_col_values(map, field, zcol);*/
+
+ 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 (zcol == NULL)
+ *rz = Points->z[0];
+ /* TODO: 2D */
+ else {
+ *rz = *z_attr;
+ z_attr++;
+ }
+
+ /* Find extends */
+ if (ctrl == 0) {
+ point->r_min = point->r_max = triple(*rx, *ry, *rz);
+ }
+ else {
+ point->r_min = triple(MIN(*point->r_min, *rx), MIN(*(point->r_min+1), *ry), MIN(*(point->r_min+2), *rz));
+ point->r_max = triple(MAX(*point->r_max, *rx), MAX(*(point->r_max+1), *ry), MAX(*(point->r_max+2), *rz));
+ }
+ if (ctrl < point->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."));
+ }
+} // end extern "C"
Added: grass-addons/grass7/vector/v.nnstat/hull.h
===================================================================
--- grass-addons/grass7/vector/v.nnstat/hull.h (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/hull.h 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,8 @@
+extern "C" {
+#include <grass/vector.h>
+}
+
+#define ALLOC_CHUNK 256
+
+struct cat_list *parse_filter_options(struct Map_info *, int , char *, char *);
+
Added: grass-addons/grass7/vector/v.nnstat/local_proto.h
===================================================================
--- grass-addons/grass7/vector/v.nnstat/local_proto.h (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/local_proto.h 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,100 @@
+#ifndef __LOCAL_PROTO__
+#define __LOCAL_PROTO__
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include <pcl/point_cloud.h>
+#include <pcl/kdtree/impl/kdtree_flann.hpp>
+
+extern "C" {
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/config.h>
+#include <grass/gis.h>
+#include <grass/gmath.h>
+#include <grass/la.h>
+#include <grass/glocale.h>
+} // end extern "C"
+
+#include <iostream>
+#include <vector>
+
+#ifndef MAX
+#define MIN(a,b) ((a<b) ? a : b)
+#define MAX(a,b) ((a>b) ? a : b)
+#endif
+
+#ifndef PI
+#define PI M_PI
+#endif
+
+struct points
+{
+ int n;
+ double *r; // xyz
+ double *r_min; // min coords
+ double *r_max; // max coords
+};
+
+struct nna_par
+{
+ int i3; // TRUE = 3D NNA, FALSE = 2D NNA (user sets by flag)
+ char v3; // TRUE = 3D layer, FALSE = 2D layer (according to Vect_is_3d())
+};
+
+typedef char *string;
+union dat
+{
+ int *i_dat;
+ double *f_dat;
+ string *s_dat;
+};
+
+struct columns
+{
+ int *lyr; // archaeological layer
+ double *z; // z coordinate
+ string *desc; // description of artifact
+};
+
+struct convex
+{
+ int n; // number of hull vertices
+ int n_faces; // number of hull faces
+ int *hull;
+ double *coord; // coordinates of faces
+ double *faces; // coordinates of vertices
+};
+
+struct nearest
+{
+ double A;
+ double rho;
+ double rA;
+ double rE;
+ double R;
+ double c;
+};
+
+extern "C" {
+ void *get_col_values(struct Map_info *, int, const char *, union dat *);
+ void read_points(struct Map_info *, int, struct nna_par *, const char *, union dat *, struct points *);
+ double *triple(double, double, double);
+
+ double bearing(double, double, double, double);
+ double distance(double, double);
+ int convexHull(struct points *, struct convex *);
+ int make3DHull(struct points *, struct convex *);
+ void convexHull3d(struct points *, struct convex *);
+ double MBR(struct points *);
+ double MBB(struct points *);
+ int nn_average_distance_real(struct points *, struct nearest *);
+ void density(struct points *, struct nna_par *, const char *, struct nearest *);
+ void nn_average_distance_expected(struct nna_par *, struct nearest *);
+ void nn_results(struct points *, struct nearest *);
+ void nn_statistics(struct points *, struct nna_par *, struct nearest *);
+} // end extern "C"
+#endif
Added: grass-addons/grass7/vector/v.nnstat/macros.h
===================================================================
--- grass-addons/grass7/vector/v.nnstat/macros.h (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/macros.h 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,39 @@
+
+/*====================================================================
+ macros.h
+
+ macros used to access data structures and perform quick tests.
+
+ ====================================================================*/
+
+/* general-purpose macros */
+#define SWAP(t,x,y) { t = x; x = y; y = t; }
+
+#define NEW(p,type) if ((p=(type *) malloc (sizeof(type))) == NULL) {\
+ printf ("Out of Memory!\n");\
+ exit(0);\
+ }
+
+#define FREE(p) if (p) { free ((char *) p); p = NULL; }
+
+
+#define ADD( head, p ) if ( head ) { \
+ p->next = head; \
+ p->prev = head->prev; \
+ head->prev = p; \
+ p->prev->next = p; \
+ } \
+ else { \
+ head = p; \
+ head->next = head->prev = p; \
+ }
+
+#define DELETE( head, p ) if ( head ) { \
+ if ( head == head->next ) \
+ head = NULL; \
+ else if ( p == head ) \
+ head = head->next; \
+ p->next->prev = p->prev; \
+ p->prev->next = p->next; \
+ FREE( p ); \
+ }
Added: grass-addons/grass7/vector/v.nnstat/main.cpp
===================================================================
--- grass-addons/grass7/vector/v.nnstat/main.cpp (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/main.cpp 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,113 @@
+
+/****************************************************************
+ *
+ * MODULE: v.nnstat
+ *
+ * AUTHOR(S): Eva Stopková
+ * functionsin files mbr.cpp and mbb.cpp are mostly taken
+ * from the module
+ * v.hull (Aime, A., Neteler, M., Ducke, B., Landa, M.)
+ *
+ * PURPOSE: Module indicates clusters, separations or random distribution of point set in 2D or 3D space.
+ *
+ * COPYRIGHT: (C) 2013-2014 by Eva Stopková and the GRASS Development Team
+ *
+ * This program is free software under the
+ * GNU General Public License (>=v2).
+ * Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ **************************************************************/
+#include "local_proto.h"
+
+extern "C" {
+ int main(int argc, char *argv[])
+ {
+ /* Vector layer and module */
+ struct Map_info map; /* input vector map */
+ struct GModule *module;
+ struct nna_par xD; /* 2D or 3D NNA */
+ struct points pnts;
+ union dat ifs; /* integer/float/string value */
+
+ struct nearest nna;
+ int field, pass;
+
+ struct {
+ struct Option *map, *A, *type, *field, *lyr, *desc, *zcol; /* A - area (minimum enclosing rectangle or specified by user) */
+ } opt;
+
+ struct {
+ struct Flag *d23;
+ } flg;
+
+ /* Module creation */
+ module = G_define_module();
+ G_add_keyword(_("Vector"));
+ G_add_keyword(_("Nearest Neighbour Analysis"));
+ module->description = _("Indicates clusters, separations or random distribution of point set in 2D or 3D space.");
+
+ /* Setting options */
+ opt.map = G_define_standard_option(G_OPT_V_INPUT); /* vector input layer */
+ opt.map->label = _("Name of input vector map");
+
+ flg.d23 = G_define_flag(); /* to process 2D or 3D Nearest Neighbour Analysis */
+ flg.d23->key = '2';
+ flg.d23->description =
+ _("Force 2D NNA even if input is 3D");
+
+ opt.A = G_define_option();
+ opt.A->key = "A";
+ opt.A->type = TYPE_DOUBLE;
+ opt.A->description = _("2D: Area. If not specified, area of Minimum Enclosing Rectangle will be used.\n3D: Volume. If not specified, volume of Minimum Enclosing Box will be used.");
+ opt.A->required = NO;
+
+ opt.field = G_define_standard_option(G_OPT_V_FIELD);
+
+ opt.zcol = G_define_standard_option(G_OPT_DB_COLUMN);
+ opt.zcol->key = "zcol";
+ opt.zcol->required = NO;
+ opt.zcol->guisection = _("Fields");
+ opt.zcol->description = _("Column with z coordinate (set for 2D vectors only)");
+
+ G_gisinit(argv[0]);
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* get parameters from the parser */
+ field = opt.field->answer ? atoi(opt.field->answer) : -1;
+
+ /* open input vector map */
+ Vect_set_open_level(2);
+
+ if (0 > Vect_open_old2(&map, opt.map->answer, "", opt.field->answer))
+ G_fatal_error(_("Unable to open vector map <%s>"), opt.map->answer);
+ Vect_set_error_handler_io(&map, NULL);
+
+ /* perform 2D or 3D NNA ? */
+ xD.v3 = Vect_is_3d(&map);
+ xD.i3 = (flg.d23->answer || xD.v3 == FALSE) ? FALSE : TRUE;
+
+ if (xD.i3 == TRUE) { /* 3D */
+ if (xD.v3 == FALSE) { /* 2D input */
+ if (opt.zcol->answer == NULL)
+ G_fatal_error(_("To process 3D Nearest Neighbour Analysis based on 2D input, "
+ "please set attribute column containing z coordinates "
+ "or switch to 2D NNA."));
+ }
+ }
+
+ read_points(&map, field, &xD, opt.zcol->answer, &ifs, &pnts);
+
+ /* Nearest Neighbour Analysis */
+ pass = nn_average_distance_real(&pnts, &nna); // Average distance (AD) between NN in pointset
+ density(&pnts, &xD, opt.A->answer, &nna); // Number of points per area/volume unit
+ nn_average_distance_expected(&xD, &nna); // Expected AD of NN in a randomly distributed pointset
+ nna.R = nna.rA / nna.rE; // Ratio of real and expected AD
+ nn_statistics(&pnts, &xD, &nna); // Student's t-test of significance of the mean
+
+ Vect_close(&map);
+ exit(EXIT_SUCCESS);
+ }
+}
Added: grass-addons/grass7/vector/v.nnstat/mbb.cpp
===================================================================
--- grass-addons/grass7/vector/v.nnstat/mbb.cpp (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/mbb.cpp 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,750 @@
+#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.
+ *******************************
+*/
+
+extern "C" {
+ /* 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);
+ }
+
+
+ /*-------------------------------------------------------------------*/
+ 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();
+ }
+
+ /*---------------------------------------------------------------------
+ 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);
+ }
+
+ /*---------------------------------------------------------------------
+ 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);
+ }
+
+ /*---------------------------------------------------------------------
+ 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);
+ }
+
+ /*---------------------------------------------------------------------
+ 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"));
+ }
+ }
+} // end of extern "C"
+
+ /* ----------------------------
+ * 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/mbb.h
===================================================================
--- grass-addons/grass7/vector/v.nnstat/mbb.h (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/mbb.h 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,84 @@
+#include "local_proto.h"
+
+/*
+ These functions are mostly taken from the module v.hull (Aime, A., Neteler, M., Ducke, B., Landa, M.)
+*/
+
+/* Define Boolean type - not necessary in C++*/
+/*typedef enum
+ { BFALSE, BTRUE } bool;
+
+/* Define vertex indices. */
+#define X 0
+#define Y 1
+#define Z 2
+
+/* Define structures for vertices, edges and faces */
+typedef struct tVertexStructure tsVertex;
+typedef tsVertex *tVertex;
+
+typedef struct tEdgeStructure tsEdge;
+typedef tsEdge *tEdge;
+
+typedef struct tFaceStructure tsFace;
+typedef tsFace *tFace;
+
+struct tVertexStructure
+{
+ double v[3];
+ int vnum;
+ tEdge duplicate; /* pointer to incident cone edge (or NULL) */
+ bool onhull; /* T iff point on hull. */
+ bool mark; /* T iff point already processed. */
+ tVertex next, prev;
+};
+
+struct tEdgeStructure
+{
+ tFace adjface[2];
+ tVertex endpts[2];
+ tFace newface; /* pointer to incident cone face. */
+ bool del; /* T iff edge should be delete. */
+ tEdge next, prev;
+};
+
+struct tFaceStructure
+{
+ tEdge edge[3];
+ tVertex vertex[3];
+ bool visible; /* T iff face visible from new point. */
+ tFace next, prev;
+};
+
+/* Define flags */
+#define ONHULL TRUE
+#define REMOVED TRUE
+#define VISIBLE TRUE
+#define PROCESSED TRUE
+
+/* Global variable definitions */
+tVertex vertices = NULL;
+tEdge edges = NULL;
+tFace faces = NULL;
+
+/* Function declarations */
+extern "C" {
+tVertex MakeNullVertex(void);
+void ReadVertices(struct points *);
+void writeVertices(struct Map_info *Map);
+void write_coord_faces(struct points *, struct convex *);
+int DoubleTriangle(void);
+int ConstructHull(void);
+bool AddOne(tVertex p);
+int VolumeSign(tFace f, tVertex p);
+tFace MakeConeFace(tEdge e, tVertex p);
+void MakeCcw(tFace f, tEdge e, tVertex p);
+tEdge MakeNullEdge(void);
+tFace MakeNullFace(void);
+tFace MakeFace(tVertex v0, tVertex v1, tVertex v2, tFace f);
+void CleanUp(void);
+void CleanEdges(void);
+void CleanFaces(void);
+void CleanVertices(void);
+bool Collinear(tVertex a, tVertex b, tVertex c);
+}
Added: grass-addons/grass7/vector/v.nnstat/mbr.cpp
===================================================================
--- grass-addons/grass7/vector/v.nnstat/mbr.cpp (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/mbr.cpp 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,203 @@
+#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 */
+extern "C" {
+ 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;
+ }
+} // end extern "C"
+
+/* 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.cpp
===================================================================
--- grass-addons/grass7/vector/v.nnstat/nearest_neighbour.cpp (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/nearest_neighbour.cpp 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,130 @@
+#include "local_proto.h"
+using namespace std;
+using namespace pcl;
+
+int nn_average_distance_real(struct points *pnts, struct nearest *nna)
+{
+ int i, n = pnts->n, k=2, ind0 = 0;
+ double sum_r;
+
+ double *r;
+ pcl::PointXYZ *pcl_r;
+
+ r = (double *) G_malloc(n*3 * sizeof(double));
+ memcpy(r, pnts->r, n*3 * sizeof(double));
+
+ /* create new vector of 3D points (use template PointCloud from namespace pcl) */
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts (new pcl::PointCloud<pcl::PointXYZ>);
+ pcl_pnts->width = n;
+ pcl_pnts->height = 1;
+ pcl_pnts->points.resize (pcl_pnts->width * pcl_pnts->height);
+
+ pcl_r = &pcl_pnts->points[0];
+
+ for (i=0; i<n; i++) {
+ pcl_r->x = *r;
+ pcl_r->y = *(r+1);
+ pcl_r->z = *(r+2);
+ pcl_r++;
+ r += 3;
+ }
+
+ pcl::KdTreeFLANN<pcl::PointXYZ> kd_tree; /* create kd-tree */
+ kd_tree.setInputCloud(pcl_pnts);
+
+ std::vector<int> ind(k);
+ std::vector<float> sqDist(k);
+
+ double dx, dy, dz, dist;
+ double *ri, *r0, *rk;
+
+ /* Average distance to NN */
+ sum_r = 0.;
+ ri = &pnts->r[0];
+ r0 = &pnts->r[0];
+ for (size_t i=0; i < pcl_pnts->points.size (); i++) {
+ if (kd_tree.nearestKSearch(pcl_pnts->points[i], k, ind, sqDist) > 0) {
+ /* sqDist not used because double precision is needed (in PCL not supported yet) */
+ rk = r0 + 3*ind[1];
+ //dist = distance(ri, rk);
+ dx = *rk - *ri;
+ dy = *(rk+1) - *(ri+1);
+ dz = *(rk+2) - *(ri+2);
+ dist = sqrt(dx*dx + dy*dy + dz*dz);
+ sum_r += dist;
+ } // end nearestKSearch
+ ri+=3;
+ } // end i
+ nna->rA = sum_r / pcl_pnts->points.size ();
+ int pass=1;
+ return pass;
+}
+
+extern "C" {
+ void density(struct points *pnts, struct nna_par *xD, const char *Acol, struct nearest *nna)
+ {
+ if (xD->i3 == TRUE) { /* 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) { /* 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 */
+ }
+
+ void nn_average_distance_expected(struct nna_par *xD, struct nearest *nna)
+ {
+ if (xD->i3 == TRUE) /* 3D */
+ nna->rE = 0.55396/pow(nna->rho,1./3.); /* according to Chandrasekhar */
+ if (xD->i3 == FALSE) /* 2D */
+ 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 nearest *nna)
+ {
+ G_message(_("\nNearest Neighbour Analysis results\n"));
+ G_message(_("Number of points .......... %d\nArea/Volume .......... %f\nDensity of points ............ %f"), pnts->n, nna->A, nna->rho);
+ G_message(_("Average distance between the nearest neighbours ............ %f m\nAverage expected distance between the nearest neighbours ... %f m\n Ratio rA/rE ... %f\n\nResults of two-tailed test of the mean\nNull hypothesis: Point set is randomly distributed within the region.\nStandard variate of the normal curve> c = %f"), nna->rA, nna->rE, nna->R, 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.\nNull 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.\nNull 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. %f"), nna->c);
+ }
+
+ 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 */
+ disp_rE = 0.347407742479038/pow(nna->rho,2./3.); /* disperzia, odvodene */
+ var_rE = 0.0405357524727094/pow(nna->rho,2./3.); /* variancia, odvodene */
+ }
+ if (xD->i3 == FALSE) { /* 2D */
+ disp_rE = 1./(nna->rho*PI); /* disperzia */
+ var_rE = (4.0 - PI)/(4.0 * PI * pnts->n * nna->rho); /* variancia */
+ }
+
+ 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, nna);
+ }
+}
Added: grass-addons/grass7/vector/v.nnstat/utils.cpp
===================================================================
--- grass-addons/grass7/vector/v.nnstat/utils.cpp (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/utils.cpp 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1,50 @@
+#include "local_proto.h"
+
+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;
+}
+
+/*double distance(double *r0, double *r1)
+{
+ double dx, dy, dz, dist;
+ dx = *r1 - *r0;
+ dy = *(r1+1) - *(r0+1);
+ dz = *(r1+2) - *(r0+2);
+ dist = sqrt(dx*dx + dy*dy + dz*dz);
+
+ return dist;
+ }*/
Added: grass-addons/grass7/vector/v.nnstat/v.nnstat.html
===================================================================
--- grass-addons/grass7/vector/v.nnstat/v.nnstat.html (rev 0)
+++ grass-addons/grass7/vector/v.nnstat/v.nnstat.html 2014-10-21 13:01:49 UTC (rev 62314)
@@ -0,0 +1 @@
+TO DO
More information about the grass-commit
mailing list