Morphological Filtering in 3D to produce skeletons

Question

Context

As a follow up of this question and that answer, I would like to identify the special lines separating 3D watersheds. These are useful in the context of astronomy to identify the filaments of the so called cosmic web. In practice, this involves identifying the lines at the intersections of theses surfaces

which are the boundaries of the watersheds of a Gaussian random field

My purpose it to draw something like this:

Given that mathematica 10.0.1 now deals with 3D WatershedComponents I would like to apply some kind of MorphologicalTransform on 3D cubes in order to identify the intersections.

Problem

Say I have a cube such as

dat={{{1, 1, 0, 4, 4, 4}, {1, 1, 0, 4, 4, 4}, {0, 0, 0, 0, 4, 4}, {2, 2, 2, 0, 0, 0},
{2, 2, 2, 0, 3, 3}, {2, 2, 2, 0, 3, 3}}, {{1, 1, 0, 4, 4, 4}, 
{0, 0, 0, 0, 4, 4}, {0, 0, 0, 0, 4, 4}, {0, 0, 0, 0, 0, 0}, 
{2, 2, 2, 0, 3, 3}, {2, 2, 2, 0, 3, 3}}, {{0, 0, 0, 0, 4, 4}, 
{0, 0, 0, 0, 4, 4}, {5, 5, 5, 0, 4, 4}, {0, 0, 0, 0, 0, 0}, 
{2, 2, 2, 0, 3, 3}, {2, 2, 2, 0, 3, 3}}, {{5, 5, 5, 0, 4, 4}, 
{5, 5, 5, 0, 4, 4}, {5, 5, 5, 0, 4, 4}, {0, 0, 0, 0, 0, 0}, 
{0, 0, 0, 0, 3, 3}, {0, 0, 0, 0, 3, 3}}, {{5, 5, 5, 0, 4, 4}, 
{5, 5, 5, 0, 4, 4}, {5, 5, 5, 0, 4, 4}, {5, 5, 5, 0, 0, 0}, 
{5, 5, 5, 0, 3, 3}, {5, 5, 5, 0, 3, 3}}, {{5, 5, 5, 0, 4, 4}, 
{5, 5, 5, 0, 4, 4}, {5, 5, 5, 0, 4, 4}, {5, 5, 5, 0, 0, 0}, 
{5, 5, 5, 0, 3, 3}, {5, 5, 5, 0, 3, 3}}};

(which is the result of the 3D watershed provided by mathematica).

Question

I am interested in identifying efficiently voxels which value is zero AND which is surrounded by at least 3 voxels with different values (zero excluded).

Any suggestions would be welcome!

Answer 1

With the understanding that the criterion for including a point in the 3D matrix is that it has 3 unique domains in any of the 26 points surrounding a zero value in the watershed here is a simple way to extract the data. In the case where zero values in the watershed are in minority it may be the fastest approach to getting a list of points that fulfill the criterion.

Simply:

intersections[ws_] := Module[{unique, dims},
  dims = Dimensions[ws];
  Reap[
    Do[
      If[ws[[i, j, k]] == 0,
        unique = (Flatten[
          ws[[i - 1 ;; i + 1, j - 1 ;; j + 1, k - 1 ;; k + 1]]] // 
          Union)~Drop~1;
        If[Length[unique] >= 3, Sow[{i, j, k, unique}]]
      ], 
      {i, 2, dims[[1]] - 1}, 
      {j, 2, dims[[2]] - 1}, 
      {k, 2, dims[[3]] - 1}
    ]
  ][[2, 1]]
]

I know - hardly elegant, but it crunches your 64^3 array in a half second on my laptop.

Then the skeleton of a cube is simply given by

  skl[cube_] := Module[{list = intersections[cube]},
             SparseArray[#[[1 ;; 3]] -> 1 & /@ list] ]

Let us try it on a Gaussian random field

  u = GaussianRandomField[n = 32*6, 3, Function[k, 1/k Exp[-1/2 k^2]]] //Chop;
  dat = WatershedComponents@Image3D[u];
  dat2 = WatershedComponents@Image3D[-u];

to produce the set of critical lines connecting peaks and minima to saddles:

  ImageAdd[Image3D[Normal[skl[#]]] & /@ {dat, dat2}] // Rasterize