6
$\begingroup$

I want to find all minimal site percolation paths from the left side to the right side, small demo grid.

A related question is here on finding bond percolation paths instead of minimal site percolation paths. A site is a neighbour to other site vertically or horizontally but not diagonally.

I want to use the FindShortestPath[], but I need to somehow get a GridGraph working: the basic idea is that unoccupied sites are removed but the vertex removal does not work as expected

enter image description here

where I try to model the lattice as a grid graph, but I cannot understand why a vertex deletion results to a path graph, this puzzle moved here.

How can I find the minimal site percolation paths?

SeedRandom[10801];
dimension = 100;
coDimension = 30;
percProbability = 0.7;
myData = Table[
   Table[Boole[RandomReal[] < percProbability], {i, dimension}], {j, 
    coDimension}];
myData // MatrixPlot

enter image description here

$\endgroup$
  • $\begingroup$ Some people here would know how to convert this to a graph… I think using Sparsify?? You should make your question not an image though… Did you see mathematica.stackexchange.com/a/5156/1089 $\endgroup$ – chris Jul 9 '16 at 18:31
  • $\begingroup$ @chris the image-based solution is computationally demanding on my comp, taking forever to compute :/ $\endgroup$ – hhh Jul 9 '16 at 18:36
8
$\begingroup$

Here's a solution using MorphologicalGraph[]:

SeedRandom[10801];
dimension = 100;
coDimension = 30;
percProbability = 0.7;
myData = Table[Table[Boole[RandomReal[] < percProbability], {i, dimension}],
          {j, coDimension}];
img = Binarize@Image@myData;

Now all you need to do is use FindShortestPath[]. For example, the shortest path from top-left to bottom-right corner:

g = MorphologicalGraph[img]
HighlightGraph[g, PathGraph[FindShortestPath[g, 1, Max@VertexList[g]]]]

enter image description here

However, MorphologicalGraph[] includes corner neighbours, which we don't want in site percolation on a square lattice. A bit of digging turned up Image`MorphologicalOperationsDump`oMorphologicalGraph[] as the function behind MorphologicalGraph. By adapting this function to ignore the thinning operation and only use corner neighbours, you can get the appropriate graph. The adapted function is at the bottom of the post.

First, let's get the left and right vertices out from your data:

