In Matlab: make mask from volume's outline

Question

I've been struggling for a while to find a good solution to the problem described below. I want to avoid for loops but feel my Matlab skills are inadequate to do otherwise. I have a 3D position matrix of size 329x230x105. This defines a 3d volume of 10000x7000x3100 meters. Most of the matix elements are zero except for a subvolume:

I need to construct a mask of the same size as my original matrix which is divided in large submatices, each defining a regular subvolume of 1000x1000x1000 meters and assign 1 to all the elements in a submatrix which contains at least one on-element (non-zero) from my original matrix. Viewed in XY :

So the end result is a 3D mask (figure below in XY) where all elements within marked cells (red) have value 1 and elements outside are set to 0 :

Please note that I am not interested in the volumes bounding box or the convex hull vertices.

Many thanks in advance.

Added information, answer to @grantnz:

well, I am not yet sure if the following code does work in all situations, but here is what I do (and it takes nearly 10 seconds on my laptop):

 % 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

Answer 1

I came up with the following solution:

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) ;

The code assumes that your matrix size is evenly divisible by the submatrix size. If this is not the case, truncate or pad your original matrix accordingly.

After the div variable specified the size of the submatrix, the cell-variables are calculated. They contain for each index of the original space, which index of the new space is there. These are expanded with meshgrid. Now, for every coordinate in the space, (x, y, z) the coordinated of the new space are (cmeshx(x, y, z), cmeshy(x, y, z), cmeshz(x, y, z))

The transformation is then easy done with logical indexing. Now, there are (x, y, z) tuples in the vol(xyz) vars. To put these into the new space, linear indexes are calculated and applied.