24
$\begingroup$

How can we calculate the volume of a 3D object using the new-in-10 computational geometry functions?

For simple objects this works:

shuttle = ExampleData[{"Geometry3D", "SpaceShuttle"}]

Mathematica graphics

Volume@BoundaryDiscretizeGraphics[shuttle]
(* 55.5217 *)

But BoundaryDiscretizeRegion tends to fail on more complex geometries:

{deimos = ExampleData[{"Geometry3D", "Deimos"}], 
 phobos = ExampleData[{"Geometry3D", "Phobos"}]}

Mathematica graphics

BoundaryDiscretizeGraphics[deimos]

BoundaryMeshRegion::binsect: The boundary curves self-intersect or cross each other in BoundaryMeshRegion[{{-0.473472,5.76096,0.000999842},{-0.32066,5.70981,0.233476},{-0.0747657,5.65315,0.231889},<<46>>,{-2.88191,5.74069,0.489521},<<9360>>},{{},{},{Polygon[{{1,2,3},<<49>>,<<18766>>}]}}]. >>

In general, what is a good way to proceed when we need to treat an object as a solid 3D volume (not surface!) for various purposes, such as triangulation, solving a PDE inside it, calculating a volume, or using it as an integration region?

Here are two more practical test cases:

bunny = ExampleData[{"Geometry3D", "StanfordBunny"}]
horse = ExampleData[{"Geometry3D", "Horse"}]

The Stanford bunny is especially relevant because originally it was created from voxel data (a common application).

Improve this question
$\endgroup$
9
  • $\begingroup$ That is an unanswered question, hopefully not. See my questions here and here $\endgroup$
    – user31001
    Commented Oct 4, 2015 at 12:35
  • $\begingroup$ I was going to suggest a Montecarlo method, but then I saw you want to use computational geometry procedures... $\endgroup$
    – Peltio
    Commented Oct 4, 2015 at 14:57
  • $\begingroup$ @Peltio With the sorts of Monte Carlo volume computations I'm familiar with, you need to be able to tell if any point is inside or outside of the object. How would you do this here? If you can do that, that alone seems valuable enough to post. $\endgroup$
    – Szabolcs
    Commented Oct 4, 2015 at 14:59
  • 3
    $\begingroup$ Off the top of my head: slice the solid surface into intersection with parallel planes and then use one of those "is point inside a 2D polygon" test to see if it's inside. (@Szabolcs - you could not have seen my comment, so your question is more than justified). An even easier way: 3Dprint the darn satellites and put them into water to measure the displaced volume :-) $\endgroup$
    – Peltio
    Commented Oct 4, 2015 at 15:02
  • 2
    $\begingroup$ In computer graphics there are various techniques for turning a bad triangle mesh into a volumetric representation; see Jacobson et al. (2013) and the references therein. If Mathematica doesn't have a good way to do it, I would export the vertex positions and normals to a point cloud, use MeshLab to create a clean mesh, and import it back into Mathematica. $\endgroup$
    – user484
    Commented Oct 5, 2015 at 1:11

4 Answers 4

Reset to default
11
$\begingroup$

Method 1:

Not perfect solution in method 1:

deimos = ExampleData[{"Geometry3D", "Deimos"}, "VertexData"];
demon = ConvexHullMesh[deimos];
Volume@BoundaryMeshRegion[demon]

Volume = 1067.07

Or:

deimos2 = ExampleData[{"Geometry3D", "Deimos"}, "VertexData"];
demon2 = DelaunayMesh[deimos2];
RegionMeasure[demon2]

Volume = 1067.07

Method 2:

EDITED!

Better solution in method 2:

deimos = ExampleData[{"Geometry3D", "Deimos"}];
Export["Test.stl", deimos];
(* Importing and editing in Blender 3D.Exporting to Test2.stl file *)
B = Import["Test2.stl", "BoundaryMeshRegion"];
Volume[B]

Deimos Volume = 1023.75

phobos = ExampleData[{"Geometry3D", "Phobos"}];
Export["TestP.stl", phobos];
(* Importing and editing in Blender 3D.Exporting to Test2P.stl file *)
Bp = Import["TestP2.stl", "BoundaryMeshRegion"];
Volume[Bp]

Phobos Volume = 5778.21

Deimos volume is 999.78 $km^3$, Phobos volume is 5783.61 $km^3$ in Wikipedia's.

Test2.stl file and Test2P.stl file

You can also export to other files such as: *.off, *.ply,*.obj because they have a "BoundaryMeshRegion" $->$ geometry representation elements.

