[GRASS-SVN] r72355 - grass/trunk/lib/vector/Vlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Mar 12 11:46:32 PDT 2018
Author: mmetz
Date: 2018-03-12 11:46:32 -0700 (Mon, 12 Mar 2018)
New Revision: 72355
Modified:
grass/trunk/lib/vector/Vlib/intersect2.c
Log:
Vlib: fix for self-intersections in Vect_line_intersection2(): snap self-intersection only once
Modified: grass/trunk/lib/vector/Vlib/intersect2.c
===================================================================
--- grass/trunk/lib/vector/Vlib/intersect2.c 2018-03-12 11:24:32 UTC (rev 72354)
+++ grass/trunk/lib/vector/Vlib/intersect2.c 2018-03-12 18:46:32 UTC (rev 72355)
@@ -85,6 +85,8 @@
#if 0
static int ident(double x1, double y1, double x2, double y2, double thresh);
#endif
+static int snap_cross(int asegment, double *adistance, int bsegment,
+ double *bdistance, double *xc, double *yc);
static int cross_seg(int i, int j, int b);
static int find_cross(int i, int j, int b);
@@ -199,11 +201,78 @@
/* shared by Vect_line_intersection, Vect_line_check_intersection, cross_seg, find_cross */
static struct line_pnts *APnts, *BPnts, *ABPnts[2], *IPnts;
+/* Snap breaks to nearest vertices within RE threshold */
+/* Calculate distances along segments */
+static int snap_cross(int asegment, double *adistance, int bsegment,
+ double *bdistance, double *xc, double *yc)
+{
+ int seg;
+ double x, y;
+ double dist, curdist;
+
+ /* 1. of A seg */
+ seg = asegment;
+ curdist = dist2(*xc, *yc, APnts->x[seg], APnts->y[seg]);
+ x = APnts->x[seg];
+ y = APnts->y[seg];
+
+ *adistance = curdist;
+
+ /* 2. of A seg */
+ dist = dist2(*xc, *yc, APnts->x[seg + 1], APnts->y[seg + 1]);
+ if (dist < curdist) {
+ curdist = dist;
+ x = APnts->x[seg + 1];
+ y = APnts->y[seg + 1];
+ }
+
+ /* 1. of B seg */
+ seg = bsegment;
+ dist = dist2(*xc, *yc, BPnts->x[seg], BPnts->y[seg]);
+ *bdistance = dist;
+
+ if (dist < curdist) {
+ curdist = dist;
+ x = BPnts->x[seg];
+ y = BPnts->y[seg];
+ }
+ /* 2. of B seg */
+ dist = dist2(*xc, *yc, BPnts->x[seg + 1], BPnts->y[seg + 1]);
+ if (dist < curdist) {
+ curdist = dist;
+ x = BPnts->x[seg + 1];
+ y = BPnts->y[seg + 1];
+ }
+
+ /* the threshold should not be too small, otherwise we get
+ * too many tiny new segments
+ * the threshold should not be too large, otherwise we might
+ * introduce new crossings
+ * the smallest difference representable with
+ * single precision floating point works well with pathological input
+ * regular input is not affected */
+ if (curdist < d_ulp(x, y)) { /* was rethresh * rethresh */
+ *xc = x;
+ *yc = y;
+
+ /* Update distances along segments */
+ seg = asegment;
+ *adistance = dist2(*xc, *yc, APnts->x[seg], APnts->y[seg]);
+ seg = bsegment;
+ *bdistance = dist2(*xc, *yc, BPnts->x[seg], BPnts->y[seg]);
+
+ return 1;
+ }
+
+ return 0;
+}
+
/* break segments */
static int cross_seg(int i, int j, int b)
{
double x1, y1, z1, x2, y2, z2;
double y1min, y1max, y2min, y2max;
+ double adist, bdist;
int ret;
y1min = APnts->y[i];
@@ -251,9 +320,11 @@
G_debug(2, " -> %d x %d: intersection type = %d", i, j, ret);
if (ret == 1) { /* one intersection on segment A */
G_debug(3, " in %f, %f ", x1, y1);
- add_cross(i, 0.0, j, 0.0, x1, y1);
+ /* snap intersection only once */
+ snap_cross(i, &adist, j, &bdist, &x1, &y1);
+ add_cross(i, adist, j, bdist, x1, y1);
if (APnts == BPnts)
- add_cross(j, 0.0, i, 0.0, x1, y1);
+ add_cross(j, bdist, i, adist, x1, y1);
}
else if (ret == 2 || ret == 3 || ret == 4 || ret == 5) {
/* partial overlap; a broken in one, b broken in one
@@ -261,12 +332,14 @@
* or b contains a; b is broken in 2 points (but 1 may be end)
* or identical */
G_debug(3, " in %f, %f; %f, %f", x1, y1, x2, y2);
- add_cross(i, 0.0, j, 0.0, x1, y1);
- add_cross(i, 0.0, j, 0.0, x2, y2);
- if (APnts == BPnts) {
- add_cross(j, 0.0, i, 0.0, x1, y1);
- add_cross(j, 0.0, i, 0.0, x2, y2);
- }
+ snap_cross(i, &adist, j, &bdist, &x1, &y1);
+ add_cross(i, adist, j, bdist, x1, y1);
+ if (APnts == BPnts)
+ add_cross(j, bdist, i, adist, x1, y1);
+ snap_cross(i, &adist, j, &bdist, &x2, &y2);
+ add_cross(i, adist, j, bdist, x2, y2);
+ if (APnts == BPnts)
+ add_cross(j, bdist, i, adist, x2, y2);
}
}
return 1; /* keep going */
@@ -616,8 +689,7 @@
{
int i, j, k, l, nl, last_seg, seg, last;
int n_alive_cross;
- double dist, curdist, last_x, last_y, last_z;
- double x, y;
+ double dist, last_x, last_y, last_z;
struct line_pnts **XLines, *Points;
struct line_pnts *Points1, *Points2; /* first, second points */
int seg1, seg2, vert1, vert2;
@@ -817,69 +889,6 @@
return 0;
}
- /* Snap breaks to nearest vertices within RE threshold */
- /* Calculate distances along segments */
- for (i = 0; i < n_cross; i++) {
-
- /* 1. of A seg */
- seg = cross[i].segment[0];
- curdist =
- dist2(cross[i].x, cross[i].y, APnts->x[seg], APnts->y[seg]);
- x = APnts->x[seg];
- y = APnts->y[seg];
-
- cross[i].distance[0] = curdist;
-
- /* 2. of A seg */
- dist =
- dist2(cross[i].x, cross[i].y, APnts->x[seg + 1],
- APnts->y[seg + 1]);
- if (dist < curdist) {
- curdist = dist;
- x = APnts->x[seg + 1];
- y = APnts->y[seg + 1];
- }
-
- /* 1. of B seg */
- seg = cross[i].segment[1];
- dist =
- dist2(cross[i].x, cross[i].y, BPnts->x[seg], BPnts->y[seg]);
- cross[i].distance[1] = dist;
-
- if (dist < curdist) {
- curdist = dist;
- x = BPnts->x[seg];
- y = BPnts->y[seg];
- }
- /* 2. of B seg */
- dist = dist2(cross[i].x, cross[i].y, BPnts->x[seg + 1], BPnts->y[seg + 1]);
- if (dist < curdist) {
- curdist = dist;
- x = BPnts->x[seg + 1];
- y = BPnts->y[seg + 1];
- }
-
- /* the threshold should not be too small, otherwise we get
- * too many tiny new segments
- * the threshold should not be too large, otherwise we might
- * introduce new crossings
- * the smallest difference representable with
- * single precision floating point works well with pathological input
- * regular input is not affected */
- if (curdist < d_ulp(x, y)) { /* was rethresh * rethresh */
- cross[i].x = x;
- cross[i].y = y;
-
- /* Update distances along segments */
- seg = cross[i].segment[0];
- cross[i].distance[0] =
- dist2(APnts->x[seg], APnts->y[seg], cross[i].x, cross[i].y);
- seg = cross[i].segment[1];
- cross[i].distance[1] =
- dist2(BPnts->x[seg], BPnts->y[seg], cross[i].x, cross[i].y);
- }
- }
-
/* l = 1 ~ line A, l = 2 ~ line B */
nl = 3;
if (same)
More information about the grass-commit
mailing list