Volume integration under surface fitting

7 次查看(过去 30 天)
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:
Any help will be appreciated. Bruno

采纳的回答

Teja Muppirala
Teja Muppirala 2013-2-18
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.

更多回答(1 个)

Shashank Prasanna
Shashank Prasanna 2013-2-16
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:
Also your third file for download is not z.mat but data.mat, hence I couldn't try this out myself.
  1 个评论
Bruno
Bruno 2013-2-16
编辑:Bruno 2013-2-16
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

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Surface and Mesh Plots 的更多信息

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by