MATLAB: Volume integration under surface fitting

integrationMATLABvolume

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

Best Answer

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.
Related Question