Vlib: fix Vect_segment_intersection()

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

lib/vector/Vlib/intersect.c

@@ -23,8 +23,8 @@
    1 - intersect at one point
    <pre>
    \  /    \  /  \  /
-   \/      \/    \/
-   /\             \
+    \/      \/    \/
+    /\             \
    /  \             \
    2 - partial overlap         ( \/                      )
    ------      a          (    distance < threshold )
@@ -113,8 +113,8 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 			      double *y2, double *z2, int with_z)
 {
     static int first_3d = 1;
-    double d, d1, d2, r1, r2, dtol, t;
-    int switched = 0;
+    double d, d1, d2, r1, dtol, t;
+    int switched;
 
     /* TODO: Works for points ? */
 
@@ -128,7 +128,48 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	first_3d = 0;
     }
 
-    /* Check identical lines */
+    /*  'Sort' lines by x, y */
+    switched = 0;
+    if (bx2 < bx1)
+	switched = 1;
+    else if (bx2 == bx1) {
+	if (by2 < by1)
+	    switched = 1;
+    }
+
+    if (switched) {
+	t = bx1;
+	bx1 = bx2;
+	bx2 = t;
+	t = by1;
+	by1 = by2;
+	by2 = t;
+	t = bz1;
+	bz1 = bz2;
+	bz2 = t;
+    }
+
+    switched = 0;
+    if (ax2 < ax1)
+	switched = 1;
+    else if (ax2 == ax1) {
+	if (ay2 < ay1)
+	    switched = 1;
+    }
+
+    if (switched) {
+	t = ax1;
+	ax1 = ax2;
+	ax2 = t;
+	t = ay1;
+	ay1 = ay2;
+	ay2 = t;
+	t = az1;
+	az1 = az2;
+	az2 = t;
+    }
+
+    /* Check identical segments */
     if ((ax1 == bx1 && ay1 == by1 && ax2 == bx2 && ay2 == by2) ||
 	(ax1 == bx2 && ay1 == by2 && ax2 == bx1 && ay2 == by1)) {
 	G_debug(2, " -> identical segments");
@@ -141,40 +182,19 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	return 5;
     }
 
