[GRASS-SVN] r55796 - grass/trunk/lib/vector/Vlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Apr 14 15:01:12 PDT 2013
Author: mmetz
Date: 2013-04-14 15:01:12 -0700 (Sun, 14 Apr 2013)
New Revision: 55796
Modified:
grass/trunk/lib/vector/Vlib/intersect.c
Log:
Vlib: fix Vect_segment_intersection()
Modified: grass/trunk/lib/vector/Vlib/intersect.c
===================================================================
--- grass/trunk/lib/vector/Vlib/intersect.c 2013-04-14 12:10:59 UTC (rev 55795)
+++ grass/trunk/lib/vector/Vlib/intersect.c 2013-04-14 22:01:12 UTC (rev 55796)
@@ -23,8 +23,8 @@
1 - intersect at one point
<pre>
\ / \ / \ /
- \/ \/ \/
- /\ \
+ \/ \/ \/
+ /\ \
/ \ \
2 - partial overlap ( \/ )
------ a ( distance < threshold )
@@ -113,8 +113,8 @@
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 @@
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 @@
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 @@
* 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 @@
/* 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 @@
*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 @@
*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 @@
*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 @@
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 @@
/* 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 @@
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 @@
*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 @@
*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 @@
*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 @@
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 @@
/* 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 @@
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 @@
* 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
More information about the grass-commit
mailing list