2D line intersection

Code for 2D line intersection

I’ve been struggling with finding the intersection between two lines so that the resulting point actually lies on the specified line segments. Most line intersection formulas use an infinite line model, so the point will tend to lie on an infinite extension of the segments. Which is not what I wanted at all.

First I implemented Paul Bourke’s excellent breakdown of the math. Then after much experimentation I realized I could calculate the distances from the intersection point to the end points of the two lines, and then check their combined length against the length of the lines. If the lenghts are the same then the intersection point is actually on the line segment. Quod erat demonstrandum.

The code also checks to see if the two lines are parallel using the dot product of the two vectors. Considering that I’m actually fairly inept at math I feel like I’ve scored a minor victory.

Code - lineIntersect.pde

// lineIntersect.pde
// Marius Watz - http://workshop.evolutionzone.com

// calculates valid intersection between two lines,
// so that the intersection will lie on the specified
// line segment.

Point p[];
int num;

void setup() {
  size(500,250);

  num=12;
  p=new Point[num*2];
  for(int i=0; i< num*2; i++) p[i]=new Point(0,0);
  initPt();
}

void draw() {
  background(200);

  p[num-1].set(mouseX,mouseY);

  // check intersections with all lines
  for(int j=0; j< num; j++) {
    line(p[j*2].x,p[j*2].y, p[j*2+1].x,p[j*2+1].y);
    for(int i=0; i< num; i++) if(i!=j) {
      Point pt=findIntersection(
        p[i*2],p[i*2+1], p[j*2],p[j*2+1]);
      if(pt!=null) ellipse(pt.x,pt.y, 14,14);
    }
  }

}

void initPt() {
  for(int i=0; i< num; i++) {
    if(random(100)>50) {
      p[i*2].set(20,random(20,height-20));
      p[i*2+1].set(width-20,random(20,height-20));
    }
    else {
      p[i*2].set(random(20,width-20),20);
      p[i*2+1].set(random(20,width-20),height-20);
    }
  }
}

void mousePressed() {
  initPt();
  p[num-2].set(mouseX,mouseY);
}

// calculates intersection and checks for parallel lines.
// also checks that the intersection point is actually on
// the line segment p1-p2
Point findIntersection(Point p1,Point p2,
  Point p3,Point p4) {
  float xD1,yD1,xD2,yD2,xD3,yD3;
  float dot,deg,len1,len2;
  float segmentLen1,segmentLen2;
  float ua,ub,div;

  // calculate differences
  xD1=p2.x-p1.x;
  xD2=p4.x-p3.x;
  yD1=p2.y-p1.y;
  yD2=p4.y-p3.y;
  xD3=p1.x-p3.x;
  yD3=p1.y-p3.y;  

  // calculate the lengths of the two lines
  len1=sqrt(xD1*xD1+yD1*yD1);
  len2=sqrt(xD2*xD2+yD2*yD2);

  // calculate angle between the two lines.
  dot=(xD1*xD2+yD1*yD2); // dot product
  deg=dot/(len1*len2);

  // if abs(angle)==1 then the lines are parallell,
  // so no intersection is possible
  if(abs(deg)==1) return null;

  // find intersection Pt between two lines
  Point pt=new Point(0,0);
  div=yD2*xD1-xD2*yD1;
  ua=(xD2*yD3-yD2*xD3)/div;
  ub=(xD1*yD3-yD1*xD3)/div;
  pt.x=p1.x+ua*xD1;
  pt.y=p1.y+ua*yD1;

  // calculate the combined length of the two segments
  // between Pt-p1 and Pt-p2
  xD1=pt.x-p1.x;
  xD2=pt.x-p2.x;
  yD1=pt.y-p1.y;
  yD2=pt.y-p2.y;
  segmentLen1=sqrt(xD1*xD1+yD1*yD1)+sqrt(xD2*xD2+yD2*yD2);

  // calculate the combined length of the two segments
  // between Pt-p3 and Pt-p4
  xD1=pt.x-p3.x;
  xD2=pt.x-p4.x;
  yD1=pt.y-p3.y;
  yD2=pt.y-p4.y;
  segmentLen2=sqrt(xD1*xD1+yD1*yD1)+sqrt(xD2*xD2+yD2*yD2);

  // if the lengths of both sets of segments are the same as
  // the lenghts of the two lines the point is actually
  // on the line segment.

  // if the point isn’t on the line, return null
  if(abs(len1-segmentLen1)>0.01 || abs(len2-segmentLen2)>0.01)
    return null;

  // return the valid intersection
  return pt;
}

