vmr_talairach.m

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

function hfile2 = vmr_Talairach(hfile, imeth, tal, acpc, inverse)
% VMR::Talairach  - (un-) Talairachize VMR
%
% FORMAT:       talvmr = vmr.Talairach(imeth, tal [, acpc, inverse]);
%
% Input fields:
%
%       imeth   interpolation method ('linear', 'cubic', 'spline')
%       tal     TAL object
%       acpc    optional ACPC TRF file to go from/back to native space
%       inverse perform inverse (un-Tal, un-ACPC) operation
%
% Output fields:
%
%       tvmr    (un-) Talairachized VMR
%
% Note: only works on VMRs in 1mm or 0.5mm ISOvoxel resolution

% TODO: piecewise interpolation !!

% Version:  v0.7b
% Build:    7090216
% Date:     Sep-02 2007, 4:54 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, 'vmr') || ...
   ~ischar(imeth) || ...
    isempty(imeth) || ...
   ~any(strcmpi(imeth(:)', {'cubic', 'linear', 'spline'})) || ...
    numel(tal) ~= 1 || ...
   ~isBVQXfile(tal, 'tal')
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if (bc.FileVersion > 1 && ...
    (any([bc.VoxResX, bc.VoxResY, bc.VoxResZ] ~= 1 & ...
         [bc.VoxResX, bc.VoxResY, bc.VoxResZ] ~= 0.5) || ...
     ~all([bc.VoxResY, bc.VoxResZ] == bc.VoxResX)))
    error( ...
        'BVQXfile:InvaldObject', ...
        'Method only valid for 1mm or 0.5mm ISOvoxel VMRs.' ...
    );
end
doovinv = false;
if nargin == 4 && ...
    islogical(acpc) && ...
   ~isempty(acpc)
    ovinv = acpc(1);
    doovinv = true;
end
if nargin < 4 || ...
    numel(acpc) ~= 1 || ...
   ~isBVQXfile(acpc, 'trf')
    acpc = [];
else
    acpcbc = bvqxfile_getcont(acpc.L);
    if ~strcmpi(acpcbc.DataFormat, 'Matrix') || ...
        acpcbc.TransformationType ~= 2 || ...
        acpcbc.CoordinateSystem ~= 0
        acpc = acpcbc.TFMatrix;
    else
        acpc = [];
    end
end
if nargin < 5 || ...
   ~islogical(inverse) || ...
    isempty(inverse)
    inverse = false;
else
    inverse = inverse(1);
end
if doovinv
    inverse = ovinv;
end
if inverse && ...
   ~isempty(acpc)
    try
        acpc = inv(acpc)';
    catch
        error( ...
            'BVQXfile:BadArgument', ...
            'Couldn''t invert ACPC transformation matrix.' ...
        );
    end
elseif ~isempty(acpc)
    acpc = acpc';
end
tvmrd = bc.VMRData(1);
tvmrd(1) = 0;

% create new VMR
hfile2 = BVQXfile('new:vmr');
cfile2 = bvqxfile_getcont(hfile2.L);
cfile2.VMR8bit = bc.VMR8bit;
if ~inverse
    cfile2.VoxResInTalairach = true;
end
cfile2.VoxResVerified = bc.VoxResVerified;

% grids depend on
res = bc.VoxResX;
if res == 1
    cs = 256;
else
    cs = 512;
    cfile2.VoxResX = res;
    cfile2.VoxResY = res;
    cfile2.VoxResZ = res;
end
[xgrd, ygrd] = ndgrid(0:res:255.5, 0:res:255.5);
tc = [xgrd(:), ygrd(:), zeros(numel(xgrd), 1)];
zval = 0:res:255.5;
zc = [zeros(numel(zval), 2), zval(:)];

% the simple way only works for forward or single transformation !
if ~inverse || ...
    isempty(acpc)
    tc(:, [3, 1, 2]) = acpc2tal(tc(:, [3, 1, 2]), tal, ~inverse);
    zc(:, [3, 1, 2]) = acpc2tal(zc(:, [3, 1, 2]), tal, ~inverse);
end
if ~isempty(acpc)
    tc = tc - 127.5;
    zc = zc - 127.5;
end
zval = zc(:, 3)';
tc(:, 4) = 1;

% determine progress bar capabilities
try
    if ~inverse
        if isempty(acpc)
            step = 'Talairach transforming';
        else
            step = 'ACPC / Talairach transforming';
        end
    else
        if isempty(acpc)
            step = 'Un-Talairach transforming';
        else
            step = 'Un-Talairach / Un-ACPC transforming';
        end
    end
        
    pbar = BVQXprogress;
    BVQXprogress(pbar, 'setposition', [80, 200, 640, 36]);
    BVQXprogress(pbar, 'settitle', sprintf('%s VMR...', step));
    BVQXprogress(pbar, 0, 'Setting up target VMR...', 'visible', 0, 1);
catch
    pbar = [];
end

% create sampling matrix
try
    sdt = single(bc.VMRData(:, :, :));
catch
    bvqxfile_clear(hfile2.L);
    error( ...
        'BVQXfile:OutOfMemory', ...
        'Out of memory.' ...
    );
end
ssz = size(sdt);
sdg = {};
stepwise = false;
if numel(bc.VMRData) < 4e7
    try
        sdg = cell(1, 3);
        sdg{1} = repmat(single(reshape(1:ssz(1), [ssz(1), 1])), [1, ssz(2:3)]);
        sdg{2} = repmat(single(reshape(1:ssz(2), [1, ssz(2)])), [ssz(1), 1, ssz(3)]);
        sdg{3} = repmat(single(reshape(1:ssz(3), [1, 1, ssz(3)])), ssz(1:2));
    catch
        sdg = {};
        stepwise = true;
    end
end
if stepwise
    smx = size(sdt);
end

% create output matrix
cfile2.VMRData = tvmrd;
tvmrd(1:(cs * cs), 1:cs) = tvmrd;

% get sampling offset and size
ssz = size(sdt) + 1;
try
    sof = [bc.OffsetX, bc.OffsetY, bc.OffsetZ];
catch
    sof = [0, 0, 0];
end

% iterate over z slices
if ~isempty(pbar)
    BVQXprogress(pbar, 0, 'Interpolating...', 'visible', 0, numel(zval));
end
for zc = 1:numel(zval)
    
    % apply transformation
    tc(:, 3) = zval(zc);
    if ~isempty(acpc)
        sc = (tc * acpc) + 127.5;
    else
        sc = tc;
    end
    
    % apply UnTal after ACPC
    if inverse && ...
       ~isempty(acpc)
        sc(:, [3, 1, 2]) = acpc2tal(sc(:, [3, 1, 2]), tal);
    end
    
    % good scaling / origin
    if res == 0.5
        sc = (2 .* sc + 1.5);
    else
        sc = sc + 1;
    end
    
    % offset
    if any(sof)
        sc = [sc(:, 1) - sof(1), sc(:, 2) - sof(2), sc(:, 3) - sof(3)];
    end
    
    % good samples
    p = (sc(:, 1) >= 0 & sc(:, 1) <= ssz(1) & ...
         sc(:, 2) >= 0 & sc(:, 2) <= ssz(2) & ...
         sc(:, 3) >= 0 & sc(:, 3) <= ssz(3));

    % sample VMR
    if any(p)
        
        % sample whole volume
        if ~stepwise
            tvmrd(p, zc) = round( ...
                interpn(sdg{:}, sdt, sc(p, 1), sc(p, 2), sc(p, 3), imeth));
            
        % or sample only required part
        else
            switch (imeth)
                case {'cubic'}
                    minc = max(floor(min(sc(p, :))) - 3, 1) - 1;
                    maxc = min(ceil(max(sc(p, :))) + 3, smx);
                case {'linear', 'nearest'}
                    minc = max(floor(min(sc(p, :))) - 1, 1) - 1;
                    maxc = min(ceil(max(sc(p, :))) + 1, smx);
                case {'spline'}
                    minc = max(floor(min(sc(p, :))) - 5, 1) - 1;
                    maxc = min(ceil(max(sc(p, :))) + 5, smx);
            end
            
        end
    end
    
    % progress bar
    if ~isempty(pbar)
        BVQXprogress(pbar, zc);
    end
end

% make sure to limit 8bit VMR
if strcmpi(class(tvmrd), 'uint8')
    tvmrd = min(tvmrd, uint8(225));
end

% set some more VMR fields
tvmrd = reshape(tvmrd, [cs, cs, cs]);
cfile2.VMRData = tvmrd;
cfile2.VoxResInTalairach = ~inverse;

% make sure that V16 files are FileVersion 1
if ~cfile2.VMR8bit || ...
    strcmpi(class(cfile2.VMRData), 'uint16')
    hfile2 = vmr_Update(hfile2, 'FileVersion', [], 1);
    cfile2.VMR8bit = false;
end
bvqxfile_setcont(hfile2.L, cfile2);

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