Message Boards Message Boards

4
|
12484 Views
|
23 Replies
|
34 Total Likes
View groups...
Share
Share this post:
GROUPS:

Making percolation paths visible

Posted 9 years ago

Dear Community,

I want to create and describe a two-dimensional network of percolating sticks. How could I make percolation paths visible by different colors?

Any suggestions? kind regards, Tommy

POSTED BY: J F
23 Replies

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

POSTED BY: Daniel Lichtblau

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

POSTED BY: Sander Huisman

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

POSTED BY: Daniel Lichtblau

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.

POSTED BY: Daniel Lichtblau

@ 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?

POSTED BY: Hans Dolhaine

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

POSTED BY: Sander Huisman

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

POSTED BY: Daniel Lichtblau

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

POSTED BY: Henrik Schachner

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

POSTED BY: Sander Huisman

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.

POSTED BY: Daniel Lichtblau

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?

POSTED BY: Sander Huisman

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

POSTED BY: Daniel Lichtblau

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

POSTED BY: Vitaliy Kaurov

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.

POSTED BY: Hans Dolhaine
Posted 9 years ago

Hi Henrik,

thank you very much!At first I have thought of an approach following graph theory, but this image analysis is sufficient. I thought about using Gwyddion, but I#d prefer keeping all in Mathematica.

Regards

POSTED BY: J F

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

POSTED BY: Henrik Schachner

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.

POSTED BY: Daniel Lichtblau
Posted 9 years ago

Thanks again, Henrik

POSTED BY: J F

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

POSTED BY: Hans Dolhaine

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.

POSTED BY: Daniel Lichtblau
Posted 9 years ago

Making percolation paths visible

POSTED BY: J F
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