Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a list of 24 points, in which 2 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 different 24 points with the same behavior.

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

share|improve this question
5  
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 at 11:19
    
@kirma Nice work! –  Chenminqi Jun 24 at 11:52

5 Answers 5

up vote 12 down vote accepted
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]

Mathematica graphics

Ref for finding intersection of 2 lines by determinants

line graph

determinant intersection

share|improve this answer
    
I shortened your code since there was useless parts :) Revert it if you mind :) –  Öskå Jun 24 at 10:50
    
The less lines there is, the better it is :) –  Öskå Jun 24 at 10:51
    
@Ö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 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 at 11:01
1  
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 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]

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.

share|improve this answer
2  
It can work in MMA 9.0.1 and your comment still exist. –  Chenminqi Jun 24 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 at 14:44
    
@kirma exciting !!! –  eldo Jun 24 at 14:56
    
Whoa - how did you come across this one? –  Yves Klett Jun 24 at 14:57
1  
@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 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 &]}]

enter image description here

share|improve this answer

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

 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

share|improve this answer

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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