calculate volume from iso-surface coordinates (x,y,z).
24 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentar
Fifteen12
am 21 Sep. 2023
Bearbeitet: Fifteen12
am 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
Antworten (3)
Walter Roberson
am 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
am 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 Kommentare
Walter Roberson
am 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
am 21 Sep. 2023
Bearbeitet: Bruno Luong
am 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 Kommentare
Bruno Luong
am 22 Sep. 2023
Bearbeitet: Bruno Luong
am 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)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!