Check if two line segments intersect

Hi,

I've written some code below to check if two line segments intersect and if they do to tell me where. As input I have the (x,y) coordinates of both ends of each line. It appeared to be working correctly but now in the scenario where line A (532.87,787.79)(486.34,769.85) and line B (490.89,764.018)(478.98,783.129) it says they intersect at (770.136, 487.08) when the lines don' intersect at all.

Has anyone any idea what is incorrect in the below code? Thanks in advance

Code:

dy[0] = y2 - y1;
	dx[0] = x2 - x1;
	dy[1] = y4 - y3;
	dx[1] = x4 - x3;
	
	m[0] = dy[0] / dx[0];
	m[1] = dy[1] / dx[1];
	b[0] = y1 - m[0] * x1;
	b[1] = y3 - m[1] * x3;

	if (m[0] != m[1])
	{
		//slopes not equal, compute intercept
		xint = (b[0] - b[1]) / (m[1] - m[0]);
		yint = m[1] * xint + b[1];
		
		//is intercept in both line segments?
		if ((xint <= max(x1, x2)) && (xint >= min(x1, x2)) &&
			(yint <= max(y1, y2)) && (yint >= min(y1, y2)) &&
			(xint <= max(x3, x4)) && (xint >= min(x3, x4)) &&
			(yint <= max(y3, y4)) && (yint >= min(y3, y4)))
		{
			if (xi && yi)
			{
				xi = xint;
				yi = yint;
				
				location_msg_ptr->current_latitude = xi;
				location_msg_ptr->current_longitude = yi;
			}
			
			return(location_msg_ptr);
		}
	}

	return(location_msg_ptr);
}

If the slopes are equal, or if the intercept is not in the line segments, then that code snippet returns location_msg_pointer. The caller has no way of knowing there is no intersect.

Sorry I should have mentioned that the location_msg_ptr is set to NULL at the beginning of the code and when I get the result back from the function I check if it's NULL. If it is I assume no intersection of lines and if not I assume there is and get back the x,y position of the intersection. However in the case that I outlined an intersection is found and a valid location_msg_ptr is returned (the one in the nested if statements).

Any idea?
Thanks.

Note, floating point computations in computers have inherent inaccuracy, so direct == or != comparisons are not recommended. Usually, you check whether they are within some tolerance of each other, like fabs(a - b) < eps, where eps is epsilon, your tolerance. If the magnitude of the difference between a and b is less than eps, you consider a and b "equal enough". If it's >= eps, you consider a and b not equal. For your case, something around 0.0001 or 0.0000001 might be a good eps.

You should read this: What Every Computer Scientist Should Know About Floating-Point Arithmetic.

Now, as for your algorithm, you find the slope correctly. You also appear to be checking whether the intersection is within the bounds of those segments. In the example you gave, the segments do intersect. Check out this link: FooPlot | Online graphing calculator and function plotter

Google Is your friend. I would start with the code for get_line_intersection from here:
geometry - How do you detect where two line segments intersect? - Stack Overflow
Though personally I would eliminate one of the divisions and defer the other division until after I know that there is an intersection. Just calculate and store the divisors separately and check for <= to them instead of <= 1.

Edit: Oh what the heck, here is the adapted code:

Code:

// Returns 1 if the lines intersect, otherwise 0. In addition, if the lines 
// intersect the intersection point may be stored in the floats i_x and i_y.
char get_line_intersection(float p0_x, float p0_y, float p1_x, float p1_y, 
    float p2_x, float p2_y, float p3_x, float p3_y, float *i_x, float *i_y)
{
    float s1_x, s1_y, s2_x, s2_y, sn, tn, sd, td, t;
    s1_x = p1_x - p0_x;     s1_y = p1_y - p0_y;
    s2_x = p3_x - p2_x;     s2_y = p3_y - p2_y;

    sn = -s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y);
    sd = -s2_x * s1_y + s1_x * s2_y;
    tn =  s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x);
    td = -s2_x * s1_y + s1_x * s2_y;

    if (sn >= 0 && sn <= sd && tn >= 0 && tn <= td)
    {
        // Collision detected
        t = tn / td;
        if (i_x != NULL)
            *i_x = p0_x + (tn * s1_x);
        if (i_y != NULL)
            *i_y = p0_y + (tn * s1_y);
        return 1;
    }

    return 0; // No collision
}

This tends to eliminate problems with division by zero as well as give a speedup. I think this about matches the exact code I've used in the past.
Up to you to adapt that to your own code.

While the biggest names of the forum seems to have given an answer, I think that I must quote this

Moreover, our methods use only additions,
subtractions, multiplications, and comparisons. We need neither division
nor trigonometric functions, both of which can be computationally expensive and
prone to problems with round-off error. For example, the “straightforward” method
of determining whether two segments intersect—compute the line equation of the
form y D mx C b for each segment (m is the slope and b is the y-intercept),
find the point of intersection of the lines, and check whether this point is on both
segments—uses division to find the point of intersection. When the segments are
nearly parallel, this method is very sensitive to the precision of the division operation
on real computers. The method in this section, which avoids division, is much
more accurate.

This comes from the third edition of Cormen's book. Cormen's book with title introduction to algorithms I guess. A book that I recommend everybody to read. If you do not own the book, just contact me.

So, in a chapter about line segments he describes how easily and naturally we can use cross products to determine things about line segments. One of those things is if two line segments intersect. The reason we use the cross products is the obvious one. (I have to attend a party in 10 minutes, so I do not have the time to analyze, but everything is described in the book).

The pseudocode:

Code:

