Volume integration under surface fitting

I have 3 set's of data: x, y and z. Fitting a surface with fit command works well, but i can't find a way to correctly evaluate the volume beneath the surface and above xy plane. I'm getting NaN and warnings: Warning: Non-finite result. The integration was unsuccessful. Singularity likely.

The code i'm using is the following:

load('data.mat')
x=data(:,1);
y=data(:,2);
z=data(:,3);
fitobject = fit([x,y],z, 'cubicinterp');
image
plot(fitobject)
xlabel('Proj_X');
ylabel('Proj_Y');
zlabel('Intensity');
disp('volume')
a=0;
b=40.34;
c=1.634;
d=80.05;
volume_under_fit = quad2d(fitobject,a,b,c,d)

The data sets are in the following file:

https://dl.dropbox.com/u/8101396/data.mat

Any help will be appreciated. Bruno

Answer 1

[Teja Muppirala]

Since your points only cover that parabolic looking region in X-Y, how is MATLAB supposed to know what the values of the function are outside that region?

It is clear that you will need to make some assumptions if you want to calculate that integral.

I can think of a few ideas:

1. Your data is more or less equal 0.043 with some oscillations, but the oscillations are relatively small in comparison to this constant term (moreover, they seem to be symmetric and kind of cancel out as well). So to a very good approximation, treat your data as constant, and the area under the curve is going to simply be the volume of the box. Like this:

mean(z)*(b-a)*(d-c) 

ans =

136.8444

2. Even though you don't have data points out there, it seems your data kind of extends out with the same behavior in the x-direction. So maybe you could make a copy of your data extended out this way, and then do the fitting.

load('data.mat')
x=data(:,1);
y=data(:,2);
z=data(:,3);

x = [x; x+max(x)]; %< --- Extend the points out periodically
y = [y; y]; %< --- Extend the points out periodically
z = [z; z]; %< --- Extend the points out periodically

fitobject = fit([x,y],z, 'cubicinterp');
plot(fitobject)
xlabel('Proj_X');
ylabel('Proj_Y');
zlabel('Intensity');
disp('volume')
a=0;
b=40.34;
c=1.634;
d=80.05;
volume_under_fit = quad2d(fitobject,a,b,c,d)

This gives you 138.1160, which is very close to the previous answer, and given the amount of noise in your data to begin with, I'm not sure you can say that one is better than another.

3. From looking at the data, you might try and model using a formula, such as a constant plus some sinusoidal terms.

z(x,y) = A + B*cos(W1*x + P1) + C*cos(W2*y + P2).

Or something like that, and then integrate that formula. But again, given the characteristics of your data to begin with, I don't think you'd get anything that you could say with certainty is more accurate than the above two answers.

Answer 2

[Shashank Prasanna]

Bruno, quad2d performs a double integral, which must only return the surface integral over the defined surface for the limits you specify.

If you want to compute the volume you should use triplequad:

http://www.mathworks.com/help/matlab/ref/triplequad.html

Also your third file for download is not z.mat but data.mat, hence I couldn't try this out myself.

[Bruno]

Thanks for trying to help.

I've figured what my problem is, quad2d (and maybe all others functions to perform volume integration) needs that the fitted surface covers all the domain of integration.But with my set of points that doesn't happen. Changing the limits of integration can't be done either (i think).

Using different integration limits, the volume can calculated (on those limits). Example:

volume_under_fit = quad2d(fitobject,0,20,30,50)

Matlab surely "knows" what is the volume beneath that fitted surface, the problem is how to "ask" him.

I've stored all datasets on data.mat. Please try again running the code, maybe you can find some alternative.

Bruno