Checking if lines intersect and if so return the c

Question

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't intersect at all.

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

double                      dy[2], dx[2], m[2], b[2];
double                      xint, yint, xi, yi;
WsqT_Location_Message       *location_msg_ptr = OPC_NIL;

FIN (intersect (<args>));

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 = (WsqT_Location_Message*)op_prg_mem_alloc(sizeof(WsqT_Location_Message));

            location_msg_ptr->current_latitude = xi;
            location_msg_ptr->current_longitude = yi;
        }

        FRET(location_msg_ptr);
    }
}

FRET(location_msg_ptr);
}

Answer 1

There is an absolutely great and simple theory about lines and their intersections that is based on adding an extra dimensions to your points and lines. In this theory a line can be created from two points with one line of code and the point of line intersection can be calculated with one line of code. Moreover, points at the Infinity and lines at the Infinity can be represented with real numbers.

You probably heard about homogeneous representation when a point [x, y] is represented as [x, y, 1] and the line ax+by+c=0 is represented as [a, b, c]? The transitioning to Cartesian coordinates for a general homogeneous representation of a point [x, y, w] is [x/w, y/w]. This little trick makes all the difference including representation of lines at infinity (e.g. [1, 0, 0]) and making line representation look similar to point one. This introduces a GREAT symmetry into formulas for numerous line/point manipulation and is an absolute MUST to use in programming. For example,

It is very easy to find line intersections through vector product

p = l1xl2

A line can be created from two points is a similar way:

l=p1xp2

In the code of OpenCV it it just:

line = p1.cross(p2);
p = line1.cross(line2);

Note that there are no marginal cases (such as division by zero or parallel lines) to be concerned with here. My point is, I suggest to rewrite your code to take advantage of this elegant theory about lines and points.

Finally, if you don't use openCV, you can use a 3D point class and create your own cross product function similar to this one:

template<typename _Tp> inline Point3_<_Tp> Point3_<_Tp>::cross(const Point3_<_Tp>& pt) const
{
    return Point3_<_Tp>(y*pt.z - z*pt.y, z*pt.x - x*pt.z, x*pt.y - y*pt.x);
}