bvcoordconv.m

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

function c = bvcoordconv(c, ct, bbox)
% bvcoordconv  - converting coordinates for BV
%
% FORMAT:       c = bvcoordconv(c, ct, bbox)
%
% Input fields:
%
%       c           Cx3 coordinates (or Cx1 / 1xC 1-based indices)
%       ct          conversion type, one of
%                   'bvc2tal'  convert BV internal volume coords to TAL
%                   'bvi2tal'  convert BV internal system coords to TAL
%                   'bvs2tal'  convert BV external system coords to TAL
%                   'bvx2tal'  convert BV internal volume indices to TAL
%                   'tal2bvc'  convert TAL to BV internal volume coords
%                   'tal2bvi'  convert TAL to BV internal system coords
%                   'tal2bvs'  convert TAL to BV external system coords
%                   'tal2bvx'  convert TAL to BV internal volume indices
%       bbox        1x1 bounding box struct with fields
%        .BBox      hardcoded bounding box
%        .DimXYZ    size of volume to sample (needed for indexing)
%        .FCube     framing cube (needed to decide about center)
%        .ResXYZ    resolution of voxels (in BV internal notation)
%
% Output fields:
%
%       c           converted coordinates
%
% For further explanations see 'dbtype bvcoordconv' !

% Version:  v0.7b
% Build:    7090512
% Date:     Sep-05 2006, 12:26 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% Note: the documentation in this file has been influenced by
%       Graham Wideman's very good notes about coordinate systems; see
%
%       http://www.grahamwideman.com/gw/brain/orientation/orientterms.htm
%  and  http://www.grahamwideman.com/gw/brain/analyze/formatdoc.htm
%
% Abstract:
%
% Coordinate systems (or their differences) commonly lead to
% misinterpretations of spatial references and/or wrongly visualized
% data (sliced images). This can only be solved by knowing about the
% "good" interpretation of data as found in the several data formats
%
% BVQX files (different for FMR/STC, VMR/VTC/VDW/GLM/VMP/...)
% SPM/Analyze/NIftI
% DICOM (scanner space)
%
% There are a few different ways to describe the orientation and
% layout of coordinate systems. In this file, the "from"-naming scheme
% is used, which gives the three letters (e.g. RPI for right-to-left
% in posterior-to-anterior in inferior-to-superior order) to describe
% how data is ordered in a chunk of data which then also determines the
% order and orientation of the coordinate system axes. So, each letter
% denotes the "side" with the **LOWER** coordinate indices (so that for
% RPI -X is right, +X is left, etc.). Other guides use the reverse
% nomenclature, so the same system could be described as LAS! Note
% that in this file, the "from"-scheme is used henceforth.
%
% Also, a yet unresolved question is whether coordinates are to be
% interpreted as zero-space or unit-space covering elements (so, is
% a voxel at (0;0;0) covering the area *around* (0;0;0) or from (0;0;0)
% to the next unit (centered around (0.5;0.5;0.5) for a 1mm dataset)?
%
% In other words, if a dataset with the center at (0;0;0) is X-mirrored,
% what happens to the sample at (0;0;0)? If it remains "in place", the
% space for this voxel is from (-0.5;-0.5;-0.5) to (0.5;0.5;0.5).
% Otherwise, the data would, after the transformation, be in (-1;0;0)!
% BrainVoyager assumes voxels to be sampled from the naming coordinate,
% so that a value at (0;0;0) would indeed be resampled to (-1;0;0) after
% mirroring. This is in accordance to the fact that BrainVoyager always
% suggests using even-numbered resolutions (unlike SPM, for which odd
% numbers in volume sizes are the standard case!)
%
% As a simplification, this last question has not yet been dealt with,
% but will hopefully be covered soon.
%
% Explanation of working mechanism:
%
% The TAL coordinate system (LPI) works (all in 1mm resolution)
%
% - 1st dimension:       left (-127) to right    (128) 
% - 2nd dimension:  posterior (-127) to anterior (128)
% - 3rd dimension:   inferior (-127) to superior (128)
%
% with a center coordinate at (0;0;0)
%
% The BV system coordinates (RAS, external, as presented to the user!) have
% the same general axes logic, but with reversed directions. To
% slighlty complicate matters, the interpretation (center) depends on
% the framing cube of the data (in older versions fixed to 256):
%
% - 1st dimension:    right (0) to left      (FCube - 1)
% - 2nd dimension: anterior (0) to posterior (FCube - 1)
% - 3rd dimension: superior (0) to inferior  (FCube - 1)
%
% with the center at (FCube / 2;FCube / 2;FCube / 2), usually at
% (128;128;128) which is still fixed for non VMR/VMP files!
%
% Please note: any rotations to coordinates are transformed around
% the center coordinate - 0.5 (so the virtual point at the center of
% the data!), assuming that the voxel center is at (0.5;0.5;0.5) !!
%
% The *second* BV system coordinates (ASR, internal ordering only,
% rarely shown to users!) on top permute the axes order as [2, 3, 1]:
%
% - 1st dimension:  anterior (0) to posterior (FCube - 1)
% - 2nd dimension:  superior (0) to inferior  (FCube - 1)
% - 3rd dimension:     right (0) to left      (FCube - 1)
%
% This leads to a sagittal slicing, and all three standard sliced
% images in BVQX have a growing indexing from (screen-space!)
% left-to-right in top-to-bottom (SAG: X,Y; COR: Z,Y; TRA: Z,X)
%
% The BV internal volume coordinates (0-based) then represent the
% indexing information into a 3-D array, so that from the BV internal
% system coordinates the Offset must be subtracted and the result
% devided by the resolution (finally adding 1 for ML's indexing).
%
% The BV volume indexing, eventually, represents this information
% in a simple number using sub2ind(DataDimXYZ, c).
%
% Note: this notation differs from the BrainVoyager QX coordinates in
%       one (more or less crucial) point:
%
% the toolbox always represents system coordinates in a 1x1x1 ISO-voxel
% resolution, while BrainVoyager does not do so for hi-res VMRs (all other
% filetypes assume an underlying 1x1x1 resolution of the hosting object!)
%
% NB: BrainVoyager also allows Neurological convention of the data,
%     leading to a LAS/ASL external/internal coordinate system. However,
%     this is only denoted in the VMR file (boolean flag telling the
%     program how to interpret all data on top), so all other files are
%     unaware of the actual orientation! Visualizing any VTC/GLM/VMP/SRF
%     file hence depends on loading the (a) correctly oriented VMR.
%     Due to this fact, it is *strongly suggested* not to use the
%     Neurological convention to allow simple data exchange between
%     users. (status as of BVQX 1.9.9)
%
%     It is still not fully clear for datasets in neurological convention
%     whether coregistration *MUST* be performed on a radiological VMR
%     to get to good results later...
%
% Other formats:
%
% SPM/Analyze:
%
% The Analyze format has (like BVQX file formats) "grown" over time.
% One of the changes is the "flipping" that has to be guessed for older
% files that do not contain this information in the Analyze file header.
% SPM assumes (for files where this information is not available) that
% data comes in a flipped orientation (see spm_defaults.m).
%
% The general logic for data storage follows a (from) RPI scheme, while the
% display is in "neurological convention" (so, in the GUI, left IS left!)
% hence, data that is stored in the radiological convention must be
% accessed in a "flipped" manner. This is denoted in the Analyze header
% either
%
% - by using a negatively oriented direction vector in the Quaternion
%   and/or AffineTransX matrix (header with nii/ni1/n+1 magic token)
% - by setting the first (otherwise unused) PixSpacing to a negative
%   value (intermediate solution)
% - by using the corresponding DataHist.Orientation value (3) for old
%   Analyze format (header without nii/ni1/n+1 magic token)
%
% If neither of these conditions is met, the toolbox's default is yet
% to assume flipped data ! To disable this behavior, use
%
% BVQXfile(0, 'config', 'hdr', 'assumeflipped', false);
%
% Also, note that, other than BrainVoyager, Analyze uses the logic that
% each voxel coordinate is in the center of the space measured, so that
% the data from the voxel around coordinate (0;0;0) occupies space from
% (-0.5;-0.5;-0.5) to (0.5;0.5;0.5), and a rotation is performed around
% (0;0;0) leaving the central voxel at it's position!

% argument check
if nargin < 3 || ...
   ~isa(c, 'double') || ...
    ndims(c) > 2 || ...
   ~ischar(ct) || ...
   ~any(strcmpi(ct(:)', ...
        {'bvc2tal', 'bvi2tal', 'bvs2tal', 'bvx2tal', ...
         'tal2bvc', 'tal2bvi', 'tal2bvs', 'tal2bvx'})) || ...
   (~strcmpi(ct(1:3), 'bvx') && ...
    ~any(size(c) == 3))
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
ct = lower(ct(:)');
cf = ct(1:3);
ci = ct(5:7);
if size(c, 1) == 3 && ...
    size(c, 2) ~= 3
    c = c';
end
bbx = bbox.BBox(1, 1:3);
xyz = bbox.DimXYZ;
fcb = bbox.FCube;
res = bbox.ResXYZ;

% which processing direction
switch (ct)
    case {'bvc2tal', 'bvi2tal', 'bvs2tal', 'bvx2tal'}
        total = true;
    case {'tal2bvc', 'tal2bvi', 'tal2bvs', 'tal2bvx'}
        total = false;
end

% in TAL direction
if total
    
    % starting at indexing
    if strcmp(cf, 'bvx')
        
        % convert to internal volume
        [cc{1:3}] = ind2sub(xyz, c(:));
        c = [cc{1}, cc{2}, cc{3}];
        cf = 'bvc';
    end
    
    % coordinates in internal volume coords but further conversion
    if strcmp(cf, 'bvc') && ...
       ~strcmp(ci, 'bvc')
   
        % convert to internal system coords
        c = [res(1) * (c(:, 1) - 1) + bbx(1), ...
             res(2) * (c(:, 2) - 1) + bbx(2), ...
             res(3) * (c(:, 3) - 1) + bbx(3)];
        cf = 'bvi';
    end
    
    % coordinates in internal system coords but further conversion
    if strcmp(cf, 'bvi') && ...
       ~strcmp(ci, 'bvi')
   
        % convert to external system coords 
        c = [c(:, 3), c(:, 1), c(:, 2)];
        cf = 'bvs';
    end
    
    % coordinates in external system coords but further conversion
    if strcmp(cf, 'bvs') && ...
       ~strcmp(ci, 'bvs')
   
        % convert to TAL space
        c = bbox.FCube / 2 - c;
    end
    
% coming from TAL direction
else
    
    % starting at TAL
    if strcmp(cf, 'tal')
        
        % convert to external system coords
        c = fcb / 2 - c;
        cf = 'bvs';
    end
    
    % coordinates in external system coords but further conversion
    if strcmp(cf, 'bvs') && ...
       ~strcmp(ci, 'bvs')
   
        % convert to internal system coords
        c = [c(:, 2), c(:, 3), c(:, 1)];
        cf = 'bvi';
    end
    
    % coordinates in internal system coords but further conversion
    if strcmp(cf, 'bvi') && ...
       ~strcmp(ci, 'bvi')
   
        % convert to internal volume coords
        c = [(c(:, 1) - bbx(1)) ./ res(1) + 1, ...
             (c(:, 2) - bbx(2)) ./ res(2) + 1, ...
             (c(:, 3) - bbx(3)) ./ res(3) + 1];
        cf = 'bvc';
    end
    
    % coordinates in internal volume coords but further conversion
    if strcmp(cf, 'bvc') && ...
       ~strcmp(ci, 'bvc')
        
        % convert to indexes
        c = sub2ind(xyz, round(c(:, 1)), round(c(:, 2)), round(c(:, 3)));
    end
end