% SKEL_STACK Skeletonize a 3D image matrix.
% (trees package)
%
% [i1 i2 i3] = skel_stack (iM, thr, options)
% ------------------------------------------
%
% inspired by algorithms described by Palagyi and Kuba. Nice piece of
% code, hopefully correctly interpreted from their papers. It involves
% numerous permutations of indices and logical operators on a
% 26-neighbourhood.
%
% Input
% -----
% - iM::matrix: contains brightness levels
% - thr::value: threshold value {DEFAULT: so that it starts with 30000
% points}
% - options::string: {DEFAULT: '-s'}
% '-m' : demo movie
% '-c' : closes stack with small window
% '-w' : waitbar
%
% Output
% ------
% - i1, i2, i3:: :coordinates within the stack of remaining points.
%
% Example
% -------
% skel_stack ([], [], '-m -c') % compare without stack closing:
% skel_stack ([], [], '-m')
%
% See also cgui_tree load_stack show_stack
% Uses (-c image processing toolbox)
%
% the TREES toolbox: edit, visualize and analyze neuronal trees
% Copyright (C) 2009 Hermann Cuntz
function [I1 I2 I3] = skel_stack (iM, thr, options)
if (nargin<1)||isempty(iM),
iM = false (250, 250, 50); iM (20 : 230, 120 : 140, 20 : 30) = 1;
iM (ceil (rand (1000, 1) * numel (iM))) = 0;
end
if (nargin<2)||isempty(thr),
minM = min (min (min (iM)));
maxM = max (max (max (iM)));
c = histc (reshape (double (iM), numel (iM), 1), [minM (maxM - minM)./99 maxM]);
cc = cumsum (flipud (c)); ic = min (find (cc > 30000));
thr = (99 - ic) * double ((maxM - minM)/99) + minM;
end
if (nargin<3)||isempty(options),
options = '-w';
end
I1 = []; I2 = []; I3 = [];
nM = ceil (numel (iM) /20000000); % memory problems and speed: cutdown in pieces
% if problems persist go down to 1000000 per piece
cM = floor (size (iM, 1) / nM);
if strfind (options, '-w'), % waitbar option: initialization
HW = waitbar (0, 'skeletonizing...');
set (HW, 'Name', '..PLEASE..WAIT..YEAH..');
end
for ward = 1 : nM,
if strfind (options, '-w'), % waitbar option: update
waitbar (ward / nM, HW);
end
if ward == nM,
MI = iM (1 + (ward - 1) * cM : end, :, :);
else
MI = iM (1 + (ward - 1) * cM : (ward) * cM, :, :);
end
if strfind (options, '-c'),
MI = imclose (MI, ones (2, 2, 2));
end
% add border:
MI = [false(1, size (MI, 2), size (MI, 3)); MI; false(1, size (MI, 2), size (MI, 3))];
MI = [false(size (MI, 1), 1, size (MI, 3)) MI false(size (MI, 1), 1, size (MI, 3))];
MI = cat (3, false (size (MI, 1), size (MI, 2)), MI, false (size (MI, 1), size (MI, 2)));
% binary threshold
MI = MI > thr;
if strfind (options, '-m'),
clf; subplot (211); imagesc (max (iM, [], 3)); colormap gray; axis equal;
set (gca, 'visible', 'off');
end
iFlo = 0; iFl = 1;
while iFlo ~= iFl
iFl = iFlo;
if strfind (options, '-m'),
subplot (212); imagesc (max (MI, [], 3));
axis equal; set (gca, 'visible', 'off');
drawnow;
end
indy = find (MI);
[i1 i2 i3] = ind2sub (size (MI), indy);
i1 = int16 (i1); i2 = int16 (i2); i3 = int16 (i3);
i1NE = repmat([i1-1 i1 i1+1], 1, 9);
i2NE = repmat([i2-1 i2-1 i2-1 i2 i2 i2 i2+1 i2+1 i2+1], 1, 3);
i3NE = [repmat(i3 - 1, 1, 9) repmat(i3, 1, 9) repmat(i3 + 1, 1, 9)];
indo = double (sub2ind (size (MI), double (i1NE), double (i2NE), double (i3NE)));
clear i1NE i2NE i3NE
NE = MI (indo); clear indo;
NEbak = NE;
P = [1 4 7 10 13 16 19 22 25 2 5 8 11 14 17 20 23 26 3 6 9 12 15 18 21 24 27;...
9 18 27 8 17 26 7 16 25 6 15 24 5 14 23 4 13 22 3 12 21 2 11 20 1 10 19;...
7 4 1 8 5 2 9 6 3 16 13 10 17 14 11 18 15 12 25 22 19 26 23 20 27 24 21;...
9 8 7 6 5 4 3 2 1 18 17 16 15 14 13 12 11 10 27 26 25 24 23 22 21 20 19;...
3 6 9 2 5 8 1 4 7 12 15 18 11 14 17 10 13 16 21 24 27 20 23 26 19 22 25;...
3 12 21 6 15 24 9 18 27 2 11 20 5 14 23 8 17 26 1 10 19 4 13 22 7 16 25;...
21 20 19 24 23 22 27 26 25 12 11 10 15 14 13 18 17 16 3 2 1 6 5 4 9 8 7;...
25 22 19 16 13 10 7 4 1 26 23 20 17 14 11 8 5 2 27 24 21 18 15 12 9 6 3;...
1 10 19 2 11 20 3 12 21 4 13 22 5 14 23 6 15 24 7 16 25 8 17 26 9 18 27;...
19 20 21 10 11 12 1 2 3 22 23 24 13 14 15 4 5 6 25 26 27 16 17 18 7 8 9;...
7 8 9 16 17 18 25 26 27 4 5 6 13 14 15 22 23 24 1 2 3 10 11 12 19 20 21];
iNE = find (sum (NE, 2) >= 3 & sum (~NE, 2) >= 5);
iF = false (length (iNE), 1);
for te = 1 : 12,
if te > 1,
NE = NEbak (:, P (te - 1, :));
end
iF1 = (sum (NE (iNE, [19 20 21 22 23 24 25 26 27]), 2) == 0) & ...
( sum (NE (iNE, [1 2 3 4 6 7 8 9 10 11 12 13 15 16 17 18]), 2) > 0) & ...
(NE (iNE, 5) == 1);
iF2 = (sum (NE (iNE, [1 4 7 10 13 16 19 22 25]), 2) == 0) & ...
( sum (NE (iNE, [2 3 5 6 8 9 11 12 17 18 20 21 23 24 26 27]), 2) > 0) & ...
(NE (iNE, 15) == 1);
iF3 = (sum (NE (iNE, [1 4 5 7 10 13 15 16 19 20 21 22 23 24 25 26 27]), 2) == 0) & ...
( sum (NE (iNE ,[2 3 8 9 11 12 17 18]), 2) > 0) & ...
(NE (iNE, 6) == 1);
iF4 = (sum (NE (iNE, [13 19 22 23 25]), 2) == 0) & ...
( sum (NE (iNE, [10 20]), 2) < 2) & ...
( sum (NE (iNE, [17 26]), 2) < 2) & ...
( sum (NE (iNE, [5 15]), 2) == 2);
iF5 = (sum (NE (iNE, [13 19 22 23]), 2) == 0) & ...
( sum (NE (iNE, [10 20]), 2) < 2) & ...
( sum (NE (iNE, [5 15 17]), 2) == 3);
iF6 = (sum (NE (iNE, [13 22 23 25]), 2) == 0) & ...
( sum (NE (iNE, [16 26]), 2) < 2) & ...
( sum (NE (iNE, [5 15 11]), 2) == 3);
iF7 = (sum (NE (iNE, [13 19 22 23]), 2) == 0) & ...
( sum (NE (iNE, [10 20]), 2) < 2) & ...
( sum (NE (iNE, [16 26]), 2) == 1) & ...
( sum (NE (iNE, [5 15 25]), 2) == 3);
iF8 = (sum (NE (iNE, [13 22 23 25]), 2) == 0) & ...
( sum (NE (iNE, [16 26]), 2) < 2) & ...
( sum (NE (iNE, [10 20]), 2) == 1) & ...
( sum (NE (iNE, [5 15 19]), 2) == 3);
iF9 = (sum (NE (iNE, [1 4 5 10 13 19 20 21 22 23 24 25 26 27]), 2) == 0) & ...
( sum (NE (iNE, [6 8]), 2) == 2);
iF10 = (sum (NE (iNE, [4 5 7 13 16 19 20 21 22 23 24 25 26 27]), 2) == 0) & ...
(sum (NE (iNE, [2 6]), 2) == 2);
iF11 = (sum (NE (iNE, [1 4 5 7 10 13 15 16 19 20 21 22 23 24 25]), 2) == 0) & ...
(sum (NE (iNE, [6 18]), 2) == 2);
iF12 = (sum (NE (iNE, [1 4 5 7 10 13 15 16 18 19 22 23 24 25 26 27]), 2) == 0) & ...
(sum (NE (iNE, [6 12]), 2) == 2);
iF = iF | iF1 | iF2 | iF3 | iF4 | iF5 | iF6 | iF7 | iF8 | iF9 | iF10 | iF11 | iF12;
end
MI (indy (iNE (iF))) = 0;
iFlo = length (iF);
end
i1 = double (i1); i2 = double (i2); i3 = double (i3);
I1 = [I1; i1+(ward - 1)*cM]; I2 = [I2; i2]; I3 = [I3; i3];
end
if strfind (options, '-w'), % waitbar option: close
close (HW);
end