## Surface area from a z-matrix

on 15 Apr 2012

### Richard Brown (view profile)

I have x, y axis values and corresponding z values. This allows me to plot the shape using surf function. Example:-
How to calculate surface area of this irregular shape?
Please i'm a novice learning this stuff, i have come along quite well since last week or two. Thank you for all your help. The other answers provided before for irregular shapes have me confused, since i don't know how to triangulate.
Additionally i also want to calculate surface area above a certain height. So i'll cut all the values below a certain height and then possibly calculate surface area above a cut-off height.
Edit:
Thanks to inputs by the community and basically whole programme written by Richard. Here's the code.
function sa(Z, cutoff)
%Credit: Richard Brown, MATLAB central forum (http://www.mathworks.com/matlabcentral/answers/)
dx=0.092; % x-axis calibration
dy=0.095; % y-axis calibration
[m, n] = size(Z);
areas = 0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(1:m-1,1:n-1))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(1:m-1,1:n-1))).^2) + ...
0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(2:m,2:n))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(2:m,2:n))).^2);
zMean = 0.25 * (Z(1:m-1,1:n-1) + Z(1:m-1,2:n) + Z(2:m,1:n-1) + Z(2:m,2:n));
areas(zMean <= cutoff) = 0;
surfaceArea = sum(areas(:));
sprintf('Total surface area is %2.4f\n', surfaceArea)
return
end

Richard Brown

### Richard Brown (view profile)

on 15 Jul 2013
I've just had it pointed out to me that there is a small mistake in this code. Assuming that x corresponds to columns, and y to rows, then the first two lines (but not the second two lines) of the areas calculation has dx and dy backwards. I've fixed it in my original answer below.

### Richard Brown (view profile)

on 16 Apr 2012
Edited by Richard Brown

### Richard Brown (view profile)

on 15 Jul 2013

OK, well how about splitting each quadrilateral into two triangles, and just summing up the areas? I'm sorry there's no way I can make this look pleasant, but ...
[m, n] = size(Z);
areas = 0.5*sqrt((dx*dy)^2 + (dy*(Z(1:m-1,2:n) - Z(1:m-1,1:n-1))).^2 + ...
(dx*(Z(2:m,1:n-1 - Z(1:m-1,1:n-1)))).^2) + ...
0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(2:m,2:n))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(2:m,2:n))).^2)
surfaceArea = sum(areas(:))
edit: Jul 15, 2013 There was a mistake, and dx and dy were backwards in the first two lines of the code. The code has been corrected now.

bobby

### bobby (view profile)

on 16 Apr 2012
That is one awesome code Richard. I can't thank you enough. Here i was searching for triangulations and TINs since last thursday. It gives me result within 3% tolerance level, which is absolutely awesome. Just a small correction in the code, it is missing a bracket.
areas = 0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(1:m-1,1:n-1))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(1:m-1,1:n-1))).^2) + ...
0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(2:m,2:n))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(2:m,2:n))).^2);
surfaceArea = sum(areas(:))
bobby

### bobby (view profile)

on 16 Apr 2012
Now how to find the surface area above a cut-off height. The present code connects the empty valleys and assumes it to a flat surface. I want to exclude them. Any intelligent ideas?
Richard Brown

### Richard Brown (view profile)

on 16 Apr 2012
For each of the areas in the 'areas' array, you could work out the mean Z value
zMean = 0.25 * (Z(1:m-1,1:n-1) + Z(1:m-1,2:n) + Z(2:m,1:n-1) + Z(2:m,2:n))
and then
areas(zMean > cutoff) = 0;
surfaceArea = sum(areas(:));
It's a bit sloppy - you should probably do this for the individual triangles, but it will certainly give you a reasonable approximation

### Walter Roberson (view profile)

on 15 Apr 2012

I do not know if you will be able to calculate the surface area of the regular shape interspersed with the irregular tops. When I look at the image with the irregular tops, it looks to me as if the tops are fairly irregular, possibly even fractal. I don't know if a meaningful surface area could be calculated: the surface area of a fractal is infinite.
For the first figure, I can think of a crude way to find the area, but I think there are simpler ways. I would need to think further about good ways to find the area. But to cross-check: are your x and y regularly spaced, or irregularly ?

