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'm trying to implement a Brillouin Zone algorithm within Mathematica, including the generation of Brillouin zones of higher order in 2D and 3D. There is a nice implementation of generating these zones in the Mathematica Guidebook for Graphics. However, this implementation uses the brute force approach in calculating line segment intersections of order $\mathcal O(n^2)$ with $n$ number of lines:

intersectionPoint[{{p1x_, p1y_}, {r1x_, r1y_}}, {{p2x_, p2y_}, {r2x_, r2y_}}, maxDist_] :=
Module[{aux},
If[PossibleZeroQ[r1y r2x - r2y r1x], Sequence @@ {},
 aux = {p1x + (r1x (p1y r2x - p2y r2x -
                 p1x r2y + p2x r2y))/(r1x r2y - r1y r2x),
           p1y + (r1y (p1y r2x - p2y r2x -
                 p1x r2y + p2x r2y))/(r1x r2y - r1y r2x)} // N;
 If[Simplify[aux.aux] <= maxDist && 
 IntervalMemberQ[
  IntervalIntersection[Interval[{p1x, p1x + r1x}], 
   Interval[{p2x, p2x + r2x}]], aux[[1]]] && 
 IntervalMemberQ[
  IntervalIntersection[Interval[{p1y, p1y + r1y}], 
   Interval[{p2y, p2y + r2y}]], aux[[2]]], aux, Sequence @@ {}]]]

I'm sure there are a lot of improvements to above code, but the fact persists that the order will be of $\mathcal O(n^2)$. Since the line intersection for the Brillouin zone algorithm for high order zones in 3D is by far the most costly step, I'm looking into smarter approaches to find intersecting line segments. The very smart algorithm by Balaban entitled An optimal algorithm for finding segments intersections achieves an order of at least $\mathcal O(n\log\,n)$. However, the algorithm involves a rather complex binary search tree implementation. Since Mathematica has very efficient search implementations within lists already natively supported, I wonder if someone has implemented the Balaban or equivalent sweep line and sweep plane algorithms within Mathematica that takes advantage of the built-in Mathematica search functions? I'm especially interested in a 2D and a 3D implementation of a line segment intersection within Mathematica.

Many thanks in advance for any help!

share|improve this question

migrated from stackoverflow.com Oct 26 '12 at 23:59

This question came from our site for professional and enthusiast programmers.

    
Did you check the link? It may help ... mathematica.stackexchange.com/questions/528/… –  s.s.o Oct 27 '12 at 0:43
add comment

2 Answers 2

up vote 16 down vote accepted

You'll be interested in the (undocumented!) functions Graphics`Mesh`IntersectQ[] (for checking the intersections) and Graphics`Mesh`FindIntersections[] (for actually finding them). As a sample:

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
            lins = Table[{Line[RandomReal[1, {2, 2}]]}, {42}];]

Graphics`Mesh`MeshInit[];

pts = FindIntersections[lins]; (* intersection points *)
Graphics[{{AbsoluteThickness[1], lins}, {Directive[Red, AbsolutePointSize[4]], Point[pts]}}]

lines with marked intersection points

share|improve this answer
    
Great info J.M.! I tested the algorithm against my brute force one and got an speedup of about 100x for 250 segments. Do you know the internals of the algorithm? –  Rainer Oct 27 '12 at 17:11
    
The approach I presented doesn't work for the 3D case, tho. That would be slightly more gnarly to do... as for the algorithm internally used, I'm not entirely sure. Maybe somebody from Wolfram can drop a comment. –  J. M. Oct 27 '12 at 17:12
2  
Now it seems after some additional searching I found a solution to the 3D problem too. Under mathworks.com/matlabcentral/newsreader/view_thread/166800 the answer by Roger Stafford seems to do at least a significant improvement over the brute-force approach in 3D. He recommends to use the 2D line sweep algorithm on the x-y and then on the x-z plane. All true 3D intersections will be the intersection of both x-y and x-z sets within a certain error limit to account for numerical errors. I'll try this out and post the result. –  Rainer Oct 27 '12 at 20:54
add comment

Is a runtime that's $O(n \log n)$ even achievable? I feel I can create a sequence of line sets with a quadratic amount of intersections:

ReIm[z_] := Through[{Re, Im}[z]]

Graphics`Mesh`MeshInit[];
numInts[n_] := numInts[n] = BlockRandom[
    Length@Graphics`Mesh`FindIntersections[
      Most@Table[Line[{-ReIm[Exp[I t]],ReIm[Exp[I t]]+RandomReal[{0,.1}]}],{t,0,2Pi,2Pi/n}]
    ]
]

Plot[x(x + 1)/2, {x,1,1001}, Epilog -> {Red, Line[{#, numInts[#]}& /@ Range[1, 1001, 100]]}]

and the relative error to $n(n+1)/2$ seems to approach zero:

relError[n_] := 1 - 2numInt[n]/(n(n + 1))

ListLinePlot[{#, relError[#]}& /@ Range[201, 1001, 100]]

share|improve this answer
add comment

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.