vmp_createsmp.m

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

function smp = vmp_CreateSMP(hfile, srf, opts)
% VMP::CreateSMP  - sample SMPs from the maps in a VMP
%
% FORMAT:       smp = vmp.CreateSMP(srf, opts)
%
% Input fields:
%
%       srf         required surface file
%       opts        1x1 struct with optional fields
%        .interp    method ('nearest', {'linear'}, 'cubic', 'spline')
%        .ipfrom    interpolate from P + n * normal vector, default: -3
%        .ipstep    interpolation stepsize, default: 1 (in normal vectors)
%        .ipto      interpolate to P + n * normal vector, default: 1
%        .mapsel    map selection, default: all maps in VMP
%        .method    method to get value ({'mean'}, 'median', 'min')
%        .recalcn   boolean, recalc normals before sampling, default: false
%
% Output fields:
%
%       smp         SMP object with as many maps as selected
%
% Note: results slightly differ from BV's results (sampling method)

% Version:  v0.7b
% Build:    7090215
% Date:     Sep-02 2007, 3: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') || ...
    numel(srf) ~= 1 || ...
   ~isBVQXfile(srf, 'srf')
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
srfs = bvqxfile_getscont(srf.L);
srfc = srfs.C;
if nargin < 3 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
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, 'ipfrom') || ...
   ~isa(opts.ipfrom, 'double') || ...
    numel(opts.ipfrom) ~= 1 || ...
    isnan(opts.ipfrom) || ...
    opts.ipfrom < -10 || ...
    opts.ipfrom > 10
    opts.ipfrom = -3;
end
if ~isfield(opts, 'ipstep') || ...
   ~isa(opts.ipstep, 'double') || ...
    numel(opts.ipstep) ~= 1 || ...
    isnan(opts.ipstep) || ...
    opts.ipstep <= 0 || ...
    opts.ipstep > 10
    opts.ipstep = 1;
end
if ~isfield(opts, 'ipto') || ...
   ~isa(opts.ipto, 'double') || ...
    numel(opts.ipto) ~= 1 || ...
    isnan(opts.ipto) || ...
    opts.ipto < opts.ipfrom || ...
    opts.ipto > 10
    opts.ipto = 1;
end
if ~isfield(opts, 'mapsel') || ...
   ~isa(opts.mapsel, 'double') || ...
    any(isinf(opts.mapsel(:)) | isnan(opts.mapsel(:)) | ...
        opts.mapsel(:) ~= fix(opts.mapsel(:)) | opts.mapsel(:) < 1)
    opts.mapsel = 1:numel(bc.Map);
else
    opts.mapsel = opts.mapsel(:)';
end
if ~isfield(opts, 'method') || ...
   ~ischar(opts.method) || ...
   ~any(strcmpi(opts.method(:)', ...
        {'mean', 'median', 'min'}))
    opts.method = 'mean';
else
    opts.method = lower(opts.interp(:))';
end
if ~isfield(opts, 'recalcn') || ...
   ~islogical(opts.recalcn) || ...
    numel(opts.recalcn) ~= 1
    opts.recalcn = false;
end
if opts.recalcn
    srf_RecalcNormals(srf);
    srfc = bqvxfile_getcont(srf.L);
end
vmpres  = bc.Resolution;
ipsamp  = opts.ipfrom:opts.ipstep:opts.ipto;
ipsnum  = numel(ipsamp);
nummaps = numel(opts.mapsel);

% get coordinates and normals
crd = srfc.VertexCoordinate;
nrm = (1 / vmpres) .* srfc.VertexNormal;
numv = size(crd, 1);

% create output
smp = BVQXfile('new:smp');
smpc = bvqxfile_getcont(smp.L);
smpc.NrOfVertices = numv;
smpc.NrOfMaps = nummaps;
smpc.NameOfOriginalSRF = srfs.F;
smpc.Map.UseValuesAboveThresh = 1;
if nummaps > 1
    smpc.Map(2:nummaps) = smpc.Map;
end

% subtract XStart, YStart, ZStart
crd = 1 + (1 / vmpres) .* [ ...
    crd(:, 1) - bc.XStart, ...
    crd(:, 2) - bc.YStart, ...
    crd(:, 3) - bc.ZStart];

