vmp_clustertable.m

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

function [cs, ctab, vmp] = vmp_ClusterTable(hfile, mapno, threshold, opts)
% VMP::StatsTable  - generate a table with clusters
%
% FORMAT:       [cs, ctab, v] = vmp.ClusterTable(mapno [, thresh [, opts]])
%
% Input fields:
%
%       mapno       map number (1 .. NrOfMaps)
%       thresh      either p/r values (0 .. 1) or t/F value (1 .. Inf)
%                   if not given or 0, uses the LowerThreshold of map
%       opts        optional settings
%        .csystem   how to report coords ('data', 'int', 'sys', {'tal'})
%        .gravcent  flag, also calculate center of gravity (default false)
%        .mcp       MCP for p values ('bonf', {'fdr'}, 'fdrln', 'none')
%        .minsize   minimum cluster size (map resolution, default by map)
%        .showneg   flag, negative values are considered (default false)
%        .showpos   flag, positive values are considered (default true)
%        .sorting   either of 'coord', 'maxstat', 'size', {'weighted'}
%        .tallabels flag, lookup closest talairach label (default false)
%        .talngm    flag, do a wider (nearest GM) search (default false)
%        .tformat   output type ('human', 'machine', {'none'})
%
% Output fields:
%
%       cs          1xC struct with properties of clusters
%       ctab        text output if either human or machine is selected
%       v           thresholded map (in VMP resolution!)
%
% Note: the weighted sorting multiplies the size * mean value, instead
%       of simply saying coord, it's OK to use coordx, coordy or coordz

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

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmp') || ...
   ~isa(mapno, 'double') || ...
    numel(mapno) ~= 1 || ...
    isinf(mapno) || ...
    isnan(mapno) || ...
    mapno < 1
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if mapno > numel(bc.Map)
    error( ...
        'BVQXfile:BadArgument', ...
        'Map number out of bounds.' ...
    );
end

% get map shortcut
mapno = fix(mapno);
map = bc.Map(mapno);
vmp = double(map.VMPData);
voff = [bc.XStart, bc.YStart, bc.ZStart];
vres = bc.Resolution;
vrs3 = vres ^ 3;

% special threshold option
if ~isa(threshold, 'double') || ...
    numel(threshold) ~= 1 || ...
    isinf(threshold) || ...
    isnan(threshold)
    threshold = 0;
end
if threshold == 0
    threshold = map.LowerThreshold;
end

