How to calculate the volume of trisurf?

Hi,I have a question.

At present, I use random points to create a sphere, and then I apply Crust (https://www.mathworks.com/matlabcentral/fileexchange/63730-surface-reconstruction-from-scattered-points-cloud) to surface it, and then I tried to calculate the volume, but the gap between the calculated volume and the actual volume is too large. I think it may be that my calculation method is wrong. I would like to ask if everyone can provide me with some opinions and methods?

Thanks.

[Matt J]

but the gap between the calculated volume and the actual volume is too large.

What is the gap that you see?

[Chen Xinzhuo]

Sphere r: 3 Actual volume calculated as 113,

But my calculated totalVolume is 0.4124.

Answer 1

[John D'Errico]

If the points are on a sphere, thus a convex object, then a convex hull is what you want, not Crust anyway. Crust will allow the result to be non-convex. And that is arguably why you get a serious under-estimate. Computing the volume, lets see, the simplest tool (that I can think of offhand) that computes a volume directly is an alphaShape. (No, it is not difficult to write. But why bother?)

X = randn(1000,3);
X = normalize(X,2,'norm');

T = alphaShape(X,inf);
V = volume(T)

V = 4.1393

The corresponding sphere volume would be

truth = 4*pi/3

truth = 4.1888

That seems reasonable for 1000 points.

abs(truth - V)/truth*100

ans = 1.1813

So within roughly 1% for the volume, for 1000 randomly sampled points on the surface of a sphere.

Odds are there is an easy way to compute the volume of the convex hull directly, but there is no need to look more deeply once I have a tool that does what I want.

[Chen Xinzhuo]

Hi John, thank you for your answer,

I have tried AlphaShape and calculated the volume, which is very similar to the actual calculation.

But for the complex point cloud AlphaShape can't perfectly fit the surface, and Crust can do this very well, so I choose Crust to find a way to calculate the volume.

[John D'Errico]

That is fine. If you want to use Crust, then use it. But then don't be surprised at the result.

[Chen Xinzhuo]

Is the way I calculated the volume correct? I'm not sure about this either.

[John D'Errico]

You appear to be using Heron's formula for the area of a triangle.

a = sqrt(sum((p(t(:,2),:) - p(t(:,1),:)).^2, 2));
b = sqrt(sum((p(t(:,3),:) - p(t(:,2),:)).^2, 2));
c = sqrt(sum((p(t(:,1),:) - p(t(:,3),:)).^2, 2));
s = (a + b + c) / 2;
area = sqrt(s .* (s - a) .* (s - b) .* (s - c));

Then for some unknown reason, you multiply it by the third point in each triangle on the surface of the object.

That will not compute the volume of a tetrahedron.

totalVolume = sum(sum(bsxfun(@times, p(t(:,3),:), area))) / 3

If I recall, you can compute the volume of a tetrahedron by ultiplying the area of a face of the tetrahedron, then multiply that by the height above that facet to the 4th vertex. Then you would divide by 6 I think.

So perhaps you were thinking to connect each facet to the point (0,0,0). But that is not the height of the tetrahedron. So I'm not sure what you were thinking there.

[John D'Errico]

For example, let me compute the volume of a convex hull, doing it the hard way. I won't use Heron's formula though.

XYZ = randn(1000,3);
XYZ = normalize(XYZ,2,'norm');

T = convhulln(XYZ);
Ctr = mean(XYZ,1); % Aan interior point in the hull. I could use (0,0,0) here too

% the cross product between vectors for two sides of a triangle, gives us
% both a vector that is normal to the facet, and has magnitude twice the
% area of the triangle.
Nrml = cross(XYZ(T(:,3),:) - XYZ(T(:,1),:),XYZ(T(:,2),:) - XYZ(T(:,1),:));

% the facet areas
Areas = sqrt(sum(Nrml.*Nrml,2));

% scale the normals to now be unit normal vectors
Nrml = Nrml./Areas;

% Now, we connect each facet to the center point, creating a set of
% tetrahedra. If the boundary surface is not a convex one, we would need to
% be careful below, using the signed areas from those cross products. I'm
% not going to bother here, since I know this is a fully convex domain.
%
% The distance to the center point, in the DIRECTION of the normal vector
% gives us the "height" of the tetrahedron.
D = sum(Nrml.*(Ctr - XYZ(T(:,1),:)),2);

% Multiply each of those distances by the area of each factet. Divide by 6.
% Add them up.
format long g
Vol = sum(Areas.*abs(D))/6

Vol = 
          4.13680491082305

And as expected, the volume is roughly

4*pi/3

ans = 
          4.18879020478639

It must be slightly less than 4*pi/3, because the convex hull lies entirely inside the sphere. Compare to the alpha shape volume. It should be the same, at least to within floating point crap.

volume(alphaShape(XYZ,inf))

ans = 
          4.13680491082305

Note that the above computation, IF you have a release older than R2016b will fail, and you must use bsxfun instead. As well, IF your domain is not a convex one, then you need to be slightly more careful, using the direction of the normal vectors more carefully.