vtc_resample.m

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

function nhfile = vtc_Resample(hfile, ifunc, ts, xs, ys, zs)
% VTC::Resample  - resample a VTC
%
% FORMAT:       newvtc = vtc.Resample(ifunc, ts, [xs, ys, zs])
%
% Input fields:
%
%       ifunc       interpolation function,
%                   'linear', 'cubic', 'nearest', 'spline'
%       ts          temporal sampling points (in MS!)
%       xs, ys, zs  spatial sampling points (BV coords, all given or none)

% Version:  v0.7b
% Build:    7083012
% Date:     Aug-30 2007, 12:58 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, 'vtc') || ...
   ~ischar(ifunc) || ...
   ~any(strcmpi(ifunc(:)', {'linear', 'cubic', 'nearest', 'spline'})) || ...
   ~isa(ts, 'double') || ...
    isempty(ts) || ...
    length(ts) ~= numel(ts) || ...
    (numel(ts) > 2 && any(abs(diff(ts, 2) > 1))) || ...
    any(isinf(ts(:)) | isnan(ts(:))) || ...
    (nargin > 5 && ( ...
     ~isa(xs, 'double') || ...
     ~isa(ys, 'double') || ...
     ~isa(zs, 'double') || ...
      length(xs) ~= numel(xs) || ...
      length(ys) ~= numel(ys) || ...
      length(zs) ~= numel(zs) || ...
      numel(xs) < 3 || ...
      numel(ys) < 3 || ...
      numel(zs) < 3 || ...
      any(xs < 0 | xs > 255) || ...
      any(ys < 0 | ys > 255) || ...
      any(zs < 0 | zs > 255) || ...
      any(diff(xs, 2)) || ...
      any(diff(ys, 2)) || ...
      any(diff(zs, 2)) || ...
      xs(1) ~= fix(xs(1)) || ...
      ys(1) ~= fix(ys(1)) || ...
      zs(1) ~= fix(zs(1)) || ...
      abs(xs(1) - xs(2)) ~= fix(abs(xs(1) - xs(2))) || ...
      abs(xs(1) - xs(2)) ~= fix(abs(ys(1) - ys(2))) || ...
      abs(xs(1) - xs(2)) ~= fix(abs(zs(1) - zs(2)))))
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
ifunc = lower(ifunc(:)');
bc = bvqxfile_getcont(hfile.L);
nv = size(bc.VTCData, 1);
if any(ts < 0 | ts > (bc.TR * (nv + 1)))
    warning( ...
        'BVQXfile:BadArgument', ...
        'Sampling points out of bounds.' ...
    );
end

% get current array size
oOffsetX = bc.XStart;
oOffsetY = bc.YStart;
oOffsetZ = bc.ZStart;
oDimX = size(bc.VTCData, 2);
oDimY = size(bc.VTCData, 3);
oDimZ = size(bc.VTCData, 4);
oRes = bc.Resolution;

% decide on spatial sampling points
if nargin < 6
    nOffsetX = oOffsetX;
    nOffsetY = oOffsetY;
    nOffsetZ = oOffsetZ;
    nRes = oRes;
    xs = single(1:oDimX);
    ys = single(1:oDimY);
    zs = single(1:oDimZ);
else
    nOffsetX = min(xs(1), xs(end));
    nOffsetY = min(ys(1), ys(end));
    nOffsetZ = min(zs(1), zs(end));
    nRes = abs(xs(1) - xs(2));
    xs = single(1 + (xs(:)' - oOffsetX) / oRes);
    ys = single(1 + (ys(:)' - oOffsetY) / oRes);
    zs = single(1 + (zs(:)' - oOffsetZ) / oRes);
end

% get newTR
if numel(ts) > 1
    nTR = round(mean(diff(ts(:))));
else
    nTR = bc.TR;
end
TRfac = nTR / bc.TR;
ts = 1 + ts(:) / bc.TR;

% try progress bar
try
    pbar = BVQXprogress;
    BVQXprogress(pbar, 'setposition', [80, 200, 640, 36]);
    BVQXprogress(pbar, 'settitle', 'Resampling VTC...');
    BVQXprogress(pbar, 0, 'Scaling original VTC...', 'visible', 0, 10);
    pvx = 9 / (numel(ys) * numel(zs));
    pvc = 1;
catch
    pbar = [];
end

% create data arrays (old in single AND new in uint16 must fit into mem!)
nVTC = uint16(0);
try
    oVTC = single(bc.VTCData(:, :, :, :));
    mVTC = mean(oVTC);
    mVTC(mVTC == 0) = 1;
    oVTC = 10000 * (oVTC ./ repmat(mVTC, [nv, 1]));
    nVTC(numel(ts), numel(xs), numel(ys), numel(zs)) = nVTC(1);
catch
    error( ...
        'BVQXfile:OutOfMemory', ...
        'Out of memory (old and new VTC must fit in memory).' ...
    );
end

% create new VTC object
if ~isempty(pbar)
    BVQXprogress(pbar, 0.5, 'Creating new object...');
end
nhfile = BVQXfile('new:vtc');
ncfile = bvqxfile_getcont(nhfile.L);
ncfile.NameOfSourceFMR = bc.NameOfSourceFMR;
ncfile.NameOfLinkedPRT = bc.NameOfLinkedPRT;
ncfile.NrOfVolumes = numel(ts);
ncfile.Resolution = nRes;
ncfile.XStart = nOffsetX;
ncfile.XEnd = nOffsetX + nRes * numel(xs);
ncfile.YStart = nOffsetY;
ncfile.YEnd = nOffsetY + nRes * numel(ys);
ncfile.ZStart = nOffsetZ;
ncfile.ZEnd = nOffsetZ + nRes * numel(zs);
ncfile.TR = nTR;
ncfile.SegmentSize = fix(bc.SegmentSize * TRfac);
ncfile.SegmentOffset = fix(bc.SegmentOffset * TRfac);

% build source matrices
if ~isempty(pbar)
    BVQXprogress(pbar, 0.55, 'Building interpolation coordinate matrices...');
end
[tf, xf, yf, zf] = ...
    ndgrid(single(1:nv), single(1:oDimX), single(1:oDimY), single(1:oDimZ));

% build sampling matrices
ts = single(ts);
nnv = numel(ts);
nvx = numel(xs);
nvt = nnv * nvx;
xs = reshape(repmat(xs(:)', [nnv, 1]), [nvt, 1]);
ts = reshape(repmat(ts(:), [nvx, 1]), [nvt, 1]);

% loop over slices dims
if ~isempty(pbar)
    BVQXprogress(pbar, 1, sprintf('Resampling with %s interpolation...', ifunc));
end
for z = 1:numel(zs)
    for y = 1:numel(ys)
        nVTC(:, :, y, z) = uint16(round(reshape(interpn(tf, xf, yf, zf, oVTC, ...
            ts, xs, repmat(ys(y), [nvt, 1]), repmat(zs(z), [nvt, 1]), ifunc), ...
            [nnv, nvx])));
        if ~isempty(pbar)
            BVQXprogress(pbar, 1 + pvc * pvx);
            pvc = pvc + 1;
        end
    end
end

% set final field
ncfile.VTCData = nVTC;
clear oVTC;

% but content into array
bvqxfile_setcont(nhfile.L, ncfile);

% remove bar
if ~isempty(pbar)
    closebar(pbar);
end