-    /*  'Sort' lines by x1, x2, y1, y2 */
-    if (bx1 < ax1)
-	switched = 1;
-    else if (bx1 == ax1) {
-	if (bx2 < ax2)
-	    switched = 1;
-	else if (bx2 == ax2) {
-	    if (by1 < ay1)
-		switched = 1;
-	    else if (by1 == ay1) {
-		if (by2 < ay2)
-		    switched = 1;	/* by2 != ay2 (would be identical */
-	    }
-	}
+    /* Check distinct (non-touching) segments */
+    if ((bx1 > ax1 && bx2 > ax1 && bx1 > ax2 && bx2 > ax2) ||
+	(bx1 < ax1 && bx2 < ax1 && bx1 < ax2 && bx2 < ax2)) {
+	/* b is to the left or right of a */
+	G_debug(2, "   -> no intersection");
+	return 0;
     }
-    if (switched) {
-	t = ax1;
-	ax1 = bx1;
-	bx1 = t;
-	t = ay1;
-	ay1 = by1;
-	by1 = t;
-	t = az1;
-	az1 = bz1;
-	bz1 = t;
-	t = ax2;
-	ax2 = bx2;
-	bx2 = t;
-	t = ay2;
-	ay2 = by2;
-	by2 = t;
-	t = az2;
-	az2 = bz2;
-	bz2 = t;
+
+    if ((by1 > ay1 && by2 > ay1 && by1 > ay2 && by2 > ay2) ||
+	(by1 < ay1 && by2 < ay1 && by1 < ay2 && by2 < ay2)) {
+	/* b is above or below a */
+	G_debug(2, "   -> no intersection");
+	return 0;
     }
 
     d = D;
@@ -188,16 +208,16 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
      *       Can it be a problem to set the tolerance to 0.0 ? */
     dtol = 0.0;
     if (fabs(d) > dtol) {
-	r1 = d1 / d;
-	r2 = d2 / d;
 
-	G_debug(2, " -> not parallel/collinear: r1 = %f, r2 = %f", r1, r2);
+	G_debug(2, " -> not parallel/collinear: d1 = %f, d2 = %f", d1, d2);
 
-	if (r1 < 0 || r1 > 1 || r2 < 0 || r2 > 1) {
+	if (d1 < 0 || d1 > d || d2 < 0 || d2 > d) {
 	    G_debug(2, "  -> no intersection");
 	    return 0;
 	}
 
+	r1 = d1 / d;
+
 	*x1 = ax1 + r1 * (ax2 - ax1);
 	*y1 = ay1 + r1 * (ay2 - ay1);
 	*z1 = 0;
@@ -219,18 +239,14 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
     /* Collinear vertical */
     /* original code assumed lines were not both vertical
      *  so there is a special case if they are */
-    if (ax1 == ax2 && bx1 == bx2 && ax1 == bx1) {
+    if (ax1 == ax2) {
+	if (bx1 != bx2)
+	    G_fatal_error("Vect_segment_intersection(): bx1 != bx2");
+	if (ax1 != bx1)
+	    G_fatal_error("Vect_segment_intersection(): ax1 != bx1");
+
 	G_debug(2, "  -> collinear vertical");
-	if (ay1 > ay2) {
-	    t = ay1;
-	    ay1 = ay2;
-	    ay2 = t;
-	}			/* to be sure that ay1 < ay2 */
-	if (by1 > by2) {
-	    t = by1;
-	    by1 = by2;
-	    by2 = t;
-	}			/* to be sure that by1 < by2 */
+	/* TODO: remove, already checked above */
 	if (ay1 > by2 || ay2 < by1) {
 	    G_debug(2, "   -> no intersection");
 	    return 0;
@@ -263,10 +279,8 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	    *x2 = bx2;
 	    *y2 = by2;
 	    *z2 = 0;
-	    if (!switched)
-		return 3;
 
-	    return 4;
+	    return 3;
 	}
 	/* b contains a */
 	if (ay1 >= by1 && ay2 <= by2) {
@@ -277,16 +291,14 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	    *x2 = ax2;
 	    *y2 = ay2;
 	    *z2 = 0;
-	    if (!switched)
-		return 4;
 
-	    return 3;
+	    return 4;
 	}
 
 	/* general overlap, 2 intersection points */
 	G_debug(2, "    -> partial overlap");
 	if (by1 > ay1 && by1 < ay2) {	/* b1 in a */
-	    if (!switched) {
+	    if (by2 >= ay2) {
 		*x1 = bx1;
 		*y1 = by1;
 		*z1 = 0;
@@ -294,9 +306,9 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 		*y2 = ay2;
 		*z2 = 0;
 	    }
-	    else {
-		*x1 = ax2;
-		*y1 = ay2;
+	    else {	/* by2 <= ay1 */
+		*x1 = ax1;
+		*y1 = ay1;
 		*z1 = 0;
 		*x2 = bx1;
 		*y2 = by1;
@@ -305,15 +317,15 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	    return 2;
 	}
 	if (by2 > ay1 && by2 < ay2) {	/* b2 in a */
-	    if (!switched) {
+	    if (by1 >= ay2) {
 		*x1 = bx2;
 		*y1 = by2;
 		*z1 = 0;
-		*x2 = ax1;
-		*y2 = ay1;
+		*x2 = ax2;
+		*y2 = ay2;
 		*z2 = 0;
 	    }
-	    else {
+	    else {	/* by1 <= ay1 */
 		*x1 = ax1;
 		*y1 = ay1;
 		*z1 = 0;
@@ -326,36 +338,32 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 
 	/* should not be reached */
 	G_warning(_("Vect_segment_intersection() ERROR (collinear vertical segments)"));
+	G_warning("a");
 	G_warning("%.15g %.15g", ax1, ay1);
 	G_warning("%.15g %.15g", ax2, ay2);
-	G_warning("x");
+	G_warning("b");
 	G_warning("%.15g %.15g", bx1, by1);
 	G_warning("%.15g %.15g", bx2, by2);
 
 	return 0;
     }
 
-    G_debug(2, "   -> collinear non vertical");
-
     /* Collinear non vertical */
-    if ((bx1 > ax1 && bx2 > ax1 && bx1 > ax2 && bx2 > ax2) ||
-	(bx1 < ax1 && bx2 < ax1 && bx1 < ax2 && bx2 < ax2)) {
-	G_debug(2, "   -> no intersection");
-	return 0;
-    }
+
+    G_debug(2, "   -> collinear non vertical");
 
     /* there is overlap or connected end points */
     G_debug(2, "   -> overlap/connected end points");
 
     /* end points */
-    if ((ax1 == bx1 && ay1 == by1) || (ax1 == bx2 && ay1 == by2)) {
+    if (ax1 == bx2 && ay1 == by2) {
 	*x1 = ax1;
 	*y1 = ay1;
 	*z1 = 0;
 	G_debug(2, "    -> connected by end points");
 	return 1;
     }
-    if ((ax2 == bx1 && ay2 == by1) || (ax2 == bx2 && ay2 == by2)) {
+    if (ax2 == bx1 && ay2 == by1) {
 	*x1 = ax2;
 	*y1 = ay2;
 	*z1 = 0;
@@ -363,23 +371,6 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	return 1;
     }
 
-    if (ax1 > ax2) {
-	t = ax1;
-	ax1 = ax2;
-	ax2 = t;
-	t = ay1;
-	ay1 = ay2;
-	ay2 = t;
-    }				/* to be sure that ax1 < ax2 */
-    if (bx1 > bx2) {
-	t = bx1;
-	bx1 = bx2;
-	bx2 = t;
-	t = by1;
-	by1 = by2;
-	by2 = t;
-    }				/* to be sure that bx1 < bx2 */
-
     /* a contains b */
     if (ax1 <= bx1 && ax2 >= bx2) {
 	G_debug(2, "    -> a contains b");
@@ -389,10 +380,8 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	*x2 = bx2;
 	*y2 = by2;
 	*z2 = 0;
-	if (!switched)
-	    return 3;
 
-	return 4;
+	return 3;
     }
     /* b contains a */
     if (ax1 >= bx1 && ax2 <= bx2) {
@@ -403,16 +392,14 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	*x2 = ax2;
 	*y2 = ay2;
 	*z2 = 0;
-	if (!switched)
-	    return 4;
 
-	return 3;
+	return 4;
     }
 
     /* general overlap, 2 intersection points (lines are not vertical) */
     G_debug(2, "    -> partial overlap");
     if (bx1 > ax1 && bx1 < ax2) {	/* b1 is in a */
-	if (!switched) {
+	if (bx2 >= ax2) {
 	    *x1 = bx1;
 	    *y1 = by1;
 	    *z1 = 0;
@@ -420,9 +407,9 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	    *y2 = ay2;
 	    *z2 = 0;
 	}
-	else {
-	    *x1 = ax2;
-	    *y1 = ay2;
+	else {		/* bx2 <= ax1 */
+	    *x1 = ax1;
+	    *y1 = ay1;
 	    *z1 = 0;
 	    *x2 = bx1;
 	    *y2 = by1;
@@ -431,15 +418,15 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 	return 2;
     }
     if (bx2 > ax1 && bx2 < ax2) {	/* b2 is in a */
-	if (!switched) {
+	if (bx1 >= ax2) {
 	    *x1 = bx2;
 	    *y1 = by2;
 	    *z1 = 0;
-	    *x2 = ax1;
-	    *y2 = ay1;
+	    *x2 = ax2;
+	    *y2 = ay2;
 	    *z2 = 0;
 	}
-	else {
+	else {		/* bx1 <= ax1 */
 	    *x1 = ax1;
 	    *y1 = ay1;
 	    *z1 = 0;
@@ -452,9 +439,10 @@ int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2,
 
     /* should not be reached */
     G_warning(_("Vect_segment_intersection() ERROR (collinear non vertical segments)"));
+    G_warning("a");
     G_warning("%.15g %.15g", ax1, ay1);
     G_warning("%.15g %.15g", ax2, ay2);
-    G_warning("x");
+    G_warning("b");
     G_warning("%.15g %.15g", bx1, by1);
     G_warning("%.15g %.15g", bx2, by2);
 
@@ -643,17 +631,6 @@ Vect_line_intersection(struct line_pnts *APoints,
     APnts = APoints;
     BPnts = BPoints;
 
-    /* BUG
-     *  sometimes aline is broken by bline in aline1,2
-     *  then aline1 is broken by aline2 in aline3,4, this should not happen
-     *  then aline3 or aline4 is broken by the same bline in aline5,6
-     *  then aline5 is broken by aline6 in aline7,8, this should not happen
-     *  then aline7 or aline8 is broken by the same bline in aline9,10
-     *  etc ad infinitum
-     *  this can be avoided by not breaking aline fragments again with bline
-     *  followed by Vect_clean_small_angles_at_nodes()
-     */
-
     /* RE (representation error).
      *  RE thresh above is nonsense of course, the RE threshold should be based on
      *  number of significant digits for double (IEEE-754) which is 15 or 16 and exponent. 
@@ -662,7 +639,7 @@ Vect_line_intersection(struct line_pnts *APoints,
      *  and we want precision in mm (+ 3) = 14 -> minimum rethresh may be around 0.001
      *  ?Maybe all nonsense? 
      *  Use rounding error of the unit in the least place ? 
-     *  max of x, y
+     *  max of fabs(x), fabs(y)
      *  rethresh = pow(2, log2(max) - 53) */
 
     /* Warning: This function is also used to intersect the line by itself i.e. APoints and