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

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)

Answer 1

[Bruno Luong]

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

[Sourav Sahoo]

Hi Bruno,

Thanks for the reply.

Your code gives me the total volume as -7.9424e-19 cubic units. I filtered my data to exclude the "peaks" and then feeded it back to your code so as to get the volume of the "valley" alone, and I'm getting the valley volume as -8.7681e-19 cubic units. Any idea why the magnitude turns out to be greater than the previous case, when we're reducing a portion of volume. Attached are the images for your reference.

[Bruno Luong]

"V is the volume between the plane x-y (z==0) and the surface z(x,y) from your data."

V is the (signed) integral of the whole thing.

If you filter out the positive surface, you get smaller volume.

The formula computed the integral, suming up the positive and negative. Up to you to separate negative and positive volumes.

If you have a question of YOUR result, show us what you have done with the code.

[Sourav Sahoo]

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]

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      

[Sourav Sahoo]

Hi,

Thanks a lot for your help! It works!!

Answer 2

[Walter Roberson]

https://www.mathworks.com/support/matlabcentral/answers/489219-alternatives-to-delaunay-triangulation