Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a set of random 3D data points.

How can I calculate the volume of the convex hull?

share|improve this question
    
What is your question? –  R Hall Jun 16 '12 at 21:53
    
The "area" of fig01 is not well defined. The usual definiton of "Area" applies to 2D surfaces, not to point sets –  belisarius Jun 16 '12 at 21:58
1  
For your new question, see this. –  J. M. Jun 17 '12 at 4:29

3 Answers 3

up vote 16 down vote accepted

This following code that uses TetGen we will compute the volume of the convex hull.

Needs["TetGenLink`"];
TetraMaker[pts_, surface_, TetGenString_?StringQ] := 
Module[{inInst, outInst, coords, surface1, meshElements, facets},
      inInst = TetGenCreate[];
      TetGenSetPoints[inInst, pts];
      facets = Partition[surface, 1];
      TetGenSetFacets[inInst, facets];
      outInst = TetGenTetrahedralize[inInst, TetGenString];
      coords = TetGenGetPoints[outInst];
      surface1 = TetGenGetFaces[outInst];
      meshElements = TetGenGetElements[outInst];
      {coords, surface1, meshElements}
     ];
TetrahedraVolume = Compile[{{coords, _Real, 2}, {elements, _Integer, 1}},
Block[{p},
      p = coords[[elements]];
      1/6*Abs[Det[p[[ {1, 2, 3}]] - p[[{2, 3, 4}]]]]
], RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed"
];

Some 3D data

data3D = RandomReal[{0, 10}, {65, 3}];

Compute the convex hull and then call the above function to form the tetrahedralization. Then call TetrahedraVolume to compute the volume.

{pts, surface} = TetGenConvexHull[data3D];
{coords, surface1, meshElements} = TetraMaker[pts, surface, "pqa2.8"];
Total[ TetrahedraVolume[coords, meshElements]]

505.135

You can use this to compute surface area of a triangulated 3D geometry

TriangleArea[pts_List?(Length[#] == 3 &)] := 
Norm[Cross[pts[[2]] - pts[[1]], pts[[3]] - pts[[1]]]]/2
TriangleArea[{pts[[#[[1]]]], pts[[#[[2]]]], pts[[#[[3]]]]}] & /@ surface // Total

337.121

I compute the area of the triangles separately and adding them gives me the area of the surface that defines the convex hull. enter image description here

In the above picture first you see the convex hull in black lines. The middle one shows the blue surface mesh created by TetGen during tetrahedralization. In the last one you can see the cell volumes of the tetrahedrons that discretize the volume of the convex hull in different random colors. We get the total volume by adding the volumes of these tetrahedrons.

Volume for Convex Hull of the point cloud in your data

data3D = Import["http://dl.dropbox.com/u/68983831/object.vtk", "VertexData"];
{pts, surface} = TetGenConvexHull[data3D];
{coords, surface1, meshElements} = TetraMaker[pts, surface, "pqa2.8"];
Total[ TetrahedraVolume[coords, meshElements]]

3120.05

Volume for 3D Geometry in your data

surface = Import["http://dl.dropbox.com/u/68983831/object.vtk", "PolygonData"];
{coords, surface1, meshElements} = TetraMaker[pts, surface, "pqa.8"];
Total[TetrahedraVolume[coords, meshElements]]

1622.23

enter image description here Code for Graphics

p1 = Graphics3D[{{Red, PointSize[0.03], Point[data3D]}, {Yellow, 
Opacity[.8], EdgeForm[{Thick, Black}], 
GraphicsComplex[pts, Polygon[surface]]}}, Boxed -> False, Axes -> True]; 
p2 = Graphics3D[{{Red, PointSize[0.03], 
Point[data3D]}, {EdgeForm[{Thick, Black}], FaceForm[None], 
GraphicsComplex[pts, Polygon[surface]]}, {Yellow, Opacity[.8], 
EdgeForm[Blue], GraphicsComplex[coords, Polygon[surface1]]}},
Boxed -> False, Axes -> True];
p3 = Graphics3D[{{Red, PointSize[0.03], Point[data3D]}, 
Table[With[{p = 
    RGBColor[RandomReal[{0, 1}, 3]]}, {Blend[{Lighter[Yellow, .8],
      Red, Yellow}, i/Length@meshElements], Glow[Darker[p, .5]], 
   Specularity[White, 20], Opacity[.2], EdgeForm[p], 
   GraphicsComplex[coords, 
    Polygon@Partition[meshElements[[i]], 3, 1, 1]]}], {i, 1, 
  Length@meshElements}]}, Boxed -> False, Axes -> True];
GraphicsGrid[Transpose@{{p1, p2, p3}}, Spacings -> {0, 0},ImageSize -> 600]
share|improve this answer
1  
There's a less messy way to compute surface area: Total[Flatten[Normal[GraphicsComplex[pts, Polygon[surface]]]] /. Polygon -> TriangleArea] –  J. M. Jun 17 '12 at 3:54

Taken together, PlatoManiac's answer and J.M.'s comment show how to do this. Here's code:

Needs["TetGenLink`"];
data3D = Import["http://dl.dropbox.com/u/68983831/object.vtk", "VertexData"]; 
{pts, surface} = TetGenConvexHull[data3D];
volumeContribution[{v1_, v2_, v3_}] := With[
   {n = Cross[v2 - v1, v3 - v1]}, Abs[v1.n]/2];
triangles = pts[[#]] & /@ surface;
Total[volumeContribution /@ triangles]/3

20491.

share|improve this answer
1  
I think that should be v1.n on the fifth line, not Abs[v1.n]. –  Heike Jun 17 '12 at 19:10
    
@ MarkMcClure: Thanks +1 –  mathew Jun 18 '12 at 3:32

Use RegionMeasure and ConvexHullMesh. Example:

points = RandomReal[{-1, 1}, {10, 3}]

RegionMeasure@ConvexHullMesh[points]

(* ==> 0.948199 *)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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