Compute intersection of two edges

Question

I have this big graph of edges and vertices in 2D space. The big graph is returned by a function computed in a C++ library. I am reading this graph and using it to compute all the intersections of its edges (the lines segements). I use sweep algorithm. For detecting the intersection of two edges I have though some problems. I have 4 different methods according to which I test if two edges intersect and if affirmative I compute and retain their intersection:

1. One which test if the two edges are the diagonals of a polygon. that is the coordinates of the edges of one line inserted into the equation of the other line have different signs.

2. One which computes the intersection each time and check whether the found intersection is between the endpoints of both segments.

3. One which is the code from this link implemented in C++ though.

4. One which implements the first method proposed by Jason Cohen ín this question.

"The problem reduces to this question: Do two lines from A to B and from C to D intersect? Then you can ask it four times (between the line and each of the four sides of the rectangle).

Here's the vector math for doing it. I'm assuming the line from A to B is the line in question and the line from C to D is one of the rectangle lines. My notation is that Ax is the "x-coordinate of A" and Cy is the "y-coordinate of C." And "*" means dot-product, so e.g.:

A*B = Ax*Bx + Ay*By.
E = B-A = ( Bx-Ax, By-Ay )
F = D-C = ( Dx-Cx, Dy-Cy ) 
P = ( -Ey, Ex )
h = ( (A-C) * P ) / ( F * P )

This h number is the key. If h is between 0 and 1, the lines intersect, otherwise they don't. If F*P is zero, of course you cannot make the calculation, but in this case the lines are parallel and therefore only intersect in the obvious cases. The exact point of intersection is C + F*h. If h is exactly 0 or 1 the lines touch at an end-point. You can consider this an "intersection" or not as you see fit."

For data that I created (small data with double values) I obtained good results with all the 4 implemented methods. When I use anyone of these methods implemented in C++ on the data from the big graph I get different results each time and not good results:

1. method returns much more intersections that I need (all the points are on the graph) but I get too many points.

2. I always obtain 0 intersections no matter what.

3. I get a lot more intersection than in 1.

4. For some example I get points which are not on the graph (so not even the intersection). But for some examples I get the intersection points plus some other points.

I have no idea where the problem can be. Any suggestion or any other idea on how to compute the intersection and test it moreover are ore than welcomed. thank you, madalina

Answer 1

What you need are unit tests for your 4 methods and to test them thoroughly. Especially with line-segment intersections there are lots of end-cases, such as parallel slopes, coincident end-points, fully or partially overlapping, in addition to all the usual tolerancing problems (e.g. what exactly does "equal slope" mean?).

Without breaking things down into smaller, more testable units you're going to have a tough time narrowing it down.

Answer 2

Although it's hard to say without being able to see your code, I'd suspect you are running into robustness issues with your floating point values. Have you considered using integer points or doing some kind of robustness enhancement like eliminating common bits in the floating point values?

