Skip to content

Commit

Permalink
Vlib: Return all intersection points from Vect_line_get_intersections (…
Browse files Browse the repository at this point in the history 
…#2358)

* Vlib: report all line intersections in Vect_line_get_intersections|2()

Current implementation of function Vect_line_get_intersections() reports
a single intersection even if lines cross multiple times
(a buggy implementation introduced by me in 2008 with
https://trac.osgeo.org/grass/changeset/34306)

Added tests should cover most of line intersection cases.

Also fixes bug #2344 (+ a test case)

Contains an update to documentation (from PR #2359 by metzm; fixes #2345)
  • Loading branch information
@marisn @neteler
marisn authored and neteler committed May 20, 2022
1 parent 1bd11ee commit cd19eb5
Show file tree
Hide file tree
Showing 4 changed files with 515 additions and 65 deletions.
72 changes: 47 additions & 25 deletions lib/vector/Vlib/intersect.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@
of limited number of decimal places and for different order of
coordinates, the results would be different)
(C) 2001-2009 by the GRASS Development Team
(C) 2001-2009, 2022 by the GRASS Development Team
This program is free software under the GNU General Public License
(>=v2). Read the file COPYING that comes with GRASS for details.
Expand Down Expand Up @@ -85,6 +85,8 @@ static int ident(double x1, double y1, double x2, double y2, double thresh);
#endif
static int cross_seg(int id, const struct RTree_Rect *rect, void *arg);
static int find_cross(int id, const struct RTree_Rect *rect, void *arg);
int line_check_intersection(struct line_pnts *APoints,
struct line_pnts *BPoints, int with_z);

#define D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2))
#define D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2))
Expand Down Expand Up @@ -635,6 +637,9 @@ static int cross_seg(int id, const struct RTree_Rect *rect, void *arg)
* intersection with B line. Points (Points->n_points == 1) are not
* supported.
*
* Superseded by the faster Vect_line_intersection2()
* Kept as reference implementation
*
* \param APoints first input line
* \param BPoints second input line
* \param[out] ALines array of new lines created from original A line
Expand Down Expand Up @@ -1215,6 +1220,7 @@ Vect_line_intersection(struct line_pnts *APoints,
static struct line_pnts *APnts, *BPnts, *IPnts;

static int cross_found; /* set by find_cross() */
static int report_all; /* should all crossings be reported or just first one */

/* break segments (called by rtree search) */
static int find_cross(int id, const struct RTree_Rect *rect, void *arg)
Expand All @@ -1239,40 +1245,29 @@ static int find_cross(int id, const struct RTree_Rect *rect, void *arg)
case 5:
break;
case 1:
if (0 > Vect_copy_xyz_to_pnts(IPnts, &x1, &y1, &z1, 1))
if (0 > Vect_append_point(IPnts, x1, y1, z1))
G_warning(_("Error while adding point to array. Out of memory"));
break;
case 2:
case 3:
case 4:
if (0 > Vect_copy_xyz_to_pnts(IPnts, &x1, &y1, &z1, 1))
if (0 > Vect_append_point(IPnts, x1, y1, z1))
G_warning(_("Error while adding point to array. Out of memory"));
if (0 > Vect_copy_xyz_to_pnts(IPnts, &x2, &y2, &z2, 1))
if (0 > Vect_append_point(IPnts, x2, y2, z2))
G_warning(_("Error while adding point to array. Out of memory"));
break;
}
/* add ALL (including end points and duplicates), clean later */
if (ret > 0) {
cross_found = 1;
return 0;
if (!report_all)
return 0;
}
return 1; /* keep going */
}

/*!
* \brief Check if 2 lines intersect.
*
* Points (Points->n_points == 1) are also supported.
*
* \param APoints first input line
* \param BPoints second input line
* \param with_z 3D, not supported (only if one or both are points)!
*
* \return 0 no intersection
* \return 1 intersection found
*/
int
Vect_line_check_intersection(struct line_pnts *APoints,
line_check_intersection(struct line_pnts *APoints,
struct line_pnts *BPoints, int with_z)
{
int i;
Expand Down Expand Up @@ -1300,15 +1295,15 @@ Vect_line_check_intersection(struct line_pnts *APoints,
if (APoints->n_points == 1 && BPoints->n_points == 1) {
if (APoints->x[0] == BPoints->x[0] && APoints->y[0] == BPoints->y[0]) {
if (!with_z) {
if (0 >
if (report_all && 0 >
Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0],
&APoints->y[0], NULL, 1))
G_warning(_("Error while adding point to array. Out of memory"));
return 1;
}
else {
if (APoints->z[0] == BPoints->z[0]) {
if (0 >
if (report_all && 0 >
Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0],
&APoints->y[0], &APoints->z[0],
1))
Expand All @@ -1330,7 +1325,7 @@ Vect_line_check_intersection(struct line_pnts *APoints,
NULL, NULL);

if (dist <= rethresh) {
if (0 >
if (report_all && 0 >
Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0], &APoints->y[0],
&APoints->z[0], 1))
G_warning(_("Error while adding point to array. Out of memory"));
Expand All @@ -1347,7 +1342,7 @@ Vect_line_check_intersection(struct line_pnts *APoints,
NULL, NULL);

if (dist <= rethresh) {
if (0 >
if (report_all && 0 >
Vect_copy_xyz_to_pnts(IPnts, &BPoints->x[0], &BPoints->y[0],
&BPoints->z[0], 1))
G_warning(_("Error while adding point to array. Out of memory"));
Expand Down Expand Up @@ -1428,9 +1423,9 @@ Vect_line_check_intersection(struct line_pnts *APoints,
rect.boundary[5] = APoints->z[i];
}

RTreeSearch(MyRTree, &rect, find_cross, &i); /* A segment number from 0 */
RTreeSearch(MyRTree, &rect, find_cross, &i); /* A segment number from 0 */

if (cross_found) {
if (!report_all && cross_found) {
break;
}
}
Expand All @@ -1441,11 +1436,37 @@ Vect_line_check_intersection(struct line_pnts *APoints,
return cross_found;
}

/*!
* \brief Check if 2 lines intersect.
*
* Points (Points->n_points == 1) are also supported.
*
* Superseded by the faster Vect_line_check_intersection2()
* Kept as reference implementation
*
* \param APoints first input line
* \param BPoints second input line
* \param with_z 3D, not supported (only if one or both are points)!
*
* \return 0 no intersection
* \return 1 intersection found
*/
int
Vect_line_check_intersection(struct line_pnts *APoints,
struct line_pnts *BPoints, int with_z)
{
report_all = 0;
return line_check_intersection(APoints, BPoints, with_z);
}

/*!
* \brief Get 2 lines intersection points.
*
* A wrapper around Vect_line_check_intersection() function.
*
* Superseded by the faster Vect_line_get_intersections2()
* Kept as reference implementation
*
* \param APoints first input line
* \param BPoints second input line
* \param[out] IPoints output with intersection points
Expand All @@ -1461,8 +1482,9 @@ Vect_line_get_intersections(struct line_pnts *APoints,
{
int ret;

report_all = 1;
IPnts = IPoints;
ret = Vect_line_check_intersection(APoints, BPoints, with_z);
ret = line_check_intersection(APoints, BPoints, with_z);

return ret;
}
89 changes: 52 additions & 37 deletions lib/vector/Vlib/intersect2.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@
of limited number of decimal places and for different order of
coordinates, the results would be different)
(C) 2001-2014 by the GRASS Development Team
(C) 2001-2014, 2022 by the GRASS Development Team
This program is free software under the GNU General Public License
(>=v2). Read the file COPYING that comes with GRASS for details.
Expand Down Expand Up @@ -89,7 +89,8 @@ 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);

int line_check_intersection2(struct line_pnts *APoints,
struct line_pnts *BPoints, int with_z, int all);

typedef struct
{ /* in arrays 0 - A line , 1 - B line */
Expand Down Expand Up @@ -666,7 +667,9 @@ static int boq_load(struct boq *q, struct line_pnts *Pnts,
* intersection with B line. Points (Points->n_points == 1) are not
* supported. If B line is NULL, A line is intersected with itself.
*
* simplified Bentley–Ottmann Algorithm
* simplified Bentley–Ottmann Algorithm:
* similar to Vect_line_intersection(), but faster
* additionally, self-intersections of a line are handled more efficiently
*
* \param APoints first input line
* \param BPoints second input line or NULL
Expand Down Expand Up @@ -1256,38 +1259,25 @@ static int find_cross(int i, int j, int b)
case 5:
break;
case 1:
if (0 > Vect_copy_xyz_to_pnts(IPnts, &x1, &y1, &z1, 1))
if (0 > Vect_append_point(IPnts, x1, y1, z1))
G_warning(_("Error while adding point to array. Out of memory"));
break;
case 2:
case 3:
case 4:
if (0 > Vect_copy_xyz_to_pnts(IPnts, &x1, &y1, &z1, 1))
if (0 > Vect_append_point(IPnts, x1, y1, z1))
G_warning(_("Error while adding point to array. Out of memory"));
if (0 > Vect_copy_xyz_to_pnts(IPnts, &x2, &y2, &z2, 1))
if (0 > Vect_append_point(IPnts, x2, y2, z2))
G_warning(_("Error while adding point to array. Out of memory"));
break;
}

return ret;
}

/*!
* \brief Check if 2 lines intersect.
*
* Points (Points->n_points == 1) are also supported.
*
* \param APoints first input line
* \param BPoints second input line
* \param with_z 3D, not supported (only if one or both are points)!
*
* \return 0 no intersection
* \return 1 intersection
* \return 2 end points only
*/
int
Vect_line_check_intersection2(struct line_pnts *APoints,
struct line_pnts *BPoints, int with_z)
line_check_intersection2(struct line_pnts *APoints,
struct line_pnts *BPoints, int with_z, int all)
{
double dist;
struct bound_box ABox, BBox, abbox;
Expand All @@ -1314,18 +1304,17 @@ Vect_line_check_intersection2(struct line_pnts *APoints,
if (APoints->n_points == 1 && BPoints->n_points == 1) {
if (APoints->x[0] == BPoints->x[0] && APoints->y[0] == BPoints->y[0]) {
if (!with_z) {
if (0 >
Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0],
&APoints->y[0], NULL, 1))
if (all && 0 >
Vect_append_point(IPnts, APoints->x[0],
APoints->y[0], APoints->z[0]))
G_warning(_("Error while adding point to array. Out of memory"));
return 1;
}
else {
if (APoints->z[0] == BPoints->z[0]) {
if (0 >
Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0],
&APoints->y[0], &APoints->z[0],
1))
if (all && 0 >
Vect_append_point(IPnts, APoints->x[0],
APoints->y[0], APoints->z[0]))
G_warning(_("Error while adding point to array. Out of memory"));
return 1;
}
Expand All @@ -1344,9 +1333,9 @@ Vect_line_check_intersection2(struct line_pnts *APoints,
NULL, NULL);

if (dist <= d_ulp(APoints->x[0], APoints->y[0])) {
if (0 >
Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0], &APoints->y[0],
&APoints->z[0], 1))
if (all && 0 >
Vect_append_point(IPnts, APoints->x[0], APoints->y[0],
APoints->z[0]))
G_warning(_("Error while adding point to array. Out of memory"));
return 1;
}
Expand All @@ -1361,9 +1350,9 @@ Vect_line_check_intersection2(struct line_pnts *APoints,
NULL, NULL);

if (dist <= d_ulp(BPoints->x[0], BPoints->y[0])) {
if (0 >
Vect_copy_xyz_to_pnts(IPnts, &BPoints->x[0], &BPoints->y[0],
&BPoints->z[0], 1))
if (all && 0 >
Vect_append_point(IPnts, BPoints->x[0], BPoints->y[0],
BPoints->z[0]))
G_warning(_("Error while adding point to array. Out of memory"));
return 1;
}
Expand Down Expand Up @@ -1540,7 +1529,7 @@ Vect_line_check_intersection2(struct line_pnts *APoints,
G_fatal_error("RB tree error!");
}
}
if (intersect == 1) {
if (!all && intersect == 1) {
break;
}
}
Expand All @@ -1551,10 +1540,36 @@ Vect_line_check_intersection2(struct line_pnts *APoints,
return intersect;
}

/*!
* \brief Check if 2 lines intersect.
*
* Points (Points->n_points == 1) are also supported.
*
* simplified Bentley–Ottmann Algorithm:
* similar to Vect_line_check_intersection(), but faster
*
* \param APoints first input line
* \param BPoints second input line
* \param with_z 3D, not supported (only if one or both are points)!
*
* \return 0 no intersection
* \return 1 intersection
* \return 2 end points only
*/
int
Vect_line_check_intersection2(struct line_pnts *APoints,
struct line_pnts *BPoints, int with_z)
{
return line_check_intersection2(APoints, BPoints, with_z, 0);
}

/*!
* \brief Get 2 lines intersection points.
*
* A wrapper around Vect_line_check_intersection() function.
* A wrapper around Vect_line_check_intersection2() function.
*
* simplified Bentley–Ottmann Algorithm:
* similar to Vect_line_get_intersections(), but faster
*
* \param APoints first input line
* \param BPoints second input line
Expand All @@ -1572,7 +1587,7 @@ Vect_line_get_intersections2(struct line_pnts *APoints,
int ret;

IPnts = IPoints;
ret = Vect_line_check_intersection2(APoints, BPoints, with_z);
ret = line_check_intersection2(APoints, BPoints, with_z, 1);

return ret;
}
Loading

0 comments on commit cd19eb5

Please sign in to comment.