Line Segment Intersection for Partial or Lapsed Mathematicians

24th February 2013 by Ian Martin

Line+Segment+Intersection+for+Partial+or+Lapsed+Mathematicians Image

Line Segment Intersection for Partial or Lapsed Mathematicians

24th February 2013 by Ian Martin

Introduction

A common problem in 2D geometry is the determination of whether two line segments intersect, and following a positive answer, optionally deriving the intersection point. The most common solution is the use of cross product to derive the intersection point (if any), and whilst to a mathematician (or student thereof) this is a relatively simple solution, research for this article quickly revealed at least two published incorrect code implementations, flawed by their author's failure to understand that maths. The purpose of this article is therefore to present a less elegant but much simpler solution, easier to illustrate and comprehend. The article is therefore written with the assumption of an average knowledge of school mathematics.

Lines and Line Segments

In common parlance a 2D line is the shortest distance between two points, that description actually applies to line segments. A line is of infinite length; a line segment is a portion of the line, between two given points. Unless the two lines are parallel they will by definition intersect somewhere, given their inifinte lengths; the question of intersection of line segments is therefore the concern of this article. The image below shows two line segments AB and CD

A Simpler Solution

To solve the line segment intersection problem we will use the fact the the intersection calculation is very simple if one of the lines is vertical or horizontal. We will translate and rotate the "world" so that this is the case; we imagine ourselves at point A, looking along the line towards B, and considering the "world" relative to our own position and orentation. This solution, described below in stages, also allows for some trivial (i.e. fast) checks that can determine that the line segments do not cross, before the full calculations are executed.

Stage 1 - Initial Checks

