map_createvmp.m

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

function vmp = map_CreateVMP(hfile, fmr, afs, vmpfile, meth, res, bbox, o)
% MAP::CreateVMP  - create VMP from FMR based MAP object
%
% FORMAT:       vmp = map.CreateVMP(fmr, afs [, vmpfile, meth, res, bb, o])
%
% Input fields:
%
%       fmr         FMR object matching MAP
%       afs         cell array with alignment files, {ia, fa, acpc, tal}
%       vmpfile     target VMP filename, if not given: unsaved
%       meth        interpolation, 'nearest', {'linear'}, 'cubic', 'spline'
%       res         VMP resolution, default: 1
%       bb          bounding box, default: small TAL box
%       o           optional struct with settings
%        .df        degrees of freedom, 1x1 or 1x2 DF setting
%
% Output fields:
%
%       vmp         resulting VMP object

% Version:  v0.7a
% Build:    7081222
% Date:     Aug-12 2007, 10:18 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, 'map') || ...
    numel(fmr) ~= 1 || ...
   ~isBVQXfile(fmr, 'fmr') || ...
   ~iscell(afs) || ...
    isempty(afs) || ...
    numel(afs{1}) ~= 1 || ...
   ~isBVQXfile(afs{1}, 'trf')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
maps = bvqxfile_getscont(hfile.L);
mapc = maps.C;
if ~any([0, 1, 2, 3] == mapc.Type)
    error( ...
        'BVQXfile:InvalidObject', ...
        'Unsupported map type: %d.', ...
         mapc.Type ...
    );
end
if nargin < 4 || ...
   ~ischar(vmpfile) || ...
    isempty(vmpfile)
    vmpfile = '';
else
    vmpfile = vmpfile(:)';
end
fmrs = bvqxfile_getscont(fmr.L);
fmrc = fmrs.C;
if numel(mapc.Map) ~= fmrc.NrOfSlices || ...
    numel(size(mapc.Map(1).Data)) ~= 2 || ...
    any(size(mapc.Map(1).Data) ~= [fmrc.ResolutionX, fmrc.ResolutionY])
    error( ...
        'BVQXfile:InvalidObject', ...
        'FMR and MAP dimension(s) mismatch.' ...
    );
