calculate volume from iso-surface coordinates (x,y,z).

Hello,

I have coordinates (x,y,z) of an isosurface. How can I calculate volume of that isosurface? I have attached an image of iso-surface and coordinates file here.

[Fifteen12]

Do isosurfaces necessarily have a volume? Are these completely closed surfaces? If you just need to calculate the surface area you could check out this approach (I haven't looked at it myself): https://www.mathworks.com/matlabcentral/fileexchange/25415-isosurface-area-calculation?s_tid=answers_rc2-2_p5_MLT

Answer 1

[Walter Roberson]

If you use boundary then it has an optional second output which is the volume.

However, I would not expect boundary() to be able to deal with disconnected components, so you would need to separate out the different components based on the vertices returned by isosurface().

You could potentially turn the vertices into a connection table and create a graph object from it, and use something like conncomp to determine which sections are disconnected from the others.

[Raju]

Are there any example code of that?

[Walter Roberson]

[x,y,z] = meshgrid([-3:0.25:3]); 
V = x.*exp(-x.^2 -y.^2 -z.^2);
S = isosurface(x, y, z, V, 1e-9);

s = S.faces;
v = S.vertices;
t = circshift(S.faces, [0 -1]);
G = graph(s(:), t(:));

bins = conncomp(G);
nbins = max(bins);

figure();
patch(S, 'edgecolor', 'none');
view(3)
title('isosurfaces');

figure();
volume_per_component = zeros(nbins,1);
component_boundaries = cell(nbins,1);
for K = 1 : nbins
    mask = bins == K;
    vertex_subset = v(mask,:);
    [component_boundaries{K}, volume_per_component(K)] = boundary(vertex_subset);
    trisurf(component_boundaries{K}, vertex_subset(:,1), vertex_subset(:,2), vertex_subset(:,3), 'FaceAlpha', 0.2);
    text(vertex_subset(1,1), vertex_subset(1,2), vertex_subset(1,3), "volume = " + volume_per_component(K));
    hold on
end
hold off
title('sub-volumes');

I am not sure at the moment why one of the volumes comes out as 0.... Ah, it looks like it might be planar, and the volume of a planar area would be 0.

[Raju]

Applying the code but not getting the expected isosurface. I have attached an isosurface and the correspording coordinates here so that you can figure out if the code works.

[Walter Roberson]

I do not see how those coordinates relate to isosurface? In order to have an isosurface you need connection information. Either that or you need to have a scattered interpolation type of situation where you have values at known positions and you interpolate a 3D volume over a cuboid grid and then run an isosurface over that volume. But your provided data does not have any face connection information and it also does not have any node value information -- just 3D vertices and vertex numbers. And your image at https://www.mathworks.com/matlabcentral/answers/2024232-calculate-volume-from-iso-surface-coordinates-x-y-z?s_tid=srchtitle#comment_2895787 shows that you are expecting a potentially-disconnected isosurface.

[Walter Roberson]

Imagine you have a solid 3 x 3 cube, and you have the list of corresponding vertices. 16 points per plane, 4 planes = 64 points

A---B---C---D

| | | |

E---F---G---H

| | | |

I---J---K---L

| | | |

M---N---O---P

%one surface of a 3 x 3 cube

Now imagine that you remove one of the interior cubes. For example remove the E/F/J/I cube.

How would that change the vertex table?

If you were to remove the A/B/F/E cube, then, Sure, the table would change because you would no longer have vertex A anywhere -- but if you remove the E/F/J/I cube then vertex E is still needed for A/B/F/E, vertex F is still needed for A/B/F/E and B/C/G/F and F/G/K/J ... and so on. So if you remove an interior cube, the list of vertices would not change. But clearly the volume should .

But if all you have is the list of vertices, then since the vertex list would be the same, you cannot tell that E/F/J/I cube was removed.

This demonstrates that just having a list of vertices is not enough, even in a very simple case. You must have the face connection information.

[Bruno Luong]

@Walter Roberson "I'am not sure at the moment why one of the volumes comes out as 0.... "

The text is clipped.

[Walter Roberson]

@Bruno Luong

Hah! That would do it!

Answer 2

[William Rose]

@Raju,

Find the delaunay triangulation of the 3D points with

DT=delaunay(x,y,z);

This gives a set of tetrahedrons which fill the volume. Then compute and add up the volumes of the tetrahedrons.

[William Rose]

@Raju,

For example:

% Find 100 random points on the units sphere
phi=asind(-1+2*rand(1,100)); % angle above the x-y plane (deg)
theta=360*rand(1,100);       % angle measured in the x-y plane (deg)
x=cosd(theta).*cosd(phi);
y=sind(theta).*cosd(phi);
z=sind(phi);
% Plot the points. You can spin it around to confirm that they are on the sphere.
figure;
plot3(x,y,z,'r*')
axis equal; grid on

% Find the Delaunay triangulation
DT=delaunay(x,y,z);

DT is a list of indices of the corners of tetrahedrons that fill the volume. There are various formulas for the volume of a tetrahedron when the corner coordinates are known. One formula is

where you must take the absolute value of the quantity above.

N=length(DT)

N = 277

xyz=[x;y;z;ones(1,100)]';
vol=zeros(1,N);  % allocate vector for the volumes
for i=1:N
    a=[xyz(DT(i,1),:);xyz(DT(i,2),:);xyz(DT(i,3),:);xyz(DT(i,4),:)];
    vol(i)=abs(det(a)/6);
end
volTot=sum(vol);
fprintf('Total volume=%.2f.\n',volTot)

Total volume=3.70.

We expect the volume to be a bit less than = volume of the unit sphere. If we use 1000 ponts incstead of 100, the volume should be closer to 4.19. Try it. Good luck.

[Raju]

@William Rose Thanks, for your comments. I have attached an image of an iso-surface and coordinates here to see if the delaunay function is applicable here. I want to calculate volume under this 3D iso-surface.

[Walter Roberson]

D = readmatrix('Q_100.txt');
XYZ = D(:,2:4);
[k, volume] = boundary(XYZ);
trisurf(k, XYZ(:,1), XYZ(:,2), XYZ(:,3), 'Facecolor','cyan','FaceAlpha',0.8);
axis equal;
hold on
scatter3(XYZ(:,1), XYZ(:,2), XYZ(:,3), 'r.');
hold off

volume

volume = 3.5096e-07

This tells us that if all we have available is the vertex coordinates, then we have no basis for the close curves that you believe should be there, and consequently the volume is probably not going to have much relationship to what you want.

Isosurface information must have the vertex connection information to be able to produce anything reasonable.

[Raju]

@Walter Roberson Thank you so much for your effort. I have attached a new text file here which have a vertex information associate with coordinates.

[Bruno Luong]

@Raju we need the connectivity (the faces output of isosurface command) as well, not the coordinates alone.

Do you understand what is connectivity?

[Walter Roberson]

When we have values for each point but no connectivity information for the vertices, then the only possibility is to treat the points as being scattered samplings of a continuous function, and to interpolate those scattered positions over a grid and construct isosurfaces of the result.

... It doesn't look very good.

data = readmatrix('Q.txt');
x = data(:,2);
y = data(:,3);
z = data(:,4);
q = data(:,5);
F = scatteredInterpolant(x, y, z, q);

N = 50;
[minx, maxx] = bounds(x);
[miny, maxy] = bounds(y);
[minz, maxz] = bounds(z);
[qX, qY, qZ] = meshgrid(linspace(minx, maxx, N), linspace(miny, maxy, N), linspace(minz, maxz, N));

qQ = F(qX, qY, qZ);

[minq, maxq] = bounds(qQ(:));
isolevels = linspace(minq, maxq, 6);
isolevels([1 end]) = [];

for V = isolevels
    isosurface(qX, qY, qZ, qQ, V);
end
view(3)
legend("q = " + isolevels);

figure()
h = scatter3(x, y, z, [], q);
%h.AlphaData = 0.3;
h.MarkerEdgeAlpha = 0.1;
h.MarkerFaceAlpha = 0.1;

[Raju]

@Walter Roberson I appreciate your effort. I will have a look if I can get some connection data.

Answer 3

[Bruno Luong]

Do you have connectivity face of these points coordinates?

If you use the command isosurface https://www.mathworks.com/help/matlab/ref/isosurface.html you should have. Please share the outputs faces and verts or structure s (save in matfile and attach here).

Or try this formula:

[x,y,z] = meshgrid([-1.1:0.05:1.1]); 
V = x.^2 + y.^2 + z.^2;
s = isosurface(x,y,z,V,1) % replace this command using your data

s = struct with fields:
    vertices: [7470×3 double]
       faces: [14936×3 double]

VF = permute(reshape(s.vertices(s.faces,:),[size(s.faces) 3]),[3 1 2]);
Vol = 1/6*sum(dot(cross(VF(:,:,1),VF(:,:,2),1),VF(:,:,3),1)) % close to 4/3*pi volume of the sphere of raduius 1

Vol = 4.1809

 4/3*pi

ans = 4.1888

This formula works for non-convex volume enclosed by the surface given by triangular connectivity (correctly oriented).

[Bruno Luong]

Example non convex volume

[X, Y, Z] = meshgrid(linspace(-2*pi, 2*pi, 200));
iR2 = 1./(X.^2+Y.^2+Z.^2);
C = iR2 .* (sin(X).*cos(Y) + sin(Y).*cos(Z) + sin(Z).*cos(X));

figure
isosurface(X, Y, Z, C, 0.05);
axis equal

s = isosurface(X, Y, Z, C, 0.05); % replace this command using your data

VF = permute(reshape(s.vertices(s.faces,:),[size(s.faces) 3]),[3 1 2]);
Vol = abs(1/6*sum(dot(cross(VF(:,:,1),VF(:,:,2),1),VF(:,:,3),1))) 

Vol = 137.2342

[Raju]

@Bruno Luong Thank you for your effort. Can you bit explain iR2 and C used in the code. Why are we using these particular functions?

[Bruno Luong]

You don't need to understand C and iR2 this is MY isosurface as example, as you refuse to share correctly your data (missing connectivity s.faces).

You have to use YOUR isosurface outputs. What is important is these commands

level = 0.05;
s = isosurface(X, Y, Z, C, level); % replace this command using your data
VF = permute(reshape(s.vertices(s.faces,:),[size(s.faces) 3]),[3 1 2]);
Vol = abs(1/6*sum(dot(cross(VF(:,:,1),VF(:,:,2),1),VF(:,:,3),1))) 

[Raju]

Thank you for your reply. I have attached a new text file here which has connectivity data you found missing. Please check.

[Bruno Luong]

@Raju Sigh, I still don't see any connectivity data. Can't help you more.

[X, Y, Z] = meshgrid(linspace(-2*pi, 2*pi, 200));
iR2 = 1./(X.^2+Y.^2+Z.^2);
C = iR2 .* (sin(X).*cos(Y) + sin(Y).*cos(Z) + sin(Z).*cos(X));
s = isosurface(X, Y, Z, C, 0.05); % replace this command using your data
% the connectivity mooke like this
s.faces(1:10,:),

ans = 10×3

1 2 3 1 4 2 2 4 5 4 6 5 5 6 7 6 8 7 9 10 11 9 12 10 10 12 1 10 1 3

The connectivity tells the mesh triangles of the surface connect which vertexes. As above the last line tell the 10th triangle is composed of of three vertices (#10, #1, #3)

[Raju]

@Bruno Luong Thanks for your effort. But these are the data I have right now. I will have a look if get some additional data (connection data you mean).