How to find volume under 2D data?

I have a set of 2D (x and y) data, that has an irregular shape and it's not axisymmetric. I want to find the volume under this surface (due to rotation around y-axis or if there is a better way). I am not sure how to do it and I thought about using solid of revolution but I don't think it is suitable for this case. Any idea how I can find the volume?

Answer 1

[Walter Roberson]

[idx, volume] = boundary(x(:), y(:), 1);

plot(x(idx), y(idx), 'y-*')

The plot() will illustrate the boundary that is found, and the variable volume will be the area contained within the boundary.

The area might not be the most accurate possible, but it is not clear what area would be accurate considering the gap.

A different calculation method would be to separate it into segments without gaps and use trapz() on each of those. You would then have to figure out area was implied by the gaps.

[Stefan Carson]

There is no gap in the data, the above image had an error in it. Please see the attached image for clarification.

The function trapz() will give the area under the 2D. How about the volume? is there a function or a way to calculate the volume? I thought about rotating the area about the y-axis to get the volume, but I am not certain it will be correct.

Also, the above function for boundary is giving the following error

Insufficient number of outputs from right hand side of equal sign to satisfy assignment.
Error in VTtest (line 3)
[idx, volume] = boundary(x(:), y(:), 1); 

[Walter Roberson]

Which MATLAB release are you using?

Is the data already sampled at equal intervals on both sides of 0? And if so is 0 itself one of the x locations? That is, for each y(x) is there a corresponding y(-x) recorded? If there is not, then interpolation can be done; it is just easier if it is already symmetric sampling.

[Stefan Carson]

I am using MATLAB 2020a. No, the data is not sampled at equal intervals on both sides of 0 as this is acquired from experimental data and It has a random shape. For example, sometimes the left side (-x) is wavy and the right side (+x) is similar to a semisphere. The following is the raw data without trying to making it symmetric.

[Walter Roberson]

Is the task to rotate around the peak with peak detection needed? Or is the location of the peak already known?

Note: it is my bedtime but perhaps someone else will have a suggestion while I sleep.

[Stefan Carson]

The location of the peak is already known, and I am trying to figure out how to calculate the volume either by slicing and assuming it's a circular disk or rotation about the y-axis. Thank you for your help and have a good night sleep.

[Walter Roberson]

To rotate around the axes x = location_of_peak:

Let x at the location of the peak be and y at peak be

For any one y in the range to , find the x before the peak at which that y occurs, and the x after the peak at which that y occurs, . Then rotating around at that y will occupy a circle of radius contributing an area of and the volume would be like

trapz(Y, pi*max(xplus(Y,xp)-xp, xp-xminus(Y,xp)).^2)

for a vector of Y from 0 to yp, and functions xplus and xminus .

I show xplus and xminus as having xp as a parameter because in practice they would be search or interpolation routines that are going to need to know the position of the peak so that they can extract the correct value. In practice it might be the index of xp rather than the value of xp that is more useful.

It looks to me as if your y values are monotonic on each side of the peak, so I guess

xminus = @(Y, xpindex) interp1(y(1:xpindex), x(1:xpindex), Y)

and if so then it should work for vector Y rather than having to loop. xplus would be similar,

xplus = @(Y, xpindex) interp1(y(xpindex:end), x(xpindex:end), Y)

Answer 2

[Bruno Luong]

% Generate some fake data, replace r with your X and y with your Y vector
% you need to peak either side r >= 0 or r <= 0 NOT both
r=linspace(0,150);
s=r/r(end);
y=polyval([-50 -100 160],s)

plot(r,y)

V = 2*pi*trapz(r,y.*abs(r)) % <- Volume formula

[Stefan Carson]

Thank you. If I used trapz() I get a value that is slightly similar to semisphere of a specific radius. However, if I use

V= pi*integral(@(X) yfun(X).*X,xmin,xmax)

I get a different value. But, I am not sure, which would be the correct one.

[Walter Roberson]

What is yfun here? It tends to imply that you have a formula for your values, but you do not. Is yfun

yfun = @(X) interp1(x,y,X)

with some interpolation specified such as spline ?

It looks to me as if this might be rotating around the wrong axes ?

[Bruno Luong]

Let me try to guess and fill what missing in your question

rho=2; % example of radius of the ball
yfun = @(r) sqrt(rho^2-r.^2);
r=linspace(0,rho,513);
y = yfun(r);

plot(r,y);
axis('equal')

% Half Volume of the 3-ball := { P=(x,y,z) in R^3: |P|<=rho; y >=0 } 
Vdiscrete = 2*pi*trapz(r,y.*abs(r)) 
Vnumerical = 2*pi*integral(@(r) yfun(r).*r, 0, rho)
Vexact = (4*pi/3*rho^3) / 2 % https://en.wikipedia.org/wiki/Volume_of_an_n-ball

[Stefan Carson]

Yes, Walter. I was trying to interpolate using a spline as the points are not always fitted to match a certain shape. I do not know if this approach is suitable because if I tried the code from Bruno, it will work for some cases where it can still be similar to a semisphere, but if the shape is random with peaks and valleys, then the result is not correct.

[Stefan Carson]

Sorry for not being clear and filling in the gaps in my question. I was trying to explain that if I assume the shape is a semisphere, the value of the integral or trapz() will be overestimated and it does not correspond to the correct value. This is mainly because the left and the right sides about the x-axis are not symmetric and it becomes complex trying to rotate or assume it's a disk to perform the calculation by slicing. Thank you.

[Bruno Luong]

Why reinvent the wheel?

The integration schemes are studied in long and large in the history, and people derive scheme such as trapezoidal rule for generic discrete data. If you know your function is smooth, then Simson has derive more accurate rules (for such function) without the need of doing explicit interpolation.

So just change the trapz rule by simpson

https://en.wikipedia.org/wiki/Simpson%27s_rule

But look at your data on left/right side; I think your data has more error than any scheme effect. Might be it's wrong to assume the solid is a revolution about y-axis/ You might be better doing some fitting and not interpolation, since those oscillation looks to me like measurement artefact than real.

But only you know the quality of your data. No numerical scheme can recover bad measurements.