I am interested in connecting neighboring points in 2/3D as list of segments. I am guessing this is something within the reach of graph theory, which is well implemented in mathematica.


To be more specific, let me consider connecting together points which correspond to critical lines in 2D.

Starting from a GaussianRandomField

u = GaussianRandomField[n = 128*2, 2, Function[k, 1/k Exp[-1/40 k^2]]] // Chop;

I can identify the 2D skeleton (the 'ridges' of the density field see this post or that answer) as

 skl = Map[If[# == 0, 1, 0] &, -u // Image // WatershedComponents, {2}] // SparseArray; 
 HighlightImage[ReliefImage[-u // Image, ColorFunction -> ColorData["ThermometerColors"]],
 skl // Image, "HighlightColor" -> Yellow]

Mathematica graphics

follow this link to upload an example of such skeleton.

For instance, looking at a small subset of the table I would have

 skl[[Range[15], Range[15]]] // Normal // MatrixPlot

So skl contains zeros and ones wherever that pixel lies on a ridge.

Mathematica graphics

I can then extract the coordinates of all points on the skeleton

 pos = skl // ArrayRules // Most // Map[First, #] & // Sort;

Indeed, they correspond to the same skeleton

 Reverse[#*{-1, 1}] & /@ pos // ListPlot[#, AspectRatio -> 1, Axes -> False] &

Mathematica graphics

In practice it is useful to know which points on the skeleton are connected to which.


I would like to join all neighbouring points of the skeleton as segments. Ideally I would like to tag each segment with a label keeping track of which branch of the skeleton it belongs to.


Using the function Nearest I have tried

 nf = Nearest[pos];
 tt = Map[nf[#, 4] &, pos, 1];
 tt2 = Select[tt, Norm[((# - Table[#[[1]], {4}]))] < 3 &];
 gr = Map[Line[{Rest[#], Table[#[[1]], {Length[#] - 1}]} // Transpose] &, tt2] // Graphics;

Mathematica graphics

which seems to do the trick but there might be more efficient/robust way to proceed?

  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – Mr.Wizard Oct 17 '14 at 18:17

As Oska noted in the comments, adding the option CornerNeighbors->False in arrayGraph gives the desired result for 2D. (See also this answer to a different Q/A)..

arrayGraph[mat_, opts : OptionsPattern[]] := 
 Module[{m = Module[{i = 1}, mat /. 1 :> i++], edges, vcs, v}, 
  v = ComponentMeasurements[m, "Label"][[All, 1]];
  vcs = ComponentMeasurements[m, "Centroid"];
  edges = UndirectedEdge @@@ DeleteDuplicates[Sort /@ Flatten[Thread /@ 
        ComponentMeasurements[m, "Neighbors", CornerNeighbors -> False]]];
  Graph[v, edges, VertexCoordinates -> vcs, opts]]

arrayGraph[Normal@skl, VertexSize -> .3, EdgeStyle -> Directive[Thick, Red], ImageSize -> 650]

enter image description here


Your Answer

