Message Boards Message Boards

GROUPS:

2D jackstraw-like network: percolation path

Posted 3 years ago
6840 Views
|
23 Replies
|
34 Total Likes
|

Dear Community,

As I am new to Mathematica I do not have much general experience which I need for my specific problem. It is as follows:

I want to create and describe a two-dimensional network of percolating sticks. The percolation threshold and the relative number of sticks contributing to percolation depends on the number of sticks per area unit n, the orientation of the sticks theta and the stick length L.

My approach is the following:

Creating n statistically distributed Centers:

n = 150; 
SeedRandom[1]
centers = RandomReal[{0, 1}, {n, 2}];
    ListPlot[centers, Frame -> True] 

enter image description here

Creating orientation: Angular distribution is a Normal distribution with mu = pi/2 because i want to have the orientation vertical and not horizontal and sigma = 1/2.3548*FWHM

\[Theta]=RandomVariate[NormalDistribution[\[Pi]/2, (1/2.3548)*\[Pi]/4], n]

Defining the stick length and orientation from the centers:

L = 0.25 
stickend1 = centers + Transpose[L {Cos[\[Theta]], Sin[\[Theta]]}/2];
stickend2 = centers - Transpose[L {Cos[\[Theta]], Sin[\[Theta]]}/2];

Plotting:

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

enter image description here

Now I am able to see whether or not there is a percolation path in x- and/or y-direction depending on the stick density n, the orientation theta and the stick length L. However, it is impossible to detect a percolation path for high n and low L. Therefore it would be great if Mathematica shows these sticks which contribute to the percolation in x- or in y-direction in another color, or at least if there is any percolation path in x- or in y-direction.

Maybe someone of the more experienced Mathematica users can help me out.

Thank you very much!

Tommy

Attachments:
23 Replies

You need to provide some more detail. What comprises percolability? (Is that even a word? I guess it is now.) Is it the existance of a vertical path with few or no collisions? Can the angles change during motion? If so, what are the allowed changes? Is there something that happens on collisions? Is there a "depth" wherein some number of sticks are allowed to overlap?

The above is probably not exhaustive. Before anyone can offer suggestions it is really necessary to pin down what exactly you are trying to simulate.

Thanks for your answer and comments. The concept is to create a two-dimensional network of electrically conductive sticks (silver, copper, etc.). For electrical percolation (diffusion of electrons along the entire observed area) structural percolation (a physical path consisting of sticks) is required. Hence, percolation means that there is a physical pathway along the entire network due to crossing sticks (I guess that is what you mean with collision).

Is it the existance of a vertical path with few or no collisions?

If I understand your question correctly, it is exactly the opposite. Here, percolation means that the sticks are the path.

My desired procedure is:

