[GRASS-SVN] r57984 - grass/trunk/vector/v.voronoi
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Oct 13 08:24:04 PDT 2013
Author: mmetz
Date: 2013-10-13 08:24:04 -0700 (Sun, 13 Oct 2013)
New Revision: 57984
Added:
grass/trunk/vector/v.voronoi/clean_topo.c
grass/trunk/vector/v.voronoi/skeleton.c
grass/trunk/vector/v.voronoi/v_voronoi_areas.png
grass/trunk/vector/v.voronoi/v_voronoi_points.png
grass/trunk/vector/v.voronoi/v_voronoi_skeleton.png
Removed:
grass/trunk/vector/v.voronoi/v_voronoi_delaunay.png
Modified:
grass/trunk/vector/v.voronoi/Makefile
grass/trunk/vector/v.voronoi/defs.h
grass/trunk/vector/v.voronoi/main.c
grass/trunk/vector/v.voronoi/sw_defs.h
grass/trunk/vector/v.voronoi/sw_main.c
grass/trunk/vector/v.voronoi/sw_voronoi.c
grass/trunk/vector/v.voronoi/v.voronoi.html
grass/trunk/vector/v.voronoi/vo_write.c
Log:
v.voronoi: abuse it for skeletons
Modified: grass/trunk/vector/v.voronoi/Makefile
===================================================================
--- grass/trunk/vector/v.voronoi/Makefile 2013-10-11 13:53:35 UTC (rev 57983)
+++ grass/trunk/vector/v.voronoi/Makefile 2013-10-13 15:24:04 UTC (rev 57984)
@@ -8,7 +8,7 @@
v_voronoi_OBJS = vo_main.o vo_extend.o vo_write.o $(SWEEP_OBJS)
-LIBES = $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB)
+LIBES = $(VECTORLIB) $(DIG2LIB) $(DBMILIB) $(GISLIB) $(MATHLIB)
DEPENDENCIES = $(VECTORDEP) $(DBMIDEP) $(GISDEP)
EXTRA_INC = $(VECT_INC)
EXTRA_CFLAGS = $(VECT_CFLAGS)
Added: grass/trunk/vector/v.voronoi/clean_topo.c
===================================================================
--- grass/trunk/vector/v.voronoi/clean_topo.c (rev 0)
+++ grass/trunk/vector/v.voronoi/clean_topo.c 2013-10-13 15:24:04 UTC (rev 57984)
@@ -0,0 +1,149 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include "sw_defs.h"
+#include "defs.h"
+
+int clean_topo(void)
+{
+ int type, line, nlines;
+ int area, nareas;
+ int verbose;
+ int err_boundaries, err_centr_out, err_centr_dupl, err_nocentr;
+ struct bound_box box;
+ double snap_thresh;
+
+ /* cleaning part 1: count errors */
+ G_message(_("Searching for topology errors..."));
+ verbose = G_verbose();
+ G_set_verbose(0);
+ Vect_build_partial(&Out, GV_BUILD_CENTROIDS);
+ G_set_verbose(verbose);
+ err_boundaries = err_centr_out = err_centr_dupl = err_nocentr = 0;
+ nlines = Vect_get_num_lines(&Out);
+ for (line = 1; line <= nlines; line++) {
+
+ if (!Vect_line_alive(&Out, line))
+ continue;
+
+ type = Vect_get_line_type(&Out, line);
+ if (type == GV_BOUNDARY) {
+ int left, right;
+
+ Vect_get_line_areas(&Out, line, &left, &right);
+
+ if (left == 0 || right == 0) {
+ G_debug(3, "line = %d left = %d right = %d", line,
+ left, right);
+ err_boundaries++;
+ }
+ }
+ if (type == GV_CENTROID) {
+ area = Vect_get_centroid_area(&Out, line);
+ if (area == 0)
+ err_centr_out++;
+ else if (area < 0)
+ err_centr_dupl++;
+ }
+ }
+
+ err_nocentr = 0;
+ nareas = Vect_get_num_areas(&Out);
+ for (area = 1; area <= nareas; area++) {
+ if (!Vect_area_alive(&Out, area))
+ continue;
+ line = Vect_get_area_centroid(&Out, area);
+ if (line == 0)
+ err_nocentr++;
+ }
+
+ /* cleaning part 2: snap */
+ /* TODO: adjust snapping treshold to ULP */
+ Vect_get_map_box(&Out, &box);
+ snap_thresh = fabs(box.W);
+ if (snap_thresh < fabs(box.E))
+ snap_thresh = fabs(box.E);
+ if (snap_thresh < fabs(box.N))
+ snap_thresh = fabs(box.N);
+ if (snap_thresh < fabs(box.S))
+ snap_thresh = fabs(box.S);
+ snap_thresh = d_ulp(snap_thresh);
+
+ if (err_nocentr || err_centr_dupl || err_centr_out) {
+ int nmod;
+
+ G_important_message(_("Cleaning output topology"));
+ Vect_snap_lines(&Out, GV_BOUNDARY, snap_thresh, NULL);
+ do {
+ Vect_break_lines(&Out, GV_BOUNDARY, NULL);
+ Vect_remove_duplicates(&Out, GV_BOUNDARY, NULL);
+ nmod =
+ Vect_clean_small_angles_at_nodes(&Out, GV_BOUNDARY, NULL);
+ } while (nmod > 0);
+
+ G_message(_("Removing dangles..."));
+ Vect_remove_dangles(&Out, GV_BOUNDARY, -1.0, NULL);
+ G_message(_("Removing bridges..."));
+ Vect_remove_bridges(&Out, NULL, NULL, NULL);
+
+ err_boundaries = 0;
+ nlines = Vect_get_num_lines(&Out);
+ for (line = 1; line <= nlines; line++) {
+
+ if (!Vect_line_alive(&Out, line))
+ continue;
+
+ type = Vect_get_line_type(&Out, line);
+ if (type == GV_BOUNDARY) {
+ int left, right;
+
+ Vect_get_line_areas(&Out, line, &left, &right);
+
+ if (left == 0 || right == 0) {
+ G_debug(3, "line = %d left = %d right = %d", line,
+ left, right);
+ err_boundaries++;
+ }
+ }
+ }
+ }
+ /* cleaning part 3: remove remaining incorrect boundaries */
+ if (err_boundaries) {
+ G_important_message(_("Removing incorrect boundaries from output"));
+ nlines = Vect_get_num_lines(&Out);
+ for (line = 1; line <= nlines; line++) {
+
+ if (!Vect_line_alive(&Out, line))
+ continue;
+
+ type = Vect_get_line_type(&Out, line);
+ if (type == GV_BOUNDARY) {
+ int left, right;
+
+ Vect_get_line_areas(&Out, line, &left, &right);
+
+ /* &&, not ||, no typo */
+ if (left == 0 && right == 0) {
+ G_debug(3, "line = %d left = %d right = %d", line,
+ left, right);
+ Vect_delete_line(&Out, line);
+ }
+ }
+ }
+ }
+
+ /* this is slow:
+ if (in_area) {
+ if (Type == GV_LINE)
+ G_message(_("Merging lines ..."));
+ else
+ G_message(_("Merging boundaries ..."));
+ Vect_merge_lines(&Out, Type, NULL, NULL);
+ }
+ */
+
+ return 1;
+}
Property changes on: grass/trunk/vector/v.voronoi/clean_topo.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass/trunk/vector/v.voronoi/defs.h
===================================================================
--- grass/trunk/vector/v.voronoi/defs.h 2013-10-11 13:53:35 UTC (rev 57983)
+++ grass/trunk/vector/v.voronoi/defs.h 2013-10-13 15:24:04 UTC (rev 57984)
@@ -5,4 +5,5 @@
extern int Type;
extern int Field;
extern int in_area;
+extern int skeleton;
extern double segf;
Modified: grass/trunk/vector/v.voronoi/main.c
===================================================================
--- grass/trunk/vector/v.voronoi/main.c 2013-10-11 13:53:35 UTC (rev 57983)
+++ grass/trunk/vector/v.voronoi/main.c 2013-10-13 15:24:04 UTC (rev 57984)
@@ -7,9 +7,10 @@
* Andrea Aime <aaime libero.it>
* Radim Blazek <radim.blazek gmail.com> (GRASS 5.1 v.voronoi)
* Glynn Clements <glynn gclements.plus.com>,
- * Markus Neteler <neteler itc.it>
+ * Markus Neteler <neteler itc.it>,
+ * Markus Metz
* PURPOSE: produce a Voronoi diagram using vector points
- * COPYRIGHT: (C) 1993-2006, 2001 by the GRASS Development Team
+ * COPYRIGHT: (C) 1993-2006, 2001, 2013 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
@@ -51,6 +52,8 @@
double x, y;
} COOR;
+
+
int cmp(void *a, void *b)
{
COOR *ca = (COOR *) a;
@@ -99,22 +102,20 @@
int i;
int **cats, *ncats, nfields, *fields;
struct {
- struct Flag *line, *table, *area;
- } flag;
- struct {
- struct Option *in, *out, *field, *smooth;
+ struct Option *in, *out, *field, *smooth, *thin;
} opt;
+ struct {
+ struct Flag *line, *table, *area, *skeleton;
+ } flag;
struct GModule *module;
struct line_pnts *Points;
struct line_cats *Cats;
- struct bound_box box;
int node, nnodes;
COOR *coor;
int verbose;
int ncoor, acoor;
- int line, nlines, type, ctype, area, nareas;
- int err_boundaries, err_centr_out, err_centr_dupl, err_nocentr;
- double snap_thresh;
+ int line, nlines, type, ctype;
+ double thresh;
G_gisinit(argv[0]);
@@ -137,13 +138,27 @@
opt.smooth->key = "smoothness";
opt.smooth->answer = "0.25";
opt.smooth->label = _("Factor for output smoothness");
- opt.smooth->description = _("Applies to input areas only. Smaller values produce smoother output but can cause numerical instability.");
+ opt.smooth->description = _("Applies to input areas only. "
+ "Smaller values produce smoother output but can cause numerical instability.");
+ opt.thin = G_define_option();
+ opt.thin->type = TYPE_DOUBLE;
+ opt.thin->key = "thin";
+ opt.thin->answer = "-1";
+ opt.thin->label = _("Maximum dangle length of skeletons");
+ opt.thin->description = _("Applies only to skeleton extraction. "
+ "Default = -1 will extract the center line.");
+
flag.area = G_define_flag();
flag.area->key = 'a';
flag.area->description =
_("Create Voronoi diagram for input areas");
+ flag.skeleton = G_define_flag();
+ flag.skeleton->key = 's';
+ flag.skeleton->description =
+ _("Extract skeletons for input areas");
+
flag.line = G_define_flag();
flag.line->key = 'l';
flag.line->description =
@@ -165,6 +180,11 @@
segf = 0.25;
G_warning(_("Option '%s' is too small, set to %g"), opt.smooth->key, segf);
}
+ thresh = atof(opt.thin->answer);
+
+ skeleton = flag.skeleton->answer;
+ if (skeleton)
+ Type = GV_LINE;
Points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
@@ -184,7 +204,7 @@
freeinit(&sfl, sizeof(struct Site));
G_message(_("Reading features..."));
- if (in_area)
+ if (in_area || skeleton)
readbounds();
else
readsites();
@@ -200,72 +220,79 @@
plot = 0;
debug = 0;
- G_message(_("Processing Voronoi triangulation..."));
+ G_message(_("Voronoi triangulation for %d points..."), nsites);
voronoi(nextone);
+ G_message(_("Writing edges..."));
+ vo_write();
- /* Close free ends by current region */
verbose = G_verbose();
G_set_verbose(0);
Vect_build_partial(&Out, GV_BUILD_BASE);
G_set_verbose(verbose);
- ncoor = 0;
- acoor = 100;
- coor = (COOR *) G_malloc(sizeof(COOR) * acoor);
+ if (skeleton) {
+ G_message(_("Thin skeletons ..."));
+ thin_skeleton(thresh);
+ }
+ else {
+ /* Close free ends by current region */
+ ncoor = 0;
+ acoor = 100;
+ coor = (COOR *) G_malloc(sizeof(COOR) * acoor);
- nnodes = Vect_get_num_nodes(&Out);
- for (node = 1; node <= nnodes; node++) {
- double x, y;
+ nnodes = Vect_get_num_nodes(&Out);
+ for (node = 1; node <= nnodes; node++) {
+ double x, y;
- if (Vect_get_node_n_lines(&Out, node) < 2) { /* add coordinates */
- Vect_get_node_coor(&Out, node, &x, &y, NULL);
+ if (Vect_get_node_n_lines(&Out, node) < 2) { /* add coordinates */
+ Vect_get_node_coor(&Out, node, &x, &y, NULL);
- if (ncoor == acoor - 5) { /* always space for 5 region corners */
- acoor += 100;
- coor = (COOR *) G_realloc(coor, sizeof(COOR) * acoor);
- }
+ if (ncoor == acoor - 5) { /* always space for 5 region corners */
+ acoor += 100;
+ coor = (COOR *) G_realloc(coor, sizeof(COOR) * acoor);
+ }
- coor[ncoor].x = x;
- coor[ncoor].y = y;
- ncoor++;
+ coor[ncoor].x = x;
+ coor[ncoor].y = y;
+ ncoor++;
+ }
}
- }
- /* Add region corners */
- coor[ncoor].x = Box.W;
- coor[ncoor].y = Box.S;
- ncoor++;
- coor[ncoor].x = Box.E;
- coor[ncoor].y = Box.S;
- ncoor++;
- coor[ncoor].x = Box.E;
- coor[ncoor].y = Box.N;
- ncoor++;
- coor[ncoor].x = Box.W;
- coor[ncoor].y = Box.N;
- ncoor++;
+ /* Add region corners */
+ coor[ncoor].x = Box.W;
+ coor[ncoor].y = Box.S;
+ ncoor++;
+ coor[ncoor].x = Box.E;
+ coor[ncoor].y = Box.S;
+ ncoor++;
+ coor[ncoor].x = Box.E;
+ coor[ncoor].y = Box.N;
+ ncoor++;
+ coor[ncoor].x = Box.W;
+ coor[ncoor].y = Box.N;
+ ncoor++;
- /* Sort */
- qsort(coor, ncoor, sizeof(COOR), (void *)cmp);
+ /* Sort */
+ qsort(coor, ncoor, sizeof(COOR), (void *)cmp);
- /* add last (first corner) */
- coor[ncoor].x = Box.W;
- coor[ncoor].y = Box.S;
- ncoor++;
+ /* add last (first corner) */
+ coor[ncoor].x = Box.W;
+ coor[ncoor].y = Box.S;
+ ncoor++;
- for (i = 1; i < ncoor; i++) {
- if (coor[i].x == coor[i - 1].x && coor[i].y == coor[i - 1].y)
- continue; /* duplicate */
+ for (i = 1; i < ncoor; i++) {
+ if (coor[i].x == coor[i - 1].x && coor[i].y == coor[i - 1].y)
+ continue; /* duplicate */
- Vect_reset_line(Points);
- Vect_append_point(Points, coor[i].x, coor[i].y, 0.0);
- Vect_append_point(Points, coor[i - 1].x, coor[i - 1].y, 0.0);
- Vect_write_line(&Out, Type, Points, Cats);
+ Vect_reset_line(Points);
+ Vect_append_point(Points, coor[i].x, coor[i].y, 0.0);
+ Vect_append_point(Points, coor[i - 1].x, coor[i - 1].y, 0.0);
+ Vect_write_line(&Out, Type, Points, Cats);
+ }
+
+ G_free(coor);
}
- G_free(coor);
-
- /* Copy input points as centroids */
nfields = Vect_cidx_get_num_fields(&In);
cats = (int **)G_malloc(nfields * sizeof(int *));
ncats = (int *)G_malloc(nfields * sizeof(int));
@@ -278,7 +305,8 @@
fields[i] = Vect_cidx_get_field_number(&In, i);
}
- if (flag.line->answer)
+ /* Copy input points */
+ if (Type == GV_LINE)
ctype = GV_POINT;
else
ctype = GV_CENTROID;
@@ -298,9 +326,9 @@
if (!Vect_point_in_box(Points->x[0], Points->y[0], 0.0, &Box))
continue;
- Vect_write_line(&Out, ctype, Points, Cats);
+ if (!skeleton)
+ Vect_write_line(&Out, ctype, Points, Cats);
-
for (i = 0; i < Cats->n_cats; i++) {
int f, j;
@@ -382,142 +410,17 @@
}
}
-
Vect_close(&In);
-
- /* cleaning part 1: count errors */
- G_message(_("Searching for topology errors..."));
- G_set_verbose(0);
- Vect_build_partial(&Out, GV_BUILD_CENTROIDS);
- G_set_verbose(verbose);
- err_boundaries = err_centr_out = err_centr_dupl = err_nocentr = 0;
- nlines = Vect_get_num_lines(&Out);
- for (line = 1; line <= nlines; line++) {
-
- if (!Vect_line_alive(&Out, line))
- continue;
-
- type = Vect_get_line_type(&Out, line);
- if (type == GV_BOUNDARY) {
- int left, right;
-
- Vect_get_line_areas(&Out, line, &left, &right);
-
- if (left == 0 || right == 0) {
- G_debug(3, "line = %d left = %d right = %d", line,
- left, right);
- err_boundaries++;
- }
- }
- if (type == GV_CENTROID) {
- area = Vect_get_centroid_area(&Out, line);
- if (area == 0)
- err_centr_out++;
- else if (area < 0)
- err_centr_dupl++;
- }
- }
-
- err_nocentr = 0;
- nareas = Vect_get_num_areas(&Out);
- for (area = 1; area <= nareas; area++) {
- if (!Vect_area_alive(&Out, area))
- continue;
- line = Vect_get_area_centroid(&Out, area);
- if (line == 0)
- err_nocentr++;
- }
-
- /* cleaning part 2: snap */
- /* TODO: adjust snapping treshold to ULP */
- Vect_get_map_box(&Out, &box);
- snap_thresh = fabs(box.W);
- if (snap_thresh < fabs(box.E))
- snap_thresh = fabs(box.E);
- if (snap_thresh < fabs(box.N))
- snap_thresh = fabs(box.N);
- if (snap_thresh < fabs(box.S))
- snap_thresh = fabs(box.S);
- snap_thresh = d_ulp(snap_thresh);
- if (err_nocentr || err_centr_dupl || err_centr_out) {
- int nmod;
+ if (Type == GV_BOUNDARY)
+ clean_topo();
- G_important_message(_("Cleaning output topology"));
- Vect_snap_lines(&Out, GV_BOUNDARY, snap_thresh, NULL);
- do {
- Vect_break_lines(&Out, GV_BOUNDARY, NULL);
- Vect_remove_duplicates(&Out, GV_BOUNDARY, NULL);
- nmod =
- Vect_clean_small_angles_at_nodes(&Out, GV_BOUNDARY, NULL);
- } while (nmod > 0);
-
- G_message(_("Removing dangles..."));
- Vect_remove_dangles(&Out, GV_BOUNDARY, -1.0, NULL);
- G_message(_("Removing bridges..."));
- Vect_remove_bridges(&Out, NULL, NULL, NULL);
-
- err_boundaries = 0;
- nlines = Vect_get_num_lines(&Out);
- for (line = 1; line <= nlines; line++) {
-
- if (!Vect_line_alive(&Out, line))
- continue;
-
- type = Vect_get_line_type(&Out, line);
- if (type == GV_BOUNDARY) {
- int left, right;
-
- Vect_get_line_areas(&Out, line, &left, &right);
-
- if (left == 0 || right == 0) {
- G_debug(3, "line = %d left = %d right = %d", line,
- left, right);
- err_boundaries++;
- }
- }
- }
- }
- /* cleaning part 3: remove remaining incorrect boundaries */
- if (err_boundaries) {
- G_important_message(_("Removing incorrect boundaries from output"));
- nlines = Vect_get_num_lines(&Out);
- for (line = 1; line <= nlines; line++) {
-
- if (!Vect_line_alive(&Out, line))
- continue;
-
- type = Vect_get_line_type(&Out, line);
- if (type == GV_BOUNDARY) {
- int left, right;
-
- Vect_get_line_areas(&Out, line, &left, &right);
-
- /* &&, not ||, no typo */
- if (left == 0 && right == 0) {
- G_debug(3, "line = %d left = %d right = %d", line,
- left, right);
- Vect_delete_line(&Out, line);
- }
- }
- }
- }
-
- /* this is slow:
- if (in_area) {
- if (Type == GV_LINE)
- G_message(_("Merging lines ..."));
- else
- G_message(_("Merging boundaries ..."));
- Vect_merge_lines(&Out, Type, NULL, NULL);
- }
- */
-
/* build clean topology */
Vect_build_partial(&Out, GV_BUILD_NONE);
Vect_build(&Out);
Vect_close(&Out);
G_done_msg(" ");
+
exit(EXIT_SUCCESS);
}
Added: grass/trunk/vector/v.voronoi/skeleton.c
===================================================================
--- grass/trunk/vector/v.voronoi/skeleton.c (rev 0)
+++ grass/trunk/vector/v.voronoi/skeleton.c 2013-10-13 15:24:04 UTC (rev 57984)
@@ -0,0 +1,472 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include "defs.h"
+
+static int next_dist(int line, int side, double mf)
+{
+ double d, dist, nextdist, totaldist;
+ int nextline, nlines;
+ int node;
+ static struct line_pnts *Points = NULL;
+
+ G_debug(3, "next_dist()");
+
+ if (!Points)
+ Points = Vect_new_line_struct();
+
+ Vect_read_line(&Out, Points, NULL, abs(line));
+ dist = Vect_line_length(Points);
+
+ if (line < 0)
+ Vect_get_line_nodes(&Out, -line, &node, NULL);
+ else
+ Vect_get_line_nodes(&Out, line, NULL, &node);
+
+ nlines = Vect_get_node_n_lines(&Out, node);
+
+ if (nlines == 1)
+ return 1;
+
+ nextdist = totaldist = 0.;
+ while (nlines > 1) {
+ nextline = dig_angle_next_line(&(Out.plus), -line, side, GV_LINE, NULL);
+ Vect_read_line(&Out, Points, NULL, abs(nextline));
+ d = Vect_line_length(Points);
+ nextdist += d;
+ totaldist += d;
+
+ if (nextline < 0)
+ Vect_get_line_nodes(&Out, -nextline, &node, NULL);
+ else
+ Vect_get_line_nodes(&Out, nextline, NULL, &node);
+
+ nlines = Vect_get_node_n_lines(&Out, node);
+ if (nlines > 2)
+ nextdist = 0.;
+ line = nextline;
+ }
+
+ if (totaldist > nextdist && dist > nextdist)
+ return (totaldist < mf * dist);
+
+ return (dist > nextdist);
+}
+
+static int loop_test(int line, int node, struct line_pnts *Points,
+ double l, double mf)
+{
+ int line1, line2, nextline;
+ int n1, n2, nspikes, nout;
+ double l1, minl, maxl;
+
+ if (Vect_get_node_n_lines(&Out, node) != 3)
+ return 0;
+
+ line1 = line2 = 0;
+ line1 = Vect_get_node_line(&Out, node, 0);
+ if (abs(line1) == abs(line))
+ line1 = 0;
+ line2 = Vect_get_node_line(&Out, node, 1);
+ if (abs(line2) == abs(line))
+ line2 = 0;
+ if (line1 == 0) {
+ line1 = line2;
+ line2 = 0;
+ }
+ if (line2 == 0)
+ line2 = Vect_get_node_line(&Out, node, 2);
+
+ if (abs(line1) == abs(line2))
+ return 1;
+
+ nextline = dig_angle_next_line(&(Out.plus), -line, GV_LEFT, GV_LINE, NULL);
+ line1 = nextline;
+
+ nspikes = 1;
+ nout = 0;
+ minl = 0;
+ maxl = 0;
+ do {
+ if (nextline < 0)
+ Vect_get_line_nodes(&Out, -nextline, &n1, NULL);
+ else
+ Vect_get_line_nodes(&Out, nextline, NULL, &n1);
+
+ if (Vect_get_node_n_lines(&Out, n1) == 1)
+ return 0;
+
+ if (n1 != node && Vect_get_node_n_lines(&Out, n1) > 2) {
+ nspikes++;
+ line2 = dig_angle_next_line(&(Out.plus), -nextline, GV_LEFT, GV_LINE, NULL);
+
+ if (line2 < 0)
+ Vect_get_line_nodes(&Out, -line2, &n2, NULL);
+ else
+ Vect_get_line_nodes(&Out, line2, NULL, &n2);
+
+ if (Vect_get_node_n_lines(&Out, n2) == 1) {
+ Vect_read_line(&Out, Points, NULL, abs(line2));
+ l1 = Vect_line_length(Points);
+ if (minl == 0 || minl > l1) {
+ minl = l1;
+ }
+ if (maxl < l1) {
+ maxl = l1;
+ }
+ }
+ else
+ nout++;
+ }
+
+ nextline = dig_angle_next_line(&(Out.plus), -nextline, GV_RIGHT, GV_LINE, NULL);
+
+ } while (abs(nextline) != abs(line1));
+
+ if (minl == 0)
+ minl = l;
+ if (maxl == 0)
+ maxl = mf * l;
+
+ nspikes -= nout;
+
+ if (mf > 1) {
+ return (nspikes < 3 || (nspikes == 3 && (l > minl || mf * l > maxl)));
+
+ if (l > minl)
+ return 1;
+ if (nspikes == 3 && l < minl)
+ return (mf * l > maxl);
+ else
+ return (nspikes < 3);
+ }
+
+ return (nspikes < 3 || (nspikes > 2 && l > minl));
+}
+
+static int break_loop(int line, int node, struct line_pnts *Points)
+{
+ int line1, line2, firstline, nextline;
+ int n1;
+ double l1, l2;
+
+ if (Vect_get_node_n_lines(&Out, node) != 3)
+ return 0;
+
+ line1 = line2 = 0;
+ line1 = Vect_get_node_line(&Out, node, 0);
+ if (abs(line1) == abs(line))
+ line1 = 0;
+ line2 = Vect_get_node_line(&Out, node, 1);
+ if (abs(line2) == abs(line))
+ line2 = 0;
+ if (line1 == 0) {
+ line1 = line2;
+ line2 = 0;
+ }
+ if (line2 == 0)
+ line2 = Vect_get_node_line(&Out, node, 2);
+
+ if (abs(line1) == abs(line2))
+ return 1;
+
+ nextline = dig_angle_next_line(&(Out.plus), -line, GV_LEFT, GV_LINE, NULL);
+ firstline = nextline;
+
+ do {
+ if (nextline < 0)
+ Vect_get_line_nodes(&Out, -nextline, &n1, NULL);
+ else
+ Vect_get_line_nodes(&Out, nextline, NULL, &n1);
+
+ if (Vect_get_node_n_lines(&Out, n1) == 1)
+ return 0;
+
+ nextline = dig_angle_next_line(&(Out.plus), -nextline, GV_RIGHT, GV_LINE, NULL);
+
+ } while (abs(nextline) != abs(firstline));
+
+ if (abs(nextline) != abs(firstline)) {
+ G_warning("no loop ???");
+ return 0;
+ }
+
+ Vect_read_line(&Out, Points, NULL, abs(line1));
+ l1 = Vect_line_length(Points);
+
+ Vect_read_line(&Out, Points, NULL, abs(line2));
+ l2 = Vect_line_length(Points);
+
+ if (l1 > l2)
+ Vect_delete_line(&Out, abs(line1));
+ else
+ Vect_delete_line(&Out, abs(line2));
+
+ Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
+
+ return 1;
+}
+
+static int length_test(int line, int node, struct line_pnts *Points,
+ double l, double mf)
+{
+ int line1, line2, nlines;
+ int n1;
+ double l1, l2;
+
+ if (Vect_get_node_n_lines(&Out, node) != 3)
+ return 0;
+
+ line1 = Vect_get_node_line(&Out, node, 0);
+ line2 = 0;
+ if (abs(line1) == abs(line))
+ line1 = 0;
+ line2 = Vect_get_node_line(&Out, node, 1);
+ if (abs(line2) == abs(line))
+ line2 = 0;
+ if (line1 == 0) {
+ line1 = line2;
+ line2 = 0;
+ }
+ if (line2 == 0)
+ line2 = Vect_get_node_line(&Out, node, 2);
+
+ if (line1 < 0)
+ Vect_get_line_nodes(&Out, -line1, &n1, NULL);
+ else
+ Vect_get_line_nodes(&Out, line1, NULL, &n1);
+
+ nlines = Vect_get_node_n_lines(&Out, n1);
+ if (nlines > 1)
+ return 0;
+
+ if (line2 < 0)
+ Vect_get_line_nodes(&Out, -line2, &n1, NULL);
+ else
+ Vect_get_line_nodes(&Out, line2, NULL, &n1);
+
+ nlines = Vect_get_node_n_lines(&Out, n1);
+ if (nlines > 1)
+ return 0;
+
+ Vect_read_line(&Out, Points, NULL, abs(line1));
+ l1 = Vect_line_length(Points);
+
+ Vect_read_line(&Out, Points, NULL, abs(line2));
+ l2 = Vect_line_length(Points);
+
+ if (l1 > mf * l2) {
+ if (mf * l < l1 && l < l2)
+ return 0;
+ }
+ if (l2 > mf * l1) {
+ if (mf * l < l2 && l < l1)
+ return 0;
+ }
+
+ return (mf * l > l1 || mf * l > l2);
+}
+
+/* thin the skeletons */
+int thin_skeleton(double thresh)
+{
+ int i;
+ int node, n1, n2;
+ struct line_pnts *Points;
+ struct ilist *list;
+ double l, minl, maxl;
+ int nlines, line, minline, line2;
+ int counter = 1;
+ double morphof = 1.6180339887;
+
+ Points = Vect_new_line_struct();
+ list = Vect_new_list();
+
+ if (thresh < 0)
+ morphof = 1;
+
+ Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
+
+ while (TRUE) {
+ G_verbose_message(_("Pass %d"), counter++);
+ Vect_reset_list(list);
+
+ for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
+ if (!Vect_node_alive(&Out, node))
+ continue;
+ nlines = Vect_get_node_n_lines(&Out, node);
+
+ if (nlines > 1)
+ continue;
+
+ line = Vect_get_node_line(&Out, node, 0);
+ if (line < 0)
+ Vect_get_line_nodes(&Out, -line, &n1, NULL);
+ else
+ Vect_get_line_nodes(&Out, line, NULL, &n1);
+
+ nlines = Vect_get_node_n_lines(&Out, n1);
+
+ if (nlines < 3)
+ continue;
+
+ Vect_read_line(&Out, Points, NULL, abs(line));
+ minline = line;
+ minl = Vect_line_length(Points);
+
+ if (nlines == 3) {
+ if (loop_test(line, n1, Points, minl, morphof))
+ continue;
+ if (length_test(line, n1, Points, minl, morphof))
+ continue;
+ }
+
+ maxl = 0;
+ for (i = 0; i < nlines; i++) {
+ line2 = Vect_get_node_line(&Out, n1, i);
+ if (abs(line2) == abs(minline) || abs(line2) == abs(line))
+ continue;
+ if (line2 < 0)
+ Vect_get_line_nodes(&Out, -line2, &n2, NULL);
+ else
+ Vect_get_line_nodes(&Out, line2, NULL, &n2);
+ if (Vect_get_node_n_lines(&Out, n2) > 1)
+ continue;
+ Vect_read_line(&Out, Points, NULL, abs(line2));
+ l = Vect_line_length(Points);
+ if (minl > l) {
+ minl = l;
+ minline = line2;
+ }
+ if (maxl == 0 || maxl < l) {
+ maxl = l;
+ }
+ }
+ if (thresh < 0) {
+ G_ilist_add(list, minline);
+ }
+ else {
+ if (minl < thresh)
+ G_ilist_add(list, minline);
+ }
+ }
+ if (list->n_values == 0)
+ break;
+ nlines = 0;
+ for (i = 0; i < list->n_values; i++) {
+ line = list->value[i];
+ if (Vect_line_alive(&Out, abs(line))) {
+ if (next_dist(line, GV_RIGHT, morphof))
+ continue;
+ if (next_dist(line, GV_LEFT, morphof))
+ continue;
+ Vect_delete_line(&Out, abs(line));
+ nlines++;
+ }
+ }
+ if (nlines == 0)
+ break;
+ else
+ Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
+ }
+
+ if (thresh >= 0)
+ return 0;
+
+ for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
+ if (!Vect_node_alive(&Out, node))
+ continue;
+ nlines = Vect_get_node_n_lines(&Out, node);
+
+ if (nlines > 1)
+ continue;
+
+ line = Vect_get_node_line(&Out, node, 0);
+ if (line < 0)
+ Vect_get_line_nodes(&Out, -line, &n1, NULL);
+ else
+ Vect_get_line_nodes(&Out, line, NULL, &n1);
+
+ nlines = Vect_get_node_n_lines(&Out, n1);
+
+ if (nlines == 3) {
+ break_loop(line, n1, Points);
+ }
+ }
+
+ while (TRUE) {
+ G_verbose_message(_("Pass %d"), counter++);
+ Vect_reset_list(list);
+
+ for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
+ if (!Vect_node_alive(&Out, node))
+ continue;
+ nlines = Vect_get_node_n_lines(&Out, node);
+
+ if (nlines > 1)
+ continue;
+
+ line = Vect_get_node_line(&Out, node, 0);
+ if (line < 0)
+ Vect_get_line_nodes(&Out, -line, &n1, NULL);
+ else
+ Vect_get_line_nodes(&Out, line, NULL, &n1);
+
+ nlines = Vect_get_node_n_lines(&Out, n1);
+
+ if (nlines < 3)
+ continue;
+
+ Vect_read_line(&Out, Points, NULL, abs(line));
+ minline = line;
+ minl = Vect_line_length(Points);
+
+ if (nlines == 3) {
+ if (break_loop(line, n1, Points))
+ continue;
+ }
+
+ for (i = 0; i < nlines; i++) {
+ line2 = Vect_get_node_line(&Out, n1, i);
+ if (abs(line2) == abs(minline) || abs(line2) == abs(line))
+ continue;
+ if (line2 < 0)
+ Vect_get_line_nodes(&Out, -line2, &n2, NULL);
+ else
+ Vect_get_line_nodes(&Out, line2, NULL, &n2);
+ if (Vect_get_node_n_lines(&Out, n2) > 1)
+ continue;
+ Vect_read_line(&Out, Points, NULL, abs(line2));
+ l = Vect_line_length(Points);
+ if (minl > l) {
+ minl = l;
+ minline = line2;
+ }
+ }
+ if (thresh < 0 || minl < thresh)
+ G_ilist_add(list, minline);
+ }
+ if (list->n_values == 0)
+ break;
+ nlines = 0;
+ for (i = 0; i < list->n_values; i++) {
+ line = list->value[i];
+ if (Vect_line_alive(&Out, abs(line))) {
+ if (next_dist(line, GV_RIGHT, morphof))
+ continue;
+ if (next_dist(line, GV_LEFT, morphof))
+ continue;
+ Vect_delete_line(&Out, abs(line));
+ nlines++;
+ }
+ }
+ if (nlines == 0)
+ break;
+ else
+ Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/vector/v.voronoi/skeleton.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass/trunk/vector/v.voronoi/sw_defs.h
===================================================================
--- grass/trunk/vector/v.voronoi/sw_defs.h 2013-10-11 13:53:35 UTC (rev 57983)
+++ grass/trunk/vector/v.voronoi/sw_defs.h 2013-10-13 15:24:04 UTC (rev 57984)
@@ -66,6 +66,12 @@
extern int PQcount;
extern int PQmin;
+/* clean_topo.c */
+int clean_topo(void);
+
+/* skeleton.c */
+int thin_skeleton(double);
+
/* sw_edgelist.c */
int ELinitialize(void);
struct Halfedge *HEcreate(struct Edge *, int);
@@ -119,5 +125,6 @@
double, double, double *, double *, int);
/* vo_write.c */
+int vo_write(void);
int write_ep(struct Edge *);
Modified: grass/trunk/vector/v.voronoi/sw_main.c
===================================================================
--- grass/trunk/vector/v.voronoi/sw_main.c 2013-10-11 13:53:35 UTC (rev 57983)
+++ grass/trunk/vector/v.voronoi/sw_main.c 2013-10-13 15:24:04 UTC (rev 57984)
@@ -34,6 +34,7 @@
int Type;
int Field;
int in_area;
+int skeleton;
double segf;
int nsites_alloc;
Modified: grass/trunk/vector/v.voronoi/sw_voronoi.c
===================================================================
--- grass/trunk/vector/v.voronoi/sw_voronoi.c 2013-10-11 13:53:35 UTC (rev 57983)
+++ grass/trunk/vector/v.voronoi/sw_voronoi.c 2013-10-13 15:24:04 UTC (rev 57984)
@@ -1,4 +1,5 @@
#include <stdlib.h>
+#include <grass/gis.h>
#include "sw_defs.h"
/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
@@ -14,6 +15,7 @@
int pm;
struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
struct Edge *e;
+ int counter = 0;
PQinitialize();
bottomsite = (*nextsite) ();
@@ -29,6 +31,8 @@
(newsite->coord.y == newintstar.y &&
newsite->coord.x < newintstar.x))) { /* new site is smallest */
+ G_percent(counter++, nsites, 2);
+
lbnd = ELleftbnd(&(newsite->coord));
rbnd = ELright(lbnd);
bot = rightreg(lbnd);
@@ -95,12 +99,7 @@
break;
}
- for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd)) {
- e = lbnd->ELedge;
- write_ep(e);
- }
-
- /* TODO: free memory */
+ G_percent(1, 1, 1);
return 0;
}
Modified: grass/trunk/vector/v.voronoi/v.voronoi.html
===================================================================
--- grass/trunk/vector/v.voronoi/v.voronoi.html 2013-10-11 13:53:35 UTC (rev 57983)
+++ grass/trunk/vector/v.voronoi/v.voronoi.html 2013-10-13 15:24:04 UTC (rev 57984)
@@ -1,61 +1,69 @@
<h2>DESCRIPTION</h2>
-<em>v.voronoi</em> uses an existing vector points map (<b>input</b>) to create
-a Voronoi diagram (Thiessen polygons) in a new vector map (<b>output</b>).
-<p>The bounds of the output map are limited by the current region.
-(see <em>g.region</em>)
+<em>v.voronoi</em> creates a Voronoi diagram (Thiessen polygons) from
+points or centroids.
+<p>The bounds of the output map are limited by the current region
+(see <em><a href="g.region.html">g.region</a></em>).
+
<p>
-<br>
-Voronoi diagram and Delaunay triangulation example:
-<center>
-<img src="v_voronoi_delaunay.png" border="1"><br>
-<table border="0" width="590">
-<tr><td><center>
-<i>Delaunay Triangulation (left pane), Voronoi diagram (center pane),
-and both (right pane)</i>
-</center></td></tr>
-</table>
-</center>
+The <em>-a</em> flag can be used to create a Voronoi diagram for areas.
+<p>
+The <em>-s</em> flag can be used to extract the center line of areas or
+skeletons of areas with <em>thin</em> >= 0. Smaller values for the
+<em>thin</em> option will preserve more detail, while negative values
+will extract only the center line.
<h2>NOTES</h2>
Voronoi diagrams may be used for nearest-neighbor flood filling.
Give the centroids attributes (start with
<em><a href="v.db.addcolumn.html">v.db.addcolumn</a></em>),
-then optionally convert to a raster map with
+then optionally convert the result to a raster map with
<em><a href="v.to.rast.html">v.to.rast</a></em>.
<h2>EXAMPLE</h2>
-Commands used with the Spearfish dataset to create the above figure.
+<h4>Voronoi diagram for points</h4>
+This example uses the hospitals in the North Carolina dataset.
<div class="code"><pre>
- g.region n=4927250 s=4919400 w=588650 e=594850
- d.frame -c fr=one at=0,100,0,33.3333
- d.frame -c fr=two at=0,100,33.3333,66.6667
- d.frame -c fr=three at=0,100,66.6667,100
+ g.region -p rast=elev_state_500m
+ v.voronoi in=hospitals out=hospitals_voronoi
+</pre></div>
- v.delaunay -lr in=archsites out=arch_delaunay
- d.frame -s one
- d.vect arch_delaunay
- d.vect archsites color=red fcolor=red size=5 icon=basic/circle
+Result:
+<center>
+<img src="v_voronoi_points.png" border="0"><br>
+<i>Voronoi diagram for hospitals in North Carolina</i>
+</center>
- v.voronoi -l in=archsites out=arch_voronoi
- d.frame -s two
- d.vect arch_voronoi type=line
- d.vect archsites color=red fcolor=red size=5 icon=basic/circle
-
- d.frame -s three
- d.vect arch_voronoi type=line
- d.vect arch_delaunay color=blue
- d.vect archsites color=red fcolor=red size=5 icon=basic/circle
+<h4>Voronoi diagram for areas</h4>
+This example uses urban areas in the North Carolina dataset.
+<div class="code"><pre>
+ g.region -p n=162500 s=80000 w=727000 e=846000 res=500
+ v.voronoi in=urbanarea out=urbanarea_voronoi -a
</pre></div>
+Result:
+<center>
+<img src="v_voronoi_areas.png" border="0"><br>
+<i>Voronoi diagram for urban areas in North Carolina</i>
+</center>
-<h2>BUGS</h2>
-Only attribute table of field 1 is copied.
+<h4>Skeletons and center lines of areas</h4>
+This example uses urban areas in the North Carolina dataset.
+<div class="code"><pre>
+ g.region -p n=161000 s=135500 w=768500 e=805500 res=500
+ v.voronoi in=urbanarea out=urbanarea_centerline -s
+ v.voronoi in=urbanarea out=urbanarea_skeleton -s thin=2000
+</pre></div>
+Result:
+<center>
+<img src="v_voronoi_skeleton.png" border="0"><br>
+<i>Skeleton (blue) and center line (red) for urban areas in North Carolina</i>
+</center>
<h2>REFERENCES</h2>
<em>Steve J. Fortune, (1987). A Sweepline Algorithm for
@@ -73,6 +81,7 @@
<h2>AUTHORS</h2>
James Darrell McCauley, Purdue University<br>
GRASS 5 update, improvements: <a href="mailto:aaime at libero.it">Andrea Aime</a>, Modena, Italy<br>
-GRASS 5.7 update: Radim Blazek
+GRASS 5.7 update: Radim Blazek<br>
+Markus Metz
<p><i>Last changed: $Date$</i>
Added: grass/trunk/vector/v.voronoi/v_voronoi_areas.png
===================================================================
(Binary files differ)
Property changes on: grass/trunk/vector/v.voronoi/v_voronoi_areas.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Deleted: grass/trunk/vector/v.voronoi/v_voronoi_delaunay.png
===================================================================
(Binary files differ)
Added: grass/trunk/vector/v.voronoi/v_voronoi_points.png
===================================================================
(Binary files differ)
Property changes on: grass/trunk/vector/v.voronoi/v_voronoi_points.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: grass/trunk/vector/v.voronoi/v_voronoi_skeleton.png
===================================================================
(Binary files differ)
Property changes on: grass/trunk/vector/v.voronoi/v_voronoi_skeleton.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Modified: grass/trunk/vector/v.voronoi/vo_write.c
===================================================================
--- grass/trunk/vector/v.voronoi/vo_write.c 2013-10-11 13:53:35 UTC (rev 57983)
+++ grass/trunk/vector/v.voronoi/vo_write.c 2013-10-13 15:24:04 UTC (rev 57984)
@@ -7,6 +7,56 @@
#include "sw_defs.h"
#include "defs.h"
+static int write_skeleton(struct line_pnts *Points)
+{
+ int i, area1, area2;
+ static struct line_pnts *APoints = NULL;
+ static struct line_cats *Cats = NULL;
+
+ if (!APoints) {
+ APoints = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+ }
+
+ if ((area1 = Vect_find_area(&In, Points->x[0], Points->y[0])) == 0)
+ return 0;
+ if ((area2 = Vect_find_area(&In, Points->x[1], Points->y[1])) == 0)
+ return 0;
+ if (area1 != area2)
+ return 0;
+
+ if (!Vect_get_area_centroid(&In, area1))
+ return 0;
+ Vect_get_area_points(&In, area1, APoints);
+ if (Vect_line_check_intersection(Points, APoints, 0))
+ return 0;
+
+ for (i = 0; i < Vect_get_area_num_isles(&In, area1); i++) {
+ Vect_get_isle_points(&In, Vect_get_area_isle(&In, area1, i), APoints);
+ if (Vect_line_check_intersection(Points, APoints, 0))
+ return 0;
+ }
+
+ Vect_get_area_cats(&In, area1, Cats);
+ Vect_write_line(&Out, GV_LINE, Points, Cats);
+
+ return 1;
+}
+
+int vo_write(void)
+{
+ struct Halfedge *lbnd;
+
+ for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd)) {
+ write_ep(lbnd->ELedge);
+ }
+
+ /* TODO: free memory */
+
+ return 1;
+}
+
+
int write_ep(struct Edge *e)
{
static struct line_pnts *Points = NULL;
@@ -39,7 +89,10 @@
Vect_reset_line(Points);
Vect_append_point(Points, x1, y1, 0.0);
Vect_append_point(Points, x2, y2, 0.0);
- Vect_write_line(&Out, Type, Points, Cats);
+ if (skeleton)
+ write_skeleton(Points);
+ else
+ Vect_write_line(&Out, Type, Points, Cats);
}
else {
int knownPointAtLeft = -1;
@@ -80,7 +133,10 @@
Vect_reset_line(Points);
Vect_append_point(Points, x1, y1, 0.0);
Vect_append_point(Points, x2, y2, 0.0);
- Vect_write_line(&Out, Type, Points, Cats);
+ if (skeleton)
+ write_skeleton(Points);
+ else
+ Vect_write_line(&Out, Type, Points, Cats);
}
}
More information about the grass-commit
mailing list