# Find intersection of pairs of straight lines

I have a list of 24 points, in which two consecutive points (1st and 2nd, 3rd and 4th, …) are supposed to form a line:

p1={{243.8, 77.}, {467.4, 12.}, {291.8, 130.}, {476., 210.5}, {103.2,
327.}, {245.2, 110.5}, {47.4, 343.}, {87.4, 108.5}, {371.,
506.5}, {384.6, 277.}, {264.6, 525.5}, {353.8, 294.5}, {113.2,
484.5}, {296., 304.5}, {459.6, 604.5}, {320.2, 466.5}, {288.2,
630.5}, {199.6, 446.5}, {138.8, 615.5}, {81.8, 410.}, {232.4,
795.}, {461.8, 727.}, {27.4, 671.5}, {206.8, 763.5}};


I also have another list of 24 points with the same property:

p2={{356.8, 32.}, {363.2, 120.}, {346., 245.}, {393.8, 158.}, {163.8,
211.5}, {230.2, 250.}, {54.6, 225.}, {139.6, 220.}, {366.,
394.5}, {451.8, 372.}, {241., 398.}, {321., 411.5}, {163.2,
347.}, {213.2, 406.5}, {332.4, 596.5}, {402.4, 528.5}, {176.,
585.5}, {256., 530.5}, {38.2, 553.}, {122.4, 507.}, {345.2,
774.5}, {345.2, 688.}, {104.6, 728.}, {161.8, 647.}};


My goal is to find the intersections between a line in p1 and a line in p2 as shown in the graph below. I really don't know how to start with this, and worse, a line on p1 does not match up with the intersecting line in p2 in terms of order in each list. This could be observed by different colors of 2 intersecting lines, and makes it harder for element-by-element manipulation*. How can I solve this?

Join[Partition[p1, 2], Partition[p2, 2]] // ListLinePlot


*I found out that this is thankfully not true, as seen below when line i in p1 and line i in p2 are plotted together, and also by Öskå in a comment to eldo's answer.

Row@Table[
ListLinePlot[{Partition[p1, 2][[i]], Partition[p2, 2][[i]]}], {i, 1,
Length@Partition[p1, 2]}]


