I have modified my question to a simple version. I have vertex clusters, polygons and a list of center points.

Attached pts and polygons of 3d tube are:

 pts =Import["https://dl.dropboxusercontent.com/u/68983831/npts.txt","Table"]; 
 polys= Import["https://dl.dropboxusercontent.com/u/68983831/npolys.txt","Table"];

These are approximate center points:

 cpts = Import["https://dl.dropboxusercontent.com/u/68983831/cpts.txt","Table"];

 Show[Graphics3D[{Opacity[.1], EdgeForm[None], GraphicsComplex[pts, Polygon@polys]},
  Boxed -> False], ListPointPlot3D[cpts]]

enter image description here

How to measure cross-sectional area, major & minor axis length as move along the center points? or may be without center points?

  • 3
    $\begingroup$ Here's an outline of a method. For a given center point do the following. (1) Find several (at least 5) boundary points. (2) Find the cylinder that best fits these points. (3) Find cross-sectional area. If you want to constrain the center and perhaps direction vector of the axis then you can use fewer points and you'll want to constrain the optimization problem for best-fitting cylinder. $\endgroup$ Commented Jul 23, 2014 at 21:54

3 Answers 3


Because when in Europe... Monte-Carlo!

This should work with any shaped tube as long you are given the center points, and a certain resolution, in your case, should be the size of the polygon that make up the surface, or for example the minimum distance between point along the surface. Let's say that distance is dis.

First we make a path along the inside of the tube this is not very important, but i might be necessary to in case the tube bends.

