function [orient,cano_rot]=cosmo_fmri_orientation(ds)
% get orientation of a dataset
%
% [orient,cano_rot]=cosmo_fmri_orientation(ds)
%
% Inputs:
%    ds             fmri dataset struct;
%
% Output:
%    orient         three letter string indicating the orientation of this
%                   dataset
%    cano_rot       canonical rotation matrix (relative to LPI)
%
% Example:
%     ds=cosmo_synthetic_dataset();
%     [orient,cano_rot]=cosmo_fmri_orientation(ds);
%     disp(orient)
%     %|| LPI
%     disp(cano_rot)
%     %|| 1 0 0 0
%     %|| 0 1 0 0
%     %|| 0 0 1 0
%     %|| 0 0 0 1
%
% Notes:
%   - there are 3!*3^2 valid orientations, these are:
%         'SAR'  'SAL'  'SPR'  'SPL'  'IAR'  'IAL'  'IPR'  'IPL'
%         'SRA'  'SLA'  'SRP'  'SLP'  'IRA'  'ILA'  'IRP'  'ILP'
%         'ASR'  'ASL'  'PSR'  'PSL'  'AIR'  'AIL'  'PIR'  'PIL'
%         'ARS'  'ALS'  'PRS'  'PLS'  'ARI'  'ALI'  'PRI'  'PLI'
%         'RAS'  'LAS'  'RPS'  'LPS'  'RAI'  'LAI'  'RPI'  'LPI'
%         'RSA'  'LSA'  'RSP'  'LSP'  'RIA'  'LIA'  'RIP'  'LIP'
%     For example, 'LPI' (used in Talairach/MNI) means that
%       * the first dimension goes from left to right
%       * the second dimension goes from posterior to anterior
%       * the third dimension goes from inferior to superior
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    mx=get_affine_matrix(ds);

    rot_orig=mx(1:3,1:3);

    % normalize the rows
    rot_row_norms=sqrt(sum(rot_orig.^2,2));
    rot=bsxfun(@rdivide,rot_orig,rot_row_norms);

    cano_rot=get_canonical_matrix(rot);

    orient=get_orientation(cano_rot);

function cano_rot=get_canonical_matrix(rot)
    % find the three major axes of the rotation matrix
    visited=false(3);
    cano_rot=zeros(4);
    cano_rot(4,4)=1;

    for dim=1:3
        max_v=0;
        for row=1:3
            for col=1:3
                v=abs(rot(row,col));
                if ~visited(row,col) && v>max_v
                    max_row=row;
                    max_col=col;
                    max_v=v;
                end
            end
        end
        if max_v==0
            error('Illegal matrix');
        end
        visited(max_row,:)=true;
        visited(:,max_col)=true;

        cano_rot(max_row,max_col)=sign(rot(max_row,max_col));
    end

    assert(all(visited(:)));
    assert(all(sum(cano_rot~=0,1)==1));
    assert(all(sum(cano_rot~=0,2)==1));

function orient=get_orientation(cano_rot)
    labs=['LR';'PA';'IS'];

    orient='   ';
    for k=1:3
        col=find(cano_rot(1:3,k)~=0);
        assert(numel(col)==1);
        v=cano_rot(col,k);
        orient(k)=labs(col,1+(v<0));
    end


function mx=get_affine_matrix(ds)
    if cosmo_isfield(ds,'a.vol.mat')
        mx=ds.a.vol.mat;
    elseif isnumeric(ds) && size(ds,1)>=3 && size(ds,2)>=3
        mx=ds;
    else
        error('cannot find affine transformation matrix');
    end
