How to calculate volume of convex hull and volume of a 3D object

Question

I have a set of random 3D data points.

How can I calculate the volume of the convex hull?

Answer 1

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.

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

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]

Answer 2

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.

Answer 3

Use RegionMeasure and ConvexHullMesh. Example:

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

RegionMeasure@ConvexHullMesh[points]

(* ==> 0.948199 *)