% options check
if nargin < 4 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'csystem') || ...
   ~ischar(opts.csystem) || ...
    isempty(opts.csystem) || ...
   ~any(strmpci(opts.csystem(:)', {'data', 'int', 'sys', 'tal'}))
    opts.csystem = 'tal';
else
    opts.csystem = lower(opts.csystem(:)');
end
if ~isfield(opts, 'gravcent') || ...
   ~islogical(opts.gravcent) || ...
    numel(opts.gravcent) ~= 1
    opts.gravcent = false;
end
if ~isfield(opts, 'mcp') || ...
   ~ischar(opts.mcp) || ...
    isempty(opts.mcp) || ...
   ~any(strcmpi(opts.mcp(:)', {'bonf', 'fdr', 'fdrln', 'none'}))
    opts.mcp = 'fdr';
else
    opts.mcp = lower(opts.mcp(:)');
end
if ~isfield(opts, 'minsize') || ...
   ~isa(opts.minsize, 'double') || ...
    numel(opts.minsize) ~= 1 || ...
    isinf(opts.minsize) || ...
    isnan(opts.minsize) || ...
    opts.minsize < 1 || ...
    opts.minsize > (numel(vmp) / 2)
    opts.minsize = map.ClusterSize;
else
    opts.minsize = fix(opts.minsize);
end
if ~isfield(opts, 'showneg') || ...
   ~islogical(opts.showneg) || ...
    numel(opts.showneg) ~= 1
    opts.showneg = false;
end
if ~isfield(opts, 'showpos') || ...
   ~islogical(opts.showpos) || ...
    numel(opts.showpos) ~= 1
    opts.showpos = true;
end
if ~isfield(opts, 'sorting') || ...
   ~ischar(opts.sorting) || ...
    isempty(opts.sorting) || ...
   ~any(strcmpi(opts.sorting, {'coord', 'coordx', 'coordy', 'coordz', ...
        'maxstat', 'size', 'weighted'}))
    opts.sorting = 'weighted';
else
    opts.sorting = lower(opts.sorting(:)');
end
if ~isfield(opts, 'tallabels') || ...
   ~islogical(opts.tallabels) || ...
    numel(opts.tallabels) ~= 1
    opts.tallabels = false;
end
if ~isfield(opts, 'talngm') || ...
   ~islogical(opts.talngm) || ...
    numel(opts.talngm) ~= 1
    opts.talngm = false;
end
if opts.talngm
    opts.tallabels = true;
end
if ~isfield(opts, 'tformat') || ...
   ~ischar(opts.tformat) || ...
    isempty(opts.tformat) || ...
   ~any(strcmpi(opts.tformat, {'human', 'machine', 'none'}))
    opts.tformat = 'none';
else
    opts.tformat = lower(opts.tformat(:)');
end

% some initial checks on threshold
if threshold < 0
    if nargin < 4
        opts.showneg = true;
        opts.showpos = false;
    end
    threshold = -threshold;
end
if threshold < 0.0001
    opts.mcp = 'none';
end

% what to show in map
vmp(isnan(vmp)) = 0;
if ~opts.showneg
    vmp(vmp < 0) = 0;
end
if ~opts.showpos
    vmp(vmp > 0) = 0;
end
if opts.showneg && ...
    opts.showpos
    showboth = true;
else
    showboth = false;
end
if ~any(vmp(:) ~= 0)
    error( ...
        'BVQXfile:InternalError', ...
        'Empty map or bad selection.' ...
    );
end

% generate cs first
cs = struct;
cs.ClusterSize = 0;
cs.ClusterSizeAnat = 0;
cs.CoordinateList = [0, 0, 0];
cs.MeanCoordinate = [0, 0, 0];
cs.MapValues = [];
cs.MaxMapValue = 0;
cs.MinMapValue = 0;
cs.MeanMapValue = 0;
cs.MeanPValue = 0;
cs.MeanPValueFDR = 0;
cs.MeanPValueBonf = 0;
if opts.gravcent
    cs.CenterOfGravity = [0, 0, 0];
end
if opts.tallabels
    cs.TalLookup = {{[], [], [], {'*', '*', '*', '*', '*'}}};
    cs.TalLabel = '*';
end
ctab = '';

% check the threshold
mtfactfdr = 0;
mtfactbonf = map.BonferroniValue;
if threshold < 1
    
    % against map type
    switch (map.Type)
        
        % t-map
        case {1}
            
            % depends on MCP method
            switch (opts.mcp)
                
                % Bonferroni correction
                case {'bonf'}
                    mthresh = custom_tinv( ...
                        (1 - threshold) / mtfactbonf, map.DF1);
                    
                % FDR correction (only on selected part)
                case {'fdr', 'fdrln'}
                    fdrval = fdr_thresholds( ...
                        1 - custom_tcdf(abs(vmp(vmp(:) ~= 0)), map.DF1), ...
                        threshold, strcmp(opts.mcp, 'fdrln'));
                    mthresh = custom_tinv(1 - fdrval, map.DF1);
                    mthresh = mthresh(end);
                    mtfactfdr = threshold / fdrval;
                    
                % no correction
                case {'none'}
                    mthresh = custom_tinv(1 - threshold, map.DF1);
            end
            
        % r/CC map
        case {2, 3}
            
            % make sure map only contains R values
            vmp = sign(vmp) .* (abs(vmp) - floor(abs(vmp)));
            mthresh = threshold;
            
        % F-Map
        case {4}
            
            % depends on MCP method
            switch (opts.mcp)
                
                % Bonferroni correction
                case {'bonf'}
                    mthresh = custom_finv( ...
                        (1 - threshold) / mtfactbonf, map.DF1, map.DF2);
                    
                % FDR correction
                case {'fdr', 'fdrln'}
                    fdrval = fdr_thresholds( ...
                        1 - custom_fcdf(abs(vmp(vmp(:) ~= 0)), map.DF1, map.DF2), ...
                        threshold, strcmp(opts.mcp, 'fdrln'));
                    mthresh = custom_finv(1 - fdrval, map.DF1, map.DF2);
                    mthresh = mthresh(end);
                    mtfactfdr = threshold / fdrval;
                    
                % no correction
                case {'none'}
                    mthresh = custom_finv(1 - threshold, map.DF1, map.DF2);
            end
            
        % PSC/z_ica/TH map
        case {11, 12, 13}
            mthresh = threshold;
            
        % else
        otherwise
            error( ...
                'BVQXfile:InvalidObject', ...
                'The map type of this VMP is not yet supported.' ...
            );
    end
   
% otherwise use it directly
else
    mthresh = threshold;
end

% threshold map
vmp(abs(vmp) < mthresh) = 0;

% build cluster lists
[ccn{1:2}] = clustercoords(vmp < 0, 'edge');
[ccp{1:2}] = clustercoords(vmp > 0, 'edge');
cl = [ccn{2}, ccp{2}];

% do some checks/calculations on the clusters
msz = opts.minsize;
vsz = size(vmp);