SEGMENTS-INTERSECT.p1; p2; p3; p4/
1 d1 D DIRECTION.p3; p4; p1/
2 d2 D DIRECTION.p3; p4; p2/
3 d3 D DIRECTION.p1; p2; p3/
4 d4 D DIRECTION.p1; p2; p4/
5 if ..d1 > 0 and d2 < 0/ or .d1 < 0 and d2 > 0// and
..d3 > 0 and d4 < 0/ or .d3 < 0 and d4 > 0//
6 return TRUE
7 elseif d1 == 0 and ON-SEGMENT.p3; p4; p1/
8 return TRUE
9 elseif d2 == 0 and ON-SEGMENT.p3; p4; p2/
10 return TRUE
11 elseif d3 == 0 and ON-SEGMENT.p1; p2; p3/
12 return TRUE
13 elseif d4 == 0 and ON-SEGMENT.p1; p2; p4/
14 return TRUE
15 else return FALSE
DIRECTION.pi; pj; pk/
1 return .pk  pi / 	 .pj  pi /
ON-SEGMENT.pi; pj; pk/
1 if min.xi; xj /  xk  max.xi; xj / and min.yi; yj /  yk  max.yi; yj /
2 return TRUE
3 else return FALSE

//the pseudocode is not readable. Contact me to provide you the info, or I will upload tomorrow some screen shots Must get ready for the party now :P

Must be a good book then if he also tries to avoid divisions.

As you would know, the cross-product is a 3D only thing. When using this algorithm in 3D it detemines not necessarily that the lines intersect, but rather finds the points of the closest approach of either line to the other.
When doing this in strictly 2D, the equivalent thing is called the "perp dot product". It just so happens that the code already shown actually does use the perp dot product. So in all liklyhood, it is the same method.

Further early out tests can bone done by noticing that the calculations for tn and td are not required until it is known that sn is between 0 and sd. Even then, until you know that sn >= 0, you don't need to calculate sd, and similar for td.
This gives rise to:

Code:

char get_line_intersection(float p0_x, float p0_y, float p1_x, float p1_y, 
float p2_x, float p2_y, float p3_x, float p3_y, float *i_x, float *i_y)
{
    float s1_x, s1_y, s2_x, s2_y, sn, tn, sd, td, t;
    s1_x = p1_x - p0_x; s1_y = p1_y - p0_y;
    s2_x = p3_x - p2_x; s2_y = p3_y - p2_y;

    sn = -s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y);
    if (sn >= 0)
    {
        sd = -s2_x * s1_y + s1_x * s2_y;
        if (sn <= sd)
        {
            tn = s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x);
            if (tn >= 0)
            {
                td = -s2_x * s1_y + s1_x * s2_y;
                if (tn <= td)
                {
                    // Collision detected
                    t = tn / td;
                    if (i_x != NULL)
                        *i_x = p0_x + (tn * s1_x);
                    if (i_y != NULL)
                        *i_y = p0_y + (tn * s1_y);
                    return 1;
                }
            }
        }
    }
    return 0; // No collision
}

Which if you modify from using deeply-nested statement to using early returns, becomes this:

Code:

char get_line_intersection(float p0_x, float p0_y, float p1_x, float p1_y, 
    float p2_x, float p2_y, float p3_x, float p3_y, float *i_x, float *i_y)
{
    float s1_x, s1_y, s2_x, s2_y, sn, tn, sd, td, t;
    s1_x = p1_x - p0_x; s1_y = p1_y - p0_y;
    s2_x = p3_x - p2_x; s2_y = p3_y - p2_y;

    sn = -s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y);
    if (sn < 0)
        return 0; // No collision
    sd = -s2_x * s1_y + s1_x * s2_y;
    if (sn > sd)
        return 0; // No collision
    tn = s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x);
    if (tn < 0)
        return 0; // No collision
    td = -s2_x * s1_y + s1_x * s2_y;
    if (tn > td)
        return 0; // No collision

    // Collision detected
    t = tn / td;
    if (i_x != NULL)
        *i_x = p0_x + (tn * s1_x);
    if (i_y != NULL)
        *i_y = p0_y + (tn * s1_y);
    return 1;
}

Look much closer to your pseudocode yet?
Hmm, wasn't that fun

Sorry, there was a slight bug in the above, the last two occurrences of tn were supposed to be t. In any case I also noticed that some calculations were repeated and those are now factored out. Here's the fixed code:

Code:

int get_line_intersection(float p0_x, float p0_y, float p1_x, float p1_y, 
    float p2_x, float p2_y, float p3_x, float p3_y, float *i_x, float *i_y)
{
    float s02_x, s02_y, s10_x, s10_y, s32_x, s32_y, s_numer, t_numer, denom, t;
    s10_x = p1_x - p0_x;
    s10_y = p1_y - p0_y;
    s02_x = p0_x - p2_x;
    s02_y = p0_y - p2_y;

    s_numer = s10_x * s02_y - s10_y * s02_x;
    if (s_numer < 0)
        return 0; // No collision

    s32_x = p3_x - p2_x;
    s32_y = p3_y - p2_y;
    t_numer = s32_x * s02_y - s32_y * s02_x;
    if (t_numer < 0)
        return 0; // No collision

    denom = s10_x * s32_y - s32_x * s10_y;
    if (s_numer > denom || t_numer > denom)
        return 0; // No collision

    // Collision detected
    t = t_numer / denom;
    if (i_x != NULL)
        *i_x = p0_x + (t * s10_x);
    if (i_y != NULL)
        *i_y = p0_y + (t * s10_y);

    return 1;
}

The division can be avoided if the intersection point isn't required, but I'll leave that to others.