function cfr = hdr_CoordinateFrame(hfile)
% HDR::CoordinateFrame  - get coordinate frame from Analyze image
%
% FORMAT:       cframe = hdr.CoordinateFrame
%
% No input fields
%
% Output fields:
%
%       cframe      struct with fields
%        .DimX/Y/Z  spatial dimension of image
%        .DimT      number of volumes for multi-volume image, default 1
%        .Dimensions combined list
%        .ResX/Y/Z  spatial resolution (in TAL/MNI/real world mm)
%        .Resolution combined list
%        .Slice1Center  center coordinate of first slice
%        .SliceNCenter  center coordinate of last slice
%        .SystemOrigin  voxel around (0;0;0)
%        .RowDir    equivalent to AffineTransX(1:3)
%        .ColDir    equivalent to AffineTransY(1:3)
%        .SlcDir    equivalent to AffineTransZ(1:3)
%        .IsRadiological  flag whether coordinate system is left-handed
%        .Trf       4x4 quaternion matrix, so that
%                   cframe.Trf * voxel := mm
%          - and -  inv(cframe.Trf) * mm := voxel

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

% global config needed
global bvqxfile_config;

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'hdr')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
sc = bvqxfile_getscont(hfile.L);
bc = sc.C;
cfr = struct;

% check for additional MAT file
afile = sc.F;
exten = afile(end-2:end);
if any(exten == upper(exten))
    mfile = [afile(1:end-4) '.MAT'];
else
    mfile = [afile(1:end-4) '.mat'];
end
mfilec = [];
mfilel = false;
if exist(mfile, 'file') == 2
    try
        mfilem = load(mfile);
        mfilec = mfilem.mat;
        mfilel = true;
    catch
        % do nothing
    end
end

% get dimensions
cfr.DimX = bc.ImgDim.Dim(2);
cfr.DimY = bc.ImgDim.Dim(3);
if bc.ImgDim.Dim(1) > 2
    cfr.DimZ = bc.ImgDim.Dim(4);
    resz = bc.ImgDim.PixSpacing(4);
else
    cfr.DimZ = 1;
    resz = 1;
end
if bc.ImgDim.Dim(1) > 3
    cfr.DimT = bc.ImgDim.Dim(5);
else
    cfr.DimT = 1;
end
cfr.Dimensions = [cfr.DimX, cfr.DimY, cfr.DimZ, cfr.DimT];

% get resolution
cfr.ResX = bc.ImgDim.PixSpacing(2);
cfr.ResY = bc.ImgDim.PixSpacing(3);
cfr.ResZ = resz;

% set full resolution element
cfr.Resolution = [cfr.ResX, cfr.ResY, cfr.ResZ];

% compute 1-based transformation matrix (as used in SPM)
if ~isempty(mfilec)
    % not needed if mat file was loaded 
elseif ~all( ...
    [bc.DataHist.NIftI1.AffineTransX(end), ...
     bc.DataHist.NIftI1.AffineTransY(end), ...
     bc.DataHist.NIftI1.AffineTransZ(end)] == 0)
    mfilec = [ ...
        bc.DataHist.NIftI1.AffineTransX; ...
        bc.DataHist.NIftI1.AffineTransY; ...
        bc.DataHist.NIftI1.AffineTransZ; ...
        0, 0, 0, 1];
    mfilec(1:3, 4) = mfilec(1:3, 4) - mfilec(1:3, 1:3) * [1; 1; 1];
    mfilel = true;
elseif ~all(bc.DataHist.OriginSPM(1:3) == 0)
    mfilec = [ ...
        cfr.ResX,     0   ,     0   , -cfr.ResX * bc.DataHist.OriginSPM(1); ...
            0   , cfr.ResY,     0   , -cfr.ResY * bc.DataHist.OriginSPM(2); ...
            0   ,     0   , cfr.ResZ, -cfr.ResZ * bc.DataHist.OriginSPM(3); ...
            0   ,     0   ,     0   ,     1];
