[GRASS-SVN] r65550 - grass-addons/grass7/vector/v.nnstat
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jul 9 01:29:00 PDT 2015
Author: evas
Date: 2015-07-09 01:29:00 -0700 (Thu, 09 Jul 2015)
New Revision: 65550
Modified:
grass-addons/grass7/vector/v.nnstat/getval.c
grass-addons/grass7/vector/v.nnstat/main.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: indented
Modified: grass-addons/grass7/vector/v.nnstat/getval.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/getval.c 2015-07-08 18:54:02 UTC (rev 65549)
+++ grass-addons/grass7/vector/v.nnstat/getval.c 2015-07-09 08:29:00 UTC (rev 65550)
@@ -4,158 +4,167 @@
/* 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;
+ int i, nrec, ctype;
+ struct field_info *Fi;
- dbCatValArray cvarr;
- dbDriver *Driver;
+ dbCatValArray cvarr;
+ dbDriver *Driver;
- double *values, *vals;
+ double *values, *vals;
- db_CatValArray_init(&cvarr); /* array of categories and values initialised */
+ 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);
+ 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 */
+ /* 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 */
+ /* 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);
+ /* 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;
- values = (double *) G_malloc(cvarr.n_values * sizeof(double));
- vals = &values[0];
+ db_close_database_shutdown_driver(Driver);
- for (i = 0; i < cvarr.n_values; i++) {
- if (ctype == DB_C_TYPE_INT) {
- *vals = (double) cvarr.value[i].val.i;
+ 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++;
}
- 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;
+ 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)
+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;
+ 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
+ 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));
+ 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];
+ rx = &pnts->r[0];
+ ry = &pnts->r[1];
+ rz = &pnts->r[2];
- pnts->R_tree = create_spatial_index(xD); // create spatial index (R-tree)
+ 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"));
+ /* 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];
}
- if (type == -2) {
- break;
+ else {
+ xD->zcol = NULL;
}
- if (type != GV_POINT) {
- nskipped++;
- continue;
- }
+ nskipped = ctrl = pass = 0;
- *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];
- }
+ 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++;
}
- 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);
+ if (nskipped > 0) {
+ G_warning(_("%d features skipped, only points are accepted"),
+ nskipped);
}
- 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);
- Vect_destroy_line_struct(Points);
+ if (ctrl <= 0) {
+ G_fatal_error(_("Unable to read coordinates of point layer"));
+ }
- 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];
- // 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);
+ 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..."));
+ G_message(_("Input coordinates have been read..."));
- return;
+ return;
}
Modified: grass-addons/grass7/vector/v.nnstat/main.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/main.c 2015-07-08 18:54:02 UTC (rev 65549)
+++ grass-addons/grass7/vector/v.nnstat/main.c 2015-07-09 08:29:00 UTC (rev 65550)
@@ -22,90 +22,95 @@
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 to be performed
- struct points pnts; // points: coordinates, number etc.
+ /* Vector layer and module */
+ struct Map_info map; // input vector map
+ struct GModule *module;
+ struct nna_par xD; // 2D or 3D NNA to be performed
+ struct points pnts; // points: coordinates, number etc.
- struct nearest nna; // structure to save results
- int field;
+ struct nearest nna; // structure to save results
+ int field;
- struct {
- struct Option *map, *A, *type, *field, *lyr, *desc, *zcol; /* A - area (minimum enclosing rectangle or specified by user) */
- } opt;
+ 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;
+ 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
+ /* 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.");
- 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");
+ /* Setting options */
+ opt.map = G_define_standard_option(G_OPT_V_INPUT); // vector input layer
- opt.A = G_define_option();
- opt.A->key = "area";
- 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;
+ 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.field = G_define_standard_option(G_OPT_V_FIELD);
+ opt.A = G_define_option();
+ opt.A->key = "area";
+ 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.zcol = G_define_standard_option(G_OPT_DB_COLUMN);
- opt.zcol->key = "zcolumn";
- opt.zcol->required = NO;
- opt.zcol->guisection = _("Fields");
- opt.zcol->description = _("Column with z coordinate (set for 2D vectors only if 3D NNA is required to be performed)");
+ opt.field = G_define_standard_option(G_OPT_V_FIELD);
- G_gisinit(argv[0]);
+ opt.zcol = G_define_standard_option(G_OPT_DB_COLUMN);
+ opt.zcol->key = "zcolumn";
+ opt.zcol->required = NO;
+ opt.zcol->guisection = _("Fields");
+ opt.zcol->description =
+ _("Column with z coordinate (set for 2D vectors only if 3D NNA is required to be performed)");
- if (G_parser(argc, argv)) {
- exit(EXIT_FAILURE);
- }
+ G_gisinit(argv[0]);
- /* get parameters from the parser */
- field = opt.field->answer ? atoi(opt.field->answer) : -1;
-
- /* open input vector map */
- Vect_set_open_level(2);
+ if (G_parser(argc, argv)) {
+ exit(EXIT_FAILURE);
+ }
- 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); // test if vector layer is 2D or 3D
- xD.i3 = (flg.d23->answer || (xD.v3 == FALSE && opt.zcol->answer == NULL)) ? FALSE : TRUE; // if -2 flag or input layer is 2D (without zcolumn), perform 2D NNA, otherwise 3D
-
- /* Warnings */
- if (flg.d23->answer && opt.zcol->answer) { // -2 flag (2D) vs. zcolumn (3D)
- G_warning(_("Flag -2 has higher priority than zcolumn considered in 2D layer. 2D Nearest Neighbour Analysis will be performed instead of 3D. If you wish to perform 3D Nearest Neighbour Analysis, please remove flag -2."));
- }
+ /* get parameters from the parser */
+ field = opt.field->answer ? atoi(opt.field->answer) : -1;
- if (xD.v3 == TRUE && opt.zcol->answer) {
- G_warning(_("Input layer <%s> is 3D - it was not necessary to set up attribute column. 3D Nearest Neighbour Analysis is being performed..."), opt.map->answer);
- }
+ /* open input vector map */
+ Vect_set_open_level(2);
- read_points(&map, field, &xD, opt.zcol->answer, &pnts); // read coordinates of the points
+ 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);
- /* Nearest Neighbour Analysis */
- nn_average_distance_real(&xD, &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);
-}
+ /* perform 2D or 3D NNA ? */
+ xD.v3 = Vect_is_3d(&map); // test if vector layer is 2D or 3D
+ xD.i3 = (flg.d23->answer || (xD.v3 == FALSE && opt.zcol->answer == NULL)) ? FALSE : TRUE; // if -2 flag or input layer is 2D (without zcolumn), perform 2D NNA, otherwise 3D
+
+ /* Warnings */
+ if (flg.d23->answer && opt.zcol->answer) { // -2 flag (2D) vs. zcolumn (3D)
+ G_warning(_("Flag -2 has higher priority than zcolumn considered in 2D layer. 2D Nearest Neighbour Analysis will be performed instead of 3D. If you wish to perform 3D Nearest Neighbour Analysis, please remove flag -2."));
+ }
+
+ if (xD.v3 == TRUE && opt.zcol->answer) {
+ G_warning(_("Input layer <%s> is 3D - it was not necessary to set up attribute column. 3D Nearest Neighbour Analysis is being performed..."),
+ opt.map->answer);
+ }
+
+ read_points(&map, field, &xD, opt.zcol->answer, &pnts); // read coordinates of the points
+
+ /* Nearest Neighbour Analysis */
+ nn_average_distance_real(&xD, &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);
+}
Modified: grass-addons/grass7/vector/v.nnstat/mbb.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/mbb.c 2015-07-08 18:54:02 UTC (rev 65549)
+++ grass-addons/grass7/vector/v.nnstat/mbb.c 2015-07-09 08:29:00 UTC (rev 65550)
@@ -16,86 +16,89 @@
/* 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. */
+ 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);
+ 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);
+ 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);
+ v = vertices;
+ do {
+ tv = v;
+ v = v->next;
+ DELETE(vertices, tv);
+ }
+ while (v != vertices);
- FREE(te);
- FREE(tf);
- FREE(tv);
+ FREE(te);
+ FREE(tf);
+ FREE(tv);
- DELETE(edges, e);
- DELETE(faces, f);
- DELETE(vertices, v);
+ DELETE(edges, e);
+ DELETE(faces, f);
+ DELETE(vertices, v);
- FREE(edges);
- FREE(faces);
- FREE(vertices);
+ FREE(edges);
+ FREE(faces);
+ FREE(vertices);
- return;
+ return;
}
/*-------------------------------------------------------------------*/
int make3DHull(struct points *pnts, struct convex *hull)
{
- int error;
+ int error;
- ReadVertices(pnts);
+ 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.");
- }
+ 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();
+ ConstructHull();
- write_coord_faces(pnts, hull);
+ write_coord_faces(pnts, hull);
- freeMem();
+ freeMem();
- return (0);
-}
+ return (0);
+}
/*---------------------------------------------------------------------
MakeNullVertex: Makes a vertex, nulls out fields.
---------------------------------------------------------------------*/
tVertex MakeNullVertex(void)
{
- tVertex v;
+ tVertex v;
- NEW(v, tsVertex); /* If memory not allocated => out of memory*/
- v->duplicate = NULL;
- v->onhull = !ONHULL;
- v->mark = !PROCESSED;
- ADD(vertices, 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;
+ return v;
}
/*---------------------------------------------------------------------
@@ -106,24 +109,25 @@
---------------------------------------------------------------------*/
void ReadVertices(struct points *pnts)
{
- tVertex v;
- int vnum = 0;
- int i;
+ tVertex v;
+ int vnum = 0;
+ int i;
- double *r;
- r = &pnts->r[0];
+ double *r;
- 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);
+ 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);
}
/*---------------------------------------------------------------------
@@ -131,61 +135,63 @@
---------------------------------------------------------------------*/
void write_coord_faces(struct points *pnts, struct convex *hull)
{
- int n = pnts->n;
+ 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];
+ /* Pointers to vertices, edges, faces. */
+ tVertex v;
+ tFace f;
+ int nv = 0, nf = 0;
+ double *hc, *hf;
- v = vertices;
- f = faces;
+ 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];
- do {
- /* Write vertex coordinates */
- *hc = ((double)(v->v[X]));
- *(hc+1) = ((double)(v->v[Y]));
- *(hc+2) = ((double)(v->v[Z]));
+ v = vertices;
+ f = faces;
- 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]));
+ do {
+ /* Write vertex coordinates */
+ *hc = ((double)(v->v[X]));
+ *(hc + 1) = ((double)(v->v[Y]));
+ *(hc + 2) = ((double)(v->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]));
+ v = v->next;
+ hc += 3;
+ nv += 3;
- *(hf+6) = ((double)(f->vertex[2]->v[X]));
- *(hf+7) = ((double)(f->vertex[2]->v[Y]));
- *(hf+8) = ((double)(f->vertex[2]->v[Z]));
+ }
+ while (v != vertices);
- f = f->next;
- hf += 9;
- nf += 9;
-
- } while (f!=faces);
+ 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]));
- /* 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));
+ *(hf + 3) = ((double)(f->vertex[1]->v[X]));
+ *(hf + 4) = ((double)(f->vertex[1]->v[Y]));
+ *(hf + 5) = ((double)(f->vertex[1]->v[Z]));
- fflush(stdout);
+ *(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);
}
@@ -205,53 +211,53 @@
int DoubleTriangle(void)
{
- tVertex v0, v1, v2, v3;
- tFace f0, f1 = NULL;
- long int vol;
+ 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);
+ /* 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;
+ v1 = v0->next;
+ v2 = v1->next;
- /* Mark the vertices as processed. */
- v0->mark = PROCESSED;
- v1->mark = PROCESSED;
- v2->mark = PROCESSED;
+ /* 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);
+ /* 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;
+ /* 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);
+ /* 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);
}
- vol = VolumeSign(f0, v3);
- }
- /* Insure that v3 will be the first added. */
- vertices = v3;
+ /* Insure that v3 will be the first added. */
+ vertices = v3;
- return (0);
+ return (0);
}
@@ -261,43 +267,45 @@
---------------------------------------------------------------------*/
int ConstructHull(void)
{
- tVertex v, vnext;
- bool changed; /* T if addition changes hull; not used. */
- int i;
- int numVertices;
+ tVertex v, vnext;
+ bool changed; /* T if addition changes hull; not used. */
+ int i;
+ int numVertices;
- G_message(_("Constructing 3D hull..."));
+ 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 = vertices;
+ i = 0;
+ do {
+ vnext = v->next;
+ v = vnext;
+ i++;
}
- v = vnext;
- i++;
+ while (v != vertices);
+ numVertices = i;
- G_percent(i, numVertices, 1);
+ v = vertices;
+ i = 0;
+ do {
+ vnext = v->next;
+ if (!v->mark) {
+ v->mark = PROCESSED;
+ changed = AddOne(v);
+ CleanUp();
+ }
+ v = vnext;
+ i++;
- } while (v != vertices);
+ G_percent(i, numVertices, 1);
- fflush(stdout);
+ }
+ while (v != vertices);
- return numVertices;
+ fflush(stdout);
+ return numVertices;
+
}
/*---------------------------------------------------------------------
@@ -309,44 +317,46 @@
---------------------------------------------------------------------*/
bool AddOne(tVertex p)
{
- tFace f;
- tEdge e, temp;
- long int vol;
- bool vis = FALSE;
+ tFace f;
+ tEdge e, temp;
+ long int vol;
+ bool vis = FALSE;
- /* Mark faces visible from p. */
- f = faces;
- do {
- vol = VolumeSign(f, p);
+ /* Mark faces visible from p. */
+ f = faces;
+ do {
+ vol = VolumeSign(f, p);
- if (vol < 0) {
- f->visible = VISIBLE;
- vis = TRUE;
+ if (vol < 0) {
+ f->visible = VISIBLE;
+ vis = TRUE;
+ }
+ f = f->next;
}
- f = f->next;
- } while (f != faces);
+ while (f != faces);
- /* If no faces are visible from p, then p is inside the hull. */
- if (!vis) {
- p->onhull = !ONHULL;
- return FALSE;
- }
+ /* 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;
+ /* 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;
}
/*---------------------------------------------------------------------
@@ -358,30 +368,29 @@
---------------------------------------------------------------------*/
int VolumeSign(tFace f, tVertex p)
{
- double vol;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
+ 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];
+ 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);
+ 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;
+ /* The volume should be an integer. */
+ if (vol > 0.0)
+ return 1;
+ else if (vol < -0.0)
+ return -1;
+ else
+ return 0;
}
@@ -392,38 +401,38 @@
---------------------------------------------------------------------*/
tFace MakeConeFace(tEdge e, tVertex p)
{
- tEdge new_edge[2];
- tFace new_face;
- int i, j;
+ 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 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);
+ /* 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;
- }
+ /* 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;
+ return new_face;
}
/*---------------------------------------------------------------------
@@ -433,33 +442,33 @@
---------------------------------------------------------------------*/
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 */
+ 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];
+ 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. */
+ /* 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;
+ f->vertex[2] = p;
}
/*---------------------------------------------------------------------
@@ -468,14 +477,14 @@
---------------------------------------------------------------------*/
tEdge MakeNullEdge(void)
{
- tEdge e;
+ 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;
+ 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;
}
/*--------------------------------------------------------------------
@@ -485,17 +494,17 @@
---------------------------------------------------------------------*/
tFace MakeNullFace(void)
{
- tFace f;
- int i;
+ 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;
+ 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;
}
/*---------------------------------------------------------------------
@@ -504,40 +513,40 @@
---------------------------------------------------------------------*/
tFace MakeFace(tVertex v0, tVertex v1, tVertex v2, tFace fold)
{
- tFace f;
- tEdge e0, e1, e2;
+ 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 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;
+ /* 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;
+ /* Link edges to face. */
+ e0->adjface[0] = e1->adjface[0] = e2->adjface[0] = f;
- return f;
+ return f;
}
/*---------------------------------------------------------------------
@@ -547,11 +556,11 @@
---------------------------------------------------------------------*/
void CleanUp(void)
{
- CleanEdges();
- CleanFaces();
- CleanVertices();
+ CleanEdges();
+ CleanFaces();
+ CleanVertices();
- return;
+ return;
}
/*---------------------------------------------------------------------
@@ -561,40 +570,42 @@
---------------------------------------------------------------------*/
void CleanEdges(void)
{
- tEdge e; /* Primary index into edge list. */
- tEdge t; /* Temporary edge pointer. */
+ 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;
+ /* 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;
}
- e = e->next;
- } while (e != edges);
+ 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);
+ /* Delete any edges marked for deletion. */
+ while (edges && edges->del) {
+ e = edges;
+ DELETE(edges, e);
}
- else
- e = e->next;
- } while (e != edges);
+ e = edges->next;
+ do {
+ if (e->del) {
+ t = e;
+ e = e->next;
+ DELETE(edges, t);
+ }
+ else
+ e = e->next;
+ }
+ while (e != edges);
- return;
+ return;
}
/*---------------------------------------------------------------------
@@ -602,26 +613,27 @@
---------------------------------------------------------------------*/
void CleanFaces(void)
{
- tFace f; /* Primary pointer into face list. */
- tFace t; /* Temporary pointer, for deleting. */
+ 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);
+ while (faces && faces->visible) {
+ f = faces;
+ DELETE(faces, f);
}
- else
- f = f->next;
- } while (f != faces);
+ f = faces->next;
+ do {
+ if (f->visible) {
+ t = f;
+ f = f->next;
+ DELETE(faces, t);
+ }
+ else
+ f = f->next;
+ }
+ while (f != faces);
- return;
+ return;
}
/*---------------------------------------------------------------------
@@ -631,42 +643,44 @@
---------------------------------------------------------------------*/
void CleanVertices(void)
{
- tEdge e;
- tVertex v, t;
+ 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);
+ /* 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) {
+ /* 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;
- 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);
+ do {
+ v->duplicate = NULL;
+ v->onhull = !ONHULL;
+ v = v->next;
+ }
+ while (v != vertices);
- /* Reset flags. */
- v = vertices;
- do {
- v->duplicate = NULL;
- v->onhull = !ONHULL;
- v = v->next;
- } while (v != vertices);
-
- return;
+ return;
}
/*---------------------------------------------------------------------
@@ -675,25 +689,25 @@
---------------------------------------------------------------------*/
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;
+ 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;
+ int error;
- error = make3DHull(pnts, hull); /* make 3D hull */
- if (error < 0) {
- G_fatal_error(_("Simple planar hulls not implemented yet"));
- }
+ error = make3DHull(pnts, hull); /* make 3D hull */
+ if (error < 0) {
+ G_fatal_error(_("Simple planar hulls not implemented yet"));
+ }
- return;
+ return;
}
/* ----------------------------
@@ -701,60 +715,77 @@
*/
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));
+ 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;
- 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 */
+ r_min = (double *)G_malloc(3 * sizeof(double));
+ r_max = (double *)G_malloc(3 * sizeof(double));
- 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;
+ double *hull_trans; /* Coordinates of hull vertices transformed into coordinate system with axes parallel to hull's edges */
- /* 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
+ hull_trans = (double *)malloc(3 * sizeof(double));
+ ht = &hull_trans[0];
- 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
+ struct convex hull;
- return V_min;
-}
+ 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;
+}
Modified: grass-addons/grass7/vector/v.nnstat/mbr.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/mbr.c 2015-07-08 18:54:02 UTC (rev 65549)
+++ grass-addons/grass7/vector/v.nnstat/mbr.c 2015-07-09 08:29:00 UTC (rev 65550)
@@ -23,129 +23,128 @@
int rightTurn(double (*P)[3], int i, int j, int k)
{
- double a, b, c, d;
+ 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 */
+ /* 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;
+ 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;
- }
+ 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;
+ int pointIdx, upPoints, loPoints;
+ int i, *upHull, *loHull;
+ int n = pnts->n;
- double *ro, (*ri)[3], (*r)[3]; /* r = [xyz] */
- r = (double (*)[3]) G_malloc(n * 3 * sizeof(double));
- ri = &r[0];
- ro = &pnts->r[0];
+ double *ro, (*ri)[3], (*r)[3]; /* r = [xyz] */
- 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);
+ r = (double (*)[3])G_malloc(n * 3 * sizeof(double));
+ ri = &r[0];
+ ro = &pnts->r[0];
- 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--;
+ for (i = 0; i < n; i++) {
+ (*ri)[0] = *ro;
+ (*ri)[1] = *(ro + 1);
+ (*ri)[2] = *(ro + 2);
+ ri++;
+ ro += 3;
}
- }
- /* 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--;
+ /* 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--;
+ }
}
- }
- 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));
+ /* 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;
- /* Obtain coordinates of hull vertices */
- hull->coord = (double *) G_malloc((hull->n + 1) * 3 * sizeof(double)); /* 1st = last pnt */
+ G_debug(3, "numPoints:%d loPoints:%d upPoints:%d", n, loPoints, upPoints);
- int *hh, *hh0;
- double *hc, (*r0)[3];
- hc = &hull->coord[0];
- hh = &hull->hull[0];
- ri = &r[0];
+ /* reclaim uneeded memory */
+ hull->hull = (int *)G_realloc(hull->hull, (hull->n + 1) * sizeof(int));
- for (i=0; i <= hull->n; i++) {
- if (i < hull->n) {
- *hc = (*(ri+*hh))[0];
- *(hc+1) = (*(ri+*hh))[1];
- *(hc+2) = (*(ri+*hh))[2];
- hc += 3;
- hh++;
+ /* 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];
+ ri = &r[0];
+
+ for (i = 0; i <= hull->n; i++) {
+ if (i < hull->n) {
+ *hc = (*(ri + *hh))[0];
+ *(hc + 1) = (*(ri + *hh))[1];
+ *(hc + 2) = (*(ri + *hh))[2];
+ 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];
+ }
}
- 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;
+ return hull->n;
}
/* ----------------------------
@@ -153,60 +152,67 @@
* ----------------------------*/
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));
+ 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;
- 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];
+ r_min = (double *)G_malloc(3 * sizeof(double));
+ r_max = (double *)G_malloc(3 * sizeof(double));
- struct convex hull;
- hull.n = convexHull(pnts, &hull);
- hc = &hull.coord[0];
+ double *hull_trans, *ht; /* Coordinates of hull vertices transformed into coordinate system with axes parallel to hull's edges */
- 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);
+ hull_trans = (double *)G_malloc(3 * sizeof(double));
+ ht = &hull_trans[0];
- 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);
+ struct convex hull;
- /* 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)));
- break;
- }
- hc_k += 3;
- } // end k
+ hull.n = convexHull(pnts, &hull);
+ hc = &hull.coord[0];
- 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
+ S_min = (pnts->r_max[0] - pnts->r_min[0]) * (pnts->r_max[1] - pnts->r_min[1]); /* Area of extent */
- G_free(r_min);
- G_free(r_max);
+ 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);
- return S_min;
+ 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)));
+ break;
+ }
+ 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
+
+ G_free(r_min);
+ G_free(r_max);
+
+ return S_min;
}
Modified: grass-addons/grass7/vector/v.nnstat/nearest_neighbour.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/nearest_neighbour.c 2015-07-08 18:54:02 UTC (rev 65549)
+++ grass-addons/grass7/vector/v.nnstat/nearest_neighbour.c 2015-07-09 08:29:00 UTC (rev 65550)
@@ -1,57 +1,61 @@
#include "local_proto.h"
-void nn_average_distance_real(struct nna_par *xD, struct points *pnts, struct nearest *nna)
+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
+ 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;
+ 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
+ G_message(_("Computing average distance between nearest neighbors..."));
+ sum_r = 0.;
- return;
+ 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)
+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 == 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"));
+ 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;
+ 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
- }
+ 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
+ }
}
/* ---------------------------------
@@ -60,70 +64,75 @@
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);
+ 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"));
- }
+ /* 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 (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(" ");
+ 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"));
+ }
- return;
+ 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)
+void
+nn_statistics(struct points *pnts, struct nna_par *xD, struct nearest *nna)
{
- double disp_rE, var_rE, sig_rE;
+ 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)
- }
+ 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
+ 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
+ nn_results(pnts, xD, nna); // write the result to the command line
- return;
+ return;
}
Modified: grass-addons/grass7/vector/v.nnstat/spatial_index.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/spatial_index.c 2015-07-08 18:54:02 UTC (rev 65549)
+++ grass-addons/grass7/vector/v.nnstat/spatial_index.c 2015-07-09 08:29:00 UTC (rev 65550)
@@ -2,145 +2,143 @@
struct RTree *create_spatial_index(struct nna_par *xD)
{
- struct RTree *R_tree; // spatial index
- struct RTree_Rect *rect; // rectangle
+ 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
- }
+ 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);
+ 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
+ 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
+ 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 == 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
- }
+ 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
+ RTreeInsertRect(rect, i + 1, pnts->R_tree); // insert rectangle to the spatial index
+ RTreeFreeRect(rect); // free rectangle memory
- return;
+ return;
}
/* find neighbors in cubic (square) surroundings */
-struct ilist *spatial_search(int i3, int i, struct points *pnts, double max_dist)
+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];
+ // local variables
+ struct RTree *R_tree = pnts->R_tree; // spatial index
+ double *r; // pointer to vector of the coordinates
- struct ilist *list;
- struct RTree_Rect *search; // search rectangle
+ r = &pnts->r[3 * i];
- list = G_new_ilist(); // create new list
- search = RTreeAllocRect(R_tree); // allocate new rectangle
+ struct ilist *list;
+ struct RTree_Rect *search; // search 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
- }
+ list = G_new_ilist(); // create new list
+ search = RTreeAllocRect(R_tree); // allocate new 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);
+ 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
+ }
- return list;
+ 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)
+ // 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
+ 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
+ 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);
+ // check spherical (circular) surrounding
+ list = spatial_search(i3, i, pnts, sqrt(2.) * max_dist);
- return list;
+ 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)
+ int j;
+ int n_vals = list->n_values; // # of selected points
+ int *value; // pointer to indices of selected points
- 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
+ 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)
- 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
+ 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
- if (i3 == TRUE) { // 3D NNA:
- *d += dz*dz; // use also the difference in z
+ 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
}
-
- 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)
+ 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
+ G_free(sqDist); // free memory of the vector of squared distances
+ G_free_ilist(list); // free list of selected points
- return dist;
+ return dist;
}
Modified: grass-addons/grass7/vector/v.nnstat/utils.c
===================================================================
--- grass-addons/grass7/vector/v.nnstat/utils.c 2015-07-08 18:54:02 UTC (rev 65549)
+++ grass-addons/grass7/vector/v.nnstat/utils.c 2015-07-09 08:29:00 UTC (rev 65550)
@@ -3,63 +3,65 @@
/* make triple of the coordinates */
double *triple(double x, double y, double z)
{
- double *t;
- t = (double *) G_malloc(3 * sizeof(double));
+ double *t;
- *t = x;
- *(t+1) = y;
- *(t+2) = z;
-
- return 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;
- }
+ double dy, dx, us;
- 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;
- }
+ dy = y1 - y0;
+ dx = x1 - x0;
+ if (dy == 0. && dx == 0.) {
+ return -9999;
+ }
- return us;
+ 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;
+ 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;
- }
+ 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