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?
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 ]
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$
Commented
Jun 11, 2022 at 16:48
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.
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}$$
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]
(* 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;
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
]
];
Anyway, this should be faster by tons than the alternatives.
Lines and returns True or False? Whether you method is fastest or not, the comparison may be a benefit.
$\endgroup$
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
Lines and returns True or False? Whether you method is fastest or not, the comparison may be a benefit.
$\endgroup$
I understand I'm late to the party. However, there's a way to do it several times faster than Lukasz algo. Namely, intersect line segments all at once using Graphics`Mesh`FindIntersections and then, post-select intersecting segments, from intersection points data:
xngSegmentBelong =
Compile[{{pts1, _Real, 1}, {pts2, _Real, 1}, {pt, _Real, 1}},
Module[{x0, y0, x1, y1, x2, y2, det},
{x0, y0} = pt;
{x1, y1} = pts1;
{x2, y2} = pts2;
det = (x2 - x1)*(y0 - y1) - (y2 - y1)*(x0 - x1);
Abs[det] < 1*^-8 && Min[{x1, x2}] <= x0 <= Max[{x1, x2}] &&
Min[{y1, y2}] <= y0 <= Max[{y1, y2}]
], RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed"]
poolIndexes = Subsets[Range@Length@pool, {2}];
n7 = "Anton";
{t7, blank} = RepeatedTiming[
xngs = Graphics`Mesh`FindIntersections[pool];
pt1list = Flatten[ConstantArray[pool[[All, 1]][[All, 1]], Length@xngs], 1];
pt2list = Flatten[ConstantArray[pool[[All, 1]][[All, 2]], Length@xngs], 1];
xngslist = Flatten[Transpose[ConstantArray[xngs, Length@pool]], 1];
xsb = xngSegmentBelong[pt1list, pt2list, xngslist];
crossingSegments =
Partition[Flatten@Position[xsb, True], 2] -
Range[0, Length@xngs - 1]*Length@pool;
asTrue = AssociationThread[crossingSegments -> True];
asFalse = AssociationThread[
Complement[poolIndexes, crossingSegments] -> False];
r7 = poolIndexes /. asTrue /. asFalse;, tmax];
r0 == r1 == r2 == r3 == r4 == r5 == r6 == r7 (*True*)
Transpose[{{n0, n1, n2, n3, n4, n5, n6, n7}, {t0, t1, t2, t3, t4, t5,
t6, t7}}] // Grid (*below*)
Here's the absolutely fastest way. It's about 20 times faster than my previous answer:
interQanton =
Compile[{{v1, _Real, 2}, {v2, _Real, 2}},
Module[{x1, x2, x3, x4, y1, y2, y3, y4, p1, p2, p3, p4, xx, yy, tst,
n1, n2, d1, d2, t1, t2, t3, t4, t5, w, t6, t7, xo, yo, xz, yz,
out, t8},
t2 = t3 = t4 = t5 = t6 = t7 = t8 = w = xz = yz = False;
xo = yo = 999.; xx = yy = 0.;
{p1, p2} = v1;
{p3, p4} = v2;
{x1, y1} = p1;
{x2, y2} = p2;
{x3, y3} = p3;
{x4, y4} = p4;
out = False;
If[p1 == p3 || p1 == p4, t8 = True; out = True;];
If[p2 == p3 || p2 == p4, t8 = True; out = True;];
If[! t8,
n1 = (x1 (x4 (-y2 + y3) + x3 (y2 - y4)) +
x2 (x4 (y1 - y3) + x3 (-y1 + y4)));
d1 = (-(x3 - x4) (y1 - y2) + (x1 - x2) (y3 - y4));
n2 = -((-x2 y1 + x1 y2) (y3 - y4) + (y1 - y2) (x4 y3 -
x3 y4));
d2 = ((x3 - x4) (y1 - y2) + (-x1 + x2) (y3 - y4));
If[Sign[Abs[d1] - 1*^-6] == -1, d1 = d2 = 1;
n1 = Min[{x1, x2}] - 1; w = True];
xx = n1/d1;
yy = n2/d2;
];
If[w,
If[(y3 - y1) != 0 && (y2 - y1) != 0,
t2 = ((x3 - x1)/(y3 - y1)) == ((x2 - x1)/(y2 - y1)); {xx,
yy} = {x3, y3};];
If[(y4 - y1) != 0 && (y2 - y1) != 0,
t3 = ((x4 - x1)/(y4 - y1)) == ((x2 - x1)/(y2 - y1)); {xx,
yy} = {x4, y4};];
If[(y1 - y3) != 0 && (y4 - y3) != 0,
t4 = ((x1 - x3)/(y1 - y3)) == ((x4 - x3)/(y4 - y3)); {xx,
yy} = {x1, y1};];
If[ (y2 - y3) != 0 && (y4 - y3) != 0,
t5 = ((x2 - x3)/(y2 - y3)) == ((x4 - x3)/(y4 - y3)); {xx,
yy} = {x2, y2};];
If[y1 == y2 == y3 == y4,
If[Min[{x1, x2}] <= x3 <= Max[{x1, x2}], xo = x3; xz = True;];
If[Min[{x1, x2}] <= x4 <= Max[{x1, x2}], xo = x4; xz = True;];
If[xz, t6 = True; {xx, yy} = {xo, y1};]];
If[x1 == x2 == x3 == x4,
If[Min[{y1, y2}] <= y3 <= Max[{y1, y2}], yo = y3; yz = True;];
If[Min[{y1, y2}] <= y4 <= Max[{y1, y2}], yo = y4; yz = True;];
If[yz, t7 = True; {xx, yy} = {x1, yo};]];
];
If[! t8,
t1 =
Min[{x1, x2}] <= xx <= Max[{x1, x2}] &&
Min[{y1, y2}] <= yy <= Max[{y1, y2}] &&
Min[{x3, x4}] <= xx <= Max[{x3, x4}] &&
Min[{y3, y4}] <= yy <= Max[{y3, y4}];
out = If[t1 || t2 || t3 || t4 || t5 || t6 || t7, True, False];
];
out
], RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed",
CompilationTarget -> "C"]
sgs = pool[[All, 1]];
sets = Subsets[Range@Length@sgs, {2}];
n8 = "Anton1";
{s1, s2} = {sets[[All, 1]], sets[[All, 2]]};
{t8, r8} = RepeatedTiming[interQanton[sgs[[s1]], sgs[[s2]]], tmax];
r0 == r1 == r2 == r3 == r4 == r5 == r6 == r7 == r8(*True*)
Transpose[{{n0, n1, n2, n3, n4, n5, n6, n7, n8}, {t0, t1, t2, t3, t4,
t5, t6, t7, t8}}] // Grid (*below*)
For larger segment sets this algo is also the fastest, I've checked with nlines = 900; for example.
Lineobjects in the plane (2D) only, or also in 3D or even nD? $\endgroup$