vtc_meanstd.m

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

function hfile2 = vtc_MeanStd(hfile, vsel)
% VTC::MeanStd  - calculate Mean and Std over a partition of volumes
%
% FORMAT:       msvmp = vtc.MeanStd([vsel])
%
% Input fields:
%
%       vsel        volume selection, default all volumes
%                   can also be a cell array of selections
%
% Output fields:
%
%       msvmp       VMP with mean and PSC total signal variance image(s)

% Version:  v0.7b
% Build:    7083011
% Date:     Aug-30 2007, 11:57 AM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vtc')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end

% get object reference
bc = bvqxfile_getcont(hfile.L);
nv = bc.NrOfVolumes;
vsz = size(bc.VTCData);

% check vsel
if nargin < 2 || ...
   (~isa(vsel, 'double') && ...
    ~iscell(vsel)) || ...
    isempty(vsel)
    vsel = 1:nv;
end
if isa(vsel, 'double')
    vsel = {vsel};
end
vsel = vsel(:)';
for cc = 1:numel(vsel)
    if ~isa(vsel{cc}, 'double') || ...
        isempty(vsel{cc}) || ...
        any(isinf(vsel{cc}(:)) | isnan(vsel{cc}(:)) | vsel{cc}(:) < 1 | ...
            vsel{cc}(:) > nv | vsel{cc}(:) ~= round(vsel{cc}(:))) || ...
        numel(vsel{cc}(:)) ~= numel(unique(vsel{cc}(:)))
        error( ...
            'BVQXfile:BadArgument', ...
            'Invalid volume selection.' ...
        );
    end
    vsel{cc} = unique(vsel{cc}(:)');
end

% create output files
bbox = aft_BoundingBox(hfile);
hfile2 = newnatresvmp(bbox.BBox, bc.Resolution);
bc2 = bvqxfile_getcont(hfile2.L);
bc2.Map.Type = 1;
bc2.Map.DF2 = 0;
bc2.Map(2:2 * numel(vsel)) = bc2.Map(1);

% get timecourse as single
tc = single(bc.VTCData(:, :, :, :));

% generate mean and std maps
for cc = 1:numel(vsel)
    mc = 2 * cc - 1;
    bc2.Map(mc).Name = sprintf('mean(vol selection %d)', cc);
    mapv = squeeze(mean(tc(vsel{cc}, :, :, :)));
    mapv(isnan(mapv) | isinf(mapv)) = 0;
    mval = mean(mapv(mapv ~= 0));
    mstd = std(mapv(mapv ~= 0));
    bc2.Map(mc).LowerThreshold = max(eps, mval - 0.5 * mstd);
    bc2.Map(mc).UpperThreshold = mval + 1.5 * mstd;
    bc2.Map(mc).DF1 = numel(vsel{cc}) - 1;
    bc2.Map(mc).VMPData = mapv;
    mc = mc + 1;
    bc2.Map(mc).Name = sprintf('std(vol selection %d)', cc);
    mapv(mapv == 0) = Inf;
    mapv = (squeeze(std(tc(vsel{cc}, :, :, :))) ./ mapv) .* ...
        (mapv > bc2.Map(mc - 1).LowerThreshold);
    mapv(isnan(mapv) | isinf(mapv)) = 0;
    mval = mean(mapv(mapv ~= 0));
    mstd = std(mapv(mapv ~= 0));
    bc2.Map(mc).LowerThreshold = max(eps, 100 * (mval - 2 * mstd));
    bc2.Map(mc).UpperThreshold = 100 * (mval + 3 * mstd);
    bc2.Map(mc).DF1 = numel(vsel{cc}) - 1;
    bc2.Map(mc).VMPData = 100 * mapv;
end
bc2.NrOfMaps = mc;

% put into object
bvqxfile_setcont(hfile2.L, bc2);