bobby

### bobby (view profile)

on 16 Apr 2012
hahaha Sean, i don't know if you are wanting me to work out the approach or commenting on Richard's script. I meant that i'm getting 20% error with Richard's script, which i think is significant. I'm in no way being facetious and putting down gentlemen giving me suggestions.
Sean de Wolski

### Sean de Wolski (view profile)

on 16 Apr 2012
I wasn't commenting on Richard's script at all, but rather the fractal nature of surface area (and length) of irregular shapes. I apologize if that came across the wrong way.
bobby

### bobby (view profile)

on 16 Apr 2012
Ah, the good old mandelbrot. Yes, I understand your point now. That is why i edited my post as i did understood that it will be nearly impossible. But the surface area of the surface that i'm after, given in the post, i'm sure a solution exists. Thanks to you i came across a very neat paper by Mandelbrot.

### Richard Brown (view profile)

on 16 Apr 2012

If your data is smooth enough (assuming that's what you're after), then there is a really quick way to work out approximately the surface area, using the normals that Matlab computes when plotting the surface. To make it simple I'll assume you have a uniform grid with spacings dx and dy.
dA = dx * dy;
h = surf( ... )
N = get(h, 'VertexNormals');
N_z = squeeze(N(:, :, 3));
% Normalise it
N_z = N_z ./ sqrt(sum(N.^2, 3));
Area = dA * sum(1./N_z(:));

Richard Brown

### Richard Brown (view profile)

on 16 Apr 2012
You *need* to smooth this first though, and the solution will only be approximate (probably to about 2sf). Otherwise you'd need to calculate the surface area of an interpolant, which would be a bit more effort intensive
bobby

### bobby (view profile)

on 16 Apr 2012
Thank you Richard for your interest and guidance. That is a very neat way of getting surface area. Yes, the grids are uniformly spaced. I am after a quite accurate solution here, since between different surfaces there can be a very small difference and i need to know that. This script is way off right now, way more than 2sigma. So i was looking at some TIN algorithms, i will create triangulated meshes and find individual triangle area. But even after spending weekend on this, i am nearly nowhere.
Richard Brown

### Richard Brown (view profile)

on 16 Apr 2012
Hence my comment about smoothing - if you apply that straight to the mesh you've got there, all bets are off! Anyway, see my next answer, it might be more in line with what you're after.

### bobby (view profile)

on 16 Apr 2012

function surfacearea(Z, cutoff)
dx=0.092;
dy=0.095;
[nrow, ncol] = size(Z);
for m=1:nrow
for n=1:ncol
zMean = mean(Z(1:m-1,1:n-1) + Z(1:m-1,2:n) + Z(2:m,1:n-1) + Z(2:m,2:n));
if(zMean >= cutoff)
areas = 0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(1:m-1,1:n-1))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(1:m-1,1:n-1))).^2) + ...
0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(2:m,2:n))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(2:m,2:n))).^2);
end
end
end
surfaceArea = sum(areas(:));
return
end
This is how the final code looks like with a cut-off limit included. Any inputs into speeding up this code? I have 1024*1280 matrix.

Richard Brown

### Richard Brown (view profile)

on 17 Apr 2012
I think you may have slightly misunderstood some of my earlier answers - you don't need to have any loops at all. All you should need to do is:
areas = 0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(1:m-1,1:n-1))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(1:m-1,1:n-1))).^2) + ...
0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(2:m,2:n))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(2:m,2:n))).^2);
zMean = 0.25 * (Z(1:m-1,1:n-1) + Z(1:m-1,2:n) + Z(2:m,1:n-1) + Z(2:m,2:n))
areas(zMean > cutoff) = 0;
surfaceArea = sum(areas(:));
And that should run pretty fast
bobby

### bobby (view profile)

on 17 Apr 2012
Oh yes. This does not require loops. Let me run it again. Thanks Richard.
bobby

### bobby (view profile)

on 17 Apr 2012
Yes, the code runs fast, only thing is there should be < sign. Thanks a lot Richard for helping me through this ordeal. I'll definitely shop again. :)