Connect neighbouring points as list of segments in 2 D

Context

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.

Problem

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] 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, Range]] // Normal // MatrixPlot

So skl contains zeros and ones wherever that pixel lies on a ridge. 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] & In practice it is useful to know which points on the skeleton are connected to which.

Question

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.

Attempt

Using the function Nearest I have tried

nf = Nearest[pos];
tt = Map[nf[#, 4] &, pos, 1];
tt2 = Select[tt, Norm[((# - Table[#[], {4}]))] < 3 &];
gr = Map[Line[{Rest[#], Table[#[], {Length[#] - 1}]} // Transpose] &, tt2] // Graphics; which seems to do the trick but there might be more efficient/robust way to proceed?

• Comments are not for extended discussion; this conversation has been moved to chat. – 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)..

ClearAll[arrayGraph];
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] 