I'm import Test.stl file to BLENDER 3D, and unchecked "Scene Unit".Click Import STL.In "Edit Mode" Click->Mesh->Vertices->Smooth Vertex. Smooth Vertex's algorithm it's not perfect ,you have to manually edit the vertex's.Click "Object Mode"and File->Export STL->Test2.stl file.

Errors in the mesh grid enter image description hereenter image description here

An HorseImp.stl file example:

 horse = ExampleData[{"Geometry3D", "Horse"}];
 Export["Horse.stl", horse];
 (* Importing and editing in Blender 3D.Exporting to HorseImp.stl file *)
 hor = Import["HorseImp.stl", "BoundaryMeshRegion"]
 {RegionCentroid[hor], Volume[hor], Area[RegionBoundary[hor]]}

Calculate Centroid,Volume,Surface Area:

{{0.536307, 0.0310618, -0.0876894}, 0.000256564, 0.0343025}

Then You can for example discretizes the region.

BoundaryDiscretizeRegion[hor, AccuracyGoal -> 0.1, 
MeshCellStyle -> {{1, All} -> Opacity[0.4, Orange]}]

enter image description here

Improve this answer
$\endgroup$
4
  • 8
    $\begingroup$ The big problem with these two approaches is that both of them measure the volume of the convex hull, not the volume of the object. For the horse and bunny it would not even be remotely close. $\endgroup$
    – Szabolcs
    Commented Oct 4, 2015 at 13:58
  • $\begingroup$ @Szabolcs. It applies only to the first method. $\endgroup$
    – Mariusz Iwaniuk
    Commented Oct 6, 2015 at 8:39
  • $\begingroup$ Thank you! Indeed the second method is accurate, but it involves having to fix the mesh manually. My actual meshes are quite badly behaved and I'm still playing with various software (Blender or others) to try to automate the process. $\endgroup$
    – Szabolcs
    Commented Oct 6, 2015 at 9:09
  • $\begingroup$ @Szabolcs Largely because of your first comment to this post I decided to try answer the question. $\endgroup$
    – Anton Antonov
    Commented Oct 7, 2015 at 0:19
10
$\begingroup$

A very simple way of calculating the volume of a non-convex polyhedron appears in equation 15 of the article by Lien and Kajiya, "A Symbolic Method for Calculating the Integral Properties of Arbitrary Nonconvex Polyhedra", Computer Graphics & Applications, Oct 1984. If all faces are triangular, this does the trick:

deimos = ExampleData[{"Geometry3D", "Deimos"}, "PolygonObjects"];
Total[Abs[deimos /. Polygon -> Det]]/6
(* 1025.7561 *)

Here, deimos is a list of triangular Polygons. The argument of each is a three-element list of 3D Cartesian vertex coordinates. In the second line, Mathematica's Polygon[] function is replaced with Det[] to evaluate the determinant of the $3\times 3$ matrix formed by the triangle's vertex coordinates. Finally, the total of the absolute values of all such determinants is divided by six.

For Phobos, the corresponding result is 5790.5116.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ That should be Abs[Total[...]], I'm afraid (see eq. 25). Otherwise you will get the wrong answer for polyhedra that do not contain the origin. $\endgroup$
    – user484
    Commented Oct 6, 2015 at 10:48
6
$\begingroup$

Here is another solution of the first, specific question about the volumes of Phobos and Deimos. The idea is to use the polygons provided by example data:

polygonObjs = 
  ExampleData[{"Geometry3D", "Deimos"}, "PolygonObjects"];
polygonObjs // Length

(* Out[74]= 18816 *)

Graphics3D[RandomSample[polygonObjs, 2000], Axes -> True]

enter image description here

Suppose we have shifted those polygons so all their points have positive z-coordinates. Then we can estimate the Deimos volume by summing the volumes of the truncated prisms made by these polygons and their projections on the xy plane while taking into account which polygons are on top and which are under.

First, we have to figure out which polygons are on top an which are under. We can use the "VertexNormals" data for this.

polygonData = ExampleData[{"Geometry3D", "Deimos"}, "PolygonData"];
polygonData // Length

(* 18816 *)

vNormals = ExampleData[{"Geometry3D", "Deimos"}, "VertexNormals"];
vNormals // Length

(* 9410 *)

