Nearest "sorts" the sticks, but because it is 2D-data it does some quadtree / KDtree technique to reduce the search time from being linear time ( $n$) to logarithmic time ( $\log(n)$) (roughly speaking).

George, Regarding use and credit of other's work, I guess there is this standard advice from mathematics.

Regarding the code, I should emphasize that most of it came right from the @George Varnavides post (all except for the part that didn't...)

George, this is very nice, and I gave it my vote. That said, there are very good reasons to make the intersection code efficient. You partly did that with your Compile version. I will illustrate why one also wants binning for examples that are larger than yours but not hugely so (say an order of magnitude each way). Specifically, I will use n=5000 and L=.05.

n = 5000;
SeedRandom[111111]
centers = RandomReal[{0, 1}, {n, 2}];
\[Theta] = 
  RandomVariate[NormalDistribution[\[Pi]/2, (1/2.3548)*\[Pi]/4], n];
L = 0.05;
stickend1 = centers + Transpose[L {Cos[\[Theta]], Sin[\[Theta]]}/2];
stickend2 = centers - Transpose[L {Cos[\[Theta]], Sin[\[Theta]]}/2];
sticks = Transpose[{stickend1, stickend2}];

lines = Map[Line, Transpose[{stickend1, stickend2}]];
Graphics[lines, Frame -> True]

I modify the compiled segment intersection code to return a 0-1 value.

lineIntComp2 = 
  Compile[{{u, _Real, 2}, {v, _Real, 2}}, 
   Block[{a = u[[1]], b = u[[2]], c = v[[1]], d = v[[2]], pt, 
     detBottom}, 
    detBottom = 
     u[[2, 2]] (v[[1, 1]] - v[[2, 1]]) + 
      u[[1, 2]] (-v[[1, 1]] + v[[2, 1]]) + (u[[1, 1]] - 
         u[[2, 1]]) (v[[1, 2]] - v[[2, 2]]);
    If[Chop[detBottom] == 0, False, 
     pt = ({(-u[[1, 2]] u[[2, 1]] + u[[1, 1]] u[[2, 2]]) (v[[1, 1]] - 
             v[[2, 1]]) + (u[[1, 1]] - 
             u[[2, 1]]) (v[[1, 2]] v[[2, 1]] - 
             v[[1, 1]] v[[2, 2]]), (-u[[1, 2]] u[[2, 1]] + 
             u[[1, 1]] u[[2, 2]]) (v[[1, 2]] - 
             v[[2, 2]]) + (u[[1, 2]] - 
             u[[2, 2]]) (v[[1, 2]] v[[2, 1]] - v[[1, 1]] v[[2, 2]])}/
        detBottom);
     Boole[
      Sign[b - a] == Sign[b - pt] && Sign[a - b] == Sign[a - pt] && 
       Sign[d - c] == Sign[d - pt] && Sign[c - d] == Sign[c - pt]]]]];

Here is the intersection matrix computation using a NearestFunction`.

nf = Nearest[centers -> Range[Length[centers]]];
candidates[k_] := Rest[nf[centers[[k]], {Infinity, L}]]
AbsoluteTiming[neighborcandidates = Table[candidates[k], {k, n}];]
AbsoluteTiming[
 mymat = SparseArray[
    Flatten[Table[{j, k} -> 
       lineIntComp2[sticks[[j]], sticks[[k]]], {j, n}, {k, 
       neighborcandidates[[j]]}]]];]

(* Out[273]= {0.028073, Null}

Out[274]= {0.739792, Null} *)

An added benefit is that the adjacency matrix is kept as a SparseArray, hence less memory usage.

Here in contrast is the full-blown Outer method.

full = Sort /@ Transpose[{stickend1, stickend2}];
AbsoluteTiming[matComp = Outer[lineIntComp, full, full, 1];]

(* Out[285]= {62.015713, Null} *)

We check that we get the same result either way:

linkedComp[mat_] := Map[Boole[# != {1. 10^6, 1. 10^6}] &, mat, {2}]
linkedComp[matComp] === Normal[mymat]

(* Out[251]= True *)

I will retain the method you chose for determining top/bottom/left/right (one might argue that it should be the stick ends rather than centers that determine this; I guess which is better depends on the physics being simulated).

bottomMost = First@Ordering[centers[[All, 2]], 1];
topMost = First@Ordering[centers[[All, 2]], -1];
leftMost = First@Ordering[centers[[All, 1]], 1];
rightMost = First@Ordering[centers[[All, 1]], -1];

The picture is not quite as expected (dots instead of segments). I do not know why this happens but it seems to correlate with a large value for n. In any case it should give a clear indication that things are working.

With[{vert = FindShortestPath[ag, bottomMost, topMost]}, 
 HighlightGraph[ag, vert, ImageSize -> 700, 
  VertexShape -> 
   Thread[vert -> (Graphics[{Red, Thick, #}] & /@ lines[[vert]])]]]

(It will probably work better if done without the factorial...). The scaling is somewhere in the vicinity of what you propose. Certainly the m*n part is there. The last factor is probably in the vicinity of log(n) though honestly I'm not sure what exactly is the behavior. There has to be a log(m) component because we maintain a binary heap of size m. But I suspect there is an additive component here that is bigger, likely in the ballpark of log(n).

What Sander proposes, which I agree is the efficient way to proceed, can be done in a fairly simple manner without explicit binning. This assumes the stick length is fixed to some common value (as is the case in the running example set). Simply create a NearestFunction based on the stick centers. This can be done as below.

nf = Nearest[centers]

Then overlaps for stick k need only be checked against all sticks that have a center within length units of the center of stick k. Code would look like this.

candidates[k_] := Rest[nf[centers[[k]], {Infinity,L}]]

It is now quite fast to find, for every stick, a set of candidate overlapping sticks.

Timing[neighbors = Table[candidates[k], {k, n}];]

(* Out[146]= {0.000913, Null} *)

In this example we find there are around 25 candidates for each stick. Much smaller than the full 149 other sticks.

Map[Length, neighbors]

(* Out[150]= {21, 26, 24, 33, 28, 28, 35, 14, 22, 31, 34, 27, 24, 28, \
29, 20, 12, 19, 25, 27, 23, 26, 18, 28, 30, 11, 23, 33, 26, 23, 31, \
26, 31, 33, 14, 33, 20, 22, 23, 20, 34, 31, 28, 5, 26, 20, 24, 29, \
25, 20, 29, 26, 30, 26, 22, 11, 17, 30, 27, 16, 23, 24, 14, 21, 21, \
35, 21, 24, 20, 33, 28, 22, 34, 32, 33, 29, 24, 24, 13, 32, 31, 21, \
26, 12, 27, 24, 31, 8, 20, 15, 11, 23, 33, 29, 28, 15, 29, 32, 22, \
20, 10, 23, 33, 23, 14, 33, 15, 27, 24, 28, 30, 18, 30, 26, 22, 20, \
22, 22, 35, 34, 13, 15, 27, 19, 28, 19, 33, 28, 10, 28, 24, 17, 31, \
23, 36, 19, 17, 19, 16, 19, 11, 28, 23, 27, 32, 14, 24, 23, 19, 28} *)

Possibly this could be further refined to take into account orientations when computing neighbors. I doubt it will be worth the effort unless one goes into significantly larger simulations.

Dear Vitaliy: Excellent solution - as always!

Dear Tommy: You can convert a graphic to an image (with a specific resolution) using:

Image[lineGraphics, ImageResolution -> 200]

Regards -- Henrik

I will give you just a crude idea using WatershedComponents. That function can find the "ridge" between so called "catchment basins". You have to read up on that to understand. Taht "ridge" is your percolation path. It can also be for example a path in a maze. Alternatively, for general education for percolation coding, take a look at Finding a percolation path, Wolfram Demos, and NKS Book.

So let's go with WatershedComponents now. Your code a bit adjusted for around a percolation threshold:

L = 0.1;
n = 900;
SeedRandom[1]
centers = RandomReal[{0, 1}, {n, 2}];
? = RandomVariate[NormalDistribution[?/2, (1/2.3548)*?/4], n];
stickend1 = centers + Transpose[L {Cos[?], Sin[?]}/2];
stickend2 = centers - Transpose[L {Cos[?], Sin[?]}/2];
lines = Map[Line, Transpose[{stickend1, stickend2}]];
i = Binarize[ColorNegate[Graphics[{Thickness[.0], lines}, PlotRangePadding -> -.02]]]

Now you can find all the ridges, - those are 0-value components spanning the outer boundaries of the catchment basins:

Colorize[MorphologicalComponents[
  WatershedComponents[i, CornerNeighbors -> True] /. {0 -> 1, x_Integer /; x > 0 -> 0}]]

For percolating cluster you select the largest one:

path = SelectComponents[MorphologicalComponents[
    WatershedComponents[i, CornerNeighbors -> True] /. {0 -> 1, x_Integer /; x > 0 -> 0}], 
   "Count", -1] // Image

And you can highlight it on the original image too:

HighlightImage[i, path, Method -> {"Boundary", 2}]

A plausible next step is to form a graph wherein sticks are vertices and edges are based on incidence (that is, overlap) relations. Designate sticks at the top as initial vertical vertices, ones at the bottom as terminal, do similarly for those at left and right. Then graph/network metods can be used to determine existence or number of paths connecting initials to terminals.

Obviously this would take some coding, I just wanted to outline an approach one might pursue following up on this nice use of morphological component separation.

Hi Thomas,

one could use mathematical morphology:

lineGraphics = Graphics[lines, Frame -> False];
ColorNegate @ Colorize[MorphologicalComponents[ColorNegate @ lineGraphics]]

which gives:

Maybe this might serve as a start.

Regards -- Henrik

Making percolation paths visible