vmp_clustertable.m

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

vsx = vsz(1);
vsy = vsz(2);
vxy = vsx * vsy;
vxo = vsx + vxy;
for c = numel(cl):-1:1
    
    % apply cluster thresholding
    if size(cl{c}, 1) < msz
        cl(c) = [];
    else
        lcl = cl{c};
        cl{c}(:, 4) = vmp(lcl(:, 1) + vsx * lcl(:, 2) + vxy * lcl(:, 3) - vxo);
        
        % switch on maptype
        switch (map.Type)
            
            % t-Map
            case {1}
                if showboth
                    cl{c}(:, 5) = 2 - 2 * custom_tcdf(abs(cl{c}(:, 4)), map.DF1);
                elseif opts.showpos
                    cl{c}(:, 5) = 1 - custom_tcdf(cl{c}(:, 4), map.DF1);
                else
                    cl{c}(:, 5) = 1 - custom_tcdf(-cl{c}(:, 4), map.DF1);
                end
                
            % F-Map
            case {4}
                cl{c}(:, 5) = 1 - custom_fcdf(cl{c}(:, 4), map.DF1, map.DF2);
                
            % in other cases
            otherwise
                cl{c}(:, 5) = zeros(size(cl{c}, 1), 1);
        end
    end
end

% fill output structure
ncls = numel(cl);
sortcls = zeros(1, ncls);
for c = 1:ncls
    
    % cluster size
    csz = size(cl{c}, 1);
    cs(c).ClusterSize = csz;
    cs(c).ClusterSizeAnat = csz * vrs3;
    
    % coordinates depend on system
    lcl = cl{c}(:, 1:3);
    tcl = 128 - (vres .* (lcl(end:-1:1, [3, 1, 2]) - 1) + ...
        repmat(voff([3, 1, 2]), [csz, 1]));
    tcm = mean(tcl);
    switch (opts.csystem)
        case {'data'}
            cs(c).CoordinateList = lcl;
        case {'int'}
            cs(c).CoordinateList = vres .* (lcl - 1) + repmat(voff, [csz, 1]) - 1;
        case {'sys'}
            cs(c).CoordinateList = vres .* (lcl(:, [3, 1, 2]) - 1) + ...
                repmat(voff([3, 1, 2]), [csz, 1]);
        case {'tal'}
            cs(c).CoordinateList = tcl;
    end
    cs(c).MeanCoordinate = mean(cs(c).CoordinateList);
    
    % map values
    clv = cl{c}(:, 4);
    cs(c).MapValues = clv;
    cs(c).MeanMapValue = mean(clv);
    if cs(c).MeanMapValue > 0
        cs(c).MaxMapValue = max(clv);
        cs(c).MinMapValue = min(clv);
    else
        cs(c).MaxMapValue = min(clv);
        cs(c).MinMapValue = max(clv);
    end
    
    % p values if possible
    if any([1, 4] == map.Type)
        cs(c).MeanPValue = mean(cl{c}(:, 5));
    else
        cs(c).MeanPValue = 0;
    end
    
    % FDR values only if factor is known
    if mtfactfdr > 0
        cs(c).MeanPValueFDR = cs(c).MeanPValue * mtfactfdr;
    else
        cs(c).MeanPValueFDR = 0;
    end
    
    % Bonferroni MCP
    cs(c).MeanPValueBonf = cs(c).MeanPValue * mtfactbonf;
    
    % center of gravity required?
    if opts.gravcent
        cog = mean(lcl .* repmat(clv, [1, 3])) ./ mean(clv);
        tcg = 129 - (vres .* cog([3, 1, 2]) + voff([3, 1, 2]));
        % again depends on coordinate system
        switch (opts.csystem)
            case {'data'}
                cs(c).CenterOfGravity = round(cog);
            case {'int'}
                cs(c).CenterOfGravity = vres .* cog + voff;
            case {'sys'}
                cs(c).CenterOfGravity = vres .* cog([3, 1, 2]) + voff([3, 1, 2]);
            case {'tal'}
                cs(c).CenterOfGravity = tcg;
        end
    end
    
    % perform Tal Lookup?
    if opts.tallabels
        
        % on center of gravity if done, otherwise on mean
        if opts.gravcent
            tcm = tcg;
        end
        
        % perform NGM lookup
        if opts.talngm
            tlp = tdlocal2(-5, tcm(1), tcm(2), tcm(3), [10,10]);
            
        % or standard (label) lookup
        else
            tlp = tdlocal2(-2, tcm(1), tcm(2), tcm(3));
            if ischar(tlp)
                tlp = splittocell(tlp, ',');
            end
            if numel(tlp) ~= 5
                tlp = '*';
            end
        end
        
        % non empty cell
        if iscell(tlp) && ...
           ~isempty(tlp)
       
            % pack result of label into NGM form
            if ~iscell(tlp{1})
                tlp = {{tcm, [0, 0, 0], 0, tlp}};
            end
            
        % make substitute for bad results
        else
            tlp = {{[], [], [], {'*', '*', '*', '*', '*'}}};
        end
        for tc = 1:numel(tlp)
            if ~iscell(tlp{tc}{4})
                tlp{tc}{4} = splittocell(tlp{tc}{4}, ',');
            end
        end
        cs(c).TalLookup = tlp;
        tld = tlp{1}{3};
        tlp = tlp{1}{4};
        taltext = '';
        if ~isempty(tlp{1}) && ...
           ~strcmp(tlp{1}, '*')
            talhemi = splittocell(tlp{1}, ' ');
            if isempty(talhemi)
                talhemi = '';
            else
                talhemi = [talhemi{1}(1) 'H '];
            end
            if ~isempty(tlp{5}) && ...
               ~strcmp(tlp{5}, '*')
                taltext = [talhemi tlp{5}];
            elseif ~isempty(tlp{3}) && ...
               ~strcmp(tlp{3}, '*')
                taltext = [talhemi tlp{3}];
            else
                taltext = [talhemi tlp{2}];
            end
        end
        if ~isempty(tld) && ...
            tld > 0
            taltext = sprintf('%s (dist: %3.1fmm)', taltext, tld);
        end
        cs(c).TalLabel = taltext;
    end
    
    % sorting for later
    switch (opts.sorting)
        
        % coordinate
        case {'coord', 'coordx', 'coordy', 'coordz'}
            
            % depends on direction
            switch (opts.sorting(end))
                case {'y'}
                    sortcls(c) = tcm(2);
                case {'z'}
                    sortcls(c) = tcm(3);
                otherwise
                    sortcls(c) = tcm(1);
            end
            
        % maxstat
        case {'maxstat'}
            
            % use maxmapvalue
            sortcls(c) = abs(cs(c).MaxMapValue);
            
        % size
        case {'size'}
            
            % use cluster size
            sortcls(c) = csz;
            
        % weighted
        case {'weighted'}
            
            % use size * mean value
            sortcls(c) = csz * abs(cs(c).MeanMapValue);
    end
