Finding all minimal paths in site percolation?

Question

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

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

Answer 1

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

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

---

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

Answer 2

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

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

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]

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.

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