test vol coordinates

function test_suite = test_vol_coordinates
% tests for cosmo_vol_coordinates
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;


function test_vol_coordinates_()
    ds=cosmo_synthetic_dataset();

    [ns,nf]=size(ds.samples);

    % sample some features
    rp=randperm(nf);
    nsel=round(nf/2);
    rp=rp(1:nsel);
    ds_sel=cosmo_slice(ds,rp,2);

    % voxel to world space
    xyz=cosmo_vol_coordinates(ds_sel);
    assertEqual(size(ds_sel.samples),[ns nsel]);

    xyz_=cosmo_vol_coordinates(ds_sel,1:nsel);
    assertEqual(xyz, xyz_);

    selsel=randperm(nsel);
    xyz__=cosmo_vol_coordinates(ds_sel,selsel);
    assertEqual(xyz(:,selsel), xyz__);

    % check coordinate transform
    ijk1=[ds_sel.fa.i; ds_sel.fa.j; ds_sel.fa.k; ones(1,nsel)];
    xyz1=ds.a.vol.mat*ijk1;
    assertElementsAlmostEqual(xyz1(1:3,:),xyz_);

    % double the dataset
    ds_sel2=cosmo_stack({ds_sel,ds_sel},2);
    xyz2=cosmo_vol_coordinates(ds_sel2);
    assertEqual(xyz2, [xyz xyz]);

    % world to voxel space
    fa_indices=cosmo_vol_coordinates(ds_sel,xyz(:,selsel));
    assertVectorsAlmostEqual(fa_indices,selsel);

    % duplicate features cannot be matched back
    assertExceptionThrown(@()cosmo_vol_coordinates(ds_sel2,...
                                            xyz(:,selsel)),'');


function test_vol_coordinates_afni()
% test coordinates using AFNI's 3dmerge
    if cosmo_skip_test_if_no_external('afni_bin')
        return;
    end
    % make an fMRI dataset
    dim_sizes=[60 50 30];
    dim_values=cellfun(@(x) 1:x,num2cell(dim_sizes),'UniformOutput',false);
    dim_labels={'i','j','k'};
    vox_sizes=[3 3 3];
    origin=-dim_sizes.*vox_sizes*.4; % slightly off-center
    nsamples=2;

    data4d=zeros([nsamples,dim_sizes]);

    ds=cosmo_flatten(data4d,dim_labels,dim_values);
    ds.sa.labels=cellfun(@(x) sprintf('sample%d',x),...
                                    num2cell(1:nsamples)',...
                                    'UniformOutput',false);
    ds.sa.stats=repmat({'Ftest(10,1)'},nsamples,1);
    ds.a.vol.mat=[diag(vox_sizes),origin';[0 0 0 1]];
    ds.a.vol.dim=dim_sizes;


    ds.sa.labels=ds.sa.labels(:,1);

    center_distance=20;
    offsets=[3 7 11]; % some offsets in x, y, z directions
    xyz=bsxfun(@plus, cosmo_cartprod(repmat({[-1,1]*center_distance},1,3)),...
                        offsets);
    ncoord=size(xyz,1);

    fa_indices=cosmo_vol_coordinates(ds,xyz');
    xyz=cosmo_vol_coordinates(ds,fa_indices); % recompute
    ds.samples(:,fa_indices)=repmat(1:ncoord,nsamples,1); % 1 to 8

    % make temporary directory
    tmp_dir=cosmo_make_temp_filename();
    tmp_dir_cleaner=onCleanup(@()remove_dir_helper(tmp_dir));

    mkdir(tmp_dir);

    % generate AFNI files
    ext_postfix={'.nii',{''};
                 '+orig',{'.BRIK','.BRIK.gz','.HEAD'}};
    for j=1:size(ext_postfix,1)
        ext=ext_postfix{j,1};
        postfixes=ext_postfix{j,2};

        get_fn=@(x) fullfile(tmp_dir,[x ext]);
        get_to_delete=@(x,exts) cellfun(@(y) sprintf('%s%s',...
                            get_fn(x),y),postfixes,'UniformOutput',false);

        base_fn=get_fn('base');
        cosmo_map2fmri(ds,base_fn);

        orients={'lpi','rai','asr'};

        for k=1:numel(orients)
            orient=orients{k};
            fn=get_fn(orient);
            cmd=sprintf('3dresample -overwrite -orient %s -prefix %s -input %s',...
                            orient, fn, base_fn);
            r=unix(cmd);
            assert(r==0);
            cmd=sprintf(['3dclust -orient LPI -1clip .5 0 0 %s''[0]'' '...
                        '2>/dev/null | grep --invert-match ''#'''], fn);
            [r,v]=unix(cmd);
            assert(r==0);

            to_delete=get_to_delete(orient,postfixes);
            cmd=sprintf('rm -f %s', cosmo_strjoin(to_delete,' '));
            unused=unix(cmd);

            s=sscanf(v,'%f');
            table=reshape(s,16,[]);

            % sort by voxel value
            [unused,i]=sort(table(11,:));
            table=table(:,i);

            line1=ones(1,ncoord);
            line0=0*line1;
            line_inc=1:ncoord;
            expected_table=[line1; xyz; kron(xyz,[1 1]'); ...
                                line_inc; line0; line_inc; xyz];
            assertEqual(table, expected_table)
        end
    end

function test_test_vol_coordinates_exceptions()
    aet=@(varargin)assertExceptionThrown(@()...
                    cosmo_vol_coordinates(varargin{:}),'');

    ds=cosmo_synthetic_dataset();

    % not 1xP or 3xP
    aet(ds, ones(2));
    aet(ds, [1;2]);

    % not in fdim
    ds_bad=ds;
    ds_bad=cosmo_dim_transpose(ds_bad,'i',1);
    aet(ds_bad,1);

    % size mismatch
    ds_bad=ds;
    ds_bad.a.vol.dim=ds_bad.a.vol.dim+1;
    aet(ds_bad,1);

    % no matrix
    ds_bad=ds;
    ds_bad.a.vol=rmfield(ds_bad.a.vol,'mat');
    aet(ds_bad,1);

    % not 4x4 matrix
    ds_bad=ds;
    ds_bad.a.vol.mat=eye(5);
    aet(ds_bad,1);



function remove_dir_helper(tmp_dir)
    if cosmo_wtf('is_octave')
        rmdir_state=confirm_recursive_rmdir();
        state_resetter=onCleanup(@()confirm_recursive_rmdir(rmdir_state));
        confirm_recursive_rmdir(false,'local');
    end
    rmdir(tmp_dir,'s');