clustercoords.m

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

function [cc, cl] = clustercoords(vol, meth)
% clustercoords  - get list of clusters in volume
%
% FORMAT:       [cc, cl] = clustercoords(vol [, meth])
%
% Input fields:
%
%       vol         XxYxZ logical array, true for indices to be clustered
%       meth        either of 'face', 'edge', or 'vertex' (default: face)
%
% Output fields:
%
%       cc          cluster coords, Cx4 list of coordinates, where
%                   the 4th column is a cluster index within 1...N
%       cl          1xN cell array with Px3 coordinates
%
% Note: the vertex based method is not yet implemented

% Version:  v0.7b
% Build:    7092510
% Date:     Sep-25 2007, 10:13 AM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   ~islogical(vol) || ...
    numel(size(vol)) ~= 3 || ...
    any(size(vol) < 2)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
st = size(vol);
sx = st(1);
sy = st(2);
sz = st(3);
if nargin < 2 || ...
   ~ischar(meth) || ...
    isempty(meth)
    meth = 'f';
else
    meth = lower(meth(1));
end
if ~any('efv' == meth)
    meth = 'f';
end

% find matching indices and their number
cf = find(vol(:));
cs = numel(cf);

% build some intermediate arrays => cluster list
ci = int32(0);
ci(sx, sy, sz) = 0;
ci(cf) = cf;

% coordinate list
[xi, yi, zi] = ind2sub(st, cf);
cc = [xi(:), yi(:), zi(:), zeros(cs, 1)];

% cluster voxel lists
cl = {[]};
cl(1:sx, 1:sy, 1:sz) = cl;
for k = 1:cs
    cl{cf(k)} = cf(k);
end

% switch on method
switch (meth)
    
    % edge connectivity
    case {'e'}
        
        % iterate over found voxels
        for k = 1:cs
            
            % get coordinates
            cxi = xi(k);
            cyi = yi(k);
            czi = zi(k);
            cxp = cxi + 1;
            cyp = cyi + 1;
            czp = czi + 1;
            
            % get initial cluster number
            cn = ci(cxi, cyi, czi);
            
            % check only three connecting faces in positive direction
            % start with x, if x + 1 is also to cluster
            if cxi < sx && ...
                vol(cxp, cyi, czi)
            
                % get other cluster number
                ncn = ci(cxp, cyi, czi);
                
                % if mismatch, unify
                if cn ~= ncn
                    [si, nc, nl] = unitecluster(cl, cn, ncn);
                    ci(si) = nc;
                    cl{nc} = nl;
                    cn = nc;
                end
            end
            
            % then do the same for y and z dirs and their combinations
            if cyi < sy && ...
                vol(cxi, cyp, czi)
                ncn = ci(cxi, cyp, czi);
                if cn ~= ncn
                    [si, nc, nl] = unitecluster(cl, cn, ncn);
                    ci(si) = nc;
                    cl{nc} = nl;
                    cn = nc;
                end
            end
            if czi < sz && ...
                vol(cxi, cyi, czp)
                ncn = ci(cxi, cyi, czp);
                if cn ~= ncn
                    [si, nc, nl] = unitecluster(cl, cn, ncn);
                    ci(si) = nc;
                    cl{nc} = nl;
                    cn = nc;
                end
            end
            if cxi < sx && ...
                cyi < sy && ...
                vol(cxp, cyp, czi)
                ncn = ci(cxp, cyp, czi);
                if cn ~= ncn
                    [si, nc, nl] = unitecluster(cl, cn, ncn);
                    ci(si) = nc;
                    cl{nc} = nl;
                    cn = nc;
                end
            end
            if cxi < sx && ...
                czi < sz && ...
                vol(cxp, cyi, czp)
                ncn = ci(cxp, cyi, czp);
                if cn ~= ncn
                    [si, nc, nl] = unitecluster(cl, cn, ncn);
                    ci(si) = nc;
                    cl{nc} = nl;
                    cn = nc;
                end
            end
            if cyi < sy && ...
                czi < sz && ...
                vol(cxi, cyp, czp)
                ncn = ci(cxi, cyp, czp);
                if cn ~= ncn
                    [si, nc, nl] = unitecluster(cl, cn, ncn);
                    ci(si) = nc;
                    cl{nc} = nl;
                end
            end
        end
        
    % face connectivity
    case {'f'}
        
        % iterate over found voxels
        for k = 1:cs
            
            % get coordinates
            cxi = xi(k);
            cyi = yi(k);
            czi = zi(k);
            cxp = cxi + 1;
            cyp = cyi + 1;
            czp = czi + 1;
            
            % get initial cluster number
            cn = ci(cxi, cyi, czi);
            
            % check only three connecting faces in positive direction
            % start with x, if x + 1 is also to cluster
            if cxi < sx && ...
                vol(cxp, cyi, czi)
            
                % get cluster numbers
                ncn = ci(cxp, cyi, czi);
                
                % if mismatch, unify
                if cn ~= ncn
                    [si, nc, nl] = unitecluster(cl, cn, ncn);
                    ci(si) = nc;
                    cl{nc} = nl;
                    cn = nc;
                end
            end
            
            % then do the same for y and z dirs
            if cyi < sy && ...
                vol(cxi, cyp, czi)
                ncn = ci(cxi, cyp, czi);
                if cn ~= ncn
                    [si, nc, nl] = unitecluster(cl, cn, ncn);
                    ci(si) = nc;
                    cl{nc} = nl;
                    cn = nc;
                end
            end
            if czi < sz && ...
                vol(cxi, cyi, czp)
                ncn = ci(cxi, cyi, czp);
                if cn ~= ncn
                    [si, nc, nl] = unitecluster(cl, cn, ncn);
                    ci(si) = nc;
                    cl{nc} = nl;
                end
            end
        end
        
    % vertex connectivity
    case {'v'}
        
        % iterate over found voxels
        for k = 1:cs
        end
end

% find unique indices in ci and set in cc
[uc{1:3}] = unique(ci(cf(:)));
uci = 1:numel(uc{1});
if numel(uc{3}) > 0
    cc(:, 4) = reshape(uci(uc{3}), [numel(uc{3}), 1]);
end

% return if no cl required
if nargout < 2
    return;
end

% generate cluster list
maxc = max(cc(:, 4));
cl = cell(1, maxc);
for c = 1:maxc
    cl{c} = cc((cc(:, 4) == c), 1:3);
end

% sub function unitecluster
function [si, nc, nl] = unitecluster(cl, c1, c2)
    
    % get voxel lists
    cl1 = cl{c1};
    cl2 = cl{c2};
    
    % what cluster is higher, then switch
    if c1 < c2
        si = cl2;
        nc = c1;
        nl = [cl1, cl2];
    else
        si = cl1;
        nc = c2;
        nl = [cl2, cl1];
    end
% end of function [ci, cl] = unitecluster(cl, c1, c2)