elseif bc.NIIFileType == 0
    switch (double(bc.DataHist.Orientation))
        case {0}
            mfilec = [ ...
                cfr.ResX,     0   ,     0   , -cfr.ResX * cfr.DimX / 2; ...
                    0   , cfr.ResY,     0   , -cfr.ResY * cfr.DimY / 2; ...
                    0   ,     0   , cfr.ResZ, -cfr.ResZ * cfr.DimZ / 2; ...
                    0   ,     0   ,     0   ,      1];
        case {1}
            mfilec = [ ...
                cfr.ResX,     0   ,     0   , -cfr.ResX * cfr.DimX / 2; ...
                    0   ,     0   , cfr.ResZ, -cfr.ResZ * cfr.DimZ / 2; ...
                    0   , cfr.ResY,     0   , -cfr.ResY * cfr.DimY / 2; ...
                    0   ,     0   ,     0   ,      1];
        case {2}
            mfilec = [ ...
                    0   , cfr.ResY,     0   , -cfr.ResY * cfr.DimY / 2; ...
                    0   ,     0   , cfr.ResZ, -cfr.ResY * cfr.DimZ / 2; ...
                cfr.ResX,     0   ,     0   , -cfr.ResX * cfr.DimX / 2; ...
                    0   ,     0   ,     0   ,      1];
        case {3}
            mfilec = [ ...
                cfr.ResX,     0   ,     0   , -cfr.ResX * cfr.DimX / 2; ...
                    0   ,-cfr.ResY,     0   ,  cfr.ResY * cfr.DimY / 2; ...
                    0   ,     0   , cfr.ResZ, -cfr.ResZ * cfr.DimZ / 2; ...
                    0   ,     0   ,     0   ,      1];
        case {4}
            mfilec = [ ...
                cfr.ResX,     0   ,     0   , -cfr.ResX * cfr.DimX / 2; ...
                    0   ,     0   ,-cfr.ResZ,  cfr.ResZ * cfr.DimZ / 2; ...
                    0   , cfr.ResY,     0   , -cfr.ResY * cfr.DimY / 2; ...
                    0   ,     0   ,     0   ,      1];
        case {5}
            mfilec = [ ...
                    0   , cfr.ResY,     0   , -cfr.ResY * cfr.DimY / 2; ...
                    0   ,     0   ,-cfr.ResZ,  cfr.ResZ * cfr.DimZ / 2; ...
                cfr.ResX,     0   ,     0   , -cfr.ResX * cfr.DimX / 2; ...
                    0   ,     0   ,     0   ,      1];
        otherwise
            warning( ...
                'BVQXfile:InvalidObject', ...
                'Unknown DataHist.Orientation value: %d.', ...
                double(bc.DataHist.Orientation) ...
            );
    end
end

% last resort
if isempty(mfilec)
    mfilec = [ ...
        cfr.ResX,     0   ,     0   , -cfr.DimX / 2; ...
            0   , cfr.ResY,     0   , -cfr.DimY / 2; ...
            0   ,     0   , cfr.ResZ, -cfr.DimZ / 2; ...
            0   ,     0   ,     0   ,      1];
end

% convention guessing needed ?
if ~mfilel && ...
   (bc.ImgDim.PixSpacing(1) < 0 || ...
    (bc.ImgDim.PixSpacing(1) == 0 && ...
     bc.NIIFileType < 1 && ...
     bc.DataHist.Orientation == 0 && ...
     bvqxfile_config.type.hdr.assumeflipped))
    if cfr.ResX < 0
        cfr.ResX = -cfr.ResX;
    end
    
    % perform x-flip to get to TAL
    mfilec(1, :) = -mfilec(1, :);
end

% compute slice 1 and N center coordinates (Trf is one-based!)
sl1c = mfilec * [0.5 + cfr.DimX / 2; 0.5 + cfr.DimY / 2; 0.5; 1];
slnc = mfilec * [0.5 + cfr.DimX / 2; 0.5 + cfr.DimY / 2; 0.5 + cfr.DimZ; 1];

% get direction components
dcp = -mfilec(1:3, 1:3);
dcp = dcp ./ repmat(sqrt(sum(dcp .* dcp, 2)), [1, 3]);
dcp(4, :) = -cross(dcp(1, :), dcp(2, :));

% make settings
cfr.Slice1Center = sl1c(1:3)';
cfr.SliceNCenter = slnc(1:3)';
cfr.SystemOrigin = -mfilec(1:3, 4)';
cfr.RowDir = dcp(1, :);
cfr.ColDir = dcp(2, :);
cfr.SlcDir = dcp(3, :);
cfr.IsRadiological = (sum(dcp(3, :) .* dcp(4, :)) > 0);
cfr.Trf = mfilec;