2

我试图在 MATLAB 中可视化包含 3D 数组的数据。

该数组的维度为 20*20*40,除了一些元素为零之外,其他所有元素都为零。

我正在寻找一种在散点图中绘制这些非零点的方法,以便非零点相互连接。

这是我到目前为止所做的:

b=zeros(6,6,3);
a=[1 1 1;2 2 1;2 1 1;1 2 1;
1 1 2;1 2 2;2 1 2;2 2 2;
6 6 3;6 5 3;5 6 3;5 5 3;
6 6 2;6 5 2;5 6 2;5 5 2];
[r c]=size(a);
for i = 1:r
 b(a(i,c-2),a(i,c-1),a(i,c)) = 1;
end
[m n o]=size(b);
figure (1)
[x,y,z] = meshgrid(1:m,1:n,1:o);
scatter3(x(:),y(:),z(:),90,b(:),'filled')

所以,我所追求的是能够将八点二中的每一个连接到立方晶格。任何想法都非常感谢。 在此处输入图像描述

编辑: 非常感谢所有提供了很多帮助的专家。现在,我面临另一个记忆问题。

我真实案例的 b 矩阵是 1000*1000*2000 ,我有大小为 4,4700,000*3 的“a”矩阵。“a”矩阵的所有元素都是整数值。虽然我有高达 48GB 的​​内存可用。但是,在上面的“for 循环”中,程序弹回了“内存不足”错误。

非常感谢任何提高内存效率的想法。

4

3 回答 3

3

考虑以下(基于先前的答案):

%# adjacency matrix
adj = false(numel(b));

%# extract first cube, and connect all its points
bb = b;
bb(:,1:3,:) = false;
idx = find(bb);
adj(idx,idx) = true;

%# extract second cube, and connect all its points
bb = b;
bb(:,4:6,:) = false;
idx = find(bb);
adj(idx,idx) = true;

%# points indices
[r c] = find(adj);
p = [r c]';

%# plot
plot3(x(p), y(p), z(p), 'LineWidth',2, 'Color',[.4 .4 1], ...
    'Marker','o', 'MarkerSize',6, ...
    'MarkerFaceColor','g', 'MarkerEdgeColor','g')
axis equal, axis vis3d, grid on, view(3)
xlabel x, ylabel y, zlabel z

在此处输入图像描述


编辑:

手动连接点会更容易:

%# edges: connecting points indices
p = [
    1 2; 2 8; 8 7; 7 1;
    37 38; 38 44; 44 43; 43 37;
    1 37; 2 38; 7 43; 8 44;
    65 66; 66 72; 72 71; 71 65;
    101 102; 102 108; 108 107; 107 101;
    65 101; 66 102; 71 107; 72 108
]';

%# plot
figure
plot3(x(p), y(p), z(p), 'LineWidth',2, 'Color',[.4 .4 1], ...
    'Marker','o', 'MarkerSize',6, ...
    'MarkerFaceColor','g', 'MarkerEdgeColor','g')
axis equal, axis vis3d, grid on, view(3)
xlabel x, ylabel y, zlabel z

%# label points
labels = strtrim(cellstr( num2str((1:numel(b))','%d') ));
idx = find(b);
text(x(idx(:)), y(idx(:)), z(idx(:)), labels(idx(:)), ...
    'Color','m', ...
    'VerticalAlignment','bottom', 'HorizontalAlignment','left')

截屏


编辑#2:(归功于@tmpearce)

您可以半自动化构建边列表的过程。p只需将上述代码中手动构造的 matrx 替换为以下内容:

%# compute edges: pairs of vertex indices 
yIdx = {1:3 4:6};           %# hack to separate each cube points
p = cell(numel(yIdx),1);
for i=1:numel(yIdx)         %# for each cube
    %# find indices of vertices in this cube
    bb = b;
    bb(:,yIdx{i},:) = false;
    idx = find(bb);

    %# compute L1-distance between all pairs of vertices,
    %# and find pairs which are unit length apart
    [r,c] = find(triu(squareform(pdist([x(idx) y(idx) z(idx)],'cityblock'))==1));

    %# store the edges
    p{i} = [idx(r) idx(c)]';
end
p = cat(2,p{:});            %# merge all edges found

这个想法是,对于每个立方体,我们计算所有顶点之间的城市街区距离。“侧”边的距离为1,而对角线的距离大于或等于2

这个过程仍然假设属于每个立方体的顶点列表提供给我们,或者像我们在上面的代码中那样轻松提取......

于 2012-06-17T00:51:02.107 回答
2

我想我会回答 Amro 的回答,并展示如何自动连接沿三次矩阵的边。此答案旨在作为 Amro 的补充,我建议将其合并到该解决方案中,以便将来的观众更容易阅读。

idx从 Amro 的解决方案开始,其中包含立方体角的线性索引。下面我只用一个立方体完成了它,所以矩阵不是很大,我们可以在这里看到它们。

>> idx' %# transposed to save space in the printed output
ans =
 1     2     7     8    37    38    43    44 '

%# now get the subindex of each point in 3D
>> [sx sy sz] = ind2sub([6 6 3],idx);
%# calculate the "edge distance" between each corner (not euclidean distance)
>> dist = abs(bsxfun(@minus,sx,sx')) +...
          abs(bsxfun(@minus,sy,sy')) +...
          abs(bsxfun(@minus,sz,sz'))
dist =

 0     1     1     2     1     2     2     3
 1     0     2     1     2     1     3     2
 1     2     0     1     2     3     1     2
 2     1     1     0     3     2     2     1
 1     2     2     3     0     1     1     2
 2     1     3     2     1     0     2     1
 2     3     1     2     1     2     0     1
 3     2     2     1     2     1     1     0 '

%# find row and column index of points one edge apart
>> [r c]=find(triu(dist==1)); %# triu means don't duplicate edges

%# create 'p' matrix like in Amro's answer
>> p = [idx(r) idx(c)]'
p =
 1     1     2     7     1     2    37     7    37     8    38    43
 2     7     8     8    37    38    38    43    43    44    44    44

为了可视化,您可以像 Amro 的答案一样绘制。

于 2012-06-18T21:14:42.190 回答
0

这是你需要的吗?

b=logical(b);                                  
scatter3(x(b(:)),y(b(:)),z(b(:)),90,'filled')

在此处输入图像描述

于 2012-06-16T11:26:59.763 回答