# Checking if two line segments intersect

What is the fastest way of checking whether two line segments of the form Line[{p1,p2}] intersect at a point? I don't need to know what the point of intersection is, but only if they do intersect.

My initial approach was to simply use RegionIntersection and write something like

LineIntersectionQ[l1_, l2_] := RegionDimension@RegionIntersection[l1, l2] == 0


But this seems to be slow. Any ideas?

• Commented Sep 19, 2021 at 1:37
• Just to get definitions clear here: Do you mean to use this for Line objects in the plane (2D) only, or also in 3D or even nD? Commented Sep 19, 2021 at 8:25
• @HenrikSchumacher I asked for general dimensions, but I'd be happy with a more efficient solution for either the 2D or 3D cases. Thanks! Commented Sep 19, 2021 at 11:30

You could try this:

Clear[interQ];
interQ[ Line[{v1_?VectorQ, v2_?VectorQ}],
Line[{v3_?VectorQ, v4_?VectorQ}]] := Block[{
m = Transpose[ {v1 - v2, v4 - v3}], vL},

If[ Det[m] == 0, False,
vL = LeastSquares[ m , v4 - v2];
(* Sow[ v1*vL[[1]] + v2*(1-vL[[1]])];  *)
0 <= vL[[1]] <= 1 && 0 <= vL[[2]] <= 1]];


Edit 2022-6-10: Adding cases for higher dimensional spaces as requested.

(* Assuming v1,v2,v3, and v4 have rational or integer coordinates *)

interExactQ[
Line[{v1_?VectorQ, v2_?VectorQ}],
Line[{v3_?VectorQ, v4_?VectorQ}] ] := Block[{
m = Transpose[{v2 - v1, v3 - v4}], vL},
vL = LeastSquares[ m, v3 - v1];
(*
Print["Dist Between Lines = ",  Norm[ m . vL - (v3-v1)]];
Print["Numerical Lambdas = ", N[vL]]; *)

Norm[ m . vL - (v3 - v1)] == 0 &&
0 <= vL[[1]] <= 1 &&
0 <= vL[[2]] <= 1 ];


(* For double precision floating point numbers there are numerical issues, the code below seems to work. If the coordinates have iD digits of precision, replace 10^14 with 10^(iD-2) *)

interDoubleRealQ[
Line[{v1_?VectorQ, v2_?VectorQ}],
Line[{v3_?VectorQ, v4_?VectorQ}] ] := Block[{
m = Transpose[{v2 - v1, v3 - v4}], vL},
vL = LeastSquares[ m, v3 - v1];
(*
Print["Dist Between Lines = ",  Norm[ m . vL - (v3 - v1)]];
Print["Numerical Lambdas = ", N[vL]]; *)
vEig = Sort[ Eigenvalues[ Transpose[ m ] . m]];
Norm[ m . vL - (v3 - v1)] <=
Norm[ {v1, v2, v3, v4}]*vEig[[2]]/vEig[[1]]/10^14 &&
0 <= vL[[1]] <= 1 && 0 <= vL[[2]] <= 1 ]

• If you're absolutely concerned about numerical stability, then Sort[Eigenvalues[Transpose[ m ] . m]] is totally the wrong way to do it (see e.g. this). Instead, do vEig = SingularValueList[m] and then you can use the check Norm[m . vL - (v3 - v1)] <= Norm[{v1, v2, v3, v4}]*(vEig[[1]]/vEig[[2]])^2/1*^14 (tho you might want to go for an adjustable tolerance instead). Commented Jun 11, 2022 at 16:48

In 2D (and only there), you can use the undocumented function GraphicsMeshFindIntersections to find the point of intersection. This is often an order of magnitude faster than RegionIntersection. If there isn't any intersection point, then an empty list is returned.

plot[l1_, l2_] := Graphics[{Thick, PointSize[0.025], Darker@Green, Point @@ l1, l1, Darker@Red, l2, Point @@ l2}];

l1 = Line@RandomReal[{-1, 1}, {2, 2}];
l2 = Line@RandomReal[{-1, 1}, {2, 2}];
pt = GraphicsMeshFindIntersections[{l1, l2}];
Graphics[{l1, l2, Red, Point[pt]}]


So something like

LineIntersectionQ2[l1_, l2_] := Length[GraphicsMeshFindIntersections[{l1, l2}]] > 0


could be quite a lot faster on average.

Apparently, this finds only intersection points if they are isolated. The following examples shall illustrate this:

