calculate volume from iso-surface coordinates (x,y,z).
24 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaire
Fifteen12
le 21 Sep 2023
Modifié(e) : Fifteen12
le 21 Sep 2023
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
Réponses (3)
Walter Roberson
le 21 Sep 2023
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().
William Rose
le 21 Sep 2023
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 commentaires
Walter Roberson
le 22 Sep 2023
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;
Bruno Luong
le 21 Sep 2023
Modifié(e) : Bruno Luong
le 21 Sep 2023
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
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
4/3*pi
This formula works for non-convex volume enclosed by the surface given by triangular connectivity (correctly oriented).
6 commentaires
Bruno Luong
le 22 Sep 2023
Modifié(e) : Bruno Luong
le 22 Sep 2023
@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,:),
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)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!