function [ label, list_num, list, label_num ] = voxels_region_3d ( maxlist, ...
nx, ny, nz, ishow )
%*****************************************************************************80
%
%% voxels_region_3d() arranges a set of voxels into contiguous regions in 3D.
%
% Discussion:
%
% On input, the ISHOW array contains zero and nonzero values. The nonzero
% values are taken to be active voxels.
%
% On output, the LABEL array is a copy of ISHOW, where the zero voxels
% remain zero, and all the active voxels have been assigned a label which
% now indicates membership in a region, or group of contiguous voxels.
%
% On output, the array LIST contains information about the regions.
% The last used element of LIST is LIST_NUM.
%
% The number of elements in region LABEL_NUM is NELEM = LIST(LIST_NUM).
% The (I,J,K) indices of the last element in this region are in
% LIST(LIST_NUM-3) through LIST(LIST_NUM-1), and the first element is
% listed in LIST(LIST_NUM-3*NELEM), LIST(LIST_NUM-3*NELEM+1),
% LIST(LIST_NUM-3*NELEM+2).
%
% The number of elements in region LABEL_NUM-1 is listed in
% LIST(LIST_NUM-3*NELEM-1), and the (I,J,K) indices of these elements
% are listed there.
%
% Picture (for a 2D example):
%
% ISHOW on input:
%
% 0 2 0 0 17 0 3
% 0 0 3 0 1 0 4
% 1 0 4 8 8 0 7
% 3 0 6 45 0 0 0
% 3 17 0 5 9 2 5
%
% LABEL on output:
%
% 0 1 0 0 2 0 3
% 0 0 2 0 2 0 3
% 4 0 2 2 2 0 3
% 4 0 2 2 0 0 0
% 4 4 0 2 2 2 2
%
% Licensing:
%
% This code is distributed under the MIT license.
%
% Modified:
%
% 22 May 2005
%
% Author:
%
% John Burkardt
%
% Input:
%
% integer MAXLIST, the maximum length of the array used to
% list the elements of the regions.
%
% integer NX, NY, NZ, the number of voxels in the X, Y and Z directions.
%
% integer ISHOW(NX,NY,NZ), the only significance to
% the entries is whether they are zero or nonzero.
%
% Output:
%
% integer LABEL(NX,NY,NZ), the nonzero entries have now been revalued
% so that contiguous entries have the same value, indicating a grouping
% into a region.
%
% integer LIST_NUM, the number of entries of LIST that were used.
% However, if MAXLIST < LIST_NUM, then there was not enough space in
% LIST to store the data properly, and LIST should not be used,
% although the data in LABEL should be correct.
%
% integer LIST(MAXLIST), contains, in stack form, a list
% of the indices of the elements in each region.
%
% integer LABEL_NUM, the number of regions discovered.
%
label = ishow;
%
% Reset all nonzero entries to -1.
%
for i = 1 : nx
for j = 1 : ny
for k = 1 : nz
if ( label(i,j,k) ~= 0 )
label(i,j,k) = -1;
end
end
end
end
%
% Start the number of items in the region list at 0.
%
list_num = 0;
%
% Start the number of regions at 0.
%
label_num = 0;
%
% The stack begins empty.
%
nstack = 0;
%
% Search for an unused "ON" voxel from which we can "grow" a new region.
%
for i = 1 : nx
for j = 1 : ny
for k = 1 : nz
%
% We found a voxel that is "ON", and does not belong to any region.
%
if ( label(i,j,k) == -1 )
%
% Increase the number of regions.
%
label_num = label_num + 1;
%
% Add this voxel to the region.
%
label(i,j,k) = label_num;
%
% Add this voxel to the stack.
%
stack(nstack+1) = i;
stack(nstack+2) = j;
stack(nstack+3) = k;
stack(nstack+4) = 1;
nstack = nstack + 4;
%
% Add this voxel to the description of the region.
%
nelements = 1;
if ( list_num + 3 <= maxlist )
list(list_num+1) = i;
list(list_num+2) = j;
list(list_num+3) = k;
end
list_num = list_num + 3;
while ( true )
%
% Find all neighbors of BASE that are "ON" but unused.
% Mark them as belonging to this region, and stack their indices.
%
ibase = stack(nstack-3);
jbase = stack(nstack-2);
kbase = stack(nstack-1);
ilo = max ( ibase-1, 1 );
ihi = min ( ibase+1, nx );
jlo = max ( jbase-1, 1 );
jhi = min ( jbase+1, ny );
klo = max ( kbase-1, 1 );
khi = min ( kbase+1, nz );
nabes = 0;
for i2 = ilo : ihi
for j2 = jlo : jhi
for k2 = klo : khi
%
% We found a neighbor to our current search point, which is "ON" and unused.
%
if ( label(i2,j2,k2) == -1 )
%
% Increase the number of neighbors.
%
nabes = nabes + 1;
%
% Mark the neighbor as belonging to the region.
%
label(i2,j2,k2) = label_num;
%
% Add the neighbor to the stack.
%
stack(nstack+1) = i2;
stack(nstack+2) = j2;
stack(nstack+3) = k2;
nstack = nstack + 3;
%
% Add the neighbor to the description of the region.
%
nelements = nelements + 1;
if ( list_num + 3 <= maxlist )
list(list_num+1) = i2;
list(list_num+2) = j2;
list(list_num+3) = k2;
end
list_num = list_num + 3;
end
end
end
end
%
% If any new neighbors were found, take the last one as the basis
% for a deeper search.
%
if ( 0 < nabes )
stack(nstack+1) = nabes;
nstack = nstack + 1;
continue
end
%
% If the current search point had no new neighbors, drop it from the stack.
%
ncan = stack(nstack) - 1;
nstack = nstack - 3;
stack(nstack) = ncan;
%
% If there are still any unused candidates at this level, take the
% last one as the basis for a deeper search.
%
if ( 0 < stack(nstack) )
continue
end
%
% If there are no more unused candidates at this level, then we need
% to back up a level in the stack. If there are any candidates at
% that earlier level, then we can still do more searching.
%
nstack = nstack - 1;
if ( nstack <= 0 )
break
end
end
%
% If we have exhausted the stack, we have completed this region.
% Tag the number of elements to the end of the region description list.
%
list_num = list_num + 1;
if ( list_num <= maxlist )
list(list_num) = nelements;
end
end
end
end
end
%
% Print some warnings.
%
if ( maxlist < list_num )
fprintf ( 1, '\n' );
fprintf ( 1, 'voxels_region_3d(): Warning!\n' );
fprintf ( 1, ' MAXLIST was too small to list the regions.\n' );
fprintf ( 1, ' Do not try to use the LIST array!\n' );
fprintf ( 1, ' The LABEL data is OK, however.\n' );
end
return
end