n = n1 (let's say 150), then 100 runs with SeedRandom[1 to 100] with a defined FWHM1 (angular distribution)

For the 100 runs I would like to know IF there is a vertical path and IF there is a horizontal path created by the sticks. hence, I get the percolation probability in x- and in y-direction (how many runs with observed percolation within the 100 runs) for a particular FWHM (FWHM_1)

Then I do the same with a narrower or wider angular distribution (FWHM_2) and so on.

My aim is to find how the stick orientation affects the percolation probability. Therefore I need some more input to the mathematica notebook in order to make it "say" whether the sticks create a path or not (separately in vertical and horizontal direction). Due to the orientation of the sticks in vertical direction, there should be a higher percolation probability in vertical direction (bottom to top) than in horizontal (left to right). This effect should increase with narrow FWHM of the angular distribution.

Then I want to do the same procedure with higher and lower n (let's say 50 to 500).

I am a materials engineer who works on such systems experimentally and I have found that the conductivity is different for the two directions if the metal nanowires are aligned as I found out from micrographs. The simulations could give further explanation about this observation.

Furthermore, it would be helpful to show the sticks contributing to the path in another color than those who don't. Again, separately for the tow directions. Hence, one picture with highlighted sticks that contribute to percolation in x-direction and respectively another picture for those sticks contributing to percolation in y-direction.

Sounds complicated..... My crude proposal would be: you have a list of sticks (straight lines from A to B). Pick the first one and see, whether it is crossed by another one, then pick this one and do the same and so on, and look whether you reach the opposite end of your domain.

regards Hans

Hi Thomas,

one could use mathematical morphology:

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

which gives:

enter image description here

Maybe this might serve as a start.

Regards -- Henrik

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.

Thanks again, Henrik

I have tried these procedure and it seems to work well and agree with my experimental results. However sometimes the sticks are aligned very close but without overlap but the resulting image shows these rods having a contact. I guess that is because of the pixel density. Therefore my question: Is there a possibility for creating smaller pixels/higher pixel density?

Regards, Tommy

Hi Henrik,

thank you very much! Yes, this might serve as a start. At first I have thought of an approach following graph theory, as Hans has suggested. However, my Mathematica skills are very rudimentary and I cannot create an algorithm for that, at least not within a time that i could justify towards the rest of my time. Your image analysis (I guess that is what you mean by mathematical morphology) is a good start.

Regards,

Tommy

To proceed in the more crude way: this gives a matrix T. If T[[ i , j ]] equals 1, the lines i and j cross each other.

    intsec[Line[{A_, B_}], Line[{U_, V_}]] := Module[{},
      dd = Det[{A - B, V - U}];
      If[Abs[dd] < 0.0001, 0,
       {t1, t2} = Inverse[{A - B, V - U}].(V - B);
       If[(0 <= Abs[t1] <= 1) && (0 <= Abs[t2] <= 1), 1, 0]
       ]
      ]

T = Table[
intsec[lines[[j]], #] & /@ lines, {j, 1, n}];
T // Short[#, 15] &

OK, T is a square matrix, which is symmetrical, so we could reduce memory used by either considering the upper (lower) triangle, or a sparse array object, but I have no experience with these.

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

enter image description here

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

enter image description here

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

enter image description here

And you can highlight it on the original image too:

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

enter image description here

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

Indeed a very nice answer!

If Thomas wants to to prevent close-calls, then one has to just do it 'geometrically', this is however much more tricky and time-consuming I'm afraid...

Finding the overlaps can be done in $n^2$ time as Hans said, but it can be done faster if one bins (the centers of) the sticks in bins with a width of the length of the bins. Then only check 'overlaps' for elements in 1 bin, and it's neighbours. This should reduce the $n^2$ a lot if the density is low. If the total 'area' is $A$ then the average density is $n/A$. So in a box of size $L$ you will find roughly $n\times L^2/A$ straws. So you want to check the straws in a box + neighbours it will be at most $ n^2 \times \frac{9L^2}{A}$. So if you're area is big compare to $L^2$, binning and checking neighbours might be a lot faster (binning just takes $n$ time).

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.

Didn't think about the NearestFunction[..][x,{inf,r}] ! that is indeed also fast. Fast and easy implementation. $m\times n\times \log(n)$ scaling right? where $m$ is the average number of candidates in the vicinity?

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

Going off of the ideas of treating this as a Graph and the geometric intersections Daniel and Sanders proposed above respectively. I use the same setup as the initial question by the user in the attachment, I just change n to 250 (to get percolation in both x and y directions from edge nodes).

We find the intersections of line segments explicitly by modifying the determinants method, (as given by the Mathworld page). given both in uncompiled and compiled versions for comparison:

lineIntersections[{a_, b_}, {c_, d_}] := 
 Block[{detBottom = Det[{a - b, c - d}], pt},
  If[Chop[detBottom] == 0, Infinity, 
   pt = ((Det[{a, b}] (c - d) - Det[{c, d}] (a - b))/detBottom);
   If[Sign[b - a] == Sign[b - pt] && Sign[a - b] == Sign[a - pt] && 
     Sign[d - c] == Sign[d - pt] && Sign[c - d] == Sign[c - pt], pt, 
    0]]]

lineIntComp = 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, {10^6, 10^6}, 
     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);
     If[Sign[b - a] == Sign[b - pt] && Sign[a - b] == Sign[a - pt] && 
       Sign[d - c] == Sign[d - pt] && Sign[c - d] == Sign[c - pt], 
      pt, {10^6, 10^6}]]
    ]];

We can then apply this on all pairs using Outer (definitely inefficient but perhaps sufficient, or user can apply earlier suggestions on scaling?)

full = Sort /@ Transpose[{stickend1, stickend2}];
mat = Outer[lineIntersections, full, full, 1]; // Timing
(*Out[16]= {1.98438,Null}*)
matComp = Outer[lineIntComp, full, full, 1]; // Timing
(*Out[17]= {0.21875,Null}*)

One can then create an Adjacency Matrix, connecting only intersecting nodes (sticks) and generate the graph whereby VertexCoordinates can be specified as the lines defined earlier by the user.

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

linkedComp[matComp] == linked[mat]
ag = AdjacencyGraph[linked[mat], VertexCoordinates -> centers, 
  EdgeStyle -> White, 
  VertexShape -> 
   Thread[Range[n] -> (Graphics[{Opacity[.5], #}] & /@ lines)], 
  VertexSize -> 250]

Now we can use various graph functionality to evaluate things like percolation. Here, I just define some extreme points for visualization purposes:

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]

Vertical percolation:

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

enter image description here

Horizontal percolation:

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

enter image description here

all connected components to bottomMost stick:

With[{vert = 
   Reap[BreadthFirstScan[ag, bottomMost, {"PrevisitVertex" -> Sow}]][[
    2, 1]]},
 HighlightGraph[ag, vert, ImageSize -> 700, 
  VertexShape -> 
   Thread[vert -> (Graphics[{Red, Thick, #}] & /@ lines[[vert]])]]]

enter image description here

Hope this helps.

Best, George

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]

enter image description here

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

enter image description here

Huge difference in speed! wow! I wouldn't have thought it would be that much! Nice code!

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

Thanks Daniel, very nice addition! Wasn't comfortable enough so as to implement your binning suggestion in the code (new to Community, and not entirely familiar how one should use/credit others). The VertexShape function seems to be dependent on VertexSize for some reason and that's why you only obtain the nodes and not the sticks, I believe. I had to deduce it empirically as 250 for my n=250 case (it doesn't seem to always vary linearly with n though).

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

@ Daniel: Very cool. I gave it a vote. And I think it leaves more opportunities than the image-processing, but that may be wrong.

Could one gain additional speed if ordering the sticks according to their - say - y-coordinate?

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

Hans, I'm not really sure but I think the answer is no. As Sander points out, Nearest already has its own internal analog to sorting. Ordering by y or other values is not likely to offer much unless one can find a way to make it cooperate with calls to a NearestFunction. I don't know how to do that (which does not mean it can't be done).

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract