graphs and networks - Calculating morphometric properties

Is it possible to calculate segment length, segment diameter and branch angles for three dimensional geometry. Geometry contains vertex points and polygon surfaces as

  file = "http://i.stack.imgur.com/N6RYT.png";
  Import[file, "Graphics3D"]

  data = Import[file, "VertexData"];
  polysurface = Import[file, "PolygonData"];
  Graphics3D[Point@data]

Answer

I decided to post another answer for several reasons:

• I made up a full working contraption by using another approach


• The previous answer could be useful for others, so I prefer to leave it there


• Both answers are quite long, and having both in one post will clutter it

That said, the plan is the following:

1. Collapse the points as before, pre-clustering them over an skeleton


2. Calculate real clusters and collapse each point to the Center of Mass of its cluster


3. From the polygon's description, get the surviving lines that form a 3D skeleton for the figure


4. Find the branch points in the skeleton, and the incoming lines to each branch point



5. Calculate the angles

Limitations/ To-dos:

1. The method used is purely metric, without regarding the topological characteristics of the body. It works very well here, because the branches are not nearby parallel structures. May fail for tangled things


2. There are two iterative loops that (perhaps) will be converted to something more Mathematica-ish if I find time/momentum and a way to do it.

Note: As the code here will be interspersed with graphics and comments, you may like to download a fully functional .nb in only one step. For that, execute the following and you'll get a new notebook ready for running, downloaded from imgur:

NotebookPut@ImportString[Uncompress@FromCharacterCode@
TakeWhile[Flatten@ImageData[Import["http://i.stack.imgur.com/tP8xo.png"],"Byte"],#!=0&],"NB"]

Ok. Let's get something done:

i = Import["http://dl.dropbox.com/u/68983831/final.vtk", "GraphicsComplex"];
j = i[[1]]; (*Point coordinates*)
(*Number of Nearest neighbors to take in account for collapsing the 3D structure*)
nn = 40;

j1 = j/2; (*some initial value for the iteration stopping test*)
(*Search for a quasi fixed point by moving towards the CoM of NN*)
While[Total[Abs[j1 - j], -1] > Length@j/100,

 j1 = j;
 (j[[#]] = Mean[Nearest[j, j[[#]], nn]]) & /@ Range@Length@j]
jp = j; (*To preserve j for later use*)
(*Show our proto-clusters*)
Show[Graphics3D[{Opacity[.1], EdgeForm[None], 
                GraphicsComplex[i[[1]], i[[2]]]}], ListPointPlot3D[j]]

Now we have a nice arrangement of all our points over one skeleton, we need to proceed further and collapse each cluster to its Center of Mass, to reduce the number of points drastically and allow us to have segments linking those points. Sorry, iteration follows:

(*Now collpse every point to its corresponding CoM*)

(*maxdist is an estimation of a safe intercluster distance*)
maxdist = Mean[EuclideanDistance[j[[#]], First@Complement[
                                              Nearest[j, j[[#]], 81], 
                                              Nearest[j, j[[#]], 80]]] & /@ 
                               RandomSample[Range@Length@j, IntegerPart[Length@j/5]]]/2;
h = 0;
agg (*buckets for clusters*) = n = {};
(*Calculate clusters, FindClusters[] doesn't do any good here :( *)
While[j != {},
  h++;

  AppendTo[agg, {First@j}];
  While[j != {} && EuclideanDistance[agg[[h, 1]], 
                                     n = (Nearest[j, agg[[h, 1]], 1][[1]])] < maxdist,
   j = DeleteCases[j, n];
   AppendTo[agg[[h]], n];
   ]
  ];
(*Clusters calculated, collapse each cluster to its mean*)
rules = (Thread[Rule[#, x]] /. x -> Mean[#]) & /@ agg;
j = (jp /. Flatten[rules, 1]);

Now we have our points completely collapsed. Lets see what we found:

(*A plot showing the lines found*)
rul = Thread[Rule[Range@Length@j, j]];
vl = (i[[2]] /. rul) /. {a_List, b_List, c_List} /; a != b != c ->  Sequence[];
uniLines = Union[Sort /@ Flatten[Subsets[#, {2}] & /@ vl[[1, 1]], 1]] /. {a_, a_} -> Sequence[];
Graphics3D@Table[{Thick, Hue@RandomReal[], Line@l}, {l, uniLines}]

Nice. Now we can explore the endpoints for each segment:

Manipulate[Show[
  Graphics3D@{Opacity[.1], EdgeForm[None], GraphicsComplex[i[[1]], i[[2]]]}, 
  Graphics3D@Table[{Hue[o/Length@uniLines], Line@uniLines[[o]]}, {o, 1, Length@uniLines}], 
  Graphics3D[{PointSize[Medium], Point[Union@Flatten[uniLines, 1]]}], 
  Graphics3D[{Red, PointSize[Large], Point[uniLines[[w]]]}]
  ], {w, 1, Length@uniLines, 1}]

We are almost there, now we identify the branch points and the related lines:

(*Now get and process the bifurcation points*)

branchPt = First /@ Extract[#, Position[Length /@ #, 3]] &@(Gather@Flatten[uniLines, 1]);
(*Get the incoming lines to each point *)
branchLn = Extract[uniLines, #] &/@ ((Position[uniLines, #] & /@ branchPt)[[All,All,1 ;;1]]);

Let's see if we identified our branch points correctly:

(*Plot Bifurcations*)
Show[Graphics3D[{Opacity[.1], EdgeForm[None],GraphicsComplex[i[[1]], i[[2]]]}],     
     Graphics3D[Line /@ branchLn],
     Graphics3D[{Red, PointSize[Large], Point[branchPt]}]]

Now we get normalized vectors along the lines, outgoing from each branch point

(*Get the normalized vectors at the branch points*)
aux[x_] := Map[# - First@Commonest[Flatten[x, 1]] &, x, {2}]; 
nv = Map[Normalize, ((Select[Flatten[aux[#], 1], # != {0, 0, 0} &]) & /@ branchLn), {2}];

Finally we are going to get our desired 9 angles. As you can observe, in all cases the three segments at each branch point are almost coplanar, so the sum of the angles ought to be near 2 Pi. We now calculate the angles and verify the coplanarity observation:

(*Angles and kind of Coplanarity test*)
(VectorAngle @@@ Subsets[#, {2}]) & /@ nv
Total /@ ((VectorAngle @@@ Subsets[#, {2}]) & /@ nv)


(* Angles: {{2.0326, 2.54025, 1.70803}, 
            {1.71701, 2.4161, 2.14805}, 
            {2.14213, 1.98282, 2.158}}*)

(* Coplanarity test (Very near 2 Pi): {6.28087, 6.28115, 6.28295} *)