25
$\begingroup$

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

Join Partition line intersection

*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]}]

line pairs

lines

Improve this question
$\endgroup$
2
  • 9
    $\begingroup$ 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}] $\endgroup$
    – kirma
    Jun 24, 2014 at 11:19
  • $\begingroup$ @kirma Nice work! $\endgroup$
    – Apple
    Jun 24, 2014 at 11:52

7 Answers 7

Reset to default
28
$\begingroup$ 
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, [email protected], 
  Point /@ MapThread[LineIntersectionPoint, {p1, p2}]}, Frame -> True]

Mathematica graphics

Ref for finding intersection of 2 lines by determinants

line graph

determinant intersection

Improve this answer
$\endgroup$
10
  • $\begingroup$ @Ö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. $\endgroup$
    – eldo
    Jun 24, 2014 at 10:57
  • $\begingroup$ 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 :) $\endgroup$
    – Öskå
    Jun 24, 2014 at 11:01
  • $\begingroup$ nice with and without edit +1 $\endgroup$
    – ubpdqn
    Jun 24, 2014 at 11:02
  • $\begingroup$ 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." $\endgroup$
    – Mr.Wizard
    Jun 24, 2014 at 17:20
  • 2
    $\begingroup$ 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. $\endgroup$
    – Mr.Wizard
    Jun 24, 2014 at 21:00
18
$\begingroup$

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]

enter image description here

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.

Improve this answer
$\endgroup$
7
  • 2
    $\begingroup$ It can work in MMA 9.0.1 and your comment still exist. $\endgroup$
    – Apple
    Jun 24, 2014 at 14:40
  • $\begingroup$ @Chenminqi This is interesting. Line definitely has no documentation hinting of that. Is it undocumented beginnings of geometric region computation? $\endgroup$
    – kirma
    Jun 24, 2014 at 14:44
  • $\begingroup$ @kirma exciting !!! $\endgroup$
    – eldo
    Jun 24, 2014 at 14:56
  • $\begingroup$ Whoa - how did you come across this one? $\endgroup$
    – Yves Klett
    Jun 24, 2014 at 14:57
  • 1
    $\begingroup$ @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. :) $\endgroup$
    – kirma
    Jun 24, 2014 at 15:07
16
$\begingroup$

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 &]}]

enter image description here

Improve this answer
$\endgroup$
9
$\begingroup$

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];

Graphics`Mesh`MeshInit[];
pts = Flatten[MapThread[FindIntersections[{Line[#1], Line[#2]}] &, {p1, p2}], 1];

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

intersecting lines

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.)

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – J. M.'s lack of A.I. ♦
    Jul 9, 2015 at 16:55
5
$\begingroup$

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]]

enter image description here

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]}]]

enter image description here

Improve this answer
$\endgroup$
5
$\begingroup$

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]}]

enter image description here

Improve this answer
$\endgroup$
1
$\begingroup$

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

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.