# 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]

-
Six hours ago you asked a question. In the answers to that question you can find at least two different ways to calculate segment lengths and angles. Could you please specify why this question is different? –  belisarius Sep 21 '12 at 5:18
@belisarius: That question was for 2d image, which had center line in the image. This question have 3d geometry and it doesnt have flat plane 2d image, I think this is entirely a different problem. –  mathew Sep 21 '12 at 5:27
One possible approach is to Rasterize the 3D graphics with various rotations applied, and use the thinning and morphological operations you already have for 2D. It should be relatively easy to solve for the 3D vertex positions from a set of known rotations and corresponding 2D vertex positions. I guess the tricky part will be to identify the same vertices in multiple views. Sorry, I don't have time to write any actual code... –  Simon Woods Sep 21 '12 at 7:47
@Jay, just to let you know I had a go using my idea, and there were more complications than I realised! In particular, the number of vertices found by MorphologicalGraph varies with the rotation of the object. It's an interesting problem. –  Simon Woods Sep 22 '12 at 17:44
@SimonWoods see my (partial) answer. Perhaps inspires you some idea on how to get there –  belisarius Sep 24 '12 at 9:53

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*)
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} *)

-
does angles are in degrees? –  mathew Sep 29 '12 at 20:57
@Jay Nope. In radians –  belisarius Sep 29 '12 at 20:58
+1 Excellent stuff :-) –  Simon Woods Sep 30 '12 at 14:30
Great effort, +1 ;-) –  Vitaliy Kaurov Sep 30 '12 at 17:04

Not a full answer, but gotta something and perhaps it's useful so someone else can help to finish.

The idea is trying to get a skeleton by collapsing the points to the mass center of a group of nearest neighbors. Like this:

file = "http://dl.dropbox.com/u/68983831/final.vtk";
i = Import[file, "GraphicsComplex"];
j = i[[1]];
k = 1 ;(*A parm. In this set up we collapse immediately*)
iter = 4; (*Number of iterations*)
nn = 40; (*Nearest neighbors to take in account*)

(* Apply elastic pressure *)
For[h = 1, h <= iter, h++,
For[s = 1, s <= Length@j, s++,
j[[s]] = (Mean[Nearest[j, j[[s]], nn]]) - j[[s]] (1 - k)]]

(*Show*)
Show[Graphics3D[{Opacity[.1], EdgeForm[None],
GraphicsComplex[i[[1]], i[[2]]]}], ListPointPlot3D[j]]


Now we have a useless set of scattered points. But not so useless if you remember that they were originally part of polygons. So:

Graphics3D@GraphicsComplex[j, i[[2]]]


And now we have a nice alien creature:

I think if one can figure out which sides of the triangles to keep, and which ones to discard, then we can replace polygons by lines. After that, an intelligent way to take means on those lines can get us somewhere ...

My 2 cents.

Edit

One step further, getting rid of some polygons:

rul = Thread[Rule[Range@Length@j, j]];
vl = i[[2]] /. rul /. Polygon -> List; (*Polygons as vertices*)

herons[{a_, b_, c_}] := Module[{n, s}, (*Heron's area formula*)
n = Norm /@ {a - b, b - c, c - a};
s = Total@n/2;
Return[Sqrt[s (s - n[[1]]) (s - n[[2]]) (s - n[[3]])]]
]
(*Get rid of fat triangles*)
Graphics3D[{Polygon@Select[vl[[1, 1]], herons[#] < 1/100 &],
Opacity[.1], EdgeForm[None], GraphicsComplex[i[[1]], i[[2]]]}]


Our alien is getting thinner

Edit 2

Getting rid of the polygons, going to lines.

(*discard the big area polygons as before*)
ps = Select[vl[[1, 1]], herons[#] < 1/100 &];
(*get the lines of those polygons*)
lines = Flatten[Subsets[#, {2}] & /@ ps, 1];
(*discard the crunched lines in the skeletonization process*)
selLines = Select[lines, EuclideanDistance[Sequence @@ #] > 1/3 &];
Graphics3D[{Line /@ selLines, Opacity[.1], EdgeForm[None], GraphicsComplex[i[[1]], i[[2]]]}]


Now we have only 1500 lines to deal with.

-
plus one thanks –  mathew Sep 24 '12 at 15:30