function ds_conv=cosmo_vol_grid_convert(ds,varargin)
% convert between volumetric (fmri) and grid-based (meeg source) dataset
%
% ds_conv=cosmo_vol_grid_convert(ds[,direction])
%
% Inputs:
%   ds              dataset struct, either in fmri or meeg source form
%   direction       (optional) conversion direction, either
%                   - 'tovol' : convert to a volumetric (fmri) dataset
%                   - 'togrid': convert to a source (meeg) dataset
%                   If this option is omitted, then 'tovol' is selected if
%                   ds is an meeg source dataset, and 'togrid' otherwise
%
% Output:
%   ds_conv         dataset in volumetric or source format
%
% Example:
%     ds=cosmo_synthetic_dataset('size','normal');
%     % feature attributes are i, j, k
%     cosmo_disp(ds.fa,'edgeitems',2);
%     %|| .i
%     %||   [ 1         2  ...  2         3 ]@1x30
%     %|| .j
%     %||   [ 1         1  ...  2         2 ]@1x30
%     %|| .k
%     %||   [ 1         1  ...  5         5 ]@1x30
%     % fmri dataset has .a.vol
%     cosmo_disp(ds.a,'edgeitems',2);
%     %|| .fdim
%     %||   .labels
%     %||     { 'i'
%     %||       'j'
%     %||       'k' }
%     %||   .values
%     %||     { [ 1         2         3 ]
%     %||       [ 1         2 ]
%     %||       [ 1         2         3         4         5 ] }
%     %|| .vol
%     %||   .mat
%     %||     [ 2         0         0        -3
%     %||       0         2         0        -3
%     %||       0         0         2        -3
%     %||       0         0         0         1 ]
%     %||   .dim
%     %||     [ 3         2         5 ]
%     %||   .xform
%     %||     'scanner_anat'
%     %
%     % convert to get grid representation
%     ds_grid=cosmo_vol_grid_convert(ds,'togrid');
%     % feature attribute is pos
%     cosmo_disp(ds_grid.fa,'edgeitems',2);
%     %|| .pos
%     %||   [ 1        11  ...  20        30 ]@1x30
%     % meeg source dataset has no .a.vol
%     cosmo_disp(ds_grid.a,'edgeitems',2);
%     %|| .fdim
%     %||   .labels
%     %||     { 'pos' }
%     %||   .values
%     %||     { [ -1        -1  ...  3         3
%     %||         -1        -1  ...  1         1
%     %||         -1         1  ...  5         7 ]@3x30 }
%     %
%     % convert it back to a format identical to the input
%     ds_vol=cosmo_vol_grid_convert(ds_grid,'tovol');
%     %
%     isequal(ds.fa,ds_vol.fa)
%     %|| true
%     isequal(ds.a.vol.mat,ds_vol.a.vol.mat)
%     %|| true
%     isequal(ds.a.fdim,ds_vol.a.fdim)
%     %|| true
%
% Notes:
%   - when ds is an meeg source dataset, it is required that the positions
%     can be placed in a regular grid
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    check_input(ds);
    direction=get_direction(ds,varargin{:});

    switch direction
        case 'tovol'
            ds_conv=convert_grid2vol(ds);

        case 'togrid'
            ds_conv=convert_vol2grid(ds);

        case 'none'
            ds_conv=ds;

        otherwise
            assert(false,'this should not happen');
    end

    cosmo_check_dataset(ds_conv);


