如何计算 3 维投影的面积?
How do I calculate the area of a 3 dimensional projection?
例如,使用以下代码,我有一个坐标矩阵,其中包含 3 个立方体对象,每个对象由 8 个角定义,总共 24 个坐标。我对我的坐标应用旋转,然后删除 y 坐标以获得 x-z 平面中的投影。如何计算这些立方体在 x-z 平面中的面积,忽略间隙并考虑重叠?我试过使用 polyarea
,但这似乎不起作用。
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
polyarea
是正确的方法,但您需要在每个投影的 convex hull 上调用它。否则,您的投影中心将有点,结果不是 "simple" 多边形。
第一步是使用 convhull
(as ) 获取每个投影多边形区域的轮廓。应该注意的是,这里适合使用凸包,因为您要处理的是长方体,它们是凸物体。我认为你的第二个长方体(位于A(9:16, :)
)的坐标有误,所以我将你的代码修改为以下内容:
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;
您可以生成多边形轮廓并像这样可视化它们:
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;
剧情如下:
如果你想计算每个多边形投影的面积并将它们相加,这将非常简单:只需更改上面的循环以捕获调用 convexhull
的第二个输出并将其求和:
totalArea = 0;
for iPoly = 1:nPoly
[~, cuboidArea] = convhull(xPoly{iPoly}, zPoly{iPoly});
totalArea = totalArea+cuboidArea;
end
但是,如果您想要多边形联合的面积,则必须考虑重叠。你有几个选择。如果你有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
。首先,我们可以采用上面创建的边约束 C
在创建投影点的三角剖分时使用:
oldState = warning('off', 'all');
DT = delaunayTriangulation(R(:, [1 3]), vertcat(C{:}));
warning(oldState);
这将在约束边相交处自动创建新顶点。不幸的是,它还会对所有点的凸包进行三角剖分,填充我们不想填充的点。这是三角剖分的样子:
figure();
triplot(DT, 'Color', 'k');
axis equal;
axis off;
我们现在必须找出不需要的多余三角形并将其移除。我们可以通过找到每个三角形的质心并使用 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
来完成此操作,从而得到投影的总面积:
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));
这个例子的总面积是:
totalArea =
9.970392341143055e+03
注意:以上代码已针对任意数量的长方体进行推广。
例如,使用以下代码,我有一个坐标矩阵,其中包含 3 个立方体对象,每个对象由 8 个角定义,总共 24 个坐标。我对我的坐标应用旋转,然后删除 y 坐标以获得 x-z 平面中的投影。如何计算这些立方体在 x-z 平面中的面积,忽略间隙并考虑重叠?我试过使用 polyarea
,但这似乎不起作用。
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
polyarea
是正确的方法,但您需要在每个投影的 convex hull 上调用它。否则,您的投影中心将有点,结果不是 "simple" 多边形。
第一步是使用 convhull
(as A(9:16, :)
)的坐标有误,所以我将你的代码修改为以下内容:
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;
您可以生成多边形轮廓并像这样可视化它们:
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;
剧情如下:
如果你想计算每个多边形投影的面积并将它们相加,这将非常简单:只需更改上面的循环以捕获调用 convexhull
的第二个输出并将其求和:
totalArea = 0;
for iPoly = 1:nPoly
[~, cuboidArea] = convhull(xPoly{iPoly}, zPoly{iPoly});
totalArea = totalArea+cuboidArea;
end
但是,如果您想要多边形联合的面积,则必须考虑重叠。你有几个选择。如果你有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
。首先,我们可以采用上面创建的边约束 C
在创建投影点的三角剖分时使用:
oldState = warning('off', 'all');
DT = delaunayTriangulation(R(:, [1 3]), vertcat(C{:}));
warning(oldState);
这将在约束边相交处自动创建新顶点。不幸的是,它还会对所有点的凸包进行三角剖分,填充我们不想填充的点。这是三角剖分的样子:
figure();
triplot(DT, 'Color', 'k');
axis equal;
axis off;
我们现在必须找出不需要的多余三角形并将其移除。我们可以通过找到每个三角形的质心并使用 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
来完成此操作,从而得到投影的总面积:
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));
这个例子的总面积是:
totalArea =
9.970392341143055e+03
注意:以上代码已针对任意数量的长方体进行推广。