MATLAB Answers

How to calculate the volume enclosed by a set of XYZ points in 3D?

55 views (last 30 days)
Sourav Sahoo
Sourav Sahoo on 23 Jan 2021
Commented: Sourav Sahoo on 24 Jan 2021
Hi,
I am trying to find the volume of a region (defined by X,Y and Z coordinates) enclosed below a (Z='constant') plane. The data has peaks (positive Z) and a valley (negative Z), and the mean surface is assigned z=0. I have tried with the following piece of code, but I doubt it gives me the total volume bound by the surface against z=0, including the peak volumes as well.
[X,Y,Z] = xyzread("data.xyz");
plot3(X,Y,Z)
[fitobject, gof, output] = fit([X,Y],Z, 'biharmonicinterp');
plot(fitobject)
a = min(X);
b = max(X);
c = min(Y);
d = max(Y);
volume_under_fit = quad2d(fitobject,a,b,c,d)

  0 Comments

Sign in to comment.

Accepted Answer

Bruno Luong
Bruno Luong on 23 Jan 2021
Edited: Bruno Luong on 23 Jan 2021
V is the volume between the plane x-y (z==0) and the surface z(x,y) from your data.
If you want the volume of the data after substract the base plane surface, you need to estimated the equation by regression.
load('data.xyz')
x=data(:,1);
y=data(:,2);
z=data(:,3);
T=delaunay(x,y);
trisurf(T,x,y,z);
xy = [x,y];
a = xy(T(:,2),:)-xy(T(:,1),:);
b = xy(T(:,3),:)-xy(T(:,1),:);
V = ((a(:,1).*b(:,2)-a(:,2).*b(:,1))' * sum(z(T),2))/6

  5 Comments

Show 2 older comments
Sourav Sahoo
Sourav Sahoo on 23 Jan 2021
Hi,
I have not changed the code. Only filtered the data to exclude values >=0.
"Up to you to separate negative and positive volumes"
  • I'm stuck with this the entire day. How should I go about separating the positive (peaks) and negative volumes (valley) from the total?
Thanks in advance, Bruno.
Bruno Luong
Bruno Luong on 23 Jan 2021
load('data.xyz')
x=data(:,1);
y=data(:,2);
z=data(:,3);
T=delaunay(x,y);
xy=[x,y];
a = xy(T(:,2),:)- xy(T(:,1),:);
b = xy(T(:,3),:)- xy(T(:,1),:);
A = (a(:,1).*b(:,2)-a(:,2).*b(:,1))/2;
h = mean(z(T),2);
pos = h > 0;
sh = trisurf(T,x,y,z,'CData',double(pos.'),'EdgeColor','none');
colormap(jet)
Vt = A .* h;
Vpos = sum(Vt(pos)) % sign (+) volume of of the red part
Vneg = sum(Vt(~pos)) % sign (-) volume of of the blue part
V = Vpos+Vneg

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!