function cosmo_plot_slices(data, dim, slice_step, slice_start, slice_stop)
% Plots a set of slices from a dataset, nifti image, or 3D data array
%
% cosmo_plot_slices(data[, dim][, slice_step][, slice_start][, slice_stop])
%
% Inputs:
%  data         an fmri dataset (e.g., from cosmo_fmri_dataset), or a 3D
%               array with data. data should contain samples from a single
%               volume (sample) only.
%  dim          dimension according to which slices are plotted
%               (default: 3). Values between 1 and 3 allows using a
%               saggital, axial or coronal view; but the mapping between
%               the 3 numbers and 3 views depends on the particular
%               orientation of the data.
%  slice_step   step between slices (default: 1). If negative then
%               -slice_step indicates the number of slices
%  slice_start  the index of the first slice to plot (default: 1).
%  slice_stop   the index of the last slice to plot (default: the number of
%               slices in the dim-th dimension).
%
% Examples:
%    % plot an fMRI dataset struct with default options
%    cosmo_plot_slices(ds)
%
%    % plot an fMRI dataset struct along the second spatial dimension
%    cosmo_plot_slices(ds, 2)
%
%    % plot a random gaussian 3D array along the first dimension
%    cosmo_plot_slices(randn([40,50,20]),1)
%
%    % plot an fMRI dataset struct along the default spatial dimension
%    % every 5-th slice
%    cosmo_plot_slices(ds, [], 5)
%
%    % plot an fMRI dataset struct along the third spatial dimension
%    % with about 12 slices
%    cosmo_plot_slices(ds, 3, -12)
%
%    % plot an fMRI dataset struct along the third spatial dimension
%    % with about 12 slices, starting at slice 10 and stopping at slice 25
%    cosmo_plot_slices(ds, 3, -12, 10, 25)
%
% Notes:
%  - Using this function only really makes sense for fMRI-like data.
%  - This function does not provide a consistent orientation for slices,
%    as this depends on the voxel-to-world transformation matrix, which is
%    completely ignored in this function. Thus left-right and top-down
%    swaps can occur. Different datasets may provide different views, for
%    example dim=1 may give a saggital view if the dataset comes from one
%    program and an axial view if it comes from another program.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if nargin<2 || isempty(dim), dim=3; end
    if nargin<3 || isempty(slice_step), slice_step=-20; end
    if nargin<4 || isempty(slice_start), slice_start=1; end
    if nargin<5 || isempty(slice_stop), slice_stop=[]; end % set later

    data3D=get_data3D(data);

    % get min and max values across the entire volume
    data_lin=data3D(:);
    mn=min(data_lin);
    mx=max(data_lin);

    % shift it so that we can walk over the first dimension
    data_sh=shiftdim(data3D, dim-1);

    if isempty(slice_stop)
        slice_stop=size(data_sh,1);
    end

    if slice_step<0
        nslices=-slice_step;
        slice_step=ceil((slice_stop-slice_start+1)/nslices);
    end

    % determine which slices to show
    slice_idxs=slice_start:slice_step:slice_stop;
    nslices=numel(slice_idxs);

    plot_ratio=.8; % ratio between number of rows and colums
    nrows=ceil(sqrt(nslices)*plot_ratio);
    ncols=ceil(nslices/nrows);

    % use header depending on dim
    header_labels={'i','j','k'};

    % order of slices and whether the slice should be transposed
    xorder=[-1 -1 1];
    yorder=[-1 1 -1];
    do_transpose=[true false true];


    for k=1:nslices
        slice_idx=slice_idxs(k);
        slice=squeeze(data_sh(slice_idx,:,:));

        if xorder(dim)<0
            slice=slice(end:-1:1,:);
        end
        if yorder(dim)<0
            slice=slice(:,end:-1:1);
        end
        if do_transpose(dim)
            slice=slice';
        end

        subplot(nrows, ncols, k);
        imagesc(slice, [mn, mx]);
        title(sprintf('%s = %d', header_labels{dim}, slice_idx));
    end


function data3D=get_data3D(data)
    if isstruct(data)
        cosmo_check_dataset(data,'fmri');
        data4D=cosmo_unflatten(data);
        sz=size(data4D);
        if sz(1)>1
            error(['expected single volume data, but found %d '...
                    'volumes. To select a single volume, use '...
                    'cosmo_slice. For example, to show the %d-th '...
                    'volume from a dataset struct ds, use:\n\n   '...
                    'cosmo_plot_slices(cosmo_slice(ds,%d))\n'],...
                    sz(1),sz(1),sz(1));
        end

        if numel(sz)>4
            error('data must be 4D');
        end

        data3D=reshape(data4D, sz(2:end));
    elseif isnumeric(data) || islogical(data)
        ndim=numel(size(data));
        if ndim>3
            error('expected 3D image - did you select a single volume?');
        end
        data3D=data;
    else
        error('illegal input: expected dataset or 3D array');
    end
