function [vols avi mov] = vols2avi(vols, avi, fps, ax, vox, step, rw, rev)

if ~exist('vols', 'var') || isempty(vols)
    vols = spm_select(inf, 'image');
end
if ischar(vols)
    vols = spm_vol(vols);
end

if ~exist('avi', 'var') || isempty(avi)
    [pth fnm] = fileparts(vols(1).fname);
    avi = fullfile(pth, [fnm '.avi']);
end

if ~exist('fps', 'var') || isempty(fps)
    fps = 24;
end

if ~exist('ax', 'var') || isempty(ax)
    ax = 2;
end
if ~exist('vox', 'var') || isempty(vox)
    vox = 55;
end

if ~exist('step', 'var') || isempty(step)
    step = 1;
end

if ~exist('rw', 'var') || isempty(rw)
    rw = false;
end

if ~exist('rev', 'var') || isempty(rev)
    rev = false;
end

dim = vols(1).dim;
odim = dim; odim(ax) = [];
[I J K] = ndgrid(1:dim(1), 1:dim(2), 1:dim(3));
switch ax
    case 1, inds = find(I == vox);
    case 2, inds = find(J == vox);
    case 3, inds = find(K == vox);
end
I = I(inds);
J = J(inds);
K = K(inds);

f = figure; colormap gray
N = length(vols);
m = 0; % movie frame
for n = 1:step:N
    m = m + 1;
    slice = spm_sample_vol(vols(n), I, J, K, 1);
    slice = fliplr(rot90(reshape(slice, odim))); %%% FIXME: check for non-sagittal
    imagesc(slice);
    axis image; axis off;
    mov(m) = getframe;
end

% rewind sequence (gives smooth looping)
if rw
    for n = N-1:-step:2
        m = m + 1;
        slice = spm_sample_vol(vols(n), I, J, K, 1);
        slice = fliplr(rot90(reshape(slice, odim))); %%% FIXME: check for non-sagittal
        imagesc(slice);
        axis image; axis off;
        mov(m) = getframe;
    end
end

% reverse whole sequence (not very neatly...)
if rev
    tmp = mov;
    N = length(mov);
    for n = 1:N
        mov(n) = tmp(N+1-n);
    end
end

movie2avi(mov, avi, 'FPS',fps);
close(f);