The currently accepted answer will only work for gridded data. If your data is scattered you can use the following approach instead:
scatteredInterpolant
+ integral2
:
f = scatteredInterpolant(x(:), y(:), z(:), 'linear');
int = integral2(@(x,y) f(x,y), xmin, xmax, ymin, ymax);
This defines the linear interpolant f
of the data z(i) = f(x(i),y(i))
and uses it as an argument to integral2
. Note that ymin
and ymax
, instead of doubles, can be function handles depending on x
. So usually you will be integrating rectangles, but this could be used for integration regions a bit more complicated.
If your integration area is rather complicated or has holes, you should consider triangulating your data.
DIY using triangulation
:
Let's say your integration area is given by the triangulation trep
, which for example could be obtained by trep = delaunayTriangulation(x(:), y(:))
. If you have your values z
corresponding to z(i) = f(trep.Points(i,1), trep.Points(i,2))
, you can use the following integration routine. It computes the exact integral of the linear interpolant. This is done by evaluating the areas of all the triangles and then using these areas as weights for the midpoint(mean)-value on each triangle.
function int = integrateTriangulation(trep, z)
P = trep.Points; T = trep.ConnectivityList;
d21 = P(T(:,2),:)-P(T(:,1),:);
d31 = P(T(:,3),:)-P(T(:,1),:);
areas = abs(1/2*(d21(:,1).*d31(:,2)-d21(:,2).*d31(:,1)));
int = areas'*mean(z(T),2);