How do you check for intersection between a line segment and a line ray emanating from a point at an angle from horizontal?

Question

Given a line segment, that is two points (x1,y1) and (x2,y2), one point P(x,y) and an angle theta. How do we find if this line segment and the line ray that emanates from P at an angle theta from horizontal intersects or not? If they do intersect, how to find the point of intersection?

Trina

Answer 1

Let's label the points q = (x1, y1) and q + s = (x2, y2). Hence s = (x2 − x1, y2 − y1). Then the problem looks like this:

Let r = (cos θ, sin θ). Then any point on the ray through p is representable as p + t r (for a scalar parameter 0 ≤ t) and any point on the line segment is representable as q + u s (for a scalar parameter 0 ≤ u ≤ 1).

The two lines intersect if we can find t and u such that p + t r = q + u s:

See this answer for how to find this point (or determine that there is no such point).

Then your line segment intersects the ray if 0 ≤ t and 0 ≤ u ≤ 1.

Lesbos

Answer 2

Here is a C# code for the algorithm given in other answers:

    /// <summary>
    /// Returns the distance from the ray origin to the intersection point or null if there is no intersection.
    /// </summary>
    public double? GetRayToLineSegmentIntersection(Point rayOrigin, Vector rayDirection, Point point1, Point point2)
    {
        var v1 = rayOrigin - point1;
        var v2 = point2 - point1;
        var v3 = new Vector(-rayDirection.Y, rayDirection.X);


        var dot = v2 * v3;
        if (Math.Abs(dot) < 0.000001)
            return null;

        var t1 = Vector.CrossProduct(v2, v1) / dot;
        var t2 = (v1 * v3) / dot;

        if (t1 >= 0.0 && (t2 >= 0.0 && t2 <= 1.0))
            return t1;

        return null;
    }

Colorant

This is only one part of the ray vs line segment intersection test. You didn't consider the possibility of the ray and the line segment being parallel. In that case the dot product dot(v2, v3) will be 0 (or really close to 0). – Salema

@Thash - thank you, I have added that check. Didn't check it though, so I hope it's working correct. – Colorant

@Colorant unfortunately, this is still not enough. In this case, they can intersect if the ray origin and line segment endpoints are collinear (also note that the ray origin could actually be inside the line segment). Solving special cases in computational geometry is fun :) – Salema

I dont understand. The intersection point is a double? I would expect intersection point to be a point. – Svoboda

@SyaifulNizamYahya The double value is the distance to the intersection point from the ray origin along the ray. You can get the intersection point this way: rayOrigin + rayDirection * distance, given rayDirection is normalized (has length of 1). – Colorant

Maybe it would be better if you change the summary to like "Returns the distance between rayorigin and intersection point" – Svoboda

Thank to you, I finally got it working. See here #53893792 – Svoboda

@SyaifulNizamYahya I clarified the summary comment, thank you. – Colorant

Answer 3

Thanks Gareth for a great answer. Here is the solution implemented in Python. Feel free to remove the tests and just copy paste the actual function. I have followed the write-up of the methods that appeared here, https://rootllama.wordpress.com/2014/06/20/ray-line-segment-intersection-test-in-2d/.

import numpy as np

def magnitude(vector):
   return np.sqrt(np.dot(np.array(vector),np.array(vector)))

def norm(vector):
   return np.array(vector)/magnitude(np.array(vector))

def lineRayIntersectionPoint(rayOrigin, rayDirection, point1, point2):
    """
    >>> # Line segment
    >>> z1 = (0,0)
    >>> z2 = (10, 10)
    >>>
    >>> # Test ray 1 -- intersecting ray
    >>> r = (0, 5)
    >>> d = norm((1,0))
    >>> len(lineRayIntersectionPoint(r,d,z1,z2)) == 1
    True
    >>> # Test ray 2 -- intersecting ray
    >>> r = (5, 0)
    >>> d = norm((0,1))
    >>> len(lineRayIntersectionPoint(r,d,z1,z2)) == 1
    True
    >>> # Test ray 3 -- intersecting perpendicular ray
    >>> r0 = (0,10)
    >>> r1 = (10,0)
    >>> d = norm(np.array(r1)-np.array(r0))
    >>> len(lineRayIntersectionPoint(r0,d,z1,z2)) == 1
    True
    >>> # Test ray 4 -- intersecting perpendicular ray
    >>> r0 = (0, 10)
    >>> r1 = (10, 0)
    >>> d = norm(np.array(r0)-np.array(r1))
    >>> len(lineRayIntersectionPoint(r1,d,z1,z2)) == 1
    True
    >>> # Test ray 5 -- non intersecting anti-parallel ray
    >>> r = (-2, 0)
    >>> d = norm(np.array(z1)-np.array(z2))
    >>> len(lineRayIntersectionPoint(r,d,z1,z2)) == 0
    True
    >>> # Test ray 6 --intersecting perpendicular ray
    >>> r = (-2, 0)
    >>> d = norm(np.array(z1)-np.array(z2))
    >>> len(lineRayIntersectionPoint(r,d,z1,z2)) == 0
    True
    """
    # Convert to numpy arrays
    rayOrigin = np.array(rayOrigin, dtype=np.float)
    rayDirection = np.array(norm(rayDirection), dtype=np.float)
    point1 = np.array(point1, dtype=np.float)
    point2 = np.array(point2, dtype=np.float)

    # Ray-Line Segment Intersection Test in 2D
    # http://bit.ly/1CoxdrG
    v1 = rayOrigin - point1
    v2 = point2 - point1
    v3 = np.array([-rayDirection[1], rayDirection[0]])
    t1 = np.cross(v2, v1) / np.dot(v2, v3)
    t2 = np.dot(v1, v3) / np.dot(v2, v3)
    if t1 >= 0.0 and t2 >= 0.0 and t2 <= 1.0:
        return [rayOrigin + t1 * rayDirection]
    return []

