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.