plot[l1_, l2_] := Graphics[{Thick, PointSize[0.025],
Darker@Green, Point @@ l1, l1,
Darker@Blue, l2, Point @@ l2,
Red, Point[GraphicsMeshFindIntersections[{l1, l2}]]
},
Frame -> True];

l1 = Line[{{-1., 0}, {1., 0.}}];
l2 = Line[{{0., 1.}, {0., -1}}];
plot[l1, l2]

l1 = Line[{{-1., 0}, {1., 0.}}];
l2 = Line[{{1., 2.}, {1., 0.}}];
plot[l1, l2]

l1 = Line[{{-1., 0}, {1., 0.}}];
l2 = Line[{{0., 2.}, {0., 0.}}];
plot[l1, l2]

l1 = Line[{{0., 0.}, {2., 0.}}];
l2 = Line[{{1., 0.}, {3., 0.}}];
plot[l1, l2]

l1 = Line[{{0., 0.}, {2., 0.}}];
l2 = Line[{{2., 0.}, {4., 0.}}];
plot[l1, l2]

l1 = Line[{{0., 0.}, {2., 0.}}];
l2 = Line[{{3., 0.}, {5., 0.}}];
plot[l1, l2]

l1 = Line[{{0., 1.}, {2., 1.}}];
l2 = Line[{{1., 0.}, {3., 0.}}];
plot[l1, l2]


I am not sure whether this is what you are looking for. You say that you want to know when the line segments intersect at a point. Do you happen to mean intersect at a single point? In that case GraphicsMeshFindIntersections seems to be your friend.

Anyways, GraphicsMeshFindIntersections is still super slow. It can only test a few thousand line pairs per second. I recently wrote an C++ implementation of the Gilbert–Johnson–Keerthi distance algorithm. That one can test up to 20 million line segment pairs per second. Just to give you an idea how far one can push this.

• Would you please elaborate upon the non-intersecting parallel cases where the two line segments could be in a horizontal line, a vertical line, or at a slope. Also the intersecting parallel case, where one line segment is a 'subset' of a larger line segment. i.e. partially or fully intersecting case. Thanks in advance.
– Syed
Commented Sep 19, 2021 at 8:22
• I added a couple of examples in my answer. Feel free to test a couple of further configurations of your interest. The slopes of the lines themselfes should not matter. (Developers at Wolfram are not dumb.) Commented Sep 19, 2021 at 8:44
• My own solution doesn't work satisfactorily if there are parallel lines intersecting at multiple intersecting points, so I wanted to ask you if a better/alternate solution existed for such cases. You seem like the experienced person to whose authority I could defer to. I never implied that the developers were dumb.
– Syed
Commented Sep 19, 2021 at 9:07
• "I never implied that the developers were dumb." Sure. And I never implied that you meant to imply that! ;) Commented Sep 19, 2021 at 9:08
• We may be dumb but we're not deaf.... Commented Sep 19, 2021 at 16:21
t := RandomReal[{-10, 10}, 2]


ln1 and ln2 will be the coordinates for a Line or InfiniteLine as needed later.

ln1 = {t, t};
ln2 = {t, t};

Clear[sol, x, y]
sol = Solve[{x, y} \[Element]
InfiniteLine@ln1 \[And] {x, y} \[Element] InfiniteLine@ln2, {x,
y}]

{x, y} = {x, y} /. sol[[1]]


Let's say this gives us the following pt of intersection for infinite lines through ln1 and ln2.

{5.77478, -1.71834}

The point of intersection has to be on both line segments.

{x, y} \[Element] Line@ln1


False

{x, y} \[Element] Line@ln2


True

Graphics[{Red, Line[ln1], Red, Line[ln2], Black, Point@Midpoint@ln1,
Point@Midpoint@ln2, Blue, PointSize[0.05], Point[{x, y}]
}]


Visualization for this this particular non-intersecting case is as follows (The smaller dots are just the midpoints of the line segments):

Visualization for an intersecting case:

How about negating the output of RegionDisjoint?

Borrowing a specific example from the ref page:

SeedRandom[1];
l1 = Line[{{0, 0}, {1, 1}}];
l2 = Line[RandomReal[1, {10, 2}]];

!RegionDisjoint[l1, l2] // AbsoluteTiming

{0.002878, False}