pNormals = Map[Mean[vNormals[[#]]] &, polygonData];
pNormals // Length

(* 18816 *)

polygons = polygonObjs[[All, 1]];

Here is a sample of polygons and their corresponding vectors (note the shift in the z-axis by 10):

sInds = RandomSample[Range[Length[polygons]], 600];
Graphics3D[
 Translate[#, {0, 0, 10}] &@{ Polygon /@ polygons[[sInds]], Gray, 
   Opacity[0.3], Thin, Arrowheads[0.005], 
   MapThread[
    Arrow[{Mean[#1], #2}] &, {polygons[[sInds]], pNormals[[sInds]]}]},
  Axes -> True]

enter image description here

Next we define a function (taken from here) to find the area of a projected (2D) polygon:

PolygonArea = 
 Compile[{{v, _Real, 2}}, Block[{x, y}, {x, y} = Transpose[v];
   Abs[x.RotateLeft[y] - RotateLeft[x].y]/2]]

Now we are ready to compute the volume by multiplying the projections` areas with hight means and normals direction signs:

Dot[Map[PolygonArea[#[[All, 1 ;; 2]]] &, polygons], 
    (Mean[#[[All, 3]]] + 10) & /@ polygons*
    Sign[pNormals[[All, 3]]]]

(* 1026.06 *)

( The result seems closer to the one pointed out in Wikipedia than the solution using ConvexHullMesh.)

The computation of the volumes of the truncated prisms can be refined to be more precise -- note that I just take the means of the z-coordinates of the shifted polygons. (We can use precise formulas or Integrate.)

This idea is very similar to the mathematical formulas for using type-2 surface integrals to compute volumes, but I was too lazy to implement that.

For Phobos this solution gives 5788.04 which is very close to the one given in Wikipedia's article on Phobos, 5783.61. Using the ConvexHullMesh solution we get again a more overestimated volume:

phobos = ExampleData[{"Geometry3D", "Phobos"}, "VertexData"];
phobos = ConvexHullMesh[phobos];
Volume@BoundaryMeshRegion[phobos]

(* Out[147]= 5987.53 *)
Improve this answer
$\endgroup$
6
  • $\begingroup$ Note the error message cited in the question: "BoundaryMeshRegion::binsect: The boundary curves self-intersect or cross each other." This gives possible explanation why your method gives overrated estimate. $\endgroup$
    – Alexey Popkov
    Commented Oct 5, 2015 at 0:16
  • $\begingroup$ @AlexeyPopkov Is this a comment for the other solution? I am not seeing any messages with this solution and I am not using BoundaryMeshRegion. $\endgroup$
    – Anton Antonov
    Commented Oct 5, 2015 at 0:43
  • 1
    $\begingroup$ I mean that polygons may actually intersect with each other what can lead to overestimate of the volume. $\endgroup$
    – Alexey Popkov
    Commented Oct 5, 2015 at 0:57
  • 1
    $\begingroup$ Nice idea for computing the volume. One could simply apply the formula that arises from the divergence theorem, $V = \frac13\sum\limits_{\text{face $i$}} \mathbf p_i\cdot\mathbf n_iA_i$, where $\mathbf p_i$ is any point on the face. However, this does not solve the general problem of obtaining an actual volumetric Region "for various purposes, such as triangulation, solving a PDE inside it, ... or using it as an integration region". $\endgroup$
    – user484
    Commented Oct 5, 2015 at 0:58
  • $\begingroup$ @Rahul Right, I made those both points in the first version of the post. $\endgroup$
    – Anton Antonov
    Commented Oct 5, 2015 at 1:02
3
$\begingroup$

I will show an approach that makes use of this α shapes code. The caveat here is that one needs to tune a parameter to achieve the "best" shape. But this could be achieved programmatically, if desired, by seeing how some property (e.g. area or volume) of the object changes as the α parameter is changed slightly. A region where that property doesn't change much as long as it's not close to it's Delaunay triangulation value, should give a close approximation of the true value. Let's import some objects:

{deimos = ExampleData[{"Geometry3D", "Deimos"}], 
 phobos = ExampleData[{"Geometry3D", "Phobos"}],
 bunny = ExampleData[{"Geometry3D", "StanfordBunny"}]}

Mathematica graphics

First, we convert these images into regions using DiscretizeGraphics

reg = DiscretizeGraphics @ deimos

Mathematica graphics

We then collect the coordinates of this region and pass them into the alphaShapes function to do its job:

coord = MeshCoordinates @ reg;
Volume @ alphaShapes[coord, 4.82]

989.749

This is impressively close to the actual value ($999.78$). Note that if you increase the α parameter slightly to 4.83 you get a huge jump to $1036.32$ and it stops changing for a while . Below I show what the triangulated output from alphaShapes looks like. Note how it captures the true shape of deimos.

Mathematica graphics

For phobos, note how close we get to the actual volume ($5783.61$). In the images below I've placed side by side the discretized image (on the left) and the result from alphaShapes triangulated (on the right).

{alphaShapes[coord2, 8.69] // Volume,
alphaShapes[coord3, 0.01] // Volume}

{5782.64, 0.0000940232}

Here coord2 and coord3 were obtained the same way coord was obtained above.

Mathematica graphics Mathematica graphics

Improve this answer
$\endgroup$

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.