end
if nargin < 5 || ...
   ~ischar(meth) || ...
    isempty(meth) || ...
   ~any(strcmpi(meth(:)', {'nearest', 'linear', 'cubic', 'spline'}))
    meth = 'linear';
else
    meth = lower(meth(:)');
end
if nargin < 6 || ...
   ~isa(res, 'double') || ...
    numel(res) ~= 1 || ...
   ~any([1, 2, 3] == res)
    res = 1;
end
if nargin < 7 || ...
   ~isa(bbox, 'double') || ...
    numel(size(bbox)) ~= 2 || ...
    any(size(bbox) ~= [2, 3]) || ...
    any(isinf(bbox(:)) | isnan(bbox(:)) | bbox(:) < 1 | bbox(:) > 254 | ...
        bbox(:) ~= fix(bbox(:))) || ...
    any((bbox(2, :) - res) <= bbox(1, :))
    bbox = [59, 57, 52; 197, 231, 172];
end
if nargin < 8 || ...
   ~isstruct(o) || ...
    numel(o) ~= 1
    o = struct;
end
if ~isfield(o, 'df') || ...
   ~isa(o.df, 'double') || ...
    numel(o.df) > 2 || ...
    any(isnan(o.df) | isinf(o.df) | o.df < 1 | o.df ~= fix(o.df))
    o.df = [];
end
if res ~= 1
    bbox(2, :) = bbox(1, :) + res .* round(diff(bbox) ./ res);
end

% generate target coordinate list and data field
if res == 1
    c = {bbox(1, 1):bbox(2, 1), ...
         bbox(1, 2):bbox(2, 2), ...
         bbox(1, 3):bbox(2, 3)};
else
    c = {bbox(1, 1):res:(bbox(2, 1) - 1), ...
         bbox(1, 2):res:(bbox(2, 2) - 1), ...
         bbox(1, 3):res:(bbox(2, 3) - 1)};
end

% check alignment files
hasia = 0;
clearia = cell(1, 1);
for ac = numel(afs):-1:1
    if numel(afs{ac}) == 1 && ...
        isBVQXfile(afs{ac}, 'trf')
        if hasia > 0
            afs(hasia) = [];
        end
        trfc = bvqxfile_getcont(afs{ac}.L);
        if trfc.TransformationType == 1 && ...
            trfc.AlignmentStep == 1
            hasia = ac;
        end
    elseif numel(afs{ac}) ~= 1 || ...
       ~isBVQXfile(afs{ac}, 'tal')
        error( ...
            'BVQXfile:BadArgument', ...
            'Invalid alignment file object list.' ...
        );
    end
end
if hasia == 0
    warning( ...
        'BVQXfile:BadArgument', ...
        'Initial alignment missing; using identity matrix instead.' ...
    );
    afs = [{BVQXfile('new:trf')}, afs(:)'];
    clearia = afs{1};
end

% generate sample output with target coordinates
nslc = fmrc.NrOfSlices;
smpo = zeros(fmrc.ResolutionX, fmrc.ResolutionY, nslc);

% fill with map content
for mc = 1:nslc
    smpo(:, :, mc) = mapc.Map(mc).Data(:, :);
end
bfv = sum(smpo(:) ~= 0);

% sample at VMP coordinates
try
    smpd = permute(single(samplefmrspace(smpo, c, fmr, afs, meth)), [2, 3, 1]);
catch
    clearbvqxobjects(clearia);
    error( ...
        'BVQXfile:InternalError', ...
        'Error sampling MAP data.' ...
    );
end

% create and fill VMP
vmp = BVQXfile('new:vmp');
vmpc = bvqxfile_getcont(vmp.L);
ofv = vmpc.FileVersion;
if res == 1
    vmpc.NativeResolutionFile = 0;
    vmpc.FileVersion = 4;
else
    vmpc.NativeResolutionFile = 1;
    vmpc.DocumentType = 1;
    vmpc.FileVersion = 5;
end
vmpc.Resolution = res;
bvqxfile_setcont(vmp.L, vmpc);
vmp_Update(vmp, 'FileVersion', struct('type', '.', 'subs', 'FileVersion'), ofv);
vmpc = bvqxfile_getcont(vmp.L);
vmpc.NrOfMaps = 1;
vmpc.XStart = bbox(1, 2);
vmpc.XEnd = bbox(2, 2);
vmpc.YStart = bbox(1, 3);
vmpc.YEnd = bbox(2, 3);
vmpc.ZStart = bbox(1, 1);
vmpc.ZEnd = bbox(2, 1);
if res ~= 1
    vmpc.LinkedPRT = fmrc.ProtocolFile;
end
vmpc.Map.Type = mapc.Type + 1;
[mapf{1:2}] = fileparts(maps.F);
vmpc.Map.Name = sprintf('VMP from MAP: %s', mapf{2});
vmpc.Map.LowerThreshold = mapc.LowerThreshold;
vmpc.Map.UpperThreshold = mapc.UpperThreshold;
if ~isempty(mapc.NrOfLags)
    vmpc.Map.NrOfLags = mapc.NrOfLags;
    vmpc.Map.MinLag = 0;
    vmpc.Map.MaxLag = mapc.NrOfLags - 1;
    vmpc.Map.CCOverLay = 1;
end
vmpc.Map.ClusterSize = ...
    round((fmrc.InplaneResolutionX * fmrc.InplaneResolutionY * ...
    (fmrc.SliceThickness + fmrc.GapThickness) * mapc.ClusterSize) / (res ^ 3));
vmpc.Map.EnableClusterCheck = 1;

% guessing here, true for r !
vmpc.Map.DF2 = 0;
if isempty(o.df)
    vmpc.Map.DF1 = fmrc.NrOfVolumes - 8;
    if vmpc.Map.Type == 4
        vmpc.Map.DF2 = 1;
    end

% df are provided
else
    vmpc.Map.DF1 = o.df(1);
    if numel(o.df) > 1
        vmpc.Map.DF2 = o.df(2);
    end
end

vmpc.Map.BonferroniValue = bfv;
if res ~= 1
    vmpc.Map.FDRThresholds = [0, 1e5, 1e5];
    if any([1, 4] == vmpc.Map.Type)
        mapv = smpo(~isnan(smpo(:)) & ~isinf(smpo(:)) & (smpo(:) ~= 0));
        fdrt = [0.1, 0.05, 0.04, 0.03, 0.02, 0.01, 0.005, 0.001];
        switch (vmpc.Map.Type)
            case {1}  % t-score
                vmpc.Map.FDRThresholds = [fdrt(:), ...
                    custom_tinv(1 - fdr_thresholds( ...
                    1 - custom_tcdf(abs(double(mapv)), ...
                    vmpc.Map.DF1), fdrt(:) / 2, true), vmpc.Map.DF1)];
            case {4}  % F-score
                vmpc.Map.FDRThresholds = [fdrt(:), ...
                    custom_finv(1 - fdr_thresholds( ...
                    1 - custom_fcdf(abs(double(mapv)), ...
                    vmpc.Map.DF1, vmpc.Map.DF2), fdrt(:), true), ...
                    vmpc.Map.DF1, vmpc.Map.DF2)];
        end
    end
end
vmpc.Map.NrOfFDRThresholds = size(vmpc.Map.FDRThresholds, 1);
vmpc.Map.VMPData = smpd;
bvqxfile_setcont(vmp.L, vmpc);

% try saving
try
    if ~isempty(vmpfile)
        aft_SaveAs(vmp, vmpfile);
    end
catch
    warning( ...
        'BVQXfile:InternalError', ...
        'Error saving VMP file: ''%s''.', ...
        lasterr ...
    );
end