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

Question

Raju 2023-9-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2024232-calculate-volume-from-iso-surface-coordinates-x-y-z

评论： Raju 2023-9-24

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Fifteen12 2023-9-21

编辑：Fifteen12 2023-9-21

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 2023-9-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2024232-calculate-volume-from-iso-surface-coordinates-x-y-z#answer_1315497

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Walter Roberson 2023-9-21

编辑：Walter Roberson 2023-9-21

在 MATLAB Online 中打开

[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.

Walter Roberson 2023-9-22

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 2023-9-22

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

The text is clipped.

Walter Roberson 2023-9-22

@Bruno Luong

Hah! That would do it!

请先登录，再进行评论。

Answer 2

William Rose 2023-9-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2024232-calculate-volume-from-iso-surface-coordinates-x-y-z#answer_1315447

@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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

William Rose 2023-9-21

在 MATLAB Online 中打开

@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.

Walter Roberson 2023-9-22

在 MATLAB Online 中打开

Q.txt

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 2023-9-24

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

请先登录，再进行评论。

Answer 3

Bruno Luong 2023-9-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2024232-calculate-volume-from-iso-surface-coordinates-x-y-z#answer_1315507

编辑：Bruno Luong 2023-9-21

在 MATLAB Online 中打开

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).

6 个评论
显示 4更早的评论隐藏 4更早的评论

Bruno Luong 2023-9-22

编辑：Bruno Luong 2023-9-22

在 MATLAB Online 中打开

@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 2023-9-24

@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).

请先登录，再进行评论。