1

以下で説明する問題の適切な解決策を見つけるのにしばらく苦労してきました。for ループを回避したいのですが、Matlab のスキルではそれ以外のことはできないと感じています。サイズ 329x230x105 の 3D 位置マトリックスがあります。これは、10000x7000x3100 メートルの 3D ボリュームを定義します。サブボリュームを除いて、ほとんどの行列要素はゼロです。

ここに画像の説明を入力

元のマトリックスと同じサイズのマスクを作成する必要があります。このマスクは大きなサブマティックに分割され、それぞれが 1000x1000x1000 メートルの通常のサブボリュームを定義し、少なくとも 1 つのオン要素 (非ゼロ) 元のマトリックスから。XY で表示: ここに画像の説明を入力

したがって、最終結果は、マークされたセル (赤) 内のすべての要素の値が 1 で、外側の要素が 0 に設定されている 3D マスク (下の XY の図) です。

ここに画像の説明を入力

ボリューム バウンディング ボックスや凸包の頂点には関心がないことに注意してください。

よろしくお願いします。

追加情報、@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 変数がサブマトリックスのサイズを指定した後、セル変数が計算されます。元のスペースのインデックスごとに、新しいスペースのどのインデックスがあるかが含まれています。これらはメッシュグリッドで展開されます。ここで、空間内のすべての座標 (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 に答える