在 MATLAB 中细分 3D 表面
Subdividing a 3D surface in MATLAB
我正在使用 MATLAB 创建一些图表来解释我的研究，但遇到了以下问题。我已经创建了一个绕 x 轴旋转的半椭圆体，但我想展示椭圆体内的一些内部结构，这些结构将绕 x 轴和 y 轴旋转。这大致如下图所示，我将椭圆体细分为四个部分（我的绘画技巧不好）。
如何使用围绕 x 轴和 y 轴进行不同旋转的平面将椭圆体细分为多个部分？我可以创建椭圆体和相交的平面，但从那里我不知道如何划分椭圆体并更改面颜色。
我在下面包含了一些基本代码以开始使用。我原以为我可以像将椭圆体切成两半那样遮住部分椭圆体坐标，但这不起作用。我猜我需要制作某种网格，但我不知道如何将相交表面和椭圆体的网格结合起来。
figure
x = 0; y = 0; z = 0;
tl = 10; tw = 4; td = 2;
% Create ellipsoid
[ex,ey,ez] = ellipsoid(x, y, z, tl, tw, td,40);
ex = ex(1:ceil(length(ez)/2),:); % Remove top half
ey = ey(1:ceil(length(ez)/2),:); % of ellipsoid
ez = ez(1:ceil(length(ez)/2),:);
% Make some planes
[ySL,zSL] = meshgrid([-10:10],[-2:0.2:2]);
xSL1 = zeros(size(ySL, 1)); % Generate z data
hSL1 = surf(xSL1,ySL,zSL);
hold on
[ySL,zSL] = meshgrid([-10:10],[-3:0.2:1]);
xSL2 = ones(size(ySL, 1)); % Generate z data
hSL2 = surf(xSL2,ySL,zSL);
% rotate(hSL,[1 0 0],5);
rotate([hSL1 hSL2],[0 1 0],-70);
hSurf1 = surf(ex,ey,ez);
set([hSurf1 hSL1 hSL2],'facecolor','blue','facealpha',.2,...
'edgecolor','none')
% Plot settings
daspect([1 1 0.3]);
hold off
view(-10,6)
非常感谢任何帮助，
这是一个有趣的问题。因此，这里有一个解决方案：1) 在椭圆内绘制三维点，以两个平面为界，以及 2) 将这些点形成一个看起来平滑的补丁表示。
我的基本想法是构建一个覆盖（根据需要密集）整个感兴趣的 X-Y-Z space 的网格。然后我做的是以参数形式 ax+by+cz+d=0 定义两个 planes。这些是要在其间绘制数据的平面。使用平面的法向量，以及关于每个平面上的一个点的信息，我可以使用点积推断 X-Y-Z 网格中的每个点是在平面上方还是下方：如果点积小于零，则点在平面上方，如果大于零，则在平面上方。
然后我使用 logical indexing 找到网格上的哪些点满足条件 1) 在平面 1 上方，2) 在平面 2 下方，以及 3) 在椭圆体内部和 -- 可选 -- 4) 在 Z=0.
最后，我绘制了椭圆体和平面，以及满足上述条件的点。
在本例中，我定义了两个平面，但我想您可以对任意数量的平面执行此操作。该代码快速而肮脏，但希望能为您指明前进的道路！
% Define the ellipse
x = 0; y = 0; z = 0;
tl = 10; tw = 4; td = 2;
% Create ellipsoid
[ex,ey,ez] = ellipsoid(x, y, z, tl, tw, td,40);
ex = ex(1:ceil(length(ez)/2),:); % Remove top half
ey = ey(1:ceil(length(ez)/2),:); % of ellipsoid
ez = ez(1:ceil(length(ez)/2),:);
% Define a 3D grid over area of interest
xx = linspace(min(ex(:)),max(ex(:)),50);
yy = linspace(min(ey(:)),max(ey(:)),50);
zz = linspace(min(ez(:)),max(ez(:)),50);
[X, Y, Z] = meshgrid(xx, yy, zz);
V = [X(:) Y(:) Z(:)]; % rearrange
% Define two planes (ax + bx + cz + d = 0)
a = [0; 0];
b = [-1; 1];
c = [-1; 1];
d = [-1; -1];
% Normal vectors of the planes
n = [a b c];
n(1,:) = n(1,:) / norm(n(1,:));
n(2,:) = n(2,:) / norm(n(2,:));
% Find a point (0,0,z) on each plane
zp = [zeros(2) (d- a * 0 - b * 0)./c];
% Define the area of interest.
% Dot product to test if grid points are below or above the planes
% We want: above plane 1, below plane 2.
V1 = sum(bsxfun(@times, bsxfun(@minus, V, -zp(1,:)), n(1,:)),2);
V2 = sum(bsxfun(@times, bsxfun(@minus, V, -zp(2,:)), n(2,:)),2);
between_planes = (V1 < 0) & (V2 < 0);
% ...and the points have to be inside the ellipsoid
in_ellipse = ((X(:) - x)/tl).^2 + ((Y(:)-y)/tw).^2 + ((Z(:)-z)/td).^2 < 1;
% Final AOI
aoi = between_planes & in_ellipse;
% Add this if you want to also have only values with Z < 0
aoi = aoi & (V(:,3) < 0);
figure;
surf(ex, ey, ez, 'facecolor','blue','facealpha',.2,...
'edgecolor','none')
所以，现在我们准备好绘制一些数据了！首先，让我们尝试绘制网格以获得一种显示单个数据点的离散可视化：
hold on;
plot3(V(aoi,1), V(aoi, 2), V(aoi, 3), 'r.')
...只是为了看看它是否有效，让我们将之前定义的平面添加到可视化中。这是您最终可能想要摆脱的东西。
% Draw the planes
[X,Y] = meshgrid(xx,yy);
Z = cell(2,1);
clr = {'r', 'g'};
for k = 1:2
Z{k} = (a(k) * X + b(k) * Y + d(k))/ (-c(k));
hs(k) = surf(xx, yy, Z{k},'facecolor',clr{k}, 'facealpha',.2, 'edgecolor','none');
drawnow;
end
view([-115 10])
legend(hs, 'Plane 1 (above)', 'Plane 2 (below');
好吧，最后，如果我们想要显示一些看起来平滑的数据，我们可以通过使用 plot 来使用 alphaShape to extract a polyhedra from those data points we have. We can visualize the polyhedra as a patch 对象。最后，我们可以为补丁中的每个顶点分配我们想要的任何颜色，通过设置 'Facecolor' 属性 设置可视化颜色的平滑变化，并通过设置 [= 删除讨厌的顶点边缘20=] 到 'none':
% Moving on to smooth representation...
shp = alphaShape(V(aoi,1), V(aoi, 2), V(aoi, 3),1);
hs = plot(shp);
% Create a colormap to use for the patch (color of each vertex)
cmap = jet(size(hs.Vertices,1));
set(hs, 'FaceVertexCData', cmap);
set(hs, 'FaceColor', 'interp'); % smooth interp coloring
set(hs, 'EdgeColor', 'none'); % Remove ugly edges
就是这样。剩下的就是用适合您的数据的任何值替换顶点颜色值。此外，如果您不希望这些数据点沿着补丁显示在显示器中，请记住提前跳过 plot3 部分。
下面的示例最终结果，也显示了两个平面。
我正在使用 MATLAB 创建一些图表来解释我的研究，但遇到了以下问题。我已经创建了一个绕 x 轴旋转的半椭圆体，但我想展示椭圆体内的一些内部结构，这些结构将绕 x 轴和 y 轴旋转。这大致如下图所示，我将椭圆体细分为四个部分（我的绘画技巧不好）。
如何使用围绕 x 轴和 y 轴进行不同旋转的平面将椭圆体细分为多个部分？我可以创建椭圆体和相交的平面，但从那里我不知道如何划分椭圆体并更改面颜色。
我在下面包含了一些基本代码以开始使用。我原以为我可以像将椭圆体切成两半那样遮住部分椭圆体坐标，但这不起作用。我猜我需要制作某种网格，但我不知道如何将相交表面和椭圆体的网格结合起来。
figure
x = 0; y = 0; z = 0;
tl = 10; tw = 4; td = 2;
% Create ellipsoid
[ex,ey,ez] = ellipsoid(x, y, z, tl, tw, td,40);
ex = ex(1:ceil(length(ez)/2),:); % Remove top half
ey = ey(1:ceil(length(ez)/2),:); % of ellipsoid
ez = ez(1:ceil(length(ez)/2),:);
% Make some planes
[ySL,zSL] = meshgrid([-10:10],[-2:0.2:2]);
xSL1 = zeros(size(ySL, 1)); % Generate z data
hSL1 = surf(xSL1,ySL,zSL);
hold on
[ySL,zSL] = meshgrid([-10:10],[-3:0.2:1]);
xSL2 = ones(size(ySL, 1)); % Generate z data
hSL2 = surf(xSL2,ySL,zSL);
% rotate(hSL,[1 0 0],5);
rotate([hSL1 hSL2],[0 1 0],-70);
hSurf1 = surf(ex,ey,ez);
set([hSurf1 hSL1 hSL2],'facecolor','blue','facealpha',.2,...
'edgecolor','none')
% Plot settings
daspect([1 1 0.3]);
hold off
view(-10,6)
非常感谢任何帮助，
这是一个有趣的问题。因此，这里有一个解决方案：1) 在椭圆内绘制三维点，以两个平面为界，以及 2) 将这些点形成一个看起来平滑的补丁表示。
我的基本想法是构建一个覆盖（根据需要密集）整个感兴趣的 X-Y-Z space 的网格。然后我做的是以参数形式 ax+by+cz+d=0 定义两个 planes。这些是要在其间绘制数据的平面。使用平面的法向量，以及关于每个平面上的一个点的信息，我可以使用点积推断 X-Y-Z 网格中的每个点是在平面上方还是下方：如果点积小于零，则点在平面上方，如果大于零，则在平面上方。
然后我使用 logical indexing 找到网格上的哪些点满足条件 1) 在平面 1 上方，2) 在平面 2 下方，以及 3) 在椭圆体内部和 -- 可选 -- 4) 在 Z=0.
最后，我绘制了椭圆体和平面，以及满足上述条件的点。
在本例中，我定义了两个平面，但我想您可以对任意数量的平面执行此操作。该代码快速而肮脏，但希望能为您指明前进的道路！
% Define the ellipse
x = 0; y = 0; z = 0;
tl = 10; tw = 4; td = 2;
% Create ellipsoid
[ex,ey,ez] = ellipsoid(x, y, z, tl, tw, td,40);
ex = ex(1:ceil(length(ez)/2),:); % Remove top half
ey = ey(1:ceil(length(ez)/2),:); % of ellipsoid
ez = ez(1:ceil(length(ez)/2),:);
% Define a 3D grid over area of interest
xx = linspace(min(ex(:)),max(ex(:)),50);
yy = linspace(min(ey(:)),max(ey(:)),50);
zz = linspace(min(ez(:)),max(ez(:)),50);
[X, Y, Z] = meshgrid(xx, yy, zz);
V = [X(:) Y(:) Z(:)]; % rearrange
% Define two planes (ax + bx + cz + d = 0)
a = [0; 0];
b = [-1; 1];
c = [-1; 1];
d = [-1; -1];
% Normal vectors of the planes
n = [a b c];
n(1,:) = n(1,:) / norm(n(1,:));
n(2,:) = n(2,:) / norm(n(2,:));
% Find a point (0,0,z) on each plane
zp = [zeros(2) (d- a * 0 - b * 0)./c];
% Define the area of interest.
% Dot product to test if grid points are below or above the planes
% We want: above plane 1, below plane 2.
V1 = sum(bsxfun(@times, bsxfun(@minus, V, -zp(1,:)), n(1,:)),2);
V2 = sum(bsxfun(@times, bsxfun(@minus, V, -zp(2,:)), n(2,:)),2);
between_planes = (V1 < 0) & (V2 < 0);
% ...and the points have to be inside the ellipsoid
in_ellipse = ((X(:) - x)/tl).^2 + ((Y(:)-y)/tw).^2 + ((Z(:)-z)/td).^2 < 1;
% Final AOI
aoi = between_planes & in_ellipse;
% Add this if you want to also have only values with Z < 0
aoi = aoi & (V(:,3) < 0);
figure;
surf(ex, ey, ez, 'facecolor','blue','facealpha',.2,...
'edgecolor','none')
所以，现在我们准备好绘制一些数据了！首先，让我们尝试绘制网格以获得一种显示单个数据点的离散可视化：
hold on;
plot3(V(aoi,1), V(aoi, 2), V(aoi, 3), 'r.')
...只是为了看看它是否有效，让我们将之前定义的平面添加到可视化中。这是您最终可能想要摆脱的东西。
% Draw the planes
[X,Y] = meshgrid(xx,yy);
Z = cell(2,1);
clr = {'r', 'g'};
for k = 1:2
Z{k} = (a(k) * X + b(k) * Y + d(k))/ (-c(k));
hs(k) = surf(xx, yy, Z{k},'facecolor',clr{k}, 'facealpha',.2, 'edgecolor','none');
drawnow;
end
view([-115 10])
legend(hs, 'Plane 1 (above)', 'Plane 2 (below');
好吧，最后，如果我们想要显示一些看起来平滑的数据，我们可以通过使用 plot 来使用 alphaShape to extract a polyhedra from those data points we have. We can visualize the polyhedra as a patch 对象。最后，我们可以为补丁中的每个顶点分配我们想要的任何颜色，通过设置 'Facecolor' 属性 设置可视化颜色的平滑变化，并通过设置 [= 删除讨厌的顶点边缘20=] 到 'none':
% Moving on to smooth representation...
shp = alphaShape(V(aoi,1), V(aoi, 2), V(aoi, 3),1);
hs = plot(shp);
% Create a colormap to use for the patch (color of each vertex)
cmap = jet(size(hs.Vertices,1));
set(hs, 'FaceVertexCData', cmap);
set(hs, 'FaceColor', 'interp'); % smooth interp coloring
set(hs, 'EdgeColor', 'none'); % Remove ugly edges
就是这样。剩下的就是用适合您的数据的任何值替换顶点颜色值。此外，如果您不希望这些数据点沿着补丁显示在显示器中，请记住提前跳过 plot3 部分。
下面的示例最终结果，也显示了两个平面。