2

I am using geom3d to model 3D geometric shapes. Specifically, I am interested to calculate a volume of the ellipsoid after slicing it with an arbitrary plane. I am using the following code (here for simplicity I will use a sphere and cut it in the center).

[x, y, z] = sphere(30);
[f, v, c] = surf2patch(x, y, z, z, 'triangles');
meshVolume(v, f)

This calculates a volume of a sphere as expected (4/3*pi).

Now, I slice it in half.

P = createPlane([0, 0, 0], [0, 0, 1]);
[v1 f1] = clipConvexPolyhedronHP(v, f, P);
drawMesh(v1, f1, 'facealpha', 0.5)

Now if I calculate the volume, I get incorrect value (significantly smaller than half of the initial volume).

 meshVolume(v1, f1)

What is the reason for the incorrect result?

Sasha
  • 5,783
  • 8
  • 33
  • 37

2 Answers2

2

(I have already answerd on Mathworks File Exchange page, but I think it can be useful to put it also here)

Using determinant without 'abs' is the correct way of computing volume of meshes. By doing that way you can handle non convex meshes, and even polyhedra with holes. But the faces of the mesh must be oriented consistently (usually outwards of the mesh). A similar topic is given here

Using the absolute value of the determinant will give correct result only for convex polyhedra. So I would use it only as a quick hack for specific application.

Regarding the original question, there was a bug in the computation of the clipped mesh, resulting in faces with inconsistent orientation, and therefore in wrong resulting volume. I have updated it so it should give better results.

(note: I am the author of the geom3d package)

Community
  • 1
  • 1
dlegland
  • 151
  • 5
0

Actually found a bug in the implementation of the meshVolume function. The function divides the volume into tetrahedra, which are 3-dimensional solids that have 4 vertices, 6 edges, and 4 triangular faces. Then the total volume of the body is the sum of the volumes of all tetrahedra. The code computes the volume by calculating the determinant of the vertices.

vols(i) = det(tetra) / 6;

It turns out that sometimes this determinant can be negative due to a different order of the rows (vertices coordinates).

The correct formula should include absolute value of the determinant.

vols(i) = abs(det(tetra) / 6);

That solves the problem. I've reported the bug on geom3d webpage.

For more details about tetrahedra see here

Sasha
  • 5,783
  • 8
  • 33
  • 37