end

% sort clusters
[sortcli{1:2}] = sort(sortcls);
cs = cs(sortcli{2}(end:-1:1));

% only print table if we must
switch (opts.tformat)
    
    % human readable table
    case {'human'}
        
        % start with header
        switch (map.Type)
            case {1}
                thtype = 't-Map';
                thdegf = sprintf('%d', map.DF1);
            case {2}
                thtype = 'correlation Map';
                thdegf = sprintf('%d', map.DF1);
            case {3}
                thtype = 'CC Map';
                thdegf = sprintf('%d', map.DF1);
            case {4}
                thtype = 'F-Map';
                thdegf = sprintf('%d, %d', map.DF1, map.DF1);
            case {11}
                thtype = 'PSC Map';
                thdegf = 'n/a';
            case {12}
                thtype = 'z(ica) Map';
                thdegf = sprintf('%d', map.DF1);
            case {13}
                thtype = 'Thickness Map';
                thdegf = 'n/a';
            otherwise
                thtype = 'unknown';
        end
        thead = cell(1, 9);
        thead{1} = sprintf('   Clustertable of map:   "%s"', map.Name);
        thead{2} = sprintf('           Type of map:   %s', thtype);
        thead{3} = sprintf('    Degrees of freedom:   %s', thdegf);
        thead{4} = sprintf('   Cluster k-threshold:   %d mm^3', opts.minsize * vrs3);
        thead{5} = sprintf(' Applied map threshold:   %.5f', mthresh);
        if mtfactfdr > 0
            thead{6} = sprintf('  FDR threshold factor:   %-6.1f', mtfactfdr);
        else
            thead{6} = '';
        end
        thead{7} = '';
        if mtfactfdr > 0
            thead{8} = '   x    y    z  |    k  |       max |  mean p  |  mean q  |  p bonf';
        else
            thead{8} = '   x    y    z  |    k  |       max |  mean p  |  p bonf';
        end
        if opts.tallabels
            thead{8} = [thead{8} ' | TAL label'];
        end
        thead{9} = '----------------------------------------------------------------------';
        thead = gluetostring(thead, char(10), 1);

        % continue with rows
        trows = cell(1, ncls);
        fdrtext = '';
        taltext = '';
        for c = 1:ncls
            if mtfactfdr > 0
                fdrtext = sprintf(' %1.5f  |', cs(c).MeanPValueFDR);
            end
            if opts.tallabels
                taltext = [' | ' cs(c).TalLabel];
            end
            if opts.gravcent
                crd = round(cs(c).CenterOfGravity);
            else
                crd = round(cs(c).MeanCoordinate);
            end
            trows{c} = sprintf('%4d %4d %4d  | %4d  | %9.6f | %1.5f  |%s %1.5f%s', ...
                crd(1), crd(2), crd(3), cs(c).ClusterSize, cs(c).MaxMapValue, ...
                cs(c).MeanPValue, fdrtext, ...
                min(cs(c).MeanPValue * mtfactbonf, 1), taltext);
        end
        
        % compile table
        ctab = sprintf('%s%s\n', thead, gluetostring(trows, char(10), 1));
        
    % machine readable table
    case {'machine'}
end