|
|
@@ -1,273 +1,37 @@
|
|
|
|
|
|
-/******************************************************************************
|
|
|
-* hull.c <s.hull>
|
|
|
-* Creates the convex hull surrounding a sites list
|
|
|
-
|
|
|
-* @Copyright Andrea Aime <aaime@libero.it>
|
|
|
-* 23 Sept. 2001
|
|
|
-* Updated 19 Dec 2003, Markus Neteler to 5.7
|
|
|
-* Last updated 16 jan 2007, Benjamin Ducke to support 3D hull creation
|
|
|
-
|
|
|
-* This file is part of GRASS GIS. It is free software. You can
|
|
|
-* redistribute it and/or modify it under the terms of
|
|
|
-* the GNU General Public License as published by the Free Software
|
|
|
-* Foundation; either version 2 of the License, or (at your option)
|
|
|
-* any later version.
|
|
|
-
|
|
|
-* This program is distributed in the hope that it will be useful,
|
|
|
-* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
-* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
-* GNU General Public License for more details.
|
|
|
-******************************************************************************/
|
|
|
-
|
|
|
+/****************************************************************
|
|
|
+ *
|
|
|
+ * MODULE: v.hull
|
|
|
+ *
|
|
|
+ * AUTHOR(S): Andrea Aime <aaime@libero.it>
|
|
|
+ * Updated 19 Dec 2003, Markus Neteler to 5.7
|
|
|
+ * Last updated 16 jan 2007, Benjamin Ducke to support 3D hull creation
|
|
|
+ * OGR support by Martin Landa <landa.martin gmail.com> (2009)
|
|
|
+ *
|
|
|
+ * PURPOSE: Creates the convex hull surrounding a vector points.
|
|
|
+ *
|
|
|
+ * COPYRIGHT: (C) 2001-2009 by 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 <stdlib.h>
|
|
|
#include <stdio.h>
|
|
|
#include <assert.h>
|
|
|
+
|
|
|
#include <grass/gis.h>
|
|
|
#include <grass/vector.h>
|
|
|
#include <grass/glocale.h>
|
|
|
|
|
|
-#include "chull.h"
|
|
|
-
|
|
|
-struct Point
|
|
|
-{
|
|
|
- double x;
|
|
|
- double y;
|
|
|
- double z;
|
|
|
-};
|
|
|
-
|
|
|
-int rightTurn(struct Point *P, int i, int j, int k)
|
|
|
-{
|
|
|
- double a, b, c, d;
|
|
|
-
|
|
|
- a = P[i].x - P[j].x;
|
|
|
- b = P[i].y - P[j].y;
|
|
|
- c = P[k].x - P[j].x;
|
|
|
- d = P[k].y - P[j].y;
|
|
|
- return a * d - b * c < 0;
|
|
|
-}
|
|
|
-
|
|
|
-int cmpPoints(const void *v1, const void *v2)
|
|
|
-{
|
|
|
- struct Point *p1, *p2;
|
|
|
-
|
|
|
- p1 = (struct Point *)v1;
|
|
|
- p2 = (struct Point *)v2;
|
|
|
- if (p1->x > p2->x)
|
|
|
- return 1;
|
|
|
- else if (p1->x < p2->x)
|
|
|
- return -1;
|
|
|
- else
|
|
|
- return 0;
|
|
|
-}
|
|
|
-
|
|
|
-int convexHull(struct Point *P, const int numPoints, int **hull)
|
|
|
-{
|
|
|
- int pointIdx, upPoints, loPoints;
|
|
|
- int *upHull, *loHull;
|
|
|
-
|
|
|
- /* sort points in ascending x order */
|
|
|
- qsort(P, numPoints, sizeof(struct Point), cmpPoints);
|
|
|
-
|
|
|
- *hull = (int *)G_malloc(numPoints * 2 * sizeof(int));
|
|
|
-
|
|
|
- /* compute upper hull */
|
|
|
- upHull = *hull;
|
|
|
- upHull[0] = 0;
|
|
|
- upHull[1] = 1;
|
|
|
- upPoints = 1;
|
|
|
- for (pointIdx = 2; pointIdx < numPoints; pointIdx++) {
|
|
|
- upPoints++;
|
|
|
- upHull[upPoints] = pointIdx;
|
|
|
- while (upPoints > 1 &&
|
|
|
- !rightTurn(P, upHull[upPoints], upHull[upPoints - 1],
|
|
|
- upHull[upPoints - 2])
|
|
|
- ) {
|
|
|
- upHull[upPoints - 1] = upHull[upPoints];
|
|
|
- upPoints--;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- /*
|
|
|
- printf("upPoints: %d\n", upPoints);
|
|
|
- for(pointIdx = 0; pointIdx <= upPoints; pointIdx ++)
|
|
|
- printf("%d ", upHull[pointIdx]);
|
|
|
- printf("\n");
|
|
|
- */
|
|
|
-
|
|
|
- /* compute lower hull, overwrite last point of upper hull */
|
|
|
- loHull = &(upHull[upPoints]);
|
|
|
- loHull[0] = numPoints - 1;
|
|
|
- loHull[1] = numPoints - 2;
|
|
|
- loPoints = 1;
|
|
|
- for (pointIdx = numPoints - 3; pointIdx >= 0; pointIdx--) {
|
|
|
- loPoints++;
|
|
|
- loHull[loPoints] = pointIdx;
|
|
|
- while (loPoints > 1 &&
|
|
|
- !rightTurn(P, loHull[loPoints], loHull[loPoints - 1],
|
|
|
- loHull[loPoints - 2])
|
|
|
- ) {
|
|
|
- loHull[loPoints - 1] = loHull[loPoints];
|
|
|
- loPoints--;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- G_debug(3, "numPoints:%d loPoints:%d upPoints:%d",
|
|
|
- numPoints, loPoints, upPoints);
|
|
|
- /*
|
|
|
- printf("loPoints: %d\n", loPoints);
|
|
|
- for(pointIdx = 0; pointIdx <= loPoints; pointIdx ++)
|
|
|
- printf("%d ", loHull[pointIdx]);
|
|
|
- printf("\n");
|
|
|
- */
|
|
|
-
|
|
|
- /* reclaim uneeded memory */
|
|
|
- *hull = (int *)G_realloc(*hull, (loPoints + upPoints) * sizeof(int));
|
|
|
- return loPoints + upPoints;
|
|
|
-}
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-void convexHull3d(struct Point *P, const int numPoints, struct Map_info *Map)
|
|
|
-{
|
|
|
-
|
|
|
- int error;
|
|
|
- int i;
|
|
|
- double *px;
|
|
|
- double *py;
|
|
|
- double *pz;
|
|
|
-
|
|
|
- px = G_malloc(sizeof(double) * numPoints);
|
|
|
- py = G_malloc(sizeof(double) * numPoints);
|
|
|
- pz = G_malloc(sizeof(double) * numPoints);
|
|
|
-
|
|
|
- for (i = 0; i < numPoints; i++) {
|
|
|
- px[i] = (P)[i].x;
|
|
|
- py[i] = (P)[i].y;
|
|
|
- pz[i] = (P)[i].z;
|
|
|
- }
|
|
|
-
|
|
|
- /* make 3D hull */
|
|
|
- error = make3DHull(px, py, pz, numPoints, Map);
|
|
|
- if (error < 0) {
|
|
|
- G_fatal_error("Simple planar hulls not implemented yet");
|
|
|
- }
|
|
|
-
|
|
|
- G_free(px);
|
|
|
- G_free(py);
|
|
|
- G_free(pz);
|
|
|
-
|
|
|
-}
|
|
|
-
|
|
|
-#define ALLOC_CHUNK 256
|
|
|
-int loadSiteCoordinates(struct Map_info *Map, struct Point **points, int all,
|
|
|
- struct Cell_head *window)
|
|
|
-{
|
|
|
- int pointIdx = 0;
|
|
|
- struct line_pnts *sites;
|
|
|
- struct line_cats *cats;
|
|
|
- struct bound_box box;
|
|
|
- int cat, type;
|
|
|
-
|
|
|
- sites = Vect_new_line_struct();
|
|
|
- cats = Vect_new_cats_struct();
|
|
|
-
|
|
|
- *points = NULL;
|
|
|
-
|
|
|
- /* copy window to box */
|
|
|
- Vect_region_box(window, &box);
|
|
|
-
|
|
|
- while ((type = Vect_read_next_line(Map, sites, cats)) > -1) {
|
|
|
-
|
|
|
- if (type != GV_POINT)
|
|
|
- continue;
|
|
|
-
|
|
|
- Vect_cat_get(cats, 1, &cat);
|
|
|
-
|
|
|
- G_debug(4, "Point: %f|%f|%f|#%d", sites->x[0], sites->y[0],
|
|
|
- sites->z[0], cat);
|
|
|
-
|
|
|
- if (all ||
|
|
|
- Vect_point_in_box(sites->x[0], sites->y[0], sites->z[0], &box)) {
|
|
|
-
|
|
|
- G_debug(4, "Point in the box");
|
|
|
-
|
|
|
- if ((pointIdx % ALLOC_CHUNK) == 0)
|
|
|
- *points =
|
|
|
- (struct Point *)G_realloc(*points,
|
|
|
- (pointIdx +
|
|
|
- ALLOC_CHUNK) *
|
|
|
- sizeof(struct Point));
|
|
|
-
|
|
|
- (*points)[pointIdx].x = sites->x[0];
|
|
|
- (*points)[pointIdx].y = sites->y[0];
|
|
|
- (*points)[pointIdx].z = sites->z[0];
|
|
|
- pointIdx++;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- if (pointIdx > 0)
|
|
|
- *points =
|
|
|
- (struct Point *)G_realloc(*points,
|
|
|
- (pointIdx + 1) * sizeof(struct Point));
|
|
|
- return pointIdx;
|
|
|
-}
|
|
|
-
|
|
|
-/*
|
|
|
- * Outputs the points that comprises the convex hull as a single closed line
|
|
|
- * and the hull baricenter as the label points (as it is a linear combination
|
|
|
- * of points on the hull is guaranteed to be inside the hull, follow from the
|
|
|
- * definition of convex polygon)
|
|
|
- */
|
|
|
-int outputHull(struct Map_info *Map, struct Point *P, int *hull,
|
|
|
- int numPoints)
|
|
|
-{
|
|
|
- struct line_pnts *Points;
|
|
|
- struct line_cats *Cats;
|
|
|
- double *tmpx, *tmpy;
|
|
|
- int i, pointIdx;
|
|
|
- double xc, yc;
|
|
|
-
|
|
|
- tmpx = (double *)G_malloc((numPoints + 1) * sizeof(double));
|
|
|
- tmpy = (double *)G_malloc((numPoints + 1) * sizeof(double));
|
|
|
-
|
|
|
- xc = yc = 0;
|
|
|
- for (i = 0; i < numPoints; i++) {
|
|
|
- pointIdx = hull[i];
|
|
|
- tmpx[i] = P[pointIdx].x;
|
|
|
- tmpy[i] = P[pointIdx].y;
|
|
|
- /* average coordinates calculation... may introduce a little
|
|
|
- numerical error but guaratees that no overflow will occurr */
|
|
|
- xc = xc + tmpx[i] / numPoints;
|
|
|
- yc = yc + tmpy[i] / numPoints;
|
|
|
- }
|
|
|
- tmpx[numPoints] = P[hull[0]].x;
|
|
|
- tmpy[numPoints] = P[hull[0]].y;
|
|
|
-
|
|
|
- Points = Vect_new_line_struct();
|
|
|
- Cats = Vect_new_cats_struct();
|
|
|
- Vect_copy_xyz_to_pnts(Points, tmpx, tmpy, 0, numPoints + 1);
|
|
|
- G_free(tmpx);
|
|
|
- G_free(tmpy);
|
|
|
-
|
|
|
- /* write out convex hull */
|
|
|
- Vect_write_line(Map, GV_BOUNDARY, Points, Cats);
|
|
|
-
|
|
|
- /* find and add centroid */
|
|
|
- Vect_reset_line(Points);
|
|
|
- Vect_append_point(Points, xc, yc, 0.0);
|
|
|
- Vect_cat_set(Cats, 1, 1);
|
|
|
- Vect_write_line(Map, GV_CENTROID, Points, Cats);
|
|
|
- Vect_destroy_line_struct(Points);
|
|
|
-
|
|
|
- return 0;
|
|
|
-}
|
|
|
+#include "hull.h"
|
|
|
|
|
|
int main(int argc, char **argv)
|
|
|
{
|
|
|
struct GModule *module;
|
|
|
- struct Option *input, *output;
|
|
|
+ struct Option *input, *output, *field;
|
|
|
struct Flag *all, *flat;
|
|
|
struct Cell_head window;
|
|
|
|
|
|
@@ -279,21 +43,20 @@ int main(int argc, char **argv)
|
|
|
int numSitePoints, numHullPoints;
|
|
|
|
|
|
int MODE2D;
|
|
|
-
|
|
|
-
|
|
|
+
|
|
|
G_gisinit(argv[0]);
|
|
|
|
|
|
module = G_define_module();
|
|
|
G_add_keyword(_("vector"));
|
|
|
G_add_keyword(_("geometry"));
|
|
|
module->description =
|
|
|
- _("Uses a GRASS vector points map to produce a convex hull vector map.");
|
|
|
+ _("Creates convex hull for vector point map.");
|
|
|
|
|
|
input = G_define_standard_option(G_OPT_V_INPUT);
|
|
|
- input->description = _("Name of input vector points map");
|
|
|
+
|
|
|
+ field = G_define_standard_option(G_OPT_V_FIELD);
|
|
|
|
|
|
output = G_define_standard_option(G_OPT_V_OUTPUT);
|
|
|
- output->description = _("Name of output vector area map");
|
|
|
|
|
|
all = G_define_flag();
|
|
|
all->key = 'a';
|
|
|
@@ -313,13 +76,13 @@ int main(int argc, char **argv)
|
|
|
Vect_check_input_output_name(input->answer, output->answer,
|
|
|
GV_FATAL_EXIT);
|
|
|
|
|
|
- /* open site file */
|
|
|
- if (Vect_open_old(&Map, sitefile, "") < 0)
|
|
|
- G_fatal_error(_("Cannot open vector map <%s>"), sitefile);
|
|
|
+ if (Vect_open_old2(&Map, sitefile, "", field->answer) < 0)
|
|
|
+ G_fatal_error(_("Unable to open vector map <%s>"), sitefile);
|
|
|
|
|
|
/* load site coordinates */
|
|
|
G_get_window(&window);
|
|
|
- numSitePoints = loadSiteCoordinates(&Map, &points, all->answer, &window);
|
|
|
+ numSitePoints = loadSiteCoordinates(&Map, &points, all->answer, &window,
|
|
|
+ Vect_get_field_number(&Map, field->answer));
|
|
|
if (numSitePoints < 0)
|
|
|
G_fatal_error(_("Error loading vector points map <%s>"), sitefile);
|
|
|
|
|
|
@@ -340,12 +103,12 @@ int main(int argc, char **argv)
|
|
|
/* create vector map */
|
|
|
if (MODE2D) {
|
|
|
if (0 > Vect_open_new(&Map, output->answer, 0)) {
|
|
|
- G_fatal_error(_("Cannot open vector map <%s>"), output->answer);
|
|
|
+ G_fatal_error(_("Unable to create vector map <%s>"), output->answer);
|
|
|
}
|
|
|
}
|
|
|
else {
|
|
|
if (0 > Vect_open_new(&Map, output->answer, 1)) {
|
|
|
- G_fatal_error(_("Cannot open vector map <%s>"), output->answer);
|
|
|
+ G_fatal_error(_("Unable to create vector map <%s>"), output->answer);
|
|
|
}
|
|
|
}
|
|
|
|
|
|
@@ -358,16 +121,12 @@ int main(int argc, char **argv)
|
|
|
|
|
|
/* output vector map */
|
|
|
outputHull(&Map, points, hull, numHullPoints);
|
|
|
-
|
|
|
}
|
|
|
else {
|
|
|
-
|
|
|
/* this does everything for the 3D hull including vector map creation */
|
|
|
convexHull3d(points, numSitePoints, &Map);
|
|
|
-
|
|
|
}
|
|
|
-
|
|
|
-
|
|
|
+
|
|
|
/* clean up and bye bye */
|
|
|
Vect_build(&Map);
|
|
|
Vect_close(&Map);
|
Martin Landa