Calculate the volume of shape descibed by 3d points.

Given a set of 3d data points, without using the Image Processing Toolbox, what is the best method for working put the volume of the object described by the points?

Answer 1

[Yu Jiang]

Hi Lucy

The volume of that object can be approximated by the volume of its convex hull. In mathematics, the convex hull or convex envelope of a set X of points in the Euclidean plane or Euclidean space is the smallest convex set that contains X. See the following link for more details

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

The Matlab function convhull can be used to find the convex hull of a given dataset and can return respectively the area or the volume of a 2D-Polygon or of a 3D-Polyaedrons. You can find more information on this function at the following address:

http://www.mathworks.com/help/matlab/ref/convhull.html

Here is one example illustrating how you can use this command to compute a volume:

clear all
N = 100;
x = linspace(0,1,N);
y = linspace(0,1,N);
Z = peaks(N)
[X,Y] = meshgrid(x,y);
[TriIdx, V] = convhull(X,Y,Z)
trisurf(TriIdx, X, Y, Z)

The variable V is the volume of the convex hull. Also, the function trisurf is used here to plot the convex hull. See the following link for more details on trisurf.

http://www.mathworks.com/help/matlab/ref/trisurf.html

-Yu

[Ali Gh]

convhull has a serious problem that it will abstract the data and for complicated shapes it really hurts ...

[Caleb Williams]

Agreed. I think you could probably do better with a sort of Monte Carlo integration method. Find the min and max of your X, Y, and Z coordinates and form a cube from those values (those are your constraints). Then scatter a bunch of points inside the cube using the "rand" function. If A% of the points fall within the 3-D solid you want to find the volume of, then the volume of the 3-D solid is A% of the volume of the cube. Does this help?

[Lateef Adewale Kareem]

I just created a function that computes the volume given the vertices.

Here are examples.

download the function Vertices2Volume from file exchange.

Vertices2Volume - File Exchange - MATLAB Central (mathworks.com)

% Volume of a cylinder of radius 2, and height 3
r = 2; h = 3;
[X, Y, Z] = cylinder([0,r,r,0], 32);
Z([1,2],:) = 0; Z([3,4],:) = h;
figure(color = 'w')
% using matlab's convex hull function
[tri, vol1] = convhull(X, Y, Z);
subplot(1,2,1)
trisurf(tri, X, Y, Z, FaceColor='cyan')
title('Convex Hull')
axis equal
% new Vertices to volume function
vol2 = Vertices2Volume(X, Y, Z);
subplot(1,2,2)
surf(X, Y, Z, FaceColor='cyan'); 
title('Actual object')
axis equal;

% analytical solution
vol3 = pi*r^2*h;
error1 = 100*abs(vol1-vol3)/vol3;
error2 = 100*abs(vol2-vol3)/vol3;

disp([vol1, vol2, vol3]);

   37.4573   37.4573   37.6991

disp([error1, error2]);

    0.6413    0.6413

% Volume of a sphere of radius 2,
r = 2;
[X, Y, Z] = sphere(32);
X = r*X; Y = r*Y; Z = r*Z;
figure(color = 'w')
% using matlab's convex hull function
[tri, vol1] = convhull(X, Y, Z);
subplot(1,2,1)
trisurf(tri, X, Y, Z, FaceColor='cyan')
title('Convex Hull')
axis equal
% new Vertices to volume function
vol2 = Vertices2Volume(X, Y, Z);
subplot(1,2,2)
surf(X, Y, Z, FaceColor='cyan'); 
title('Actual object')
axis equal;

% analytical solution
vol3 = 4*pi*r^3/3;
error1 = 100*abs(vol1-vol3)/vol3;
error2 = 100*abs(vol2-vol3)/vol3;

disp([vol1, vol2, vol3]);

   33.2153   33.2153   33.5103

disp([error1, error2]);

    0.8805    0.8805

% Volume of a Torus with minor radius 1 and major radius 5
r = 1.5; R = 4; u = linspace(0, 2*pi, 33); 
z = r*sin(u); x = R + r*cos(u); 
v = linspace(0, 2*pi, 65)'; 
X = cos(v)*x; Y = sin(v)*x;
Z = ones(size(v))*z; 
figure(color = 'w')
% using matlab's convex hull function
[tri, vol1] = convhull(X, Y, Z);
subplot(1,2,1)
trisurf(tri, X, Y, Z, FaceColor='cyan')
title('Convex Hull')
axis equal
% new Vertices to volume function
vol2 = Vertices2Volume(X, Y, Z);
subplot(1,2,2)
surf(X, Y, Z, FaceColor='cyan'); 
title('Actual object')
axis equal;

% analytical solution
vol3 = 2*pi^2*r^2*R;
error1 = 100*abs(vol1-vol3)/vol3;
error2 = 100*abs(vol2-vol3)/vol3;

disp([vol1, vol2, vol3]);

  252.6483  176.2302  177.6529

disp([error1, error2]);

   42.2146    0.8008

function vol = Vertices2Volume(X, Y, Z)
% This function computes the volume of a 3D shape given its vertices X, Y,
% Z, povided X, Y, Z are such that can be plotted using surf function. 

%% Input
 % X = X-Coordinates
 % Y = Y-Coordinates
 % Z = Z-Coordinates
%% Output
 % vol = volume exclosed by the figure

%% Copyright
% Lateef Adewale Kareem 2023.

%% Solution
    % extract size and set vol to zero;
    [M, N] = size(X); vol = 0; 
    for m = 1:M-1
        for n = 1:N-1
            % convert subscripts of each quadrilateral to indices
            IND = sub2ind([M, N],[m, m, m+1, m+1],[n, n+1, n+1, n])';
            % split into indices of 2 triangles
            j = IND([1,2,4]); k = IND([2,3,4]);
            % compute the volume of the tetrahedral formed when each triangle 
            % is joined to the origin
            vol = vol + det([X(j), Y(j), Z(j)])/6 + det([X(k), Y(k), Z(k)])/6;
        end
    end
    vol = abs(vol);
end

[jixiong fei]

can this function compute the intersection part volume between two objects?

[Lateef Adewale Kareem]

if you can provide separate surface for each objects, and a single surface for the two combined, then, you can compute the volume of the intersection using V1 + V2 - V12

V1 is the volume of the surface1, V2 is the volume of surface2 and V12 is the volume of the combined surface12