floodfill3.m

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

function vol = floodfill3(vol, sx, sy, sz, meth)
% floodfill3  - keep adjacent values from seed point
%
% FORMAT:       vol = floodfill3(vol, sx, sy, sz [, meth])
%
% Input fields:
%
%       vol         XxYxZ logical/uint8 array, true/1 voxels are linked
%       sx, sy, sz  start index (1-based)
%       meth        either of 'face', 'edge', or 'vertex' (default: face)
%
% Output fields:
%
%       vol         same as input array with only one cluster

% Version:  v0.7a
% Build:    7082500
% Date:     Aug-25 2007, 12:03 AM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 4 || ...
   (~islogical(vol) && ...
    ~strcmpi(class(vol), 'uint8')) || ...
    numel(size(vol)) ~= 3 || ...
    any(size(vol) < 2) || ...
   ~isa(sx, 'double') || ...
    numel(sx) ~= 1 || ...
    isnan(sx) || ...
   ~isa(sy, 'double') || ...
    numel(sy) ~= 1 || ...
    isnan(sy) || ...
   ~isa(sz, 'double') || ...
    numel(sz) ~= 1 || ...
    isnan(sz) || ...
    sx < 1 || ...
    sx > size(vol, 1) || ...
    sy < 1 || ...
    sy > size(vol, 2) || ...
    sz < 1 || ...
    sz > size(vol, 3) || ...
    any([sx, sy, sz] ~= fix([sx, sy, sz]))
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
if nargin < 5 || ...
   ~ischar(meth) || ...
    isempty(meth)
    meth = 'f';
else
    meth = lower(meth(1));
end
if ~any('efv' == meth)
    meth = 'f';
end
st = size(vol);
vx = st(1);
vy = st(2);
vz = st(3);

% check whether sx, sy, sz are in volume
if ~vol(sx, sy, sz)
    vol(:) = false;
    return;
end

% convert vol to uint8
vol = uint8(vol);
vol(sx, sy, sz) = 2;

% and list of directions (1-3 face, 1-9 edge, 1-13 vertices)
%  1: (  +x  ,  -x  )  2: (  +y  ,  -y  )  3: (  +z  ,  -z  ) 
%  4: ( +x+y , -x-y )  5: ( +x-y , -x+y )
%  6: ( +x+z , -x-z )  7: ( +x-z , -x+z )
%  8: ( +y+z , -y-z )  9: ( +y-z , -y+z )
% 10: (+x+y+z,-x-y-z) 11: (+x+y-z,-x-y+z)
% 12: (+x-y+z,-x+y-z) 13: (+x-y-z,-x+y+z)
dl = [...
     1,  0,  0;  0,  1,  0;  0,  0,  1; ...
     1,  1,  0;  1, -1,  0; ...
     1,  0,  1;  1,  0, -1; ...
     0,  1,  1;  0,  1, -1; ...
     1,  1,  1;  1,  1, -1; ...
     1, -1,  1;  1, -1, -1  ...
    ];
if meth == 'f'
    md = 3;
elseif meth == 'e'
    md = 9;
else
    md = 13;
end

% setup hook point arrays
hx = zeros(1, md * ceil(numel(vol) ^ 0.66));
hy = zeros(1, md * ceil(numel(vol) ^ 0.66));
hz = zeros(1, md * ceil(numel(vol) ^ 0.66));

% iterate while new points found
while ~isempty(sx)
    
    % get indices of to checking voxels
    hc = 1;
    for dc = 1:md
        dr = dl(dc, :);
        dx = dr(1);
        dy = dr(2);
        dz = dr(3);
        if dx ~= 0
            tx = [sx - dx, sx + dx];
        else
            tx = [sx, sx];
        end
        if dy ~= 0
            ty = [sy - dy, sy + dy];
        else
            ty = [sy, sy];
        end
        if dz ~= 0
            tz = [sz - dz, sz + dz];
        else
            tz = [sz, sz];
        end
        gt = find(tx > 0 & tx <= vx & ty > 0 & ty < vy & tz > 0 & tz < vz);
        gs = numel(gt);
        hx(hc:(hc + gs - 1)) = tx(gt);
        hy(hc:(hc + gs - 1)) = ty(gt);
        hz(hc:(hc + gs - 1)) = tz(gt);
        hc = hc + gs;
    end
    hc = hc - 1;
    tc = unique(sub2ind(st, hx(1:hc), hy(1:hc), hz(1:hc)));
    
    % only accept where further coordinates
    tc = tc(vol(tc) == 1);
    
    % set these to 2 as well
    vol(tc) = 2;
    
    % then continue with these
    [sx, sy, sz] = ind2sub(st, tc);
end

% convert back
vol = (vol == 2);