How do I calculate the area of a 3 dimensional projection?

Question

For example, using the following code, I have a coordinate matrix with 3 cubical objects defined by 8 corners each, for a total of 24 coordinates. I apply a rotation to my coordinates, then delete the y coordinate to obtain a projection in the x-z plane. How do I calculate the area of these cubes in the x-z plane, ignoring gaps and accounting for overlap? I have tried using polyarea, but this doesn't seem to work.

clear all
clc
A=[-100 -40 50
-100    -40 0
-120    -40 50
-120    -40 0
-100    5   0
-100    5   50
-120    5   50
-120    5   0
-100    0   52
-100    0   52
20  0   5
20  0   5
-100    50  5
-100    50  5
20  50  52
20  50  52
-30 70  53
-30 70  0
5   70  0
5   70  53
-30 120 53
-30 120 0
5   120 53
5   120 0]; %3 Buildings Coordinate Matrix
theta=60; %Angle
rota = [cosd(theta) -sind(theta) 0; sind(theta) cosd(theta) 0; 0 0 1]; %Rotation matrix
R=A*rota; %rotates the matrix
R(:,2)=[];%deletes the y column

Answer 1

polyarea is the right way to go, but you need to call it on the convex hull of each projection. If not, you will have points in the centers of your projections and the result is not a "simple" polygon.

Answer 2

The first step will be to use convhull (as yar suggests) to get an outline of each projected polygonal region. It should be noted that a convex hull is appropriate to use here since you are dealing with cuboids, which are convex objects. I think you have an error in the coordinates for your second cuboid (located in A(9:16, :)), so I modified your code to the following:

A = [-100   -40    50
     -100   -40     0
     -120   -40    50
     -120   -40     0
     -100     5     0
     -100     5    50
     -120     5    50
     -120     5     0
     -100     0    52
     -100     0     5
       20     0    52
       20     0     5
     -100    50     5
     -100    50    52
       20    50     5
       20    50    52
      -30    70    53
      -30    70     0
        5    70     0
        5    70    53
      -30   120    53
      -30   120     0
        5   120    53
        5   120     0];
theta = 60;
rota = [cosd(theta) -sind(theta) 0; sind(theta) cosd(theta) 0; 0 0 1];
R = A*rota;

And you can generate the polygonal outlines and visualize them like so:

nPerPoly = 8;
nPoly = size(R, 1)/nPerPoly;
xPoly = mat2cell(R(:, 1), nPerPoly.*ones(1, nPoly));
zPoly = mat2cell(R(:, 3), nPerPoly.*ones(1, nPoly));
C = cell(1, nPoly);
for iPoly = 1:nPoly
  P = convhull(xPoly{iPoly}, zPoly{iPoly});
  xPoly{iPoly} = xPoly{iPoly}(P);
  zPoly{iPoly} = zPoly{iPoly}(P);
  C{iPoly} = P([1:end-1; 2:end].')+nPerPoly.*(iPoly-1);  % Constrained edges, needed later
end

figure();
colorOrder = get(gca, 'ColorOrder');
nColors = size(colorOrder, 1);
for iPoly = 1:nPoly
  faceColor = colorOrder(rem(iPoly-1, nColors)+1, :);
  patch(xPoly{iPoly}, zPoly{iPoly}, faceColor, 'EdgeColor', faceColor, 'FaceAlpha', 0.6);
  hold on;
end
axis equal;
axis off;

And here's the plot:

If you wanted to calculate the area of each polygonal projection and add them up it would be very easy: just change the above loop to capture and sum the second output from the calls to convexhull:

totalArea = 0;
for iPoly = 1:nPoly
  [~, cuboidArea] = convhull(xPoly{iPoly}, zPoly{iPoly});
  totalArea = totalArea+cuboidArea;
end

However, if you want the area of the union of the polygons, you have to account for the overlap. You have a few alternatives. If you have the Mapping Toolbox then you could use the function polybool to get the outline, then use polyarea to compute its area. There are also utilities you can find on the MathWorks File Exchange (such as this and this). I'll show you another alternative here that uses delaunayTriangulation. First we can take the edge constraints C created above to use when creating a triangulation of the projected points:

oldState = warning('off', 'all');
DT = delaunayTriangulation(R(:, [1 3]), vertcat(C{:}));
warning(oldState);

This will automatically create new vertices where the constrained edges intersect. Unfortunately, it will also perform the triangulation on the convex hull of all the points, filling in spots that we don't want filled. Here's what the triangulation looks like:

figure();
triplot(DT, 'Color', 'k');
axis equal;
axis off;

We now have to identify the extra triangles we don't want and remove them. We can do this by finding the centroids of each triangle and using inpolygon to test if they are outside all 3 of our individual cuboid projections. We can then compute the areas of the remaining triangles and sum them up using polyarea, giving us the total area of the projection:

dtFaces = DT.ConnectivityList;
dtVertices = DT.Points;
meanX = mean(reshape(dtVertices(dtFaces, 1), size(dtFaces)), 2);
meanZ = mean(reshape(dtVertices(dtFaces, 2), size(dtFaces)), 2);
index = inpolygon(meanX, meanZ, xPoly{1}, zPoly{1});
for iPoly = 2:nPoly
  index = index | inpolygon(meanX, meanZ, xPoly{iPoly}, zPoly{iPoly});
end
dtFaces = dtFaces(index, :);
xUnion = reshape(dtVertices(dtFaces, 1), size(dtFaces)).';
yUnion = reshape(dtVertices(dtFaces, 2), size(dtFaces)).';
totalArea = sum(polyarea(xUnion, yUnion));

And the total area for this example is:

totalArea =

     9.970392341143055e+03

NOTE: The above code has been generalized for an arbitrary number of cuboids.