There is an open source geometry library called GEOS (http://trac.osgeo.org/geos/) which might be useful to test your data set out on. There are also algorithms in GEOS to perform snap rounding to an integer grid and eliminate common bits to help you determine if you are encountering a robustness issue. Of further note is how GEOS computes the intersection using point-line duality in homogenous space (which I can't quite tell if the dot product projection technique you describe is mathmatically equivalant to).

As an aside, my favorite solution to compute intersections in a graph of edges is to use a sweepline in conjunction with a monotone chain of edges. When you use monotone chains, you can eliminate a lot of edge intersection tests, since chains never intersect themselves. This is what GEOS does.

Answer 3

I am using a sweep line on monotone edges when I compute the intersection (I have a sort function which sort the edges inside the constructor and I test them they I well-ordered), even though I obtain as intersection some points which are vertices for some edges, even though I have a lot of tests to eliminate this cases. This is the code for the 4th method (which gives me the best results so far, but not all intersections, but at least some good intersection as compared with the others):

//4 method
double* QSweep::findIntersection(edge_t edge1,edge_t edge2) {  
point_t p1=myPoints_[edge1[0]];
point_t p2=myPoints_[edge1[1]];
point_t p3=myPoints_[edge2[0]];
point_t p4=myPoints_[edge2[1]];

double xD1,yD1,xD2,yD2,xD3,yD3,xP,yP,h,denom;
double* pt=new double[3];

// calculate differences  
xD1=p2[0]-p1[0];  
xD2=p4[0]-p3[0];  
yD1=p2[1]-p1[1];  
yD2=p4[1]-p3[1];  
xD3=p1[0]-p3[0];  
yD3=p1[1]-p3[1];    

xP=-yD1;
yP=xD1;
denom=xD2*(-yD1)+yD2*xD1;
if (denom==0) {
 return NULL;
}
else{
h=(xD3*(-yD1)+yD3*xD1)/denom;
}
std::cout<<"h is"<<h<<endl;
if ((h>0)&&(h<1)){
 pt[0]=p3[0]+xD2*h;  
 pt[1]=p3[1]+yD2*h;
 pt[2]=0.00;
}
else{
 return NULL;
}

// return the valid intersection  
return pt;

}

void QSweep:: intersection(){
nbIntersections=0;
double* vector;
for (int i=2;i<myEdges_.size()-1;i++){
 vector=findIntersection(myEdges_[i-1],myEdges_[i]);
 if (vector!=NULL){
      nbIntersections++;
   myIntersections.push_back(vector);
   swap(myEdges_[i-1], myEdges_[i]);
 }
}
}

The intersection function is always the same I just find the findIntersection function. For the 1 and 2 method I am using the following version of the findIntersection function:

double* QSweep::computeIntersection(double m1, double b1, double m2, double b2){
double* v=new double[3];

v[0]= (b2-b1)/(m1-m2);
v[1]= (m1*b2-m2*b1)/(m1-m2);
v[2]=0.00;
return v;
}

double* QSweep::findIntersection(edge_t edge1, edge_t edge2){


//equation for the support line of the first edge

double a=myPoints_[edge1[0]][0];
double b=myPoints_[edge1[0]][1];
double c=myPoints_[edge1[1]][0];
double d=myPoints_[edge1[1]][1];
double m1=getSlope(a,b,c,d);
double b1=getYIntercept(a,b,c,d);

double x=myPoints_[edge2[0]][0];
double y=myPoints_[edge2[0]][1];
double s=myPoints_[edge2[1]][0];
double t=myPoints_[edge2[1]][1];
double m2=getSlope(x,y,s,t);
double b2=getYIntercept(x,y,s,t);
double* vector;

double line2InLine1=evalEqnLineD(myPoints_[edge2[0]],edge1);
double line2InLine1New=evalEqnLineD(myPoints_[edge2[1]],edge1);
double line1InLine2=evalEqnLineD(myPoints_[edge1[0]],edge2);
double line1InLine2New=evalEqnLineD(myPoints_[edge1[1]],edge2);

if (m1==m2){
 return NULL;
}
else{

     if ((line1InLine2*line1InLine2New)<0 &&   (line2InLine1*line2InLine1New)<0){
  vector=computeIntersection(m1,b1,m2,b2);

  if ((vector[0]>a && vector[0]<c)&&(vector[0]>x && vector[0]<s)){
   return vector;
  }
  else{
   return NULL;
  }
 }
 else{
  return NULL;
 }
}
return vector;

 }

I will start again from the beginning to see what the intersection points that I want have in common. Even tought with some tests I dont even get the good intersection points but other points which are on the graph so I am really confused.

thank you for your suggestions, madalina

Answer 4

Why re-invent the wheel?

If you can handle relying on Boost, I'd check out the GTL library in the sandbox. For instructions on how to download this, read the wiki (instructions right there on the front page).

The GTL library lends itself to any implementation you have for your data structure (i.e. it's designed in the spirit of STL, where data structures and algorithms are orthogonal). The code is both fast and correct. Check it out if you can.

Answer 5

You may want to try Boost.Geometry (aka Generic Geometry Library).