I wanted to use this code to calculate the volume of a complex object (a volcanic rock) based on a surface scan. The volume returned appeared to be of the right order of magnitude (based on very gross calculations) so I'm reasonably happy. However, could the author include a little more detail in his documentation, and ideally a slightly more complicated example than just a tetrahedron to convince more cynical readers? Thank you.
I think I found a way to fix the code. On the line where Gauss' law is used:
vol = vol + abs(P(1) * N(1) * a);
The absolute value there is not correct and should be removed - when the vector field and normal vector are antiparallel the addition to the volume needs to be negative. With the current implementation you can change the volume by shifting the X coordinate. If you get rid of the absolute value everything is fine, so long as you can ensure that the normals on the surface point OUT. Although the addition of the absolute value seems intuitive, it is not how Gauss's law works.
The normals can be computed with another code on the filie exchange: computeNormalVectorTriangulation by David Gingras. David's code will give you normals that point outwards so long as the triangles are defined in a certain way... Refer to his code.
This update seems to give the right volume and remove the dependence on translation of vertices, at least in my testing.
Is there a sample dataset you can provide to verify the calculations?
I got the Area calculation to work for a unit sphere, however the Volume calculation gave me a 'Nan' result.
I too have used the Div. Theorm. to calculate volume of a triangular mesh, albeit differently. This is a clipping from my program:
Mtot=0;
rho=1;
T=(3-by-n set of your x,y,z data points)
tri=(m-by-num_elements mesh connectivity matrix)
for i=1:num_elements
% Get verticies of element
verts=tri(i,:);
% Get 3d coords of points on element
p1=T(:,verts(1));
p2=T(:,verts(2));
p3=T(:,verts(3));
A=[p1 p2 p3];
% Calculate Mass of individual tetrahedron
M(i)=rho*det(A)/6;
% Calculate total mass of volume;
Mtot=Mtot+M(i);
end
Also, I do not take the absolute value of my volume since I verify the direction of the normals to point out of the volume prior to calculation (not shown in above program)