% iterate over maps
for mc = 1:nummaps
    
    % get source map and values
    srcmap = bc.Map(opts.mapsel(mc));
    mapval = double(srcmap.VMPData);
    
    % prepare interpolation array
    ipvals = zeros(numv, ipsnum);
    
    % fill some SMP fields
    smpc.Map(mc).Type = srcmap.Type;
    smpc.Map(mc).LowerThreshold = srcmap.LowerThreshold;
    smpc.Map(mc).UpperThreshold = srcmap.UpperThreshold;
    smpc.Map(mc).UseValuesAboveThresh = srcmap.UseValuesAboveThresh;
    smpc.Map(mc).RGBLowerThreshPos = srcmap.RGBLowerThreshPos;
    smpc.Map(mc).RGBUpperThreshPos = srcmap.RGBUpperThreshPos;
    smpc.Map(mc).RGBLowerThreshNeg = srcmap.RGBLowerThreshNeg;
    smpc.Map(mc).RGBUpperThreshNeg = srcmap.RGBUpperThreshNeg;
    smpc.Map(mc).UseRGBColor = srcmap.UseRGBColor;
    smpc.Map(mc).TransColorFactor = srcmap.TransColorFactor;
    smpc.Map(mc).DF1 = srcmap.DF1;
    smpc.Map(mc).DF2 = srcmap.DF2;
    smpc.Map(mc).BonferroniValue = min(numv, srcmap.BonferroniValue);
    smpc.Map(mc).Name = srcmap.Name;
    
    % only simple approach for non-CC maps
    if srcmap.Type ~= 3
        
        % interpolate according to method
        for ipc = 1:ipsnum
            nf = ipsamp(ipc);
            if ~strcmpi(opts.interp, 'linear')
                ipvals(:, ipc) = interpn(mapval, ...
                    crd(:, 1) + nf * nrm(:, 1), ...
                    crd(:, 2) + nf * nrm(:, 2), ...
                    crd(:, 3) + nf * nrm(:, 3), opts.interp);
            else
                ipvals(:, ipc) = interpn_linnonull(mapval, ...
                    crd(:, 1) + nf * nrm(:, 1), ...
                    crd(:, 2) + nf * nrm(:, 2), ...
                    crd(:, 3) + nf * nrm(:, 3), 0, 0.125);
            end
        end
        
    % CC maps are dealt with more complicated
    else
        
        % fill some
        smpc.Map(mc).NrOfLags = srcmap.NrOfLags;
        smpc.Map(mc).MinLag = srcmap.MinLag;
        smpc.Map(mc).MaxLag = srcmap.MaxLag;
        smpc.Map(mc).CCOverlay = srcmap.CCOverlay;
        
        % extract lags
        maplag = max(0, floor(mapval));
        mapval = max(0, mapval - maplag);
        
        % prepate additional array
        iplags = zeros(numv, ipsnum);
        
        % interpolate according to method
        for ipc = 1:ipsnum
            nf = ipsamp(ipc);
            if ~strcmpi(opts.interp, 'linear')
                ipvals(:, ipc) = max(0, min(1, interpn(mapval, ...
                    crd(:, 1) + nf * nrm(:, 1), ...
                    crd(:, 2) + nf * nrm(:, 2), ...
                    crd(:, 3) + nf * nrm(:, 3), opts.interp)));
                iplags(:, ipc) = max(0, min(srcmap.MaxLag, round(interpn( ...
                    maplag, ...
                    crd(:, 1) + nf * nrm(:, 1), ...
                    crd(:, 2) + nf * nrm(:, 2), ...
                    crd(:, 3) + nf * nrm(:, 3), opts.interp))));
            else
                ipvals(:, ipc) = max(0, min(1, interpn_linnonull(mapval, ...
                    crd(:, 1) + nf * nrm(:, 1), ...
                    crd(:, 2) + nf * nrm(:, 2), ...
                    crd(:, 3) + nf * nrm(:, 3), 0, 0.125)));
                iplags(:, ipc) = max(0, min(srcmap.MaxLag, ...
                    round(interpn_linnonull(maplag, ...
                    crd(:, 1) + nf * nrm(:, 1), ...
                    crd(:, 2) + nf * nrm(:, 2), ...
                    crd(:, 3) + nf * nrm(:, 3), 0, 0.125))));
            end
        end
    end
    
    % what summary method
    switch opts.method
        
        % mean value
        case {'mean'}
            if srcmap.Type ~= 3
                smpc.Map(mc).SMPData = single(mean(ipvals, 2));
            else
                smpc.Map(mc).SMPData = single(mean(ipvals, 2)) + ...
                    1000 * single(round(mean(iplags, 2)));
            end
            
        % median
        case {'median'}
            if srcmap.Type ~= 3
                smpc.Map(mc).SMPData = single(median(ipvals, 2));
             else
                smpc.Map(mc).SMPData = single(median(ipvals, 2)) + ...
                    1000 * single(median(iplags, 2));
            end
            
        % min value
        case {'min'}
            mapsgn = sign(ipvals(:, 1));
            difsgn = (sign(min(ipvals, [], 2)) ~= sign(max(ipvals, [], 2)));
            smpc.Map(mc).SMPData = mapsgn .* min(abs(ipvals), [], 2);
            if srcmap.Type == 3
                smpc.Map(mc).SMPData = smpc.Map(mc).SMPData + ...
                    1000 * min(iplags, [], 2);
            end
            smpc.Map(mc).SMPData = single(smpc.Map(mc).SMPData);
            smpc.Map(mc).SMPData(difsgn) = 0;
    end
end

% put content in new object
bvqxfile_setcont(smp.L, smpc);