if __name__ == "__main__":
    import doctest
    doctest.testmod()

Sarawak

# Test ray 6 is incorrect. It is a duplicate of # Test ray 5. – Modiolus

Answer 4

Note: this solution works without making vector classes or defining vector multiplication/division, but is longer to implement. It also avoids division by zero errors. If you just want a block of code and don’t care about the derivation, scroll to the bottom of the post.

Let’s say we have a ray defined by x, y, and theta, and a line defined by x1, y1, x2, and y2.

First, let’s draw two rays that point from the ray’s origin to the ends of the line segment. In pseudocode, that’s

theta1 = atan2(y1-y, x1-x);
theta2 = atan2(y2-y, x2-x);

Next we check whether the ray is inside these two new rays. They all have the same origin, so we only have to check the angles.

To make this easier, let’s shift all the angles so theta1 is on the x axis, then put everything back into a range of -pi to pi. In pseudocode that’s

dtheta = theta2-theta1; //this is where theta2 ends up after shifting
ntheta = theta-theta1; //this is where the ray ends up after shifting
dtheta = atan2(sin(dtheta), cos(dtheta))
ntheta = atan2(sin(ntheta), cos(ntheta))

(Note: Taking the atan2 of the sin and cos of the angle just resets the range of the angle to within -pi and pi without changing the angle.)

Now imagine drawing a line from theta2’s new location (dtheta) to theta1’s new location (0 radians). That’s where the line segment ended up.

The only time where the ray intersects the line segment is when theta is between theta1 and theta2, which is the same as when ntheta is between dtheta and 0 radians. Here is the corresponding pseudocode:

sign(ntheta)==sign(dtheta)&&
abs(ntheta)<=abs(dtheta)

This will tell you if the two lines intersect. Here it is in one pseudocode block:

theta1=atan2(y1-y, x1-x);
theta2=atan2(y2-y, x2-x);
dtheta=theta2-theta1;
ntheta=theta-theta1;
dtheta=atan2(sin(dtheta), cos(dtheta))
ntheta=atan2(sin(ntheta), cos(ntheta))
return (sign(ntheta)==sign(dtheta)&&
abs(ntheta)<=abs(dtheta));

Now that we know the points intersect, we need to find the point of intersection. We’ll be working from a completely clean slate here, so we can ignore any code up to this part. To find the point of intersection, you can use the following system of equations and solve for xp and yp, where lb and rb are the y-intercepts of the line segment and the ray, respectively.

y1=(y2-y1)/(x2-x1)*x1+lb
yp=(y2-y1)/(x2-x1)*xp+lb
y=sin(theta)/cos(theta)*x+rb
yp=sin(theta)/cos(theta)*x+rb

This yields the following formulas for xp and yp:

xp=(y1-(y2-y1)/(x2-x1)*x1-y+sin(theta)/cos(theta)*x)/(sin(theta)/cos(theta)-(y2-y1)/(x2-x1));
yp=sin(theta)/cos(theta)*xp+y-sin(theta)/cos(theta)*x

Which can be shortened by using lm=(y2-y1)/(x2-x1) and rm=sin(theta)/cos(theta)

xp=(y1-lm*x1-y+rm*x)/(rm-lm);
yp=rm*xp+y-rm*x;

You can avoid division by zero errors by checking if either the line or the ray is vertical then replacing every x and y with each other. Here’s the corresponding pseudocode:

if(x2-x1==0||cos(theta)==0){
    let rm=cos(theta)/sin(theta);
    let lm=(x2-x1)/(y2-y1);
    yp=(x1-lm*y1-x+rm*y)/(rm-lm);
    xp=rm*yp+x-rm*y;
}else{
    let rm=sin(theta)/cos(theta);
    let lm=(y2-y1)/(x2-x1);
    xp=(y1-lm*x1-y+rm*x)/(rm-lm);
    yp=rm*xp+y-rm*x;
}

TL;DR:

bool intersects(x1, y1, x2, y2, x, y, theta){
    theta1=atan2(y1-y, x1-x);
    theta2=atan2(y2-y, x2-x);
    dtheta=theta2-theta1;
    ntheta=theta-theta1;
    dtheta=atan2(sin(dtheta), cos(dtheta))
    ntheta=atan2(sin(ntheta), cos(ntheta))
    return (sign(ntheta)==sign(dtheta)&&abs(ntheta)<=abs(dtheta));
}
point intersection(x1, y1, x2, y2, x, y, theta){
    let xp, yp;
    if(x2-x1==0||cos(theta)==0){
        let rm=cos(theta)/sin(theta);
        let lm=(x2-x1)/(y2-y1);
        yp=(x1-lm*y1-x+rm*y)/(rm-lm);
        xp=rm*yp+x-rm*y;
    }else{
        let rm=sin(theta)/cos(theta);
        let lm=(y2-y1)/(x2-x1);
        xp=(y1-lm*x1-y+rm*x)/(rm-lm);
        yp=rm*xp+y-rm*x;
    }
    return (xp, yp);
}

Sayce

Do you intend to say the two new rays point from all three rays' origin to the beginning and end of the line segment? It might also be helpful if you include the intermediate steps in getting to "This yields the following formulas for xp and yp:" – Cotswold

Could you write out more of the calculations for xx an yp, specifically the intermediate steps? – Cotswold

@EllieKesselman yes, the two new rays share an origin with the original ray, and point to the endpoints of the line segment. Because the rays extend forever, all you need to know to determine if the intersection occurs is whether or not the original ray is between the two new rays or not. – Sayce