How to use isosurface to measure volume of the cone??

3 views (last 30 days)
M.S. Khan on 30 Aug 2020
Edited: M.S. Khan on 31 Aug 2020
Using code below , i have measure volume of a cone but i want to measure it using isosuface() i.e. using voxels. Can anyone help me please. I will be very thankful.
Thanks in advance for all help.
Regards
R = 10;
h = 20;
N = 10000;
% Generate points in cone
dx = (pi*R^2*h/3/N)^(1/3);
rvec = linspace(-R,R,ceil(2*R/dx));
hvec = linspace(0,h,ceil(h/dx));
[X,Y,Z] = ndgrid(rvec,rvec,hvec);
is_in_cone = (X.^2+Y.^2) <= (R/h*(h-Z)).^2;
x = X(is_in_cone);
y = Y(is_in_cone);
z = Z(is_in_cone);
% Estimate the volume of the 3D object assuming the object is convex
tri = delaunay(x,y,z);
trisurf(tri,x,y,z)
% indices of triangle
i1 = tri(:,1);
i2 = tri(:,2);
i3 = tri(:,3);
i4 = tri(:,4);
%Volume by summing tetrahedron
v1 = [x(i1)-x(i2) y(i1)-y(i2) z(i1)-z(i2)];
v2 = [x(i1)-x(i3) y(i1)-y(i3) z(i1)-z(i3)];
v3 = [x(i1)-x(i4) y(i1)-y(i4) z(i1)-z(i4)];
A = 1/2*cross(v1,v2,2); % surface of a triangle
V = 1/3*dot(A,v3,2); % volume of a tetrahedron
format long
V = sum(abs(V))

Bruno Luong on 30 Aug 2020
Edited: Bruno Luong on 31 Aug 2020
Theoretical volume = 2094.4, Isosurface estimation = 2092.4.
Pretty good to me.
R = 10;
h = 20;
N = 1e5;
% Theoretical volume
Vol = pi*R^2*h/3
% Generate points in cone
dx = (pi*R^2*h/3/N)^(1/3);
margin = 1e-3;
rmin = -(1+margin)*R;
rmax = +(1+margin)*R;
rvec = linspace(rmin,rmax,ceil((rmax-rmin)/dx));
zmin = -0.5*margin*h;
zmax = (1+0.5*margin)*h;
hvec = linspace(zmin,zmax,ceil((zmax-zmin)/dx));
[X,Y,Z] = meshgrid(rvec,rvec,hvec);
Vcone = sqrt(X.^2+Y.^2) - abs(R/h*(h-Z));
Vcone = max(Vcone,-Z);
Vcone = max(Vcone,Z-h);
[F,V] = isosurface(X,Y,Z,Vcone,0);
% Volume estimated from isosurface
VF = permute(reshape(V(F,:),[size(F) 3]),[3 1 2]);
Vol = 1/6*sum(dot(cross(VF(:,:,1),VF(:,:,2),1),VF(:,:,3),1))
% close all
% isosurface(X,Y,Z,Vcone,0)
% axis equal
% axis([-20 20 -20 20 -10 30])

1 Comment

M.S. Khan on 31 Aug 2020
Dear Bruno, Thank you so much. you are shining star. If possible, please help me how to measure its volume layer by layer. i have added this question of layers in the new thread also. Warmest regards.

Walter Roberson on 30 Aug 2020
I absolutely positively would not use isosurface() to measure volumes. It makes the task into a Computer Vision task, requiring that you take screen captures of the surface from several different angles (6 different angles at absolute minimum), and use Stereo Vision techniques to try to reconstruct depths. It is too many unnecessary complications for too low an accuracy.
Do not do this. Find a different approach.

Show 1 older comment
Bruno Luong on 30 Aug 2020
"I absolutely positively would not use isosurface() to measure volumes"
Never said never. The iso surface return an oriented boundary mesh that can be worked out. Just need to be careful about the implicit function that returns all the relevant faces and not extra faces.
"Thanks for reply. could you please suggest any different appraoch. i will be very thankful. "
I did show you 3 other ways of computing the volume: theoretical, convhull, tetrahedron (the code you copy above)... how many methods do you need?
M.S. Khan on 31 Aug 2020
Dear Bruno, thanks for your kind reply. i want to learn and practicing different methods and techniques using matlab. i am very very thankful for your cooperation and guidance.Your contribution of knowlege to the globe of humanity is marvelous. Thanks for your time and efforts for new learners of matlab. i am very much impressed. Good wishes and prayers for you.
M.S. Khan on 31 Aug 2020
Dear Bruno, i need one more favor. Using these methods, we measured volume of whole object. Now, if we want to divide this object into layers and then we measure volume layer by layer of the object. so, the sum of volume of all layers will give us volume of the object. Is it possible? Sorry for asking again and again. i know i am taking your valuable time.