Surface area from a z-matrix

I have x, y axis values and corresponding z values. This allows me to plot the shape using surf function. Example:-

http://i.imgur.com/89KoR.png

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]

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.

Answer 1

[Richard Brown]

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]

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]

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]

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

Answer 2

[Walter Roberson]

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]

x, y are regularly spaced.

I know the regular tops with interspersed ones will be very difficult, but i think i'll somehow alienate the z-matrix for each individual pillar tops and then use the algorithm for second ones on them, and surface area of a cylinder will be easy. I know a very crude way, but since i'm dealing with fractal surface i think a trend will be good enough too.

[Walter Roberson]

Sorry, I do not recognize "alienate" in that context? "isolate" perhaps?

If you are dealing with a surface that is literally fractal, then the surface area is infinite.

[bobby]

Yes, i meant isolated and not alienated. You are right, for self-affine surfaces i cannot find surface area and hence i'm droppping the idea of finding Surface area for irregular top structures, and also edited the original post to show the same.

Any direction you can give me to solve the surface area for first surface.

[Sean de Wolski]

How long is the coast of England?

[bobby]

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]

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.

http://www.google.com/#hl=en&gs_nf=1&cp=19&gs_id=1k&xhr=t&q=how+long+is+the+coast+of+britain&fp=46b06bf633a62cc4

[bobby]

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.

Answer 3

[Richard Brown]

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]

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]

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]

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.

Answer 4

[bobby]

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]

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]

Oh yes. This does not require loops. Let me run it again. Thanks Richard.

[bobby]

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. :)