How to measure cross sectional area of Tube

Question

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

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

Answer 1

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.

ClearAll[myderivate];
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[{x}, 
      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.

dis=0.1;
SeedRandom@1;
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)

0.1838

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

0.483758

Or an area of π*r^2:

0.7352

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],
      ListPointPlot3D[newcpts], 
      Graphics3D[{ Blue, 
      Point[myTransform[randpointsofmy, newcpts[[1]], 
      mytangcurve[[1]]]]}], 
      Graphics3D[{Red, Point[Cases[
      myTransform[randpointsofmy, newcpts[[1]], mytangcurve[[1]]], 
      x_ /; isitinsideQ[myquicknormals, x]]]}]]

Answer 2

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]

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

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}}];

Monte-Carlo

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]

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

Now, we find the cross-sectional area:

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

{0.854673757}

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.

Answer 3

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 *)
Needs["ComputationalGeometry`"]
Graphics`Mesh`MeshInit[];
c = fpts[[ConvexHull[fpts]]];
PolygonArea[c]

0.727818

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

Graphics@Polygon@c

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:

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