getIndex[sites_, site_] := Position[sites, _?(# == site &)]
getLeftAndRightVertices[data_] := 
 Module[{sites, leftsites, rightsites},
  sites = Position[data, _?(# == 1 &)];
  leftsites = Select[sites, #[[2]] == 1 &];
  rightsites = Select[sites, #[[2]] == Last@Dimensions@data &];
  {Flatten[getIndex[sites, #] & /@ leftsites],
   Flatten[getIndex[sites, #] & /@ rightsites]}]

{leftvertices, rightvertices} = getLeftAndRightVertices[myData];

Now you can find the shortest path between any left vertex and any right vertex. To get the overall shortest path from left to right, you can do the following, with a warning when no path can be found.

g2 = myMorphologicalGraph[img, VertexCoordinates -> Automatic];

allpairs = Tuples[{leftvertices, rightvertices}];
allpaths = Quiet@MapThread[FindShortestPath[g2, #1, #2] &, Transpose@allpairs];
pathlengths = Length@# & /@ allpaths;

(* Workaround to deal with unconnected components which *)
(* give lengths of 0. The minimum possible path length *)
(* is of course == dimension *)
If[Max@pathlengths >= dimension,
 sortlengths = Ordering[pathlengths];
 pos = FirstPosition[pathlengths[[sortlengths]], _?(# >= dimension &)];
 shortestpair = Flatten@allpairs[[sortlengths]][[pos]];
 shortestpath = Flatten@allpaths[[sortlengths]][[pos]];
 Show[img, HighlightGraph[g2, PathGraph[shortestpath]]],
 (* Print warning *)
 Print["No connected path found"]]

(* Left = 770, Right = 1242 *)
(* Length = 128 *)

enter image description here


Code for myMorphologicalGraph[]:

myMorphologicalGraph[skeleton_, 
   opts : OptionsPattern[MorphologicalGraph]] := 
  Module[{vertices, vertexComponents, vertexCoordinates, 
    vertexCount = 0, linkComponents, linkCount, linkWeights, 
    onePixelLinks, all4all, vertex4all, vertex4vertex, vertex4link, 
    all8all, vertex8all, vertex8vertex, vertex8link, link8all, 
    link8vertex, link8link, redundantEdges, directEdges, linkedEdges, 
    loopEdges, cleanEdges, extraEdges, allEdges, 
    edgeWeights},(*vertices=ImageAdd[MorphologicalTransform[skeleton,
   "SkeletonEndPoints",Padding\[Rule]0],MorphologicalTransform[
   skeleton,"SkeletonBranchPoints",Padding\[Rule]100]];*)
   (* Use the direct binarized image *)
   vertices = skeleton;
   vertexComponents = 
    Replace[ImageData[vertices], 1 :> ++vertexCount, {2}];
   vertexCoordinates = OptionValue[VertexCoordinates];
   If[vertexCoordinates === Automatic, 
    vertexCoordinates = 
     ComponentMeasurements[vertexComponents, "Centroid"]];
   linkComponents = 
    Image`MorphologicalOperationsDump`ConstrainedMComponents[
     ImageSubtract[skeleton, vertices], 
     Dilation[vertices, CrossMatrix[1]]];
   linkCount = Max[linkComponents];
   linkComponents = 
    Replace[linkComponents, Except[0, n_] :> n + vertexCount, {2}];
   linkWeights = 
    Dispatch[ComponentMeasurements[linkComponents, "Count"]];
   all4all = 
    ComponentMeasurements[vertexComponents + linkComponents, 
     "Neighbors", CornerNeighbors -> False];
   vertex4all = Select[all4all, First[#1] <= vertexCount &];
   vertex4vertex = 
    DeleteCases[vertex4all, _?(#1 > vertexCount &), {3}];
   vertex4link = DeleteCases[vertex4all, _?(#1 <= vertexCount &), {3}];
   (* Set corner neighbours to False here *)
   all8all = 
    ComponentMeasurements[vertexComponents + linkComponents, 
     "Neighbors", CornerNeighbors -> False];
   vertex8all = Select[all8all, First[#1] <= vertexCount &];
   vertex8vertex = 
    DeleteCases[vertex8all, _?(#1 > vertexCount &), {3}];
   vertex8link = DeleteCases[vertex8all, _?(#1 <= vertexCount &), {3}];
   link8all = Select[all8all, First[#1] > vertexCount &];
   link8vertex = DeleteCases[link8all, _?(#1 > vertexCount &), {3}];
   link8link = DeleteCases[link8all, _?(#1 <= vertexCount &), {3}];
   redundantEdges = 
    Image`MorphologicalOperationsDump`sortEdges[
     DeleteCases[
       Image`MorphologicalOperationsDump`growEdges[
        Image`MorphologicalOperationsDump`growEdges[
         Image`MorphologicalOperationsDump`toEdges[vertex4vertex], 
         vertex4link], link8vertex], 
       edge[_, v_, _, v_] | edge[v_, _, _, v_]][[All, {1, 3, 4}]]];
   directEdges = 
    Select[Image`MorphologicalOperationsDump`toEdges[vertex8vertex], 
     OrderedQ];
   linkedEdges = 
    Image`MorphologicalOperationsDump`growEdges[
     Image`MorphologicalOperationsDump`toEdges[vertex8link], 
     link8vertex];
   loopEdges = 
    Cases[Tally[linkedEdges, #1[[2]] === #2[[2]] &], {e_, 1} -> e];
   loopEdges = 
    Pick[loopEdges, 
     Thread[Replace[loopEdges[[All, 2]], linkWeights, {1}] > 2]];
   linkedEdges = Select[linkedEdges, First[#1] < Last[#1] &];
   cleanEdges = 
    Join[Complement[linkedEdges, redundantEdges], loopEdges];
   onePixelLinks = 
    Union[Flatten[vertex4link[[All, 2]]]] \[Intersection] 
     Cases[Tally[Flatten[linkComponents]], {_, 1}][[All, 1]];
   extraEdges = 
    Image`MorphologicalOperationsDump`growEdges[
     Image`MorphologicalOperationsDump`growEdges[link8vertex, 
      Select[Cases[
        Image`MorphologicalOperationsDump`growEdges[
         edge[onePixelLinks], link8link], 
        edge[_, Alternatives @@ onePixelLinks]], OrderedQ]], 
     link8vertex];
   extraEdges = 
    Image`MorphologicalOperationsDump`sortEdges[
     DeleteCases[extraEdges, 
      Alternatives @@ 
       Cases[extraEdges, 
        edge[v_, l1_, l2_, v_] -> edge[_, l1, l2, _]]]];
   allEdges = 
    Apply[UndirectedEdge, 
     directEdges \[Union] cleanEdges[[All, {1, -1}]] \[Union] 
      extraEdges[[All, {1, -1}]], {1}];
   edgeWeights = OptionValue[EdgeWeight];
   If[edgeWeights === Automatic, 
    edgeWeights = 
     Join[Thread[Apply[UndirectedEdge, directEdges, {1}] -> 1], 
      Sort[Thread[
        Apply[UndirectedEdge, cleanEdges[[All, {1, -1}]], {1}] -> 
         1 + (cleanEdges[[All, 2]] /. linkWeights)]], 
      Thread[Apply[UndirectedEdge, extraEdges[[All, {1, -1}]], {1}] ->
         3]]; edgeWeights = 
     Replace[allEdges, Dispatch[edgeWeights], {1}]];
   Graph[allEdges, VertexCoordinates -> vertexCoordinates, 
    EdgeWeight -> edgeWeights, 
    Sequence @@ 
     FilterRules[opts, 
      DeleteCases[
       Options[MorphologicalGraph], (VertexCoordinates -> _) | \
(EdgeWeight -> _)]]]];
|improve this answer|||||
$\endgroup$
  • $\begingroup$ Did you try to solve the shortest path between all left vertices and all right vertices? For me, the computer freezes in trying to do it. $\endgroup$ – hhh Jul 9 '16 at 21:20
  • $\begingroup$ @hhh did you try the update to see if your computer can solve it in a suitable amount of time? It's pretty quick on my PC. $\endgroup$ – dr.blochwave Jul 10 '16 at 11:07
  • 1
    $\begingroup$ Yes it works, nice one +1! $\endgroup$ – hhh Jul 10 '16 at 11:12
  • $\begingroup$ Hmmm, I tried this with different lattice sizes and for some reason, it is unable to find the minimal path. Have you tried this with different lattice sizes? $\endgroup$ – hhh Jul 10 '16 at 12:04
  • 1
    $\begingroup$ @hhh issue fixed I think, now issues a warning if no path found. $\endgroup$ – dr.blochwave Jul 10 '16 at 12:19
3
$\begingroup$

Solution based on the GridGraph

SeedRandom[10801];
dimension = 20;
coDimension = 30;
percProbability = 0.7;    
deleteMe = 
  Pick[Table[i, {i, dimension*coDimension}], 
   Table[RandomReal[] > percProbability, {i, 
     dimension*coDimension}]];
G = GridGraph[{dimension, coDimension}, VertexLabels -> "Name", 
   ImagePadding -> 30];
G = SetProperty[G, VertexCoordinates -> GraphEmbedding[G]];
H = VertexDelete[G, deleteMe]
FindShortestPath[H, 1, 600]
HighlightGraph[H, PathGraph[%]]

enter image description here

that finds the shortest path from the site 1 to the site 600. And next I want to find the shortest path from the left side to the right side

rightSide = Complement[Table[i, {i, 581, 600}], deleteMe]
Table[FindShortestPath[H, 1, i], {i, rightSide}]
Table[Length[FindShortestPath[H, 1, i]], {i, rightSide}]
shortest = Table[FindShortestPath[H, 1, i], {i, rightSide}][[4]];
HighlightGraph[H, PathGraph[shortest]]

enter image description here

that is the shortest path of length 47 between the vertex 1 and the left side. Next I need to do this over each site on the right side to find the shortest path between the right side and the left side

paths = Table[
   Table[FindShortestPath[H, j, i], {i, rightSide}], {j, leftSide}];
pathLengths = Table[
  Table[Length[paths[[k]][[h]]], {k, Length[paths]}],
  {h, Length[paths[[]][[1]]]}]
Histogram[pathLengths, 50]
pathLengths // Max
FindShortestPath[H, rightSide[[17]], leftSide[[13]]]
HighlightGraph[H, PathGraph[%]]
Pick[pathLengths, pathLengths // Positive]

enter image description here

where the zero entries tell me that there are sites from which no path to the other side. So we need to look for positive entries only to find the minimal paths and we found it to be of 33 size.

enter image description here

The largest pairwise minimal path length is 53 and it turns out to be between 1-596.

FindShortestPath[H, rightSide[[13]], leftSide[[1]]]
HighlightGraph[H, PathGraph[%]]

enter image description here

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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