This answer provides results of timing tests of the other answers, at least the ones that provided compatible MMA code. The test methods are shown below. With a smaller number indicating a faster method of determining whether the line segments intersect, here is a summary of the latest results: $$\begin{array}{cc} \text{Original} & 11.5 \\ \text{Henrik} & 3.019 \\ \text{Greg} & 5.3 \\ \text{Syed} & 9.990 \\ \text{irchans} & 0.16 \\ \text{Eric} & 6.99 \\ \text{Lukasz} & 0.12 \\ \end{array}$$

### Test Procedure

First, generate a test set of line pairs from a pool of candidate lines.

ClearAll["Global*"]
SeedRandom[123]

nlines = 90;
pool = Line /@ RandomReal[{-100, 100}, {nlines, 2, 2}];
set1 = Subsets[pool, {2}];
tmax = 10  (* time in seconds for repeated timings *);

\$Version  (*12.1.0 for Linux x86 (64-bit) (March 14, 2020)*)
Graphics[pool, Frame -> True, ImageSize -> Small, AspectRatio -> 1]


### Test functions

(*  Original function *)
ClearAll[LineIntersectionQ]
LineIntersectionQ[l1_, l2_] :=
RegionDimension@RegionIntersection[l1, l2] == 0

n0 = "Original";
{t0, r0} = RepeatedTiming[LineIntersectionQ @@@ set1, tmax];

ClearAll[hsFunc]
hsFunc[l1_, l2_] :=
Length[GraphicsMeshFindIntersections[{l1, l2}]] > 0

n1 = "Henrik";
{t1, r1} = RepeatedTiming[hsFunc @@@ set1, tmax];

ClearAll[ghFunc]
ghFunc[l1_, l2_] := Not@RegionDisjoint[l1, l2]

n2 = "Greg";
{t2, r2} = RepeatedTiming[ghFunc @@@ set1, tmax];

ClearAll[x, y, sFunc]
sFunc[ln1_,
ln2_] := ({x, y} \[Element] ln1 \[And] {x, y} \[Element] ln2) /.
First@
Solve[{x, y} \[Element]
InfiniteLine @@ ln1 \[And] {x, y} \[Element]
InfiniteLine @@ ln2, {x, y}]

n3 = "Syed";
{t3, r3} = RepeatedTiming[sFunc @@@ set1, tmax];

ClearAll[interQ];
interQ[Line[{v1_?VectorQ, v2_?VectorQ}],
Line[{v3_?VectorQ, v4_?VectorQ}]] :=
Block[{m = Transpose[{v1 - v2, v4 - v3}], vL},
If[Det[m] == 0, False, vL = LeastSquares[m, v4 - v2];
0 <= vL[[1]] <= 1 && 0 <= vL[[2]] <= 1]]

n4 = "irchans";
{t4, r4} = RepeatedTiming[interQ @@@ set1, tmax];

ClearAll[intersect]
intersect[Line[{a_, b_}], Line[{c_, d_}]] :=
Module[{cap},
cap = Solve[{a t + b (1 - t) == c s + d (1 - s), 0 <= t <= 1,
0 <= s <= 1}, {s, t}, Reals];
And[Length[cap] == 1, NumberQ[s /. cap[[1]]], NumberQ[t /. cap[[1]]]]]

n5 = "Eric";
{t5, r5} = RepeatedTiming[intersect @@@ set1, tmax] // Quiet;

ClearAll[interQ6]
interQ6[Line@{u1_?VectorQ, u2_?VectorQ},
Line@{v1_?VectorQ, v2_?VectorQ}] := Block[{},
And[Det[{u2 - u1, v1 - u2}] Det[{u2 - u1, v2 - u2}] < 0,
Det[{v2 - v1, u1 - v2}] Det[{v2 - v1, u2 - v2}] < 0]]

n6 = "Lukasz";
{t6, r6} = RepeatedTiming[interQ6 @@@ set1, tmax] // Quiet;


### Test results comparison

r0 == r1 == r2 == r3 == r4 == r5 == r6 (* True *)
Transpose[{{n0, n1, n2, n3, n4, n5, n6},
{t0, t1, t2, t3, t4, t5, t6}}] // Grid  (* above *)


For line segments in 2D, you may try also this one:

interQ[{u1_?VectorQ, u2_?VectorQ}, {v1_?VectorQ, v2_?VectorQ}] :=
Block[{},
And[
Det[{u2 - u1, v1 - u2}] Det[{u2 - u1, v2 - u2}] < 0,
Det[{v2 - v1, u1 - v2}] Det[{v2 - v1, u2 - v2}] < 0
]
];

