[GRASS-SVN] r57621 - grass/trunk/vector/v.voronoi
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Sep 8 10:36:42 PDT 2013
Author: mmetz
Date: 2013-09-08 10:36:42 -0700 (Sun, 08 Sep 2013)
New Revision: 57621
Modified:
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/vo_write.c
Log:
v.voronoi: add tesselation of areas
Modified: grass/trunk/vector/v.voronoi/defs.h
===================================================================
--- grass/trunk/vector/v.voronoi/defs.h 2013-09-08 16:52:44 UTC (rev 57620)
+++ grass/trunk/vector/v.voronoi/defs.h 2013-09-08 17:36:42 UTC (rev 57621)
@@ -3,5 +3,6 @@
extern struct bound_box Box;
extern struct Map_info In, Out;
extern int Type;
-extern int All;
extern int Field;
+extern int in_area;
+extern double segf;
Modified: grass/trunk/vector/v.voronoi/main.c
===================================================================
--- grass/trunk/vector/v.voronoi/main.c 2013-09-08 16:52:44 UTC (rev 57620)
+++ grass/trunk/vector/v.voronoi/main.c 2013-09-08 17:36:42 UTC (rev 57621)
@@ -99,10 +99,10 @@
int i;
int **cats, *ncats, nfields, *fields;
struct {
- struct Flag *line, *table;
+ struct Flag *line, *table, *area;
} flag;
struct {
- struct Option *in, *out, *field;
+ struct Option *in, *out, *field, *segf;
} opt;
struct GModule *module;
struct line_pnts *Points;
@@ -132,6 +132,18 @@
opt.out = G_define_standard_option(G_OPT_V_OUTPUT);
+ opt.segf = G_define_option();
+ opt.segf->type = TYPE_DOUBLE;
+ opt.segf->key = "segment";
+ opt.segf->answer = "0.25";
+ opt.segf->label = _("Factor to control boundary interpolation");
+ opt.segf->description = _("Applies to input areas only. Smaller values produce smoother output but can cause numerical instability.");
+
+ flag.area = G_define_flag();
+ flag.area->key = 'a';
+ flag.area->description =
+ _("Create Voronoi diagram for input areas");
+
flag.line = G_define_flag();
flag.line->key = 'l';
flag.line->description =
@@ -147,7 +159,12 @@
else
Type = GV_BOUNDARY;
- All = 0;
+ in_area = flag.area->answer;
+ segf = atof(opt.segf->answer);
+ if (segf < GRASS_EPSILON) {
+ segf = 0.25;
+ G_warning(_("Option '%s' is too small, set to %g"), opt.segf->key, segf);
+ }
Points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
@@ -167,7 +184,10 @@
freeinit(&sfl, sizeof(struct Site));
G_message(_("Reading features..."));
- readsites();
+ if (in_area)
+ readbounds();
+ else
+ readsites();
Vect_open_new(&Out, opt.out->answer, 0);
@@ -483,6 +503,16 @@
}
}
+ /* 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);
Modified: grass/trunk/vector/v.voronoi/sw_defs.h
===================================================================
--- grass/trunk/vector/v.voronoi/sw_defs.h 2013-09-08 16:52:44 UTC (rev 57620)
+++ grass/trunk/vector/v.voronoi/sw_defs.h 2013-09-08 17:36:42 UTC (rev 57621)
@@ -103,7 +103,7 @@
int scomp(const void *, const void *);
struct Site *nextone(void);
int readsites(void);
-struct Site *readone(void);
+int readbounds(void);
/* sw_memory.c */
int freeinit(struct Freelist *, int);
Modified: grass/trunk/vector/v.voronoi/sw_main.c
===================================================================
--- grass/trunk/vector/v.voronoi/sw_main.c 2013-09-08 16:52:44 UTC (rev 57620)
+++ grass/trunk/vector/v.voronoi/sw_main.c 2013-09-08 17:36:42 UTC (rev 57621)
@@ -1,6 +1,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
+#include <math.h>
#include <grass/gis.h>
#include <grass/vector.h>
#include <grass/glocale.h>
@@ -31,9 +32,12 @@
struct bound_box Box;
struct Map_info In, Out;
int Type;
-int All;
int Field;
+int in_area;
+double segf;
+int nsites_alloc;
+
/* sort sites on y, then x, coord */
int scomp(const void *v1, const void *v2)
{
@@ -103,12 +107,48 @@
}
+int addsite(double x, double y, double z, int id)
+{
+ if (nsites >= nsites_alloc) {
+ nsites_alloc += 100;
+ sites =
+ (struct Site *)G_realloc(sites,
+ (nsites_alloc) * sizeof(struct Site));
+ }
+ sites[nsites].coord.x = x;
+ sites[nsites].coord.y = y;
+ sites[nsites].coord.z = z;
+
+ sites[nsites].sitenbr = id;
+ sites[nsites].refcnt = 0;
+
+ if (nsites > 0) {
+ if (xmin > sites[nsites].coord.x)
+ xmin = sites[nsites].coord.x;
+ if (xmax < sites[nsites].coord.x)
+ xmax = sites[nsites].coord.x;
+ if (ymin > sites[nsites].coord.y)
+ ymin = sites[nsites].coord.y;
+ if (ymax < sites[nsites].coord.y)
+ ymax = sites[nsites].coord.y;
+ }
+ else {
+ xmin = xmax = sites[nsites].coord.x;
+ ymin = ymax = sites[nsites].coord.y;
+ }
+
+ nsites++;
+
+ return nsites;
+}
+
/* read all sites, sort, and compute xmin, xmax, ymin, ymax */
int readsites(void)
{
int nlines, ltype;
struct line_pnts *Points;
struct line_cats *Cats;
+ double z = 0.;
Points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
@@ -129,40 +169,16 @@
if (!(ltype & GV_POINTS))
continue;
/* G_percent(Vect_get_next_line_id(&In), nlines, 2); */
-
- if (!All) {
- if (!Vect_point_in_box(Points->x[0], Points->y[0], 0.0, &Box))
- continue;
- }
- sites[nsites].coord.x = Points->x[0];
- sites[nsites].coord.y = Points->y[0];
+ if (!Vect_point_in_box(Points->x[0], Points->y[0], 0.0, &Box))
+ continue;
+
if (mode3d) {
G_debug(3, "Points->z[0]: %f", Points->z[0]);
- sites[nsites].coord.z = Points->z[0];
+ z = Points->z[0];
}
- else
- sites[nsites].coord.z = 0.0;
- sites[nsites].sitenbr = nsites;
- sites[nsites].refcnt = 0;
-
- if (nsites > 1) {
- if (xmin > sites[nsites].coord.x)
- xmin = sites[nsites].coord.x;
- if (xmax < sites[nsites].coord.x)
- xmax = sites[nsites].coord.x;
- if (ymin > sites[nsites].coord.y)
- ymin = sites[nsites].coord.y;
- if (ymax < sites[nsites].coord.y)
- ymax = sites[nsites].coord.y;
- }
- else {
- xmin = xmax = sites[nsites].coord.x;
- ymin = ymax = sites[nsites].coord.y;
- }
-
- nsites++;
+ addsite(Points->x[0], Points->y[0], z, nsites);
}
if (nsites < 2) {
@@ -186,18 +202,278 @@
return 0;
}
-/* read one site */
-struct Site *readone(void)
+/* valid boundary: one area with centroid */
+int n_areas(int line, int *aid)
{
- struct Site *s;
+ int larea, rarea, ncentroids;
+
+ ncentroids = 0;
+ Vect_get_line_areas(&In, line, &larea, &rarea);
+
+ if (larea < 0)
+ larea = Vect_get_isle_area(&In, -larea);
+ if (larea > 0) {
+ if (Vect_get_area_centroid(&In, larea) > 0) {
+ ncentroids++;
+ *aid = larea;
+ }
+ }
+ if (rarea < 0)
+ rarea = Vect_get_isle_area(&In, -rarea);
+ if (rarea > 0) {
+ if (Vect_get_area_centroid(&In, rarea) > 0) {
+ ncentroids++;
+ *aid = rarea;
+ }
+ }
+
+ return ncentroids;
+}
- s = (struct Site *)getfree(&sfl);
- s->refcnt = 0;
- s->sitenbr = siteidx;
- siteidx++;
+/* read all boundaries, sort, and compute xmin, xmax, ymin, ymax */
+int readbounds(void)
+{
+ int line, nlines, ltype, node, nnodes;
+ struct line_pnts *Points;
+ struct line_cats *Cats;
+ int i;
+ int area_id;
+ double x, y, z, x1, y1, z1;
+ double maxdist, sdist, l, dx, dy, dz;
+ int nconnected;
+ struct ilist *linelist, *arealist;
+
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ nlines = Vect_get_num_lines(&In);
- if (scanf("%lf %lf", &(s->coord.x), &(s->coord.y)) == EOF)
- return ((struct Site *)NULL);
+ nsites = 0;
+ nsites_alloc = nlines * 2;
+ sites = (struct Site *)G_malloc(nsites_alloc * sizeof(struct Site));
- return (s);
+ Vect_set_constraint_type(&In, GV_BOUNDARY);
+ Vect_set_constraint_field(&In, Field);
+
+ l = 0;
+ maxdist = 0;
+ for (line = 1; line <= nlines; line++) {
+
+ if (!Vect_line_alive(&In, line))
+ continue;
+ ltype = Vect_get_line_type(&In, line);
+
+ if (!(ltype & GV_BOUNDARY))
+ continue;
+
+ area_id = 0;
+ if (n_areas(line, &area_id) != 1)
+ continue;
+
+ Vect_read_line(&In, Points, Cats, line);
+ Vect_line_prune(Points);
+
+ l += Vect_line_length(Points);
+ nsites += Points->n_points;
+ }
+ if (nsites)
+ maxdist = segf * l / nsites;
+ G_verbose_message("Maximum segment length: %g", maxdist);
+
+ nsites = 0;
+ z = 0.;
+ z1 = 0;
+ dz = 0;
+ for (line = 1; line <= nlines; line++) {
+
+ if (!Vect_line_alive(&In, line))
+ continue;
+ ltype = Vect_get_line_type(&In, line);
+
+ if (!(ltype & GV_BOUNDARY))
+ continue;
+
+ area_id = 0;
+ if (n_areas(line, &area_id) != 1)
+ continue;
+
+ Vect_read_line(&In, Points, Cats, line);
+ Vect_line_prune(Points);
+
+ if (nsites + Points->n_points > nsites_alloc) {
+ nsites_alloc = nsites + Points->n_points;
+ sites =
+ (struct Site *)G_realloc(sites,
+ (nsites_alloc) * sizeof(struct Site));
+ }
+
+ for (i = 0; i < Points->n_points; i++) {
+ if (!Vect_point_in_box(Points->x[i], Points->y[i], 0.0, &Box))
+ continue;
+
+ x = Points->x[i];
+ y = Points->y[i];
+ if (mode3d) {
+ G_debug(3, "Points->z[i]: %f", Points->z[i]);
+ z = Points->z[i];
+ }
+
+ if (i > 0 && i < Points->n_points - 1)
+ addsite(x, y, z, area_id);
+
+ /* densify */
+ if (maxdist > 0 && i < Points->n_points - 1) {
+ dx = Points->x[i + 1] - Points->x[i];
+ dy = Points->y[i + 1] - Points->y[i];
+ if (mode3d)
+ dz = Points->z[i + 1] - Points->z[i];
+ l = sqrt(dx * dx + dy * dy);
+
+ if (l > maxdist) {
+ int n = ceil(l / maxdist) + 0.5;
+ double step = l / n;
+
+ while (--n) {
+ sdist = (step * n) / l;
+ x1 = x + sdist * dx;
+ y1 = y + sdist * dy;
+ if (mode3d)
+ z1 = z + sdist * dz;
+
+ addsite(x1, y1, z1, area_id);
+ }
+ }
+ }
+ }
+ }
+
+ /* process nodes */
+ nnodes = Vect_get_num_nodes(&In);
+ linelist = Vect_new_list();
+ arealist = Vect_new_list();
+
+ for (node = 1; node <= nnodes; node++) {
+ Vect_get_node_coor(&In, node, &x, &y, &z);
+ if (!mode3d)
+ z = 0.;
+
+ if (!Vect_point_in_box(x, y, 0.0, &Box))
+ continue;
+
+ /* count number of connected boundaries
+ * must be >= 2 */
+ /* count number of valid boundaries */
+ nlines = Vect_get_node_n_lines(&In, node);
+ nconnected = 0;
+ Vect_reset_list(linelist);
+ Vect_reset_list(arealist);
+
+ for (i = 0; i < nlines; i++) {
+ line = Vect_get_node_line(&In, node, i);
+
+ ltype = Vect_get_line_type(&In, abs(line));
+
+ if (!(ltype & GV_BOUNDARY))
+ continue;
+
+ if (n_areas(abs(line), &area_id) == 1) {
+ Vect_list_append(linelist, line);
+ Vect_list_append(arealist, area_id);
+ nconnected++;
+ }
+ }
+
+ if (arealist->n_values == 1) {
+
+ area_id = arealist->value[0];
+ addsite(x, y, z, area_id);
+ }
+ else if (arealist->n_values > 1) {
+ /* displacement */
+ double displace;
+
+ if (fabs(x) > fabs(y))
+ displace = d_ulp(fabs(x));
+ else
+ displace = d_ulp(fabs(y));
+
+ displace *= 2;
+
+ for (i = 0; i < linelist->n_values; i++) {
+
+ line = linelist->value[i];
+
+ if (n_areas(abs(line), &area_id) != 1)
+ G_fatal_error(_("All boundaries in the list should be valid"));
+
+ Vect_read_line(&In, Points, Cats, abs(line));
+ Vect_line_prune(Points);
+
+ if (Points->n_points < 2)
+ G_fatal_error(_("Boundary is degenerate"));
+
+ if (line < 0) {
+ dx = Points->x[Points->n_points - 2] - Points->x[Points->n_points - 1];
+ dy = Points->y[Points->n_points - 2] - Points->y[Points->n_points - 1];
+ }
+ else {
+ dx = Points->x[1] - Points->x[0];
+ dy = Points->y[1] - Points->y[0];
+ }
+ l = sqrt(dx * dx + dy * dy);
+ if (displace > l)
+ displace = l;
+ }
+
+ for (i = 0; i < linelist->n_values; i++) {
+
+ line = linelist->value[i];
+
+ if (n_areas(abs(line), &area_id) != 1)
+ G_fatal_error(_("All boundaries in the list should be valid"));
+
+ Vect_read_line(&In, Points, Cats, abs(line));
+ Vect_line_prune(Points);
+
+ if (Points->n_points < 2)
+ G_fatal_error(_("Boundary is degenerate"));
+
+ if (line < 0) {
+ dx = Points->x[Points->n_points - 2] - Points->x[Points->n_points - 1];
+ dy = Points->y[Points->n_points - 2] - Points->y[Points->n_points - 1];
+ }
+ else {
+ dx = Points->x[1] - Points->x[0];
+ dy = Points->y[1] - Points->y[0];
+ }
+ l = sqrt(dx * dx + dy * dy);
+ if (l > displace * 2) {
+ sdist = displace / l;
+ x1 = x + sdist * dx;
+ y1 = y + sdist * dy;
+ z1 = 0;
+
+ addsite(x1, y1, z1, area_id);
+ }
+ }
+ }
+ }
+
+ if (nsites < 2) {
+ const char *name = Vect_get_full_name(&In);
+ Vect_close(&In);
+ G_fatal_error(_("Found %d vertices in <%s>, but at least 2 are needed"),
+ nsites, name);
+ }
+
+
+ qsort(sites, nsites, sizeof(struct Site), scomp);
+ removeDuplicates();
+
+ Vect_destroy_line_struct(Points);
+ Vect_destroy_cats_struct(Cats);
+ Vect_destroy_list(linelist);
+ Vect_destroy_list(arealist);
+
+ return 0;
}
Modified: grass/trunk/vector/v.voronoi/sw_voronoi.c
===================================================================
--- grass/trunk/vector/v.voronoi/sw_voronoi.c 2013-09-08 16:52:44 UTC (rev 57620)
+++ grass/trunk/vector/v.voronoi/sw_voronoi.c 2013-09-08 17:36:42 UTC (rev 57621)
@@ -25,8 +25,11 @@
if (!PQempty())
newintstar = PQ_min();
- if (newsite != (struct Site *)NULL && (PQempty()
- || newsite->coord.y < newintstar.y || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x))) { /* new site is smallest */
+ if (newsite != (struct Site *)NULL &&
+ (PQempty() || newsite->coord.y < newintstar.y ||
+ (newsite->coord.y == newintstar.y &&
+ newsite->coord.x < newintstar.x))) { /* new site is smallest */
+
out_site(newsite);
lbnd = ELleftbnd(&(newsite->coord));
rbnd = ELright(lbnd);
@@ -54,9 +57,8 @@
temp->coord.y == newsite->coord.y);
newsite = temp;
}
- else if (!PQempty())
+ else if (!PQempty()) {
/* intersection is smallest */
- {
lbnd = PQextractmin();
llbnd = ELleft(lbnd);
rbnd = ELright(lbnd);
Modified: grass/trunk/vector/v.voronoi/vo_write.c
===================================================================
--- grass/trunk/vector/v.voronoi/vo_write.c 2013-09-08 16:52:44 UTC (rev 57620)
+++ grass/trunk/vector/v.voronoi/vo_write.c 2013-09-08 17:36:42 UTC (rev 57621)
@@ -20,6 +20,9 @@
if (!triangulate) {
double x1, y1, x2, y2;
+
+ if (in_area && e->reg[le]->sitenbr == e->reg[re]->sitenbr)
+ return 0;
if (e->ep[le] != NULL && e->ep[re] != NULL) { /* both end defined */
x1 = e->ep[le]->coord.x;
More information about the grass-commit
mailing list