Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Volume integration under surface fitting

Asked by Bruno on 15 Feb 2013

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

0 Comments

Bruno

Products

2 Answers

Answer by Teja Muppirala on 18 Feb 2013
Accepted 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.

0 Comments

Teja Muppirala
Answer by Shashank Prasanna on 16 Feb 2013

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.

1 Comment

Bruno on 16 Feb 2013

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

Shashank Prasanna

Contact us