• Nice answer. Thanks for contributing it. I added your solution to my comparison answer. I hope you become a regular contributor to this forum. Commented Jun 12, 2022 at 4:53
• The line (AB) intersects the inner part of a non-collinear segment [CD] iff (AB × AC)(AB × AD) < 0 (here × stands for the 2D vector product — which are equal to determinants).
• Moreover, the non-collinear segments [AB] and [BC] intersect iff the line (AB) intersects [CD], and the line (CD) intersects [AB].
• In non-collinear case, the same holds when one allows intersection at endpoints — provided the length of the segments is not 0.
• If one wishes to treat the collinear case too, one needs to special-case the situation when all 4 vector products are 0.
• Similarly when the length may be 0: then the corresponding pair of vector products is 0, so one needs only to special-case this situation.

Anyway, this should be faster by tons than the alternatives.

• In fact, the cases with length 0 are going to give one of the products positive in the first two lines — unless collinear. So only the case of all the products being 0 should be special-cased. (And then one can check the projections ontp coordinate axes: min(max(A₁,B₁),max(C₁,D₁)) ≥ max(min(A₁,B₁),min(C₁,D₁)); same for the 2nd coordinate.) Commented Sep 20, 2021 at 4:34
• I would like very much to add your method to my answer comparing the various methods. Can you include a function, like the one in the original post, that takes to Lines and returns True or False? Whether you method is fastest or not, the comparison may be a benefit. Commented Sep 20, 2021 at 21:10
• I do not program Mathematica. Commented Sep 21, 2021 at 0:19

We can parametrize the first line segment on the parameter interval $$t \in [0,1]$$ and the second on the parameter interval $$s \in [0,1]$$, then use Solve to find the intersection, restricting to those two intervals. If there is no intersection, Solve returns an empty list, so the Length of the result is zero. If there is a single point of intersection, Solve returns a list containing a list of Rules giving values to the two parameters. If there is more than one point of intersection, Solve returns a list containing a list containing a single ConditionalExpression that parametricly describes the resulting segment.

intersect[Line[{a_, b_}], Line[{c_, d_}]] :=
Module[{cap},
cap = Solve[{a t + b (1 - t) == c s + d (1 - s),
0 <= t <= 1, 0 <= s <= 1}, {s, t}, Reals];
And[Length[cap] == 1,
NumberQ[s /. cap[[1]]],
NumberQ[t /. cap[[1]]]]
]

RepeatedTiming[intersect[Line[{{-1, 0}, {1, 0}}], Line[{{0, 1}, {0, -1}}]]]
(*  {0.00048, True}  *)

RepeatedTiming[intersect[Line[{{-1, 0}, {1, 0}}], Line[{{1, 2}, {1, 0}}]]]
(*  {0.00047, True}  *)

RepeatedTiming[intersect[Line[{{-1, 0}, {1, 0}}], Line[{{0, 2}, {0, 0}}]]]
(*  {0.00047, True}  *)

RepeatedTiming[intersect[Line[{{0, 0}, {2, 0}}], Line[{{1, 0}, {3, 0}}]]]
(*  {0.002, False}  *)

RepeatedTiming[intersect[Line[{{0, 0}, {2, 0}}], Line[{{2, 0}, {4, 0}}]]]
(*  {0.002, True}  *)

RepeatedTiming[intersect[Line[{{0, 0}, {2, 0}}], Line[{{3, 0}, {5, 0}}]]]
(*  {0.001, False}  *)

RepeatedTiming[intersect[Line[{{0, 1}, {2, 1}}], Line[{{1, 0}, {3, 0}}]]]
(*  {0.000037, False}  *)

R1[p, q] := m p + (1 - m) q;
R2[p, q] := n p + (1 - n) q;
Test[pts] :=
And[(0 <= m <= 1), (0 <= n <= 1)] /.
Solve[R1[pts[[1]], pts[[2]]] == R2[pts[[3]], pts[[4]]], {m, n}];
Draw[pts] :=
Graphics[{Line[{pts[[1]], pts[[2]]}], Line[{pts[[3]], pts[[4]]}]}];
{Test[#][[1]], Draw[#]} &@RandomReal[{-1, 1}, {4, 2}] // Timing


• I would like very much to add your method to my answer comparing the various methods. Can you include a function, like the one in the original post, that takes to Line`s and returns True or False? Whether you method is fastest or not, the comparison may be a benefit. Commented Sep 20, 2021 at 21:08