[GRASS-SVN] r57820 - grass/branches/releasebranch_6_4/vector/v.voronoi
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Sep 23 07:16:58 PDT 2013
Author: mmetz
Date: 2013-09-23 07:16:56 -0700 (Mon, 23 Sep 2013)
New Revision: 57820
Modified:
grass/branches/releasebranch_6_4/vector/v.voronoi/sw_defs.h
grass/branches/releasebranch_6_4/vector/v.voronoi/sw_geometry.c
grass/branches/releasebranch_6_4/vector/v.voronoi/sw_main.c
grass/branches/releasebranch_6_4/vector/v.voronoi/sw_output.c
grass/branches/releasebranch_6_4/vector/v.voronoi/vo_extend.c
grass/branches/releasebranch_6_4/vector/v.voronoi/vo_main.c
grass/branches/releasebranch_6_4/vector/v.voronoi/vo_write.c
Log:
v.voronoi: backport fixes for #957, #1682, #2019
Modified: grass/branches/releasebranch_6_4/vector/v.voronoi/sw_defs.h
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.voronoi/sw_defs.h 2013-09-23 14:10:27 UTC (rev 57819)
+++ grass/branches/releasebranch_6_4/vector/v.voronoi/sw_defs.h 2013-09-23 14:16:56 UTC (rev 57820)
@@ -110,6 +110,7 @@
int makevertex(struct Site *);
int deref(struct Site *);
int ref(struct Site *);
+double d_ulp(double);
/* sw_heap.c */
int PQinsert(struct Halfedge *, struct Site *, double);
Modified: grass/branches/releasebranch_6_4/vector/v.voronoi/sw_geometry.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.voronoi/sw_geometry.c 2013-09-23 14:10:27 UTC (rev 57819)
+++ grass/branches/releasebranch_6_4/vector/v.voronoi/sw_geometry.c 2013-09-23 14:16:56 UTC (rev 57820)
@@ -1,5 +1,6 @@
#include <stdlib.h>
#include <math.h>
+#include <grass/gis.h>
#include "sw_defs.h"
int geominit(void)
@@ -32,12 +33,22 @@
newedge->ep[0] = (struct Site *)NULL;
newedge->ep[1] = (struct Site *)NULL;
- dx = s2->coord.x - s1->coord.x;
- dy = s2->coord.y - s1->coord.y;
+ if (s1->coord.x < s2->coord.x ||
+ (s1->coord.x == s2->coord.x && s1->coord.y < s2->coord.y)) {
+ dx = s2->coord.x - s1->coord.x;
+ dy = s2->coord.y - s1->coord.y;
+ newedge->c =
+ s1->coord.x * dx + s1->coord.y * dy + (dx * dx + dy * dy) * 0.5;
+ }
+ else {
+ dx = s1->coord.x - s2->coord.x;
+ dy = s1->coord.y - s2->coord.y;
+ newedge->c =
+ s2->coord.x * dx + s2->coord.y * dy + (dx * dx + dy * dy) * 0.5;
+ }
+
adx = dx > 0 ? dx : -dx;
ady = dy > 0 ? dy : -dy;
- newedge->c =
- s1->coord.x * dx + s1->coord.y * dy + (dx * dx + dy * dy) * 0.5;
if (adx > ady) {
newedge->a = 1.0;
newedge->b = dy / dx;
@@ -55,12 +66,29 @@
return (newedge);
}
+/* single precision ULP */
+double d_ulp(double d)
+{
+ int exp;
+ if (d == 0)
+ return GRASS_EPSILON;
+
+ if (d < 0)
+ d = fabs(d);
+
+ d = frexp(d, &exp);
+ exp -= 22;
+ d = ldexp(d, exp);
+
+ return d;
+}
+
struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2)
{
struct Edge *e1, *e2, *e;
struct Halfedge *el;
- double d, xint, yint;
+ double d, dt, xint, yint;
int right_of_site;
struct Site *v;
@@ -72,9 +100,21 @@
return ((struct Site *)NULL);
d = e1->a * e2->b - e1->b * e2->a;
- if (-1.0e-10 < d && d < 1.0e-10)
+ if (fabs(e1->a * e2->b) > fabs(e1->b * e2->a)) {
+ dt = fabs(e1->a * e2->b);
+ }
+ else
+ dt = fabs(e1->b * e2->a);
+
+ if (dt != dt)
return ((struct Site *)NULL);
+ dt = d_ulp(dt);
+ G_debug(4, "dt = %g", dt);
+
+ if (-dt < d && d < dt)
+ return ((struct Site *)NULL);
+
xint = (e1->c * e2->b - e2->c * e1->b) / d;
yint = (e2->c * e1->a - e1->c * e2->a) / d;
Modified: grass/branches/releasebranch_6_4/vector/v.voronoi/sw_main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.voronoi/sw_main.c 2013-09-23 14:10:27 UTC (rev 57819)
+++ grass/branches/releasebranch_6_4/vector/v.voronoi/sw_main.c 2013-09-23 14:16:56 UTC (rev 57820)
@@ -3,6 +3,7 @@
#include <unistd.h>
#include <grass/gis.h>
#include <grass/Vect.h>
+#include <grass/glocale.h>
#include "sw_defs.h"
#include "defs.h"
@@ -41,10 +42,8 @@
void removeDuplicates()
{
int i, j;
- int foundDupe;
i = j = 1;
- foundDupe = 0;
while (i < nsites)
if (mode3d) {
if (sites[i].coord.x == sites[i - 1].coord.x &&
@@ -80,23 +79,27 @@
/* read all sites, sort, and compute xmin, xmax, ymin, ymax */
int readsites(void)
{
- int nlines, line;
+ int nlines, ltype;
struct line_pnts *Points;
-
+ struct line_cats *Cats;
+
Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ nlines = Vect_get_num_primitives(&In, GV_POINTS);
- nlines = Vect_get_num_lines(&In);
-
nsites = 0;
sites = (struct Site *)G_malloc(nlines * sizeof(struct Site));
- for (line = 1; line <= nlines; line++) {
- int type;
+ while(TRUE) {
+ ltype = Vect_read_next_line(&In, Points, Cats);
+ if(ltype == -2)
+ break;
- type = Vect_read_line(&In, Points, NULL, line);
- if (!(type & GV_POINTS))
+ 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;
@@ -131,7 +134,14 @@
nsites++;
}
- if (nsites < nlines - 1)
+ if (nsites < 2) {
+ const char *name = Vect_get_full_name(&In);
+ Vect_close(&In);
+ G_fatal_error(_("Found %d points/centroids in <%s>, but at least 2 are needed"),
+ nsites, name);
+ }
+
+ if (nsites < nlines)
sites =
(struct Site *)G_realloc(sites,
(nsites) * sizeof(struct Site));
@@ -139,6 +149,9 @@
qsort(sites, nsites, sizeof(struct Site), scomp);
removeDuplicates();
+ Vect_destroy_line_struct(Points);
+ Vect_destroy_cats_struct(Cats);
+
return 0;
}
Modified: grass/branches/releasebranch_6_4/vector/v.voronoi/sw_output.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.voronoi/sw_output.c 2013-09-23 14:10:27 UTC (rev 57819)
+++ grass/branches/releasebranch_6_4/vector/v.voronoi/sw_output.c 2013-09-23 14:16:56 UTC (rev 57820)
@@ -47,17 +47,9 @@
int out_ep(struct Edge *e)
{
- static struct line_pnts *Points = NULL;
- static struct line_cats *Cats = NULL;
-
- if (!Points) {
- Points = Vect_new_line_struct();
- Cats = Vect_new_cats_struct();
- }
-
- if (!triangulate & plot)
+ if ((!triangulate) & plot)
clip_line(e);
- if (!triangulate & !plot) {
+ if ((!triangulate) & !plot) {
/*
fprintf (stdout,"e %d", e->edgenbr);
fprintf (stdout," %d ", e->ep[le] != (struct Site *)NULL ? e->ep[le]->sitenbr : -1);
@@ -85,7 +77,7 @@
int out_site(struct Site *s)
{
- if (!triangulate & plot & !debug)
+ if ((!triangulate) & plot & (!debug))
circle(s->coord.x, s->coord.y, cradius);
/*
Modified: grass/branches/releasebranch_6_4/vector/v.voronoi/vo_extend.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.voronoi/vo_extend.c 2013-09-23 14:10:27 UTC (rev 57819)
+++ grass/branches/releasebranch_6_4/vector/v.voronoi/vo_extend.c 2013-09-23 14:16:56 UTC (rev 57820)
@@ -33,7 +33,7 @@
{
double nx, ny; /* intersection coordinates */
- if (x > w && x < e && y > s && y < n) {
+ if (x >= w && x <= e && y >= s && y <= n) {
/* vertical line? */
if (a == 0) {
*c_x = knownPointAtLeft ? e : w;
Modified: grass/branches/releasebranch_6_4/vector/v.voronoi/vo_main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.voronoi/vo_main.c 2013-09-23 14:10:27 UTC (rev 57819)
+++ grass/branches/releasebranch_6_4/vector/v.voronoi/vo_main.c 2013-09-23 14:16:56 UTC (rev 57820)
@@ -108,11 +108,14 @@
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;
G_gisinit(argv[0]);
@@ -190,7 +193,10 @@
voronoi(triangulate, nextone);
/* 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;
@@ -336,7 +342,7 @@
if (fields[i] == 0)
continue;
- G_message(_("Layer %d"), fields[i]);
+ G_debug(1, "Layer %d", fields[i]);
/* Make a list of categories */
IFi = Vect_get_field(&In, fields[i]);
@@ -369,7 +375,9 @@
Vect_close(&In);
/* cleaning part 1: count 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++) {
@@ -409,11 +417,22 @@
}
/* 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(_("Output needs topological cleaning"));
- Vect_snap_lines(&Out, GV_BOUNDARY, 1e-7, NULL);
+ Vect_snap_lines(&Out, GV_BOUNDARY, snap_thresh, NULL);
do {
Vect_break_lines(&Out, GV_BOUNDARY, NULL);
Vect_remove_duplicates(&Out, GV_BOUNDARY, NULL);
Modified: grass/branches/releasebranch_6_4/vector/v.voronoi/vo_write.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.voronoi/vo_write.c 2013-09-23 14:10:27 UTC (rev 57819)
+++ grass/branches/releasebranch_6_4/vector/v.voronoi/vo_write.c 2013-09-23 14:16:56 UTC (rev 57820)
@@ -42,23 +42,43 @@
Vect_write_line(&Out, Type, Points, Cats);
}
else {
- int knownPointAtLeft;
+ int knownPointAtLeft = -1;
if (e->ep[le] != NULL) {
x1 = e->ep[le]->coord.x;
y1 = e->ep[le]->coord.y;
knownPointAtLeft = 1;
}
- else {
+ else if (e->ep[re] != NULL) {
x1 = e->ep[re]->coord.x;
y1 = e->ep[re]->coord.y;
knownPointAtLeft = 0;
}
+ if (knownPointAtLeft == -1) {
+
+ /*
+ G_warning("Dead edge!");
+ return 1;
+ */
+ x2 = (e->reg[le]->coord.x + e->reg[re]->coord.x) / 2.0;
+ y2 = (e->reg[le]->coord.y + e->reg[re]->coord.y) / 2.0;
+ knownPointAtLeft = 0;
+ if (!extend_line(Box.S, Box.N, Box.W, Box.E,
+ e->a, e->b, e->c, x2, y2, &x1, &y1,
+ knownPointAtLeft)) {
+
+ G_warning("Undefined edge, unable to extend line");
+
+ return 0;
+ }
+ G_debug(0, "x1 = %g, y1 = %g, x2 = %g, y2 = %g", x1, y1, x2, y2);
+ knownPointAtLeft = 1;
+ }
if (extend_line(Box.S, Box.N, Box.W, Box.E,
e->a, e->b, e->c, x1, y1, &x2, &y2,
- knownPointAtLeft)
- ) {
+ knownPointAtLeft)) {
+
/* Don't write zero length */
if (x1 == x2 && y1 == y2)
return 0;
More information about the grass-commit
mailing list