averagevmrs.m

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

function avmr = averagevmrs(v, cv)
% averagevmrs  - average VMRs
%
% FORMAT:       avmr = averagevmrs(v [, centervmr])
%
% Input fields:
%
%       v           cell array with list of VMR (objects or filenames)
%       centervmr   center around middle coordinate (default false)
%
% Output fields:
%
%       avmr        average VMR with maximum bounding box

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

% argument check
if nargin < 1 || ...
   ~iscell(v) || ...
    numel(v) < 2
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
if nargin < 2 || ...
   ~islogical(cv) || ...
    isempty(cv)
    cv = false;
else
    cv = cv(1);
end
nvmr = numel(v);

% test BVQXprogress
try
    pbar = BVQXprogress;
    BVQXprogress(pbar, 'setposition', [80, 200, 640, 36]);
    BVQXprogress(pbar, 'settitle', 'Averaging VMRs...');
    BVQXprogress(pbar, 0, 'Checking VMRs...', 'visible', 0, nvmr + 3);
catch
    pbar = [];
end

% check each cell
v16auto = BVQXfile(0, 'config', 'vmr', 'v16auto');
BVQXfile(0, 'config', 'vmr', 'v16auto', false);
vmrtios = BVQXfile(0, 'transiosize', 'vmr');
BVQXfile(0, 'transiosize', 'vmr', 65536);
v = v(:);
res = zeros(nvmr, 3);
vm8 = true(1, nvmr);
clearloaded = false(1, nvmr);
for cc = 1:nvmr
    if ischar(v{cc}) && ...
        numel(v{cc}) > 4 && ...
        exist(v{cc}(:)', 'file') == 2 && ...
        any(strcmpi(v{cc}(end-3:end), {'.vmr', '.v16'}))
        try
            [v{cc}, clearloaded(cc)] = BVQXfile(v{cc}(:)');
        catch
            clearbvqxobjects(v(clearloaded));
            if ~isempty(pbar)
                closebar(pbar);
            end
            BVQXfile(0, 'config', 'vmr', 'v16auto', v16auto);
            BVQXfile(0, 'transiosize', 'vmr', vmrtios);
            error( ...
                'BVQXtools:BVQXfileError', ...
                'Error loading VMR from file ''%s'': %s.', ...
                v{cc}(:)', lasterr ...
            );
        end
    elseif numel(v{cc}) ~= 1 || ...
       ~isBVQXfile(v{cc}, 'vmr')
        clearbvqxobjects(v(clearloaded));
        if ~isempty(pbar)
            closebar(pbar);
        end
        BVQXfile(0, 'config', 'vmr', 'v16auto', v16auto);
        BVQXfile(0, 'transiosize', 'vmr', vmrtios);
        error( ...
            'BVQXtools:BadArgument', ...
            'Invalid cell (no VMR filename or object).' ...
        );
    end
    res(cc, :) = v{cc}.BoundingBox.ResXYZ;
    vm8(cc) = v{cc}.VMR8bit;
end
BVQXfile(0, 'config', 'vmr', 'v16auto', v16auto);
BVQXfile(0, 'transiosize', 'vmr', vmrtios);

% check resolutions
if any(any(diff(res) ~= 0)) || ...
    any(vm8 ~= vm8(1))
    clearbvqxobjects(v(clearloaded));
    if ~isempty(pbar)
        closebar(pbar);
    end
    error( ...
        'BVQXtools:BadArgument', ...
        'VMRs must match in resolution and bytesize.' ...
    );
end
vm8 = vm8(1);

% get max bounding box
if ~isempty(pbar)
    BVQXprogress(pbar, 1, 'Getting bounding box...');
end
bbox = minmaxbbox(v);
bbox = bbox(3:4, :);

% create array with matching size
vmrsz = diff(bbox) + 1;
if vm8
    avmrd = uint16([]);
    avmrd(vmrsz(1), vmrsz(2), vmrsz(3)) = uint16(0);
else
    avmrd = uint32([]);
    avmrd(vmrsz(1), vmrsz(2), vmrsz(3)) = uint32(0);
end

% iterate over VMRs
if ~isempty(pbar)
    BVQXprogress(pbar, 2, 'Averaging...');
end
for cc = 1:nvmr
    
    % reframe VMR and get data
    av = v{cc}.CopyObject;
    if istransio(av.VMRData)
        av.VMRData = av.VMRData(:, :, :);
    end
    av.Reframe(bbox);
    
    % add data to average
    if vm8
        avmrd = avmrd + min(uint16(av.VMRData), uint16(225));
    else
        avmrd = avmrd + av.VMRData;
    end
    if ~isempty(pbar)
        BVQXprogress(pbar, 2 + cc);
    end
    av.ClearObject;
end

% copy first VMR and then clear loaded VMRs
avmr = bless(v{1}.CopyObject, 1);
clearbvqxobjects(v(clearloaded));

% build average VMR
if max(avmrd(:)) > (225 * nvmr)
    avmrd = uint16(single(avmrd) ./ single(nvmr));
    vm8 = false;
else
    avmrd = uint8(single(avmrd) ./ single(nvmr));
    vm8 = true;
end
if ~isempty(pbar)
    BVQXprogress(pbar, 3 + nvmr);
end

% copy data into avmr
avmr.VMRData = avmrd;
avmr.VMR8bit = vm8;

% reframe around center ?
if cv
    
    % which resolution
    fc = avmr.FramingCube / 2;
    bb = avmr.BoundingBox.BBox;
    bb(2, :) = bb(2, :) + 1;
    ns = max(abs(fc - bb));
    avmr.Reframe([fc - ns; (fc - 1) + ns]);
end

% clear no longer needed settings
if avmr.FileVersion == 3
    avmr.Slice1CenterX = -127.5;
    avmr.Slice1CenterY = 0;
    avmr.Slice1CenterZ = 0;
    avmr.SliceNCenterX = 127.5;
    avmr.SliceNCenterY = 0;
    avmr.SliceNCenterZ = 0;
    avmr.RowDirX = 0;
    avmr.RowDirY = 1;
    avmr.RowDirZ = 0;
    avmr.ColDirX = 0;
    avmr.ColDirY = 0;
    avmr.ColDirZ = -1;
    avmr.NrOfPastSpatialTransformations = 0;
    avmr.Trf(:) = [];
end

% clear progress bar
if ~isempty(pbar)
    closebar(pbar);
end