test fmri orientation

function test_suite = test_fmri_orientation
    % tests for cosmo_fmri_orientation
    %
    % #   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_orientation()
    % make dataset
    ds = cosmo_synthetic_dataset('size', 'big');
    ds = cosmo_slice(ds, 1);
    ds.a.vol.mat(1:3, 4) = ds.a.vol.mat(1:3, 4) - 20;

    orig_orient = cosmo_fmri_orientation(ds);
    assertEqual(orig_orient, 'LPI');

    % get all orientation and
    orients = get_orients();
    norient = numel(orients);

    fmts = {'nii'};
    nfmt = numel(fmts);

    nperm = nfmt * norient;

    ntest = 5;
    rp = randperm(nperm);
    ncoord = 10;
    for k = 1:ntest
        p = rp(k);
        orient_idx = ceil(p / nfmt);
        fmt_idx = mod(k - 1, nfmt) + 1;

        orient = orients{orient_idx};
        fmt = fmts{fmt_idx};

        ijk = ceil(bsxfun(@times, rand(4, ncoord), [ds.a.vol.dim(:); 1]));
        orient_idx = get_feature_index(ds, ijk);
        xyz = ds.a.vol.mat * ijk;

        % verify equality with re-oriented dataset
        ds_ro = cosmo_fmri_reorient(ds, orient);  % using CoSMoMVPA
        assertEqual(cosmo_fmri_orientation(ds_ro), orient);
        assertEqual(cosmo_fmri_orientation(ds_ro.a.vol.mat), orient);

        ijk_ro = ds_ro.a.vol.mat \ xyz;
        idx_ro = get_feature_index(ds_ro, ijk_ro);
        assertElementsAlmostEqual(ds.samples(:, orient_idx), ...
                                  ds_ro.samples(:, idx_ro), ...
                                  'absolute', 1e-5);

        % verify that going back works fine
        ds2 = cosmo_fmri_reorient(ds_ro, orig_orient);
        assertEqual(ds, ds2);

        % check mapping to and from fmri dataset
        img = cosmo_map2fmri(ds, ['-' fmt]);
        ds3 = cosmo_fmri_dataset(img);
        assertEqual(ds.a, ds3.a);
    end

function test_fmri_orientation_exceptions()
    aet = @(varargin)assertExceptionThrown(@() ...
                                           cosmo_fmri_orientation(varargin{:}), '');
    aet([1 2 3]);
    aet({});
    aet(zeros(4));

function test_fmri_reorient_exceptions()
    aet = @(varargin)assertExceptionThrown(@() ...
                                           cosmo_fmri_reorient(varargin{:}), '');
    ds = cosmo_synthetic_dataset();
    aet(ds, '  RAI');
    aet(ds, '  RAA');

function test_fmri_orientation_with_afni_binary()
    if cosmo_skip_test_if_no_external('afni_bin')
        return
    end

    ds = cosmo_synthetic_dataset('size', 'big');
    fmt = '.nii';
    orients = get_orients();
    i = ceil(rand() * numel(orients));
    orient = orients{i};

    nfeatures = size(ds.samples, 2);
    xyz = cosmo_vol_coordinates(ds);

    ds_rs = afni_resample(ds, fmt, orient); % in AFNI
    assertEqual(cosmo_fmri_orientation(ds_rs), orient);

    xyz_rs = cosmo_vol_coordinates(ds_rs);

    xyz_cell = mat2cell(xyz, [1 1 1], nfeatures);
    xyz_rs_cell = mat2cell(xyz_rs, [1 1 1], nfeatures);

    mp = cosmo_align(xyz_cell, xyz_rs_cell);
    assertElementsAlmostEqual(ds.samples(:, mp), ds_rs.samples, ...
                              'absolute', 1e-5);

function idxs = get_feature_index(ds, ijk)
    n = size(ijk, 2);
    idxs = zeros(1, n);

    i = ds.fa.i;
    j = ds.fa.j;
    k = ds.fa.k;

    for m = 1:n
        msk = i == ijk(1, m) & j == ijk(2, m) & k == ijk(3, m);
        idx = find(msk);
        assert(numel(idx) == 1);
        idxs(m) = idx;
    end

function orients = get_orients()
    % get all 48 possible orientations
    order = perms(1:3);
    signs = [1 0 1 0 1 0 1 0
             1 1 0 0 1 1 0 0
             1 1 1 1 0 0 0 0];

    labs = ['LR'; 'PA'; 'IS'];

    orients = cell(48, 1);

    for k = 1:6
        for j = 1:8
            orient = '   ';
            for dim = 1:3
                row = order(k, dim);
                orient(dim) = (labs(row, signs(row, j) + 1));
            end
            orients{(k - 1) * 8 + j} = orient;
        end
    end

function [success, output] = run_afni_command(cmd, raise)
    if nargin < 2
        raise = false;
    end
    dyld_path = getenv('DYLD_LIBRARY_PATH');
    cleaner = onCleanup(@()setenv('DYLD_LIBRARY_PATH', dyld_path));
    setenv('DYLD_LIBRARY_PATH', '/usr/local/bin');

    dyld_fb_path = getenv('DYLD_FALLBACK_LIBRARY_PATH');
    cleaner2 = onCleanup(@()setenv('DYLD_FALLBACK_LIBRARY_PATH', ...
                                   dyld_fb_path));
    setenv('DYLD_FALLBACK_LIBRARY_PATH', ...
           [getenv('DYLD_FALLBACK_LIBRARY_PATH') ':/sw/lib']);

    [s, output] = unix(cmd);
    success = s == 0;

    if ~success && raise
        error(output);
    end

function ds_rs = afni_resample(ds, fmt, orient)
    tmp_pat = ['tmp_%d' fmt];
    [fn_orig, cleaner] = get_temp_filename(tmp_pat);
    cosmo_map2fmri(ds, fn_orig);

    [fn_rs, cleaner2] = get_temp_filename(tmp_pat);

    % create command
    cmd = sprintf('3dresample -rmode NN -prefix %s -inset %s -orient %s', ...
                  fn_rs, fn_orig, orient);
    [success, output] = run_afni_command(cmd);
    if ~success
        error(output);
    end

    ds_rs = cosmo_fmri_dataset(fn_rs);

function temp_path = get_temp_path()
    c = cosmo_config();
    if ~isfield(c, 'temp_path')
        error('temp_path is not set in cosmo_config');
    end
    temp_path = c.temp_path;

function [fn, cleaner] = get_temp_filename(pat)
    temp_path = get_temp_path();
    k = 0;
    while true
        fn = fullfile(temp_path, sprintf(pat, k));
        if ~exist(fn, 'file')
            cleaner = onCleanup(@()delete(fn));
            return
        end
        k = k + 1;
    end