Vlib: fix for self-intersections in Vect_line_intersection2(): snap self-intersection only once

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@72355 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/vector/Vlib/intersect2.c

@@ -85,6 +85,8 @@ static int debug_level = -1;
 #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 @@ static int ident(double x1, double y1, double x2, double y2, double thresh)
 /* 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 @@ static int cross_seg(int i, int j, int b)
 	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 @@ static int cross_seg(int i, int j, int b)
 	     *  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 @@ Vect_line_intersection2(struct line_pnts *APoints,
 {
     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 @@ Vect_line_intersection2(struct line_pnts *APoints,
 	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)