voi_combine.m

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

function [hfile, res] = voi_Combine(hfile, vspec, cmdname, cmdspec)
% VOI::Combine  - combine regions within a VOI object
%
% FORMAT:       [voi, res] = voi.Combine(vspec, cmdname [,cmdspec])
%
% Input fields:
%
%       vspec       VOI region selection, 1xN list or regexp pattern
%       cmdname     either of
%                   - 'intersect'  (BVQX'S 'a AND b')
%                   - 'overlap'    (probabilistic overlap between VOIs)
%                   - 'setdiff'    (voxels in a NOT b)
%                   - 'union'      (BVQX'S 'a OR b')
%       cmpspec     optional struct with settings
%        .bbox      2x3 mask/bounding box (in BVSystem coordinates)
%        .bothdirs  combine in both dirs (where useful, default false)
%        .cmbsubj   boolean flag, over-subject combination (default true)
%        .nostore   boolean flag, do not store combined VOIs into object,
%                   only return the VOI struct (*as first output*)
%        .thresh    relative value, when e.g. overlap is set (default 0.5)
%
% Output fields:
%
%       voi         altered object (unless nostore is set to true)
%       res         resulting VOI structures (1xN struct)

% Version:  v0.7b
% Build:    7090215
% Date:     Sep-02 2007, 3:08 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 3 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'voi') || ...
   ((~isa(vspec, 'double') || ...
     isempty(vspec) || ...
     any(isinf(vspec(:)) | isnan(vspec(:)) | vspec(:) < 0 | ...
         vspec(:) ~= fix(vspec(:)))) && ...
    (~ischar(vspec) || ...
     isempty(vspec) || ...
     ~any(vspec(:) == '+' | vspec(:) == '*'))) || ...
   ~ischar(cmdname) || ...
    isempty(cmdname) || ...
   ~any(strcmpi(cmdname(:)', {'intersect', 'overlap', 'setdiff', 'union'}))
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
vspec = vspec(:)';
if any(vspec > numel(bc.VOI))
    error( ...
        'BVQXfile:BadArgument', ...
        'Selected VOI(s) out of bounds.' ...
    );
end
cmdname = lower(cmdname(:)');
if nargin < 4 || ...
   ~isstruct(cmdspec) || ...
    numel(cmdspec) ~= 1
    cmdspec = struct;
end
if ~isfield(cmdspec, 'bbox') || ...
   ~isa(cmdspec.bbox, 'double') || ...
    numel(size(cmdspec.bbox)) ~= 2 || ...
    any(size(cmdspec.bbox) ~= [2, 3]) || ...
    any(isnan(cmdspec.bbox(:)) | cmdspec.bbox(:) < 0 | cmdspec.bbox > 255)
    cmdspec.bbox = [];
else
    cmdspec.bbox = round(cmdspec.bbox);
end
bb = cmdspec.bbox;
if ~isfield(cmdspec, 'bothdirs') || ...
   ~islogical(cmdspec.bothdirs) || ...
    numel(cmdspec.bothdirs) ~= 1
    cmdspec.bothdirs = false;
end
if ~isfield(cmdspec, 'cmbsubj') || ...
   ~islogical(cmdspec.cmbsubj) || ...
    numel(cmdspec.cmbsubj) ~= 1
    cmdspec.cmbsubj = true;
end
if ~isfield(cmdspec, 'nostore') || ...
   ~islogical(cmdspec.nostore) || ...
    numel(cmdspec.nostore) ~= 1
    cmdspec.nostore = false;
end
if ~isfield(cmdspec, 'thresh') || ...
   ~isa(cmdspec.thresh, 'double') || ...
    numel(cmdspec.thresh) ~= 1 || ...
    isnan(cmdspec.thresh) || ...
    cmdspec.thresh <= 0 || ...
    cmdspec.thresh > 1
    cmdspec.thresh = 0.5;
end

% check vspec
vst = bc.VOI;
nvoi = numel(vst);
if ischar(vspec)
    vspecname = vspec;
    vspecname(vspec == '.' | vspec == '+' | vspec == '*' | ...
              vspec == '(' | vspec == ')' | vspec == '|' | ...
              vspec == '[' | vspec == ']' | vspec == '-' | ...
              vspec == '{' | vspec == '}' | vspec == ',' | ...
              vspec == '^' | vspec == '$' | vspec == '?') = [];
    nsv = false(1, nvoi);
    for vc = 1:nvoi
        rxr = regexpi(vst(vc).Name, vspec);
        if ~isempty(rxr)
            nsv(vc) = true;
        end
    end
    vspec = find(nsv);
else
    vspecname = '';
end
nsel = numel(vspec);
if nsel < 2
    error( ...
        'BVQXfile:BadArgument', ...
        'VOI selection too small' ...
    );
end

% good new name pattern given
if isfield(bc, 'SubjectVOINamingConvention') && ...
    strcmpi(bc.SubjectVOINamingConvention, '<subj>_<voi>')
    subjvoi = true;
else
    subjvoi = false;
end
if isempty(vspecname)
    vnames = cell(nsel, 1);
    emptyname = false;
    if subjvoi
        vname1 = regexprep(vst(vspec(1)).Name, '^[^_]*_', '');
        if isempty(vname1)
            emptyname = true;
        end
        for vc = 1:nsel
            vnames{vc} = regexprep(vst(vspec(vc)).Name, '^[^_]*_', '');
            if isempty(vnames{vc})
                emptyname = true;
                break;
            end
        end
    else
        vname1 = regexprep(vst(vspec(1)).Name, '_[^_]*$', '');
        if isempty(vname1)
            emptyname = true;
        end
        for vc = 1:nsel
            vnames{vc} = regexprep(vst(vspec(vc)).Name, '_[^_]*$', '');
            if isempty(vnames{vc})
                emptyname = true;
                break;
            end
        end
    end
    if ~emptyname
        for cc = numel(vname1):-1:0
            if cc > 0 && ...
                numel(strmatch(vname1(1:cc), vnames)) == numel(vnames)
                break;
            end
        end
    else
        cc = 0;
    end
    if cc > 0
        vspecname = vname1(1:cc);
    else
        if ~isempty(vname1)
            vspecname = vname1;
        else
            vspecname = vst(vspec(1)).Name;
        end
    end
end

% get coords
vistal = (lower(bc.CoordsType(1)) == 't');
crdnum = zeros(1, nsel);
coords = cell(nsel, 1);
maxc = 1;
minc = 16777216;
for vc = 1:nsel
    c = vst(vspec(vc)).Voxels;
    if vistal
        c = 128 - c;
    end
    if ~isempty(bb)
        c = c(:, ( ...
            c(:, 1) >= bb(1, 1) & c(:, 1) <= bb(2, 1) & ...
            c(:, 2) >= bb(1, 2) & c(:, 2) <= bb(2, 2) & ...
            c(:, 3) >= bb(1, 3) & c(:, 3) <= bb(2, 3)));
    end
    c = round(c(:, 1)) + 256 .* round(c(:, 2)) + 65536 .* round(c(:, 3)) + 1;
    if any(c < 1 | c > 16777216)
        error( ...
            'BVQXfile:InvalidObject', ...
            'Invalid coordinates in VOI (out of VMR space).' ...
        );
    end
    if ~isempty(c)
        maxc = max(max(c), maxc);
        minc = min(min(c), minc);
    end
    coords{vc} = c;
    crdnum(vc) = numel(c);
end
minc = minc - 1;
if minc > 0
    for vc = 1:nsel
        coords{vc} = coords{vc} - minc;
    end
end
minc = minc - 1;

% build resulting array and intermittent arrays
res = vst(vspec(1));
res.NrOfVoxels = 0;
res.Voxels = zeros(0, 3);
msk = false(1, maxc);
msk(coords{1}) = true;
vec = false(1, maxc);

% rest depends of command
switch (cmdname)
    
    % intersect (AND)
    case {'intersect'}
        
        % perform intersection with boolean operators
        for vc = 2:nsel
            vec(:) = false;
            vec(coords{vc}) = true;
            msk = msk & vec;
        end
        res.Voxels = find(msk) + minc;
        res.NrOfVoxels = numel(res.Voxels);
        if subjvoi == cmdspec.cmbsubj
            res.Name = ['ISECT_' vspecname];
        else
            res.Name = [vspecname '_ISECT'];
        end
        
    % overlap
    case {'overlap'}
        
        % perform intersection with boolean operators
        clear vec;
        msk = uint16(msk);
        for vc = 2:nsel
            msk(coords{vc}) = msk(coords{vc}) + uint16(1);
        end
        res.Voxels = find(msk >= uint16(ceil(nsel * cmdspec.thresh))) + minc;
        res.NrOfVoxels = numel(res.Voxels);
        if subjvoi == cmdspec.cmbsubj
            res.Name = ['OVERLAP_' vspecname];
        else
            res.Name = [vspecname '_OVERLAP'];
        end
        
    % setdiff (NOT)
    case {'setdiff'}
        
        % vec not needed
        clear vec;
        
        % without both/all directions
        if ~cmdspec.bothdirs
            for vc = 2:nsel
                msk(coords{vc}) = false;
            end
            res.Voxels = find(msk) + minc;
            res.NrOfVoxels = numel(res.Voxels);
            if subjvoi == cmdspec.cmbsubj
                res.Name = [res.Name '_UNIQUE'];
            else
                res.Name = ['UNIQUE_' res.Name];
            end
            
        % otherwise perform setdiff with each as a start!
        else
            
            % grow result
            res(2:nsel) = res;
            
            % iterate over first part of combination
            for vc = 1:nsel
                
                % setup array
                msk(:) = false;
                msk(coords{vc}) = true;
                
                % iterate over others
                for svc = 1:nsel
                    if svc ~= vc
                        msk(coords{svc}) = false;
                    end
                end
                
                % fill output
                res(vc).Voxels = find(msk) + minc;
                res(vc).NrOfVoxels = numel(res(vc).Voxels);
                if subjvoi == cmdspec.cmbsubj
                    res(vc).Name = [vst(vspec(vc)).Name '_UNIQUE'];
                else
                    res(vc).Name = ['UNIQUE_' vst(vspec(vc)).Name];
                end
            end
        end
        
    % union (OR)
    case {'union'}
        
        % perform union with boolean operators
        for vc = 2:nsel
            vec(coords{vc}) = true;
            msk = msk | vec;
        end
        res.Voxels = find(msk) + minc;
        res.NrOfVoxels = numel(res.Voxels);
        if subjvoi == cmdspec.cmbsubj
            res.Name = ['UNION_' vspecname];
        else
            res.Name = [vspecname '_UNION'];
        end
        
end

% recompute coordinate system
for vc = 1:numel(res)
    c = res(vc).Voxels(:);
    c = [mod(c, 256), mod(floor(c ./ 256), 256), floor(c ./ 65536)];
    if vistal
        c = 128 - c(end:-1:1, :);
    end
    res(vc).Voxels = c;
    res(vc).Name = strrep(res(vc).Name, '__', '_');
end

% only output
if cmdspec.nostore
    hfile = res;
    return;
end

% add to VOI file
bc.VOI = [bc.VOI(:)', res(:)'];
bc.NrOfVOIs = numel(bc.VOI);
bvqxfile_setcont(hfile.L, bc);