class Point{
  float x,y;
  Point(float x, float y){
    this.x = x;
    this.y = y;
  }

  void set(float x, float y){
    this.x = x;
    this.y = y;
  }
}

There are 7 comments to "Code: 2D line intersection". You may leave your own comment.
1. guillaume, September 11th, 2007 at 12:20

Here is a faster way to do this, no sqrt:

/*
MŽthode pour dŽtecter une intersection de deux lignes
AdaptŽ d’un programme en langage C
*/

int DONT_INTERSECT = 0;
int COLLINEAR = 1;
int DO_INTERSECT = 2;

float x =0, y=0;

void setup(){
size(320,320);
fill(255,0,0);
}

void draw(){

int intersect;

background(255);

// lignes
stroke(0);

// ligne fixe
line(20,height/2, width-20, (height/2)-20);

// ligne en mouvement
line(10,10,mouseX, mouseY);

intersect = intersect(20, height/2, width, (height/2)-20, 10, 10, mouseX, mouseY);

// dessiner le point d’intersection
noStroke();
if (intersect == DO_INTERSECT) ellipse(x, y, 5, 5);

}

int intersect(float x1, float y1, float x2, float y2, float x3, float y3, float x4, float y4){

float a1, a2, b1, b2, c1, c2;
float r1, r2 , r3, r4;
float denom, offset, num;

// Compute a1, b1, c1, where line joining points 1 and 2
// is “a1 x + b1 y + c1 = 0″.
a1 = y2 - y1;
b1 = x1 - x2;
c1 = (x2 * y1) - (x1 * y2);

// Compute r3 and r4.
r3 = ((a1 * x3) + (b1 * y3) + c1);
r4 = ((a1 * x4) + (b1 * y4) + c1);

// Check signs of r3 and r4. If both point 3 and point 4 lie on
// same side of line 1, the line segments do not intersect.
if ((r3 != 0) && (r4 != 0) && same_sign(r3, r4)){
return DONT_INTERSECT;
}

// Compute a2, b2, c2
a2 = y4 - y3;
b2 = x3 - x4;
c2 = (x4 * y3) - (x3 * y4);

// Compute r1 and r2
r1 = (a2 * x1) + (b2 * y1) + c2;
r2 = (a2 * x2) + (b2 * y2) + c2;

// Check signs of r1 and r2. If both point 1 and point 2 lie
// on same side of second line segment, the line segments do
// not intersect.
if ((r1 != 0) && (r2 != 0) && (same_sign(r1, r2))){
return DONT_INTERSECT;
}

//Line segments intersect: compute intersection point.
denom = (a1 * b2) - (a2 * b1);

if (denom == 0) {
return COLLINEAR;
}

if (denom < 0){
offset = -denom / 2;
}
else {
offset = denom / 2 ;
}

// The denom/2 is to get rounding instead of truncating. It
// is added or subtracted to the numerator, depending upon the
// sign of the numerator.
num = (b1 * c2) - (b2 * c1);
if (num < 0){
x = (num - offset) / denom;
}
else {
x = (num + offset) / denom;
}

num = (a2 * c1) - (a1 * c2);
if (num < 0){
y = ( num - offset) / denom;
}
else {
y = (num + offset) / denom;
}

// lines_intersect
return DO_INTERSECT;
}

boolean same_sign(float a, float b){

return (( a * b) >= 0);
}

2. guillaume, September 11th, 2007 at 12:28

sorry apparently I can’t directly copy/paste java code in the comment…
the complete code can be found at : http://www.ecole-art-aix.fr/article425.html

3. marius watz, September 11th, 2007 at 14:57

Hi Guillaume, thanks for the code. I knew mine wouldn’t be the fastest but I was happy just to get it working. I’ll try your version now.

4. christian, December 16th, 2007 at 12:29

I tried both the algorithms and marius’ one generates a more precise result. In my case I needed more precision than performance, cheers!

5. marius watz, December 17th, 2007 at 02:36

My algorithm more precise? Now that is a surprise.

6. zack, March 10th, 2008 at 20:44

shouldn’t be surprising… since any sqrt or division leads to a loss in precision… so the order of the math in the algorithm is just as important as which operations to do for both performance and precision. Any Numerical Analysis textbook can fill you in on the details on why this happens.

Comment on this entry

You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>