• In Wolfram Programming Cloud, and supposedly future Mathematica 10, you can find intersections using geometric computation: {x, y} /. Solve[Apply[And, Element[{x, y}, Line[Partition[#, 2]]] & /@ {p1, p2}], {x, y}] – kirma Jun 24 '14 at 11:19
• @kirma Nice work! – Apple Jun 24 '14 at 11:52

p1 = Partition[{{243.8, 77.}, {467.4, 12.}, {291.8, 130.}, {476.,
210.5}, {103.2, 327.}, {245.2, 110.5}, {47.4, 343.}, {87.4,
108.5}, {371., 506.5}, {384.6, 277.}, {264.6, 525.5}, {353.8,
294.5}, {113.2, 484.5}, {296., 304.5}, {459.6, 604.5}, {320.2,
466.5}, {288.2, 630.5}, {199.6, 446.5}, {138.8, 615.5}, {81.8,
410.}, {232.4, 795.}, {461.8, 727.}, {27.4, 671.5}, {206.8,
763.5}}, 2];

p2 = Partition[{{356.8, 32.}, {363.2, 120.}, {346., 245.}, {393.8,
158.}, {163.8, 211.5}, {230.2, 250.}, {54.6, 225.}, {139.6,
220.}, {366., 394.5}, {451.8, 372.}, {241., 398.}, {321.,
411.5}, {163.2, 347.}, {213.2, 406.5}, {332.4, 596.5}, {402.4,
528.5}, {176., 585.5}, {256., 530.5}, {38.2, 553.}, {122.4,
507.}, {345.2, 774.5}, {345.2, 688.}, {104.6, 728.}, {161.8,
647.}}, 2];

LineIntersectionPoint[{a_, b_}, {c_, d_}] :=
(Det[{a, b}] (c - d) - Det[{c, d}] (a - b))/Det[{a - b, c - d}]

Graphics[{Line /@ {p1, p2}, Red, PointSize@.05,
Point /@ MapThread[LineIntersectionPoint, {p1, p2}]}, Frame -> True]


Ref for finding intersection of 2 lines by determinants

• @Öskå - thanks for your elegant edit. With the too complicated "line production" I wanted to verify for myself that the lines cross pairwise in ascending order, because the OP suggested that "line on p1 does not match up with the intersecting line in p2 in terms of order in each list." Probably he was distracted by the coloring of ListLinePlot. – eldo Jun 24 '14 at 10:57
• I guess so :) Using Manipulate[ Graphics[{Line /@ {p1, p2}, Red, Line /@ {p1[[i]], p2[[i]]}}], {i, 1, Length@p1, 1}] helps to see that intersect with each others :) – Öskå Jun 24 '14 at 11:01
• nice with and without edit +1 – ubpdqn Jun 24 '14 at 11:02
• I don't believe this is sufficient if we believe the OP when he writes "... and worse, a line on p1 does not match up with the intersecting line in p2 in terms of order in each list." – Mr.Wizard Jun 24 '14 at 17:20
• I edited your code to make it somewhat cleaner. You could also use LineIntersectionPoint[p_, q_] := (Det[p] #2 - Det[q] #)/Det[{#, #2}] & @@ Subtract @@@ {p, q} which is more terse and eliminates redundant subtraction, but it is of a significantly different style so I didn't substitute it. – Mr.Wizard Jun 24 '14 at 21:00

Turning my comment into an answer per (now deleted?) comment which requested it.

This is documented to work only in Wolfram Language at this point (specifically Wolfram Programming Cloud). Interestingly enough, it does work also with Mathematica 9.0.1., although documentation has no indication of Line or Solve supporting geometric regions.

p1 = {{243.8, 77.}, {467.4, 12.}, {291.8, 130.}, {476.,
210.5}, {103.2, 327.}, {245.2, 110.5}, {47.4, 343.}, {87.4,
108.5}, {371., 506.5}, {384.6, 277.}, {264.6, 525.5}, {353.8,
294.5}, {113.2, 484.5}, {296., 304.5}, {459.6, 604.5}, {320.2,
466.5}, {288.2, 630.5}, {199.6, 446.5}, {138.8, 615.5}, {81.8,
410.}, {232.4, 795.}, {461.8, 727.}, {27.4, 671.5}, {206.8,
763.5}};

p2 = {{356.8, 32.}, {363.2, 120.}, {346., 245.}, {393.8,
158.}, {163.8, 211.5}, {230.2, 250.}, {54.6, 225.}, {139.6,
220.}, {366., 394.5}, {451.8, 372.}, {241., 398.}, {321.,
411.5}, {163.2, 347.}, {213.2, 406.5}, {332.4, 596.5}, {402.4,
528.5}, {176., 585.5}, {256., 530.5}, {38.2, 553.}, {122.4,
507.}, {345.2, 774.5}, {345.2, 688.}, {104.6, 728.}, {161.8,
647.}};

(* Convert coordinate-lists to two collections of lines which can be used as
primitives in both in graphics and new geometric computation. *)
{lines1, lines2} = Line[Partition[#, 2]]& /@ {p1, p2};

(* Create points which belong to both geometric regions
consisting of line collections, that is any intersections. *)
points = Point[{x, y}] /. Solve[{x, y} \[Element] lines1 &&
{x, y} \[Element] lines2, {x, y}];

(* Represent all these as Graphics. *)
Graphics[{Blue, lines1, Red, lines2,
Black, PointSize[Large], points}, Frame->True]


EDIT:

You can also write above Solvein v10 as:

Solve[{x, y} \[Element] RegionIntersection[lines1, lines2], {x, y}]


This gets interesting when you consider the fact these regions can be much more than lines, for instance circles, filled regions such as disks, implicit and parametric regions, and derived regions. Also in higher dimensions, and symbolically. And they can be discretized, among other things for use of FEM in v10.

• It can work in MMA 9.0.1 and your comment still exist. – Apple Jun 24 '14 at 14:40
• @Chenminqi This is interesting. Line definitely has no documentation hinting of that. Is it undocumented beginnings of geometric region computation? – kirma Jun 24 '14 at 14:44
• @kirma exciting !!! – eldo Jun 24 '14 at 14:56
• Whoa - how did you come across this one? – Yves Klett Jun 24 '14 at 14:57
• @YvesKlett See "Geometric computation" in Wolfram Language documentation. I have been playing with it on Programming Cloud which was unveiled yesterday. Finding out it works with 9.0.1 is entirely thanks to Chenminqi. :) – kirma Jun 24 '14 at 15:07

Here is a direct vector calculation that verifies the segments (not just the infinite lines) intersect.

 segsegintersection[ lines_ ] := Module[{
md = Subtract @@ (Plus @@ # & /@ lines),
sub = Subtract @@ # & /@ lines, det},
det = -Det[sub];
If[And @@ (Abs[#] <= 1 & /@ #) ,
(Plus @@ #[[1]] - Subtract @@ #[[1]] Last@#[[2]])/2 & @
{First@lines, # }, False] &@
(Det[{#[[1]], md}]/det & /@ ( {#, Reverse@#}  &@ sub))];


in the example provided they all intersect.. but I thought it useful to included here for completeness. This is way faster than using Solve with constraints. Note @eldo's LineIntersectionPoint is faster than this by a factor of 2 if you do not need the intersection check.

  Graphics[ {Line /@ p1  , Line /@ p2 , Red, PointSize[.025],
Point@ MapThread[segsegintersection[{ #1 , #2 }] & , {p1, p2} ]}]


same plot as the others..

An example with only some intersections:

 lines = RandomReal[{-1, 1}, {20, 2, 2}];
Graphics[{Line@lines, Red, PointSize[.02],
Point@Select[ segsegintersection[#] & /@
Subsets[lines, {2}] ,  # =!= False &]}]


As I noted in this answer, there is a built-in, but undocumented function that can do this:

p1 = Partition[{{243.8, 77.}, {467.4, 12.}, {291.8, 130.}, {476., 210.5}, {103.2, 327.},
{245.2, 110.5}, {47.4, 343.}, {87.4, 108.5}, {371., 506.5}, {384.6, 277.},
{264.6, 525.5}, {353.8, 294.5}, {113.2, 484.5}, {296., 304.5},
{459.6, 604.5}, {320.2, 466.5}, {288.2, 630.5}, {199.6, 446.5},
{138.8, 615.5}, {81.8, 410.}, {232.4, 795.}, {461.8, 727.}, {27.4, 671.5},
{206.8, 763.5}}, 2];

p2 = Partition[{{356.8, 32.}, {363.2, 120.}, {346., 245.}, {393.8, 158.}, {163.8, 211.5},
{230.2, 250.}, {54.6, 225.}, {139.6, 220.}, {366., 394.5}, {451.8, 372.},
{241., 398.}, {321., 411.5}, {163.2, 347.}, {213.2, 406.5}, {332.4, 596.5},
{402.4, 528.5}, {176., 585.5}, {256., 530.5}, {38.2, 553.}, {122.4, 507.},
{345.2, 774.5}, {345.2, 688.}, {104.6, 728.}, {161.8, 647.}}, 2];

GraphicsMeshMeshInit[];
pts = Flatten[MapThread[FindIntersections[{Line[#1], Line[#2]}] &, {p1, p2}], 1];

Graphics[{{Red, Line[p1]}, {Blue, Line[p2]}, {AbsolutePointSize[5], Point[pts]}}]


The method in eldo's answer will work only in two dimensions. Additionally, due to its determinantal form (Cramer's rule), it can suffer from instability for certain special configurations. Thus, it is profitable to recast the formula used as an explicit linear least squares problem that can be dealt with by LeastSquares[], as well as reformulating it so that it makes no assumptions on the dimensions of the given endpoints. Thus:

lineIntersection[Line[{p1_, p2_}], Line[{q1_, q2_}]] :=
With[{sol = LeastSquares[Transpose[{p2 - p1, q2 - q1}], q1 - p1]},
{p1 + sol[[1]] (p2 - p1), q1 - sol[[2]] (q2 - q1)}]


This returns a list of two points. In two dimensions, the two points are identical, and this is the unique point of intersection:

pts == MapThread[lineIntersection[Line[#1], Line[#2]][[1]] &, {p1, p2}]
True


In three dimensions (and beyond!), if the two lines do intersect, the two points should be identical. If the two lines are in fact skew lines, the two points returned will be different, and these are the points of closest approach in each of the two lines.

(N.B. The previous version of this answer used a formulation equivalent to the normal equations, which can also fail in some cases. The current version should now be quite robust.)

• It should probably be noted that in fact, lineIntersection[] does not at all depend on the dimensions of the points contained within the Line[] objects given to it. – J. M. is away Jul 9 '15 at 16:55

This is not as neat as eldo but post as another (rather hamfisted) approach:

lin[p_, q_] := p + (q - p) s
linv[p_, q_, v_] := p + (q - p) v
eqn[list_, v_] := linv[##, v] & @@@ Partition[list, 2]


Visualizing:

all = Show[ParametricPlot[lin @@@ Partition[p2, 2], {s, 0, 1}],
ParametricPlot[lin @@@ Partition[p1, 2], {s, 0, 1},
PlotStyle -> Red]]


Finding points of intersection:

soln = {#1, Solve[#1 == #2, {s, t}]} & @@@
Tuples[{eqn[p1, s], eqn[p2, t]}];
pts = #1 /. #2[[1]] & @@@
Select[soln,
And[0 < #[[2, 1, 1, 2]] < 1, 0 < #[[2, 1, 2, 2]] < 1] &];


The points:

{{357.665, 43.8996}, {386.456, 171.367}, {174.78, 217.866}, {67.659,
224.232}, {377.821, 391.4}, {309.378, 409.539}, {203.788,
395.3}, {392.522, 538.096}, {244.015, 538.739}, {110.507,
513.497}, {345.2, 761.563}, {113.367, 715.586}}


Visualizing:

int = Show[all, Graphics[{PointSize[0.02], Point[pts]}]]


Since it only has 24 lines, we can use brute force to list all the possible intersections.

 test[{{pp1_, pp2_}, {pp3_, pp4_}}] :=
Module[{a, b, x1, y1, x2, y2, sol},
{x1, y1} = pp1 + a (pp2 - pp1); {x2, y2} = pp3 + b (pp4 - pp3);
sol = Solve[{x1 == x2, y1 == y2, 0 <= a <= 1, 0 <= b <= 1}, {a, b}];
If[sol != {}, pp1 + a (pp2 - pp1) /. sol, {}]
]
points = test /@ Tuples[{Partition[p1, 2], Partition[p2, 2]}] // Flatten[#, 1] &;
pic = ListLinePlot[Join[Partition[p1, 2], Partition[p2, 2]],
Axes -> None, PlotStyle -> Black, ImageSize -> {500, 300},
Epilog -> {Red, PointSize[0.015], Point[points]}]


This works for any lines, including vertical lines, but does not check for the cases where there are no or infinitely many intersections:

LineIntersectionPoint[{p1_, p2_}, {q1_, q2_}] := {\[FormalX], \[FormalY]} /.
First@Solve[
Cross[p2 - p1].({\[FormalX], \[FormalY]} - p1) ==
Cross[q2 - q1].({\[FormalX], \[FormalY]} - q1) == 0, {\[FormalX], \[FormalY]}]


or equivalently:

LineIntersectionPoint[{p1_, p2_}, {q1_, q2_}] := {\[FormalX], \[FormalY]} /.
First@FindInstance[
Cross[p2 - p1].({\[FormalX], \[FormalY]} - p1) ==
Cross[q2 - q1].({\[FormalX], \[FormalY]} - q1) == 0, {\[FormalX], \[FormalY]}]