Considering again two line segments AB and CD. If the lines are identical (i.e. A and C, and B and D are at the same position (or A and D and B and C) the line segments should be designated as crossing or not crossing, subject to the ultimate objective of the test. No further checks are needed. If both line segments are actually points (i.e. A=B and C=D), execute a simple test to determine if the the points coincide. No further checks are needed. Finally if one of the line segment ends is the same as one of the other ends (but only one, not both), the line segments would normally be designated as not crossing. This is shown below, where B and C are coincidental.

Stage 2 - Trivial Rejection

Consider the line segments below. We can clearly see that they do not cross. If we check the bounding boxes of each line segment, we can quickly in this case deduce that the line segments cannot cross, and no further checks are needed. The follow pseudo-code is used:- if (box1.xmax < box2.xmin OR box1.ymax < box2.ymin OR box1.xmin > box2.xmax OR box1.ymin > box2.ymax) then reject In the case shown below, line1's box.xmax is less than line2's box xmin, so they can never cross.

Stage 3 - Translate and Rotate Relative to A

Considering again the line segments AB and CD. Clearly they cross; we will now show how to determine whether they do, and optionally calculate the point at which they do.

In order to determine whether they cross, we will rotate and shift the "world" of AB and CD so that AB is vertical, with A being at the centre (0,0). We will view the "world" from the perspective of point A, looking at point B.

First, translate A to (0,0) and move B,C and D by the same amount:- Subtract A's X and Y coordinates from points B,C and D:-

Now we will rotate the "world" around 0,0 so that AB lies at X=0 (on the Y axis). We will rotate both AB and CD around (0,0) by the angle between AB and the Y axis (X=0); we will designate this angle a. Here's where the school maths comes in. The Pythagoras formula shows that for a right angled triangle, the square of the longest edge = the sum of the square of the other two edges. Considering the line segment AB shown below:-

If we complete the right angled triangle using point E which is (A.x,B.y), the angle a is that between AB and AE. Using Pythagoras, the length or magnitude of line segment AB is derived thus:- sqrt( (B.x-A.x)^2 + (B.y-A.y>^2) ) In order to rotate points by a given angle, we need the sine and cosine of that angle. For angle a these are:- sin a = length(EB) / length(AB) or sin = B.x / length(AB) cos a = length(EA) / length(AB) or cos = B.y / length(AB) Remember that we calculated length(AB) using Pythagoras, above. To rotate points around (0,0) by angle a, use these formulae, which you may have encountered in school maths:- x' = x cos a - y sin a y' = y cos a + x sin a When the above formula is applied to AB and CD we see the result below. (A shrewd observer will note that A (0,0) does not need to be rotated, and the new position B' of B is (0,length(AB)).) Remember that this is the same "world", merely viewed from the perspetive of A, looking towards B. We will calculate the intersection point (if any) relative to point A - the fact of intersection is not changed by the co-ordinates being relative to A.

Allowing Tolerance

Even when using high precision floating point code, slight errors must be allowed for. Furthermore, although lines and line segments are normally regarded as having zero with, your application may consider them to be "thick". For these reasons we need to add a small tolerance to some checks applied after the rotation has been applied. We will define a small value which will be used for the tolerance. For example:- const double small_value = 0.00001;

Stage 4 - Further Trivial Checks

The new position of C'D' relative to AB' allows us to easily ascertain if C'D' crosses A'B'. We can see from the above diagram that the following pseudo-code checks can be applied:- if C'.y < 0 AND D'.y < 0 they do not cross if C'.y > B'.y AND D'.y > B'.y they do not cross if C'.x < 0 AND D'.x < 0 they do not cross if C'.x > 0 AND D'.x > 0 they do not cross

Stage 5 - Check for Co-linear

If line segment C'D' is now on the Y axis (abs(C'.x) < small_value) AND abs(D'.x) < small_value) then the line segments lie on the same infinite line. The only questions that now concerns is whether they share the same space on that line. The stage 4 trivial checks have determined that A'B' and C'D' are not above or below each other, we know that if C'D' lies at or very close to X=0 it must cross A'B'.

Stage 6 - Check for Parallel

If line segment C'D' is now vertical ((abs(C'.x - D'.x) < small_value) we know that they cannot cross, and we also know that it is not co-linear, having already performed the stage 5 check.

Stage 7 - Do They Cross?

Having completed the above trivial checks, if the question is still not settled, we must calculate where C'D' crosses the Y axis (X=0). Here is some more maths, but really simple, similar triangles. Look at C'D' crossing the Y axis below at point F:-

The smaller shaded triangle is "similar" to the larger triangle formed by C'D' and (D'.x,C'.y). To compute the Y position of point F:- F.y = C'.y - (((D'.y - C'.y) * C'.x) / (D'.x - C'.x) ) Note that we already know that C'D' is not vertical, so D'.x-C'x cannot be zero - there is no need to apply a division by zero check before executing the above calculation. We now have the position of F. We simply need to check it against line A'B'. Thus, if F.y < 0 or F.y> B'.y, the line segments do not cross.

Stage 8 - Calculate the Intersection Point (optional)

If we need the intersection position, we will need to reverse the rotation and translation of stage 3, upon point F. First, rotate point F the other way, using the rotation through the angle -a, using again these formulae:- x' = x cos a - y sin a y' = y cos a + x sin a ...but since F.x=0, sin(-a) = -sin(a), and cos(-a) = cos(a), the operation becomes much simpler:- F'.x = F.y * sin a F'.y = F.y * cos a Then, translate the rotated F' by the position of A:- F'.x = F'x + A.x F'.y = F'y + A.y ... to derive the actual "world" position of the intersection point F, rather than relative to AB.

Disadvantages

Clearly this method is not as elegant as the cross product solution. The single square root required is relatively slow.

Advantages

This method much is easier to comprehend and illustrate. The trivial checks at different stages can save un-needed calculations.

Improvements

If a complex system, where many line segments must be checked against each other, a simple improvment to performance could be derived by caching the line segment lengths as they are calculated.

Adding Figures

  1. Using the line segments AB and CD in the worked example:-AB is (-5,-2) to (1,2) CD is (-8,-2) to (2,1)
  2. The boundary boxes cross, so we proceed to translate so that A is at (0,0).
  3. Subtract A from B,C,D i.e. subtract (-5,-2)
  4. AB is now (0,0) to (6,4) CD is now (-3,0) to (7,3)
  5. Calculate length(AB) using Pythagoras. Length(AB) = 7.2111025
  6. Calculate sin(a) = B.x / length(AB) = 0.83205029433784372
  7. Calculate cos(a) = B.y / length(AB) = 0.55470019622522915
  8. Rotate AB and CD around (0,0) A'B' is (0,0) to (0,7.2111025) C'D' is (-1.66410059,-2.49615) to (1.38675,7.4884)
  9. C'D' is not co-linear or vertical, and crosses X=0, so calculate where it crosses.
  10. C'D' crosses X=0 at point F.
  11. F.Y = 2.9499964 (F.X = 0, for clarity)
  12. Check if F.Y is between 0 and B'.Y or 7.2111025. It is, so the line segments do cross.
  13. Rotating F the opposite way, and translating by -A, we calculate the intersection point(-2.5454545,-0.3636363740)

Code Example

This C# example uses two classes Point2D and LineSeg2D 
public struct ConstantValue

{

    internal const double SmallValue = 0.00001;

}



public class Point2D

{

    private float X;

    private float Y;



    public Point2D(double src_x, double src_y)

    {

        this.X = (float)src_x;

        this.Y = (float)src_y;

    }



    public Point2D(float src_x, float src_y)

    {

        this.X = src_x;

        this.Y = src_y;

    }



    public float X

    {

        set

        {

            this.X = value;

        }

        get

        {

            return this.X;

        }

    }



    public float Y

    {

        set

        {

            this.Y = value;

        }

        get

        {

            return this.Y;

        }

    }

}



public class LineSeg2D

{

    private Point2D m_Start;

    private Point2D m_End;



    public LineSeg2D(Point2D start, Point2D end)

    {

        this.m_Start = startPoint;

        this.m_End = endPoint;

    }



    // returns true if intersection, and returns position

    // if a buffer is supplied (intersect_point)

    //

    public bool CheckIntersect(Point2D intersect_point, LineSeg2D other_line)

    {



        //check if points are same



        if (

            this.m_End.X == other_line.m_End.X &&

            this.m_End.Y == other_line.m_End.Y &&

            this.m_Start.X == other_line.m_Start.X &&

            this.m_Start.Y == other_line.m_Start.Y

            )

        {

            //same other_line - mark as crossing

            if (intersect_point != null)

            {

                intersect_point.X = this.m_End.X;

                intersect_point.Y = this.m_End.Y;

            }



            return (true);

        }



        if (

            this.m_End.X == other_line.m_Start.X &&

            this.m_End.Y == other_line.m_Start.Y &&

            this.m_Start.X == other_line.m_End.X &&

            this.m_Start.Y == other_line.m_End.Y

            )

        {

            //same other_line - mark as crossing

            if (intersect_point != null)

            {

                intersect_point.X = this.m_End.X;

                intersect_point.Y = this.m_End.Y;

            }

            return true;

        }



        //allow end points to coincide

        if (

            this.m_End.X == other_line.m_Start.X &&

            this.m_End.Y == other_line.m_Start.Y

            )

        {

            return false;

        }



        if (

            other_line.m_End.X == this.m_Start.X &&

            other_line.m_End.Y == this.m_Start.Y

            )

        {

            return false;

        }



        //check line1 is not a point



        if (

            this.m_End.X == this.m_Start.X &&

            this.m_End.Y == this.m_Start.Y

            )

        {

            //line1 is a point

            return false;

        }



        //check line2 is not a point

        if (

            this.m_End.X == this.m_Start.X &&

            this.m_End.Y == this.m_Start.Y

            )

        {

            //line1 is a point

            return false;

        }



        //get the magnitude of line1

        double line1_dx = this.m_Start.X - this.m_End.X;

        double line1_dy = this.m_Start.Y - this.m_End.Y;

        double line1_mag = Math.Sqrt((line1_dx * line1_dx)

        + (line1_dy * line1_dy));



        //shift line2 relative to line1 START POINT



        double line2_sx = other_line.m_End.X - this.m_End.X;

        double line2_sy = other_line.m_End.Y - this.m_End.Y;



        double line2_fx = other_line.m_Start.X - this.m_End.X;

        double line2_fy = other_line.m_Start.Y - this.m_End.Y;



        double sin = line1_dx / line1_mag;

        double cos = line1_dy / line1_mag;



        //rotate



        //x' = x cos f - y sin f

        //y' = y cos f + x sin f



        double line2a_sx = (line2_sx * cos) - (line2_sy * sin);

        double line2a_sy = (line2_sy * cos) + (line2_sx * sin);



        double line2a_fx = (line2_fx * cos) - (line2_fy * sin);

        double line2a_fy = (line2_fy * cos) + (line2_fx * sin);



        //check if co-linear

        if (Math.Abs(line2a_sx) < ConstantValue.SmallValue &&

           Math.Abs(line2a_fx) < ConstantValue.SmallValue         )         {                  //line2 lies on the Y axis                if (line2a_sy > line2a_fy)

            {

                double t = line2a_sy;

                line2a_sy = line2a_fy;

                line2a_fy = t;

            }



            if (line2a_sy > line1_mag)

            {

                //line2 is ABOVE line1

                return false;

            }

            if (line2a_fy < 0)

            {

                //line2 is BELOW line1

                return false;

            }



            //lines are co-linear and DO interesect



            if (intersect_point != null)

            {

                intersect_point.X = this.m_End.X;

                intersect_point.Y = this.m_End.Y;

            }

            return true;



        }



        //ensure line2 crosses Y axis

        if (line2a_sx < 0 && line2a_fx < 0)         {                 //both ends are BELOW x=0                 return false;         }                  if (line2a_fx > 0 && line2a_sx > 0)

        {

            //both ends ABOVE y=0

            return false;

        }



        //check if line2a is all

        if (line2a_sy < 0 && line2a_fy < 0)                      {                              return false;                     }                              //check if line2a is all >mag1

        if (line2a_sy >line1_mag && line2a_fy > line1_mag)

        {

            return false;

        }



        //check if line2a is parallel

        if (Math.Abs(line2a_sx - line2a_fx) < ConstantValue.SmallValue)

        {

            return false;

        }



        //now calculate where line2a crosses Y axis



        //using "similar trianges"

        double ry = ((line2a_fy - line2a_sy) * line2a_sx) /

             (line2a_fx - line2a_sx);



        ry = line2a_sy-ry;        //ry is where it crosses x=0;



        if (ry < 0 || ry > line1_mag)

        {

            //crosses Y axis OUTSIDE of line1

            return false;

        }



        if (intersect_point != null)

        {

            // calculate intersection point by roating (0,ry)

            // by the reverse of the original

            // rotation, and adding to line1 start point

            //

            intersect_point.X = (float)((ry * sin) + this.m_End.X);

            intersect_point.Y = (float)((ry * cos) + this.m_End.Y);

        }

        return true;

    }



}