calculate volume from iso-surface coordinates (x,y,z).
조회 수: 24 (최근 30일)
이전 댓글 표시
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
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
답변 (3개)
Walter Roberson
2023년 9월 21일
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
2023년 9월 21일
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
Walter Roberson
2023년 9월 22일
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
2023년 9월 21일
편집: Bruno Luong
2023년 9월 21일
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
Bruno Luong
2023년 9월 22일
편집: Bruno Luong
2023년 9월 22일
@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)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!