1

我一直在努力寻找一个好的解决方案来解决下面描述的问题。我想避免 for 循环，但觉得我的 Matlab 技能不足以做其他事情。我有一个大小为 329x230x105 的 3D 位置矩阵。这定义了 10000x7000x3100 米的 3d 体积。除子体积外，大多数 matix 元素为零：

在此处输入图像描述

我需要构建一个与我的原始矩阵相同大小的掩码，该掩码被划分为大型子矩阵，每个子矩阵定义一个 1000x1000x1000 米的常规子体积，并将 1 分配给包含至少一个元素的子矩阵中的所有元素（非零）来自我的原始矩阵。在 XY 中查看： 在此处输入图像描述

所以最终结果是一个 3D 蒙版（下图 XY），其中标记单元格（红色）内的所有元素都具有值 1 并且外部元素设置为 0 ：

在此处输入图像描述

请注意，我对体积边界框或凸包顶点不感兴趣。

提前谢谢了。

添加信息，回答@grantnz：

好吧，我还不确定以下代码是否适用于所有情况，但这是我所做的（在我的笔记本电脑上需要将近 10 秒）：

 % get subscripts of non-zero 3d grid cells
 [u v w]=ind2sub(size(tmpT),find(tmpT));  % tmpT is the original 230x329x105 matrix
 increm = 30.480061;        % grid cell size in meters
 U=(u-1)*increm;            % convert subscripts to position in meters
 V=(v-1)*increm;
 W=(w-1)*increm;
 U=U/1000;                  % round down to nearest 1000 meter
 V=V/1000;
 W=W/1000;
 U=floor(U);
 V=floor(V);
 W=floor(W);
 U=U*1000;
 V=V*1000;
 W=W*1000;
 U=round(U/increm);        % find subscripts of 1000 meter cells 
 V=round(V/increm);
 W=round(W/increm);
 U(U==0)=1;                % make sure no zero 
 V(V==0)=1;
 W(W==0)=1;
 IX=[U V W];               % make subscript matrix
 [q,i,j]=unique(IX,'rows');  % find unique vectors in subscript matrix
 myMask=zeros(size(tmpT)); % initiate mask matrix
 oneKM = 33;                % this many cells in 1000 meters
 for i=1:length(q)          % for each position vector, set mask (from:from+1000 meter) to 1
    myMask(q(i,1):q(i,1)+oneKM,q(i,2):q(i,2)+oneKM,q(i,3):q(i,3)+oneKM)=1;
 end
4

1 回答 1

1

我想出了以下解决方案：

clear all
clc

x=1:336;
y=1:240;
z=1:120;

[meshy, meshx, meshz] = meshgrid(y, x, z);

% myspace is just your volumentric dataset as a logical matrix
s1 = [168 120 60];
myspace = sqrt((meshx - s1(1)).^2 + (meshy - s1(2)).^2 + (meshz - s1(3)).^2 ) < 15;

% Now let's do the work
div = [6 6 6]; % size of one submatrix

cells_x = ceil(x ./ div(1));
cells_y = ceil(y ./ div(2));
cells_z = ceil(z ./ div(3));

% The order of x and y is correct
[cmeshy, cmeshx, cmeshz] = meshgrid(cells_y, cells_x, cells_z);

% Here, the transformation happens.
volx = cmeshx(myspace);
voly = cmeshy(myspace);
volz = cmeshz(myspace);

newspace = zeros(size(myspace) ./ div);
indices = sub2ind(size(newspace), volx, voly, volz);
newspace(indices) = 1;

vol3d('cdata', newspace);
view(3) ;

该代码假定您的矩阵大小可以被子矩阵大小整除。如果不是这种情况，请相应地截断或填充原始矩阵。

在 div 变量指定子矩阵的大小后，计算单元变量。它们包含原始空间的每个索引，新空间的索引在那里。这些是用meshgrid扩展的。现在，对于空间中的每个坐标，(x, y, z) 新空间的坐标是 (cmeshx(x, y, z), cmeshy(x, y, z), cmeshz(x, y, z) )

然后通过逻辑索引轻松完成转换。现在，vol(xyz) 变量中有 (x, y, z) 元组。为了将这些放入新空间，需要计算和应用线性索引。

于 2013-08-04T16:55:45.377 回答