vmp_importspmmaps.m

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

function hfile = vmp_ImportSPMMaps(hfile, maps, opts)
% VMP::ImportSPMMaps  - import SPM based maps from Analzye/NIftI files
%
% FORMAT:       vmp.ImportSPMMaps(maps [, opts])
%
% Input fields:
%
%       maps        filenames of spmT_* / spmF_* maps
%       opts        1x1 struct with optional fields
%        .interp    method ('nearest', {'linear'}, 'cubic', 'spline')
%        .maptype   1x1 or Mx1 map type, 'F', 'r', 't' (default from name)
%        .thresh    1x2 or Mx2 lower and upper thresholds (default: p=0.05)
%        .trf       optional cell array with transformation files
%
% No output fields.

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

% check arguments
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmp') || ...
    isempty(maps) || ...
   (~ischar(maps) && ...
    ~iscell(maps))
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
res = bc.Resolution;
if isempty(bc.Map)
    error( ...
        'BVQXfile:BadArgument', ...
        'The VMP must not be empty (for default settings).' ...
    );
end
if nargin < 3 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ischar(maps)
    if size(maps) == 1
        maps = {maps(:)'};
    else
        maps = cellstr(maps);
    end
end
if ~isfield(opts, 'interp') || ...
   ~ischar(opts.interp) || ...
   ~any(strcmpi(opts.interp(:)', {'cubic', 'linear', 'nearest', 'spline'}))
    opts.interp = 'linear';
else
    opts.interp = lower(opts.interp(:)');
end
if ~isfield(opts, 'maptype') || ...
   ~ischar(opts.maptype) || ...
   (numel(opts.maptype) ~= 1 && ...
    numel(opts.maptype) ~= numel(maps))
    opts.maptype = 'a';
else
    opts.maptype = lower(opts.maptype(:)');
end
if numel(opts.maptype) == 1
    opts.maptype(2:numel(maps)) = opts.maptype(1);
end
if ~isfield(opts, 'thresh') || ...
   (~isa(opts.thresh, 'double') && ...
    ~isa(opts.thresh, 'single')) || ...
    size(opts.thresh, 2) ~= 2 || ...
   (size(opts.thresh, 1) ~= 1 && ...
    size(opts.thresh, 1) ~= numel(maps)) || ...
    any(isinf(opts.thresh(:)) | isnan(opts.thresh(:)))
    opts.thresh = [0.05, 0.001];
else
    opts.thresh = abs(opts.thresh(:, 1:2, 1));
    opts.thresh = opts.thresh(:, 1:2);
end
if size(opts.thresh, 1) == 1
    opts.thresh = repmat(opts.thresh, [numel(maps), 1]);
end
if ~isfield(opts, 'trf') || ...
   ~iscell(opts.trf)
    opts.trf = {};
else
    opts.trf = opts.trf(:)';
end
hastal = [];
for oc = numel(opts.trf):-1:1
    if numel(opts.trf{oc}) ~= 1 || ...
       (~isBVQXfile(opts.trf{oc}, 'tal') && ...
        ~isBVQXfile(opts.trf{oc}, 'trf'))
        opts.trf(oc) = [];
        continue;
    end
    if isBVQXfile(opts.trf{oc}, 'tal')
        hastal = opts.trf{oc};
        opts.trf(oc) = [];
    else
        tbc = bvqxfile_getcont(opts.trf{oc}.L);
        if tbc.TransformationType ~= 2
            opts.trf(oc) = [];
        end
    end
end
for mc = numel(maps):-1:1
    if ~ischar(maps{mc}) || ...
        isempty(maps{mc}) || ...
        exist(maps{mc}(:)', 'file') ~= 2
        maps(mc) = [];
        opts.maptype(mc) = [];
        opts.thresh(mc, :) = [];
        continue;
    end
    maps{mc} = maps{mc}(:)';
    if numel(maps{mc}) > 4 && ...
        strcmpi(maps{mc}(end-3:end), '.img')
        if exist([maps{mc}(1:end-3) 'hdr'], 'file') == 2
            maps{mc} = [maps{mc}(1:end-3) 'hdr'];
        elseif exist([maps{mc}(1:end-3) 'HDR'], 'file') == 2
            maps{mc} = [maps{mc}(1:end-3) 'HDR'];
        else
            opts.maptype(mc) = [];
            opts.thresh(mc, :) = [];
            maps(mc) = [];
        end
    end
end
if isempty(maps)
    error( ...
        'BVQXfile:BadArgument', ...
        'Map file(s) not found.' ...
    );
end
nobj = numel(maps);
mobj = cell(1, nobj);
for mc = 1:nobj
    try
        mobj{mc} = BVQXfile(maps{mc});
        if numel(mobj{mc}) ~= 1 || ...
           ~isBVQXfile(mobj{mc}, 'hdr')
            error('BAD_HDR_FILE');
        end
        tbc = bvqxfile_getcont(mobj{mc}.L);
        if isempty(regexpi(tbc.DataHist.Description, ...
                '^spm (beta|contrast) - (\d+)\:\s*(.*)$')) && ...
            isempty(regexpi(tbc.DataHist.Description, ...
                '^spm\{([a-z]+)_\[([0-9][0-9\.]+)(\,[0-9][0-9\.]+)?\]\}\s*\-?\s*(.*)$'))
            error('BAD_HDR_FILE');
        end
    catch
        clearbvqxobjects(mobj);
        error( ...
            'BVQXfile:BadArgument', ...
            'Invalid map file (not an SPM result map): %s', ...
            maps{mc} ...
        );
    end
end

% build coordinate indices for sampling
if res > 1 || ...
    bc.NativeResolutionFile
    xs = 0.5 * res + (bc.XStart:res:bc.XEnd-1);
    ys = 0.5 * res + (bc.YStart:res:bc.YEnd-1);
    zs = 0.5 * res + (bc.ZStart:res:bc.ZEnd-1);
else
    xs = 0.5 * res + (bc.XStart:res:bc.XEnd);
    ys = 0.5 * res + (bc.YStart:res:bc.YEnd);
    zs = 0.5 * res + (bc.ZStart:res:bc.ZEnd);
end
xyzn = [numel(xs), numel(ys), numel(zs)];
[xs, ys, zs] = ndgrid(single(xs), single(ys), single(zs));
xs = xs(:);
ys = ys(:);
zs = zs(:);

% correct axes order in system coordinates
c = [zs, xs, ys];

% apply un-TAL ?
if ~isempty(hastal)
    c = acpc2tal(c, hastal, false);
end

% apply transformations in backwards order
for oc = numel(opts.trf):-1:1
    tbc = bvqxfile_getcont(opts.trf{oc}.L);
    c = applybvtrf(c, tbc.TFMatrix, false);
end

% to TAL system
c = 128 - c;

% add maps to structure once !
bc.Map = bc.Map(:);
onmaps = numel(bc.Map);
bc.Map(onmaps+1:onmaps+nobj) = bc.Map(onmaps);

% FDR thresholds
fdrv = [0.1, 0.05, 0.04, 0.03, 0.02, 0.01, 0.005, 0.001]';

% iterate over map objects
for oc = 1:nobj
    
    % get target position
    tobj = onmaps + oc;
    
    % get content of map
    sobj = bvqxfile_getscont(mobj{oc}.L);
    cobj = sobj.C;
    if istransio(cobj.VoxelData)
        cobj.VoxelData = resolve(cobj.VoxelData);
    end
	cobj.VoxelData = single(cobj.VoxelData);
	cobj.VoxelData(isnan(cobj.VoxelData)) = single(0);
    bvqxfile_setcont(mobj{oc}.L, cobj);
    bfv = sum(cobj.VoxelData(:) ~= 0);
    
    % parse description
    [sfn{1:2}] = fileparts(sobj.F);
    if ~isempty(sfn{1})
        [pfn{1:2}] = fileparts(sfn{1});
        sfn{2} = [pfn{2} '/' sfn{2}];
    end
    desc = cobj.DataHist.Description;
    [iti{1:3}] = regexpi(desc, '^spm (beta|contrast) - (\d+)\:\s*(.*)$');
    [dfi{1:3}] = regexpi(desc, ...
        '^spm\{([a-z]+)_\[([0-9][0-9\.]+)(\,[0-9][0-9\.]+)?\]\}\s*\-?\s*(.*)$');
    
    % maptype
    if ~any('bcfrt' == opts.maptype(oc))
        if ~isempty(iti{3})
            opts.maptype(oc) = lower(desc(iti{3}{1}(1)));
            mapname = sprintf('%s - %s %s: %s', sfn{2}, ...
                desc(iti{3}{1}(1, 1):iti{3}{1}(1, 2)), ...
                desc(iti{3}{1}(2, 1):iti{3}{1}(2, 2)), ...
                desc(iti{3}{1}(3, 1):iti{3}{1}(3, 2)));
        else
            opts.maptype(oc) = lower(desc(dfi{3}{1}(1)));
            mapname = sprintf('%s: %s', sfn{2}, ...
                desc(dfi{3}{1}(4, 1):dfi{3}{1}(4, 2)));
        end
    end
    if ~any('bcfrt' == opts.maptype(oc))
        opts.maptype(oc) = 'b';
    end
    
    % fill in map values
    bc.Map(tobj).DF2 = 0;
    switch (opts.maptype(oc))
        case {'b', 'c'}
            bc.Map(tobj).Type = 11;
        case {'f'}
            bc.Map(tobj).Type = 4;
            if dfi{3}{1}(3, 2) >= dfi{3}{1}(3, 1)
                bc.Map(tobj).DF2 = ...
                    str2double(desc(dfi{3}{1}(3, 1):dfi{3}{1}(3, 2)));
            else
                bc.Map(tobj).DF2 = 1;
            end
        case {'r'}
            bc.Map(tobj).Type = 2;
        case {'t'}
            bc.Map(tobj).Type = 1;
    end
    if any('frt' == opts.maptype(oc))
        bc.Map(tobj).DF1 = str2double(desc(dfi{3}{1}(2, 1):dfi{3}{1}(2, 2)));
        bc.Map(tobj).FDRThresholds = [fdrv, ...
            applyfdr(cobj.VoxelData(:), opts.maptype(oc), fdrv, ...
            bc.Map(tobj).DF1, bc.Map(tobj).DF2, true)];
    else
        bc.Map(tobj).DF1 = 1;
        bc.Map(tobj).FDRThresholds = [1, 1e9, 1e9];
    end
    bc.Map(tobj).NrOfFDRThresholds = size(bc.Map(tobj).FDRThresholds, 1);
    bc.Map(tobj).Name = mapname;
    bc.Map(tobj).BonferroniValue = bfv;
    bc.Map(tobj).UnknownValue = -1;
    bc.Map(tobj).VMPData = reshape(hdr_SampleCoords( ...
        mobj{oc}, c, 1, opts.interp, 'tal'), xyzn);
end

% put content in new object
bc.NrOfMaps = numel(bc.Map);
bvqxfile_setcont(hfile.L, bc);