14
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Related mathematica.stackexchange.com/q/51391/58370 $\endgroup$
    – Rohit Namjoshi
    Sep 19, 2021 at 1:37
  • $\begingroup$ 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? $\endgroup$
    – Henrik Schumacher
    Sep 19, 2021 at 8:25
  • 1
    $\begingroup$ @HenrikSchumacher I asked for general dimensions, but I'd be happy with a more efficient solution for either the 2D or 3D cases. Thanks! $\endgroup$
    – sam wolfe
    Sep 19, 2021 at 11:30

9 Answers 9

Reset to default
17
$\begingroup$

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 ]
Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ 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). $\endgroup$
    – J. M.'s lack of A.I. ♦
    Jun 11, 2022 at 16:48
  • $\begingroup$ @J.M. Cool! Let me think about this! :) $\endgroup$
    – irchans
    Jun 11, 2022 at 17:05
14
$\begingroup$

In 2D (and only there), you can use the undocumented function Graphics`Mesh`FindIntersections 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 = Graphics`Mesh`FindIntersections[{l1, l2}];
Graphics[{l1, l2, Red, Point[pt]}]

So something like

LineIntersectionQ2[l1_, l2_] := Length[Graphics`Mesh`FindIntersections[{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[Graphics`Mesh`FindIntersections[{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 Graphics`Mesh`FindIntersections seems to be your friend.

Anyways, Graphics`Mesh`FindIntersections 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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Syed
    Sep 19, 2021 at 8:22
  • $\begingroup$ 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.) $\endgroup$
    – Henrik Schumacher
    Sep 19, 2021 at 8:44
  • $\begingroup$ 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. $\endgroup$
    – Syed
    Sep 19, 2021 at 9:07
  • $\begingroup$ "I never implied that the developers were dumb." Sure. And I never implied that you meant to imply that! ;) $\endgroup$
    – Henrik Schumacher
    Sep 19, 2021 at 9:08
  • 4
    $\begingroup$ We may be dumb but we're not deaf.... $\endgroup$
    – Daniel Lichtblau
    Sep 19, 2021 at 16:21
11
$\begingroup$ 
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):

enter image description here

Visualization for an intersecting case:

enter image description here

Improve this answer
$\endgroup$
10
$\begingroup$

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}
Improve this answer
$\endgroup$
10
$\begingroup$

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]

enter image description here

Test functions

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

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

(* Henrik Schumacher's answer *)
ClearAll[hsFunc]
hsFunc[l1_, l2_] := 
 Length[Graphics`Mesh`FindIntersections[{l1, l2}]] > 0

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

(* Greg Hurst's answer *)
ClearAll[ghFunc]
ghFunc[l1_, l2_] := Not@RegionDisjoint[l1, l2]

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

(* Syed's answer *)
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];

(*  irchans' answer  *)
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];

(*  Eric Tower's answer  *)
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;

(* Lukasz's answer *)
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 *)
Improve this answer
$\endgroup$
6
$\begingroup$

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
    ]
   ];
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Nice answer. Thanks for contributing it. I added your solution to my comparison answer. I hope you become a regular contributor to this forum. $\endgroup$
    – LouisB
    Jun 12, 2022 at 4:53
4
$\begingroup$
  • 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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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.) $\endgroup$
    – Ilya Zakharevich
    Sep 20, 2021 at 4:34
  • 1
    $\begingroup$ 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. $\endgroup$
    – LouisB
    Sep 20, 2021 at 21:10
  • $\begingroup$ I do not program Mathematica. $\endgroup$
    – Ilya Zakharevich
    Sep 21, 2021 at 0:19
3
$\begingroup$

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}  *)
Improve this answer
$\endgroup$
3
$\begingroup$ 
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

enter image description here enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – LouisB
    Sep 20, 2021 at 21:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.