Area Calculation of Surf plot in MATLAB

Question

I have irregular 3D cartesian coordinates that make up the one eighth of the surface of a sphere. Thanks to Benoit_11 Answer to a previously posed question the surface can now be plotted outside of the MATLAB cftool in a normal command line script

Since this point I have been attempting to calculate the area of the surface using the following code patched together from other answers within the area. The code basically calculates the area from the vertices that produce the surface and then sums them all to produce an area.

surface = [ansx1,ansy1,ansz1];
[m,n] = size(zdata1);
area = 0;
for i = 1:m-1
      for j = 1:n-1
          v0_1 = [xdata1(i,j)     ydata1(i,j)     zdata1(i,j)    ];
          v1_1 = [xdata1(i,j+1)   ydata1(i,j+1)   zdata1(i,j+1)  ];
          v2_1 = [xdata1(i+1,j)   ydata1(i+1,j)   zdata1(i+1,j)  ];
          v3_1 = [xdata1(i+1,j+1) ydata1(i+1,j+1) zdata1(i+1,j+1)];
          a_1= v1_1 - v0_1;
          b_1 = v2_1 - v0_1;
          c_1 = v3_1 - v0_1;
          A_1 = 1/2*(norm(cross(a_1, c_1)) + norm(cross(b_1, c_1)));
          area = area + A_1;
      end
end
fprintf('\nTotal area is: %f\n\n', area);`

However the issue I am having is that the calculated surface over estimates the possible surface. This is due to the removal of NaN from the original matrix and their replacement with 0 this results in the figure 1. Figure 2 provides the only area I would like to calculate

Does anybody have a way of ignoring the zeros within the code provided to calculate the surface area of the data that generates Figure 1?

Thanks in advance

Answer 1

I think you only have to check if one of the four points of an field equal zero.

What about this:

% example surface
[X,Y,Z] = peaks(30);

% manipulate it
[lza, lzb] = size(Z);
for nza = 1:lza
   for nzb = 1:lzb
      if Z(nza,nzb) < 0
         Z(nza,nzb) = Z(nza,nzb)-1;
      else
         Z(nza,nzb) = 0;
      end
   end
end

surfc(X,Y,Z)

% start calculating the surface area
A = 0;
lX = length(X);
lY = length(Y);

for nx = 1:lX-1
   for ny = 1:lY-1

      eX = [X(ny,nx)   X(ny,nx+1)
         X(ny+1,nx) X(ny+1,nx+1)];
      eY = [Y(ny,nx)   Y(ny,nx+1)
         Y(ny+1,nx) Y(ny+1,nx+1)];
      eZ = [Z(ny,nx)   Z(ny,nx+1)
         Z(ny+1,nx) Z(ny+1,nx+1)];

      % check the field
      if eZ(1,1)==0 || eZ(1,2)==0 || eZ(2,1)==0 || eZ(2,2)==0
         continue
      end

      % take two triangles, calculate the cross product to get the surface area
      % and sum them.
      v1 = [eX(1,1) eY(1,1) eZ(1,1)];
      v2 = [eX(1,2) eY(1,2) eZ(1,2)];
      v3 = [eX(2,1) eY(2,1) eZ(2,1)];
      v4 = [eX(2,2) eY(2,2) eZ(2,2)];
      A  = A + norm(cross(v2-v1,v3-v1))/2;
      A  = A + norm(cross(v2-v4,v3-v4))/2;

   end
end