newcpts = cpts[[#]] & /@ (Last@FindShortestTour[cpts]);

We order the points so they are one after in order in a particular direction. It might break down if the bends are too tight. But it is not too important, and it might be completely unnecessary for what we will do.

Next we will make a estimate of the tangent vectors along that curve.

myderivate[x__] :=
  Mean@(Subtract @@@ (Reverse /@ Partition[x, 2, 1])/(Length[x] - 1));
mytangcurve = myderivate /@ Partition[newcpts, 3, 1, {2, 2}, {}];

Next we compute the normals of the surface. We use the same from this answer

newareavectors = 
  SortBy[Flatten[#, 1] &@
      Map[Function[{y}, {y, 
        If[#.(First@Nearest[newcpts, pts[[y]]] - pts[[y]]) >= 0, -#, #] &[
          Cross @@ ((pts[[#]] - pts[[y]]) & /@ DeleteCases[x, y])]}], x]], polys], First];

{thepoints, thenorms} = 
  Thread[({Mean[#[[1]]], Mean[#[[2]]]} & /@ (Thread /@ GatherBy[newareavectors, #[[1]] &]))];
myquicknormals = #[[2]] & /@ SortBy[Thread[{thepoints /. 
  Rule @@@ (Thread[{#, Range[Length@#]}] &
    [Sort@DeleteDuplicates@Flatten[polys]]), thenorms}], First];

First we create a box of a large enough width, but of thickness dis.

randpointsofmy = 
  Thread[{RandomReal[{-1, 1}, 10000], RandomReal[{-1, 1}, 10000], 
    RandomReal[{-dis/2., dis/2.}, 10000]}];

dis = 0.1, seems to be right for your polygons.

Now we move this box to each point in the center and check the percentage of points that are inside.

myTransform[points_, translation_, planevector_] := 
  Dot[RotationMatrix[VectorAngle[{0, 0, 1}, planevector], 
    Cross[{0, 0, 1}, planevector]], #] & /@ ( # + translation & /@ points)

and create a boolean to check whether points are inside the object

isitinsideQ[normals_, point_] := 
   And @@ ((MapIndexed[VectorAngle[#1, pts[[Last@#2]] - point] < Pi/2. &, normals]));

There are more efficient ways out there, but you got to make sure the points in the surface polygons are all oriented in the right direction (clockwise or anticlowise). This should work either way, but you need the inside points for reference to create the normals.

So we can try for one point say newcpts[[1]]

numofpt = N@Count[isitinsideQ[myquicknormals, #] & /@ 
  myTransform[randpointsofmy, newcpts[[1]], mytangcurve[[1]]], True]/(Length@randpointsofmy)

The ration of points inside versus outside is just

$$ \frac{V_{cylinder}}{V_{box}} =\frac{\pi r^2}{4}$$

since we our cases the widths of the box were 2. This leads to a radius of

r = First@Select[rr /. Solve[numofpt/1 == (Pi*rr^2)/4, rr], # > 0 &]

Or an area of π*r^2:


Now you can do that for each point. I will leave the uncertainty analysis to you.

And just to show you we are getting the right points:

 Show[Graphics3D[{Opacity[.4], EdgeForm[None], 
      GraphicsComplex[pts, Polygon@polys]}, Boxed -> True, Axes -> True],
      Graphics3D[{ Blue, 
      Point[myTransform[randpointsofmy, newcpts[[1]], 
      Graphics3D[{Red, Point[Cases[
      myTransform[randpointsofmy, newcpts[[1]], mytangcurve[[1]]], 
      x_ /; isitinsideQ[myquicknormals, x]]]}]]

Mathematica graphics

  • $\begingroup$ The volume ratios are really unnecessary, it just gives a more accurate answer. You can just measure the area of the points that are inside the shape after rotating the back to the flat coordinate system. $\endgroup$
    – lalmei
    Commented Jul 26, 2014 at 21:37
  • $\begingroup$ I edited your code quite heavily, most of all, I added a SeedRandom so anyone testing your code should have the same values and figure. If you feel unhappy with the changes, feel free to rollback. $\endgroup$
    – Öskå
    Commented Jul 27, 2014 at 10:11
  • $\begingroup$ @Öskå ok, thanks for the higher resolution image too! I am still editing it offline for now, and may change in the future. $\endgroup$
    – lalmei
    Commented Jul 28, 2014 at 10:22
  • 1
    $\begingroup$ @lalmei No problem :) Just bear in mind that it's a good habit to add seeds when using Random* so everyone has the same results and figures :) $\endgroup$
    – Öskå
    Commented Jul 28, 2014 at 10:24

Here, I'll give an approximate solution also based on the Monte-Carlo approach but using very convenient functions introduced in V10. Here it is:

pts = Import["https://dl.dropboxusercontent.com/u/68983831/npts.txt", "Table"];
polys = Import["https://dl.dropboxusercontent.com/u/68983831/npolys.txt", "Table"];
cpts = Import["https://dl.dropboxusercontent.com/u/68983831/cpts.txt", "Table"];

Here's the cylinder

gr = Graphics3D[GraphicsComplex[pts, Polygon@polys], Axes -> True]

Mathematica graphics

We tetrahedralize the cylinder using the points:

del = DelaunayMesh[pts];

We visualize it for our sanity:

HighlightMesh[del, {Style[2, Opacity[0.2], Blue], 
  Style[1, Directive[Thickness[0.0001], Orange]], 
  Style[0, Directive[PointSize[0.02], Green]]}]

Mathematica graphics

Yes, the choice of colors are awful :)

Using one of the center points for example (we use cpts[[5]]), we generate a plane that passes through the point and intersects the cylinder. It's easy to pick two other points that lie on the plane (they must have same z-axis in this case, where the length of the cylinder is along the z-axis so our plane is an x-y plane):

ptxy1 = cpts[[5]]; (* first point on the plane *)
ptxy2 = {cpts[[5, 1]], RandomReal[{0, 1}], ptxy1[[3]]};
ptxy3 = {RandomReal[{0,1}], ptxy2[[2]], ptxy1[[3]]};

Now we create two vectors that lie on the plane and take their Cross product to get a normal to the plane, from which we can get an equation for the plane.

vecA = ptxy2 - ptxy1;
vecB = ptxy3 - ptxy1;
normal = Cross[vecA, vecB];

And the equation of the plane becomes:

eqn = normal.({x, y, z} - ptxy1) == 0;

Lets create the plane Region:

plane = ImplicitRegion[eqn, {x, y, z}];

z = ptxy1[[3]];

We discretize the region:

preg = DiscretizeRegion[plane, {{-1, 1}, {-0.5, 1}, {-z, z}}];


Let's generate points on the plane:

ptsonplane = Table[{RandomReal[{-1, 1}], RandomReal[{-0.5, 1}], cpts[[5, 3]]}, {10000}];

We now find what points fall inside and outside the cylinder:

rm = RegionMember[del]

Mathematica graphics

inreg = rm[ptsonplane];

pinreg = Pick[ptsonplane, inreg] (* Points in the cylinder *)

poutreg = Complement[ptsonplane, pinreg] (* Points that fall outside the cylinder *)

Let's visualize it:

Show[{plane, Graphics3D[{Red, Point[pinreg], Blue, Point[poutreg]}], 
  HighlightMesh[DiscretizeGraphics[gr], Style[2, Opacity[0.1]]]}]

Mathematica graphics

Now, we find the cross-sectional area:

   x /. NSolve[Length[pinreg] / Length[inreg] == x / RegionMeasure[plane], x]


You can do the same for other center points.

Note: I've used Delaunay tetrahedralization here since it's a cylinder. This method will not work if you have a curved cylinder or other similar non-convex shapes. I'm still thinking of a more general approach.


This might be a bit naive but if you are approximating a cylinder, why wouldn't you simply do the following?

pts = Import["https://dl.dropboxusercontent.com/u/68983831/npts.txt", "Table"];
cpts = Import["https://dl.dropboxusercontent.com/u/68983831/cpts.txt", "Table"];
polys = Import["https://dl.dropboxusercontent.com/u/68983831/npolys.txt", "Table"];
fpts = (First@GatherBy[pts, Last])[[1 ;; All, 1 ;; 2]]; (* taking one slice *)
c = fpts[[ConvexHull[fpts]]];

I'm simply taking advantage of the fact that ConvexHull "eats up" the inner points, so I have:


Mathematica graphics

From that I can calculate the radius of the circle and extract the bottom and top points of your "tube":

r = Mean[EuclideanDistance[Mean@c, #] & /@ c];
e = With[{l = DeleteDuplicates[Last /@ pts]}, {Min@l, Max@l}];

giving me a rough approximation of your tube:

  Graphics3D[{Opacity[.5], EdgeForm[None], 
    GraphicsComplex[pts, Polygon@polys]}, Boxed -> False], 
  Graphics3D[{[email protected], Cylinder[{Append[Mean@c, First@e], Append[Mean@c, Last@e]}, r]}]]

Mathematica graphics


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.