13
$\begingroup$

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

Mathematica graphics

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

My purpose it to draw something like this:

Mathematica graphics

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!

Improve this question
$\endgroup$
1
  • $\begingroup$ Possibly related: (46631) $\endgroup$ Commented Oct 7, 2014 at 20:03

1 Answer 1

Reset to default
12
+250
$\begingroup$

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

Mathematica graphics

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks a lot! I edited slightly your answer I hope you don't mind. I would suggest you remove your other answer (?) $\endgroup$ Commented Oct 9, 2014 at 6:50
  • $\begingroup$ I like the convolution idea but now that I understand the problem formulation it doesn't make a lot of sense. Cheers -S. $\endgroup$ Commented Oct 9, 2014 at 14:08

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.