function ds_conv=convert_vol2grid(ds)
    % ensure required fields are there
    pos=cosmo_vol_coordinates(ds);
    nfeatures=size(ds.samples,2);


    [idxs,unq_pos_tr]=cosmo_index_unique(pos');
    pos_values=unq_pos_tr';
    pos_fa=zeros(1,nfeatures);
    for k=1:numel(idxs)
        pos_fa(idxs{k})=k;
    end
    assert(all(pos_fa~=0));
    % remove i, j, k fields
    ds_conv=cosmo_dim_remove(ds,{'i','j','k'});

    % remove volumetric field
    ds_conv.a=rmfield(ds_conv.a,'vol');

    % insert position field
    ds_conv=cosmo_dim_insert(ds_conv,2,0,{'pos'},...
                    {pos_values},{pos_fa},'matrix_labels',{'pos'});


function ds_conv=convert_grid2vol(ds)
    pos_label='pos';
    ijk_labels={'i';'j';'k'};

    [dim, index]=cosmo_dim_find(ds,pos_label,true);
    assert(dim==2);
    values=ds.a.fdim.values{index};

    ndim=size(values,1);
    if ndim~=3
        error('''pos'' in .a.fdim.values{%d} must be 3xN',index);
    end

    tolerance=1e-6;
    min_deltas=zeros(ndim,1);
    min_values=zeros(ndim,1);
    dim_sizes=zeros(ndim,1);
    dim_values=cell(ndim,1);

    fa=struct();
    for dim=1:ndim
        % attempt to convert grid to linear space
        vs=values(dim,:);
        unq_vs=unique(vs);

        min_value=min(unq_vs);

        deltas=unique(diff(unq_vs));
        deltas=deltas(deltas>tolerance);

        if numel(deltas)==0
            min_delta=1;
            max_ratio=0;
        else
            min_delta=min(deltas);

            ratios=(unq_vs-min_value)/min_delta;

            max_ratio=round(max(ratios));

            if numel(ratios)~=(1+max_ratio) || ...
                            max(abs(ratios-(0:max_ratio)))>tolerance
                error(['''pos'' in .a.fdim.values{%d}(%d,:) do not '...
                        'follow grid-like structure'],index,dim);
            end
        end
        % store data
        min_deltas(dim)=min_delta;
        min_values(dim)=min_value;
        dim_sizes(dim)=max_ratio+1;

        % convert position to integer in range 1:dim_sizes(dim)
        vs_idxs=1+round((vs-min_value)/min_delta);
        dim_values{dim}=1:dim_sizes(dim);

        % set feature attribute
        ijk_label=ijk_labels{dim};
        fa.(ijk_label)=(vs_idxs(ds.fa.(pos_label)));
    end

    mat=diag([min_deltas;1]);
    mat(1:3,4)=min_values-min_deltas;

    ds_conv=cosmo_dim_remove(ds,'pos');
    ds_conv=cosmo_dim_insert(ds_conv,2,0,ijk_labels,dim_values,fa);
    ds_conv.a.vol.mat=mat;
    ds_conv.a.vol.dim=dim_sizes(:)';

    if cosmo_isfield(ds_conv,'a.meeg.dim')
        ds_conv.a.meeg=rmfield(ds_conv.a.meeg,'dim');
    end



function check_input(ds)
    cosmo_check_dataset(ds);


function tf=ds_has_vol(ds)
    tf=all(cosmo_isfield(ds,{'fa.i','fa.j','fa.k','a.vol.mat'}));

function tf=ds_has_pos(ds)
    tf=cosmo_isfield(ds,'fa.pos');


function direction=get_direction(ds,varargin)
    narg=numel(varargin);

    if narg>1
        error('More than two input arguments are not supported');
    end

    has_direction=narg==1;

    if ds_has_vol(ds)
        guessed_direction='togrid';
    elseif ds_has_pos(ds)
        guessed_direction='tovol';
    else
        error(['Input dataset must either have .fa.pos '...
                    '(for MEEG source), or '...
                    'both fa.{i,j,k} and .a.vol.mat (for fMRI volume)']);
    end

    if has_direction
        direction=varargin{1};

        if ~cosmo_match({direction},{'togrid','tovol'})
            error('direction must be ''togrid'' or ''tovol''');
        end

        if ~strcmp(direction,guessed_direction)
            direction='none';
        end
    else
        direction=guessed_direction;
    end

