How to calculate the volume of a convex hull?

Question

Given a spatial curve represented by a parametric equation, is it possible in Mathematica 9 to calculate symbolically (or at least numerically) the volume of its convex hull?

Answer 1

Here is a way to do it numerically:

Reap the points from a parametric plot:

pf = {Cos[u], Sin[u] + Cos[v], Sin[v]};
data = Reap[
    ParametricPlot3D[
     Sow[pf], {u, 0, 2 \[Pi]}, {v, -\[Pi], \[Pi]}]][[2, 1]];

Get all the numeric values, pack the data and delete the duplicates:

pts = Cases[data, {_?NumericQ, _?NumericQ, _?NumericQ}];
pts = DeleteDuplicates[Developer`ToPackedArray@N[pts]];

To look at the points:

Graphics3D[Point[pts], Boxed -> False]

Load TetGenLink and compute and visualize the convex hull:

<<TetGenLink`
ch = TetGenConvexHull[pts];
Graphics3D[GraphicsComplex[ch[[1]], Polygon[ch[[2]]]], Boxed -> False]

Compute the delaunay mesh of the convex hull, write a function to compute the volume of a tetrahedron, apply it and total the volume of the tetrahedra:

dt = TetGenDelaunay[pts];
tetVol = Compile[{{coords, _Real, 2}, {incidents, _Integer, 1}}, 
   Block[{a, b, c, d},
    {a, b, c, d} = coords[[incidents]];
    1/6*Det[Transpose[{a, c, b}] - d]
    ], RuntimeAttributes -> Listable];

Total[tetVol @@ dt]
(* 12.5 *)

Sanity check: for a sphere we get:

pf = {Sin[u]*Sin[v], Cos[u]*Sin[v], Cos[v]};
4.138` - 4./3*\[Pi]*1^3
(*-0.05*)

which seams reasonable.

Update:

As suggested by Daniel Lichtblau, this also works for parametric curves:

pf = {(t^2 - 1)/(t^2 + 1), 
   2 t/(t^2 + 1), ((t^2 - 1)/(t^2 + 1))^2 - (2 t/(t^2 + 1))^2};
data = Reap[ParametricPlot3D[Sow[pf], {t, -10, 10}]][[2, 1]];
pts = Cases[data, {_?NumericQ, _?NumericQ, _?NumericQ}];
pts = DeleteDuplicates[Developer`ToPackedArray@N[pts]];
Graphics3D[Point[pts], Boxed -> False]

<< TetGenLink`
ch = TetGenConvexHull[pts];
Graphics3D[GraphicsComplex[ch[[1]], Polygon[ch[[2]]]], Boxed -> False]

And what I find most amazing is the suggestion from george2079:

shoelace3rd = Compile[{{coords, _Real, 2}, {tris, _Integer, 2}},
   -Total[Det[coords[[#]]] & /@ tris]/6
   ];
shoelace3rd @@ ch
(* 3.12 *)

Answer 2

Following ruebenko's approach, after computing the convex hull you can much faster directly get the volume without triangulating the volume as 1/6 Sum[ Det [ triangle vertices ] ]:

Do[ -Total[
    Det[{ch[[1, #[[1]]]], ch[[1, #[[2]]]], ch[[1, #[[3]]]]}] & /@ 
       ch[[2]]]/6 , {10}]
(* 12.5077 *)

This is the 3d shoelace formula:

https://stackoverflow.com/questions/14710682/algorithm-for-partially-filling-a-polygonal-mesh/14713130#14713130

http://en.wikipedia.org/wiki/Shoelace_formula

I must say I've searched pretty good for a link to derivation/proof of this useful thing and come up empty.

Ok the derivation was not so hard to work out.. Note implict here is an assumption that the vertices of each triangle are ordered consistantly CW or CCW. Also this does not require convexity, only whatever connectivity restrictions follow from the divergance theorem.

Answer 3

In Version 10, once the points have been obtained as per user21's approach, we can tetrahedralize them directly using DelaunayMesh

pf = {Cos[u], Sin[u] + Cos[v], Sin[v]};
pp = ParametricPlot3D[pf, {u, 0, 2 Pi}, {v, -Pi, Pi}]

data = Reap[ParametricPlot3D[Sow[pf], {u, 0, 2 Pi}, {v, -Pi, Pi}]][[2, 1]];
pts = Cases[data, {_?NumericQ, _?NumericQ, _?NumericQ}];

Then:

del = DelaunayMesh[pts]

And the volume is computed easily using RegionMeasure

RegionMeasure[del]

12.5077258

Answer 4

Here is another way that uses the Graphics object directly:

gr = ParametricPlot3D[{Cos[u], Sin[u] + Cos[v], Sin[v]}, {u, 0, 2 Pi}, {v, -Pi, Pi}]

We discretize the graphics using DiscretizeGraphics

mr = DiscretizeGraphics[Normal[gr /. (Lighting -> _) :> Lighting -> Automatic]]

We compute the convex hull

hull = ConvexHullMesh[MeshCoordinates[mr]]

And finally, we get the volume using RegionMeasure

RegionMeasure[hull]

12.5077258

OR Volume

Volume[hull]

12.5077258

As a bonus we can get the surface area of the ParametricPlot object:

RegionMeasure[mr, 2]

33.1870442

Same result can be obtained using Area

Area[mr]

33.1870442