function [base, flup] = make_diffs(thr, means_too, base, flup)
%make_diffs -- create difference images for selected (matched) images,
% follow-up(s) minus baseline, where baselines are selected first, and
% follow-ups are selected in the same order, repeating this order if
% multiple follow-up images are selected. String matrices or cell-strings
% can be passed in for base and flup to avoid GUI selection.
% 
% Syntax:
%  make_diffs(absolute_threshold, means_too, base, flup)
%
% Write mean images as well if means_too flag is true.
% Save outputs with follow-up names prefixed with 'diff_' ('mean_').
% Output images have NaN where either input was below absolute_threshold.
%
% Examples:
%  make_diffs; % GUI selection, no threshold, just create differences
%  make_diffs(-inf, true); % GUI selection, no threshold, diffs and means
%  make_diffs(thr, true, base, flup); no GUI, absolute threshold of thr

% check spm version:
if exist('spm_select','file') % should be true for spm5
    spm5 = 1;
    select = @(msg) spm_select(inf, 'image', msg);
    fparts = @spm_fileparts; % separates ',number', from e.g. 'blah.nii,1'
    float32 = spm_type('float32');
elseif exist('spm_get','file') % should be true for spm2
    spm5 = 0;
    select = @(msg) spm_get(inf, 'img', msg);
    fparts = @fileparts; % standard Matlab fileparts command (cf. above)
    float32 = spm_type('float');
else
    error('make_diffs:spm_missing',...
        'Failed to find spm_get/spm_select, please add SPM to Matlab path')
end
spm_defaults;

% default to not threshold 
if ~exist('thr','var') || isempty(thr)
    thr = -inf;
end

% get filenames
if ~exist('base', 'var') || isempty(base)
    base = select('Select baseline images');
elseif iscellstr(base)
    base = char(base);
end
if ~exist('flup', 'var') || isempty(flup)
    flup = select('Select follow-up images in set(s) with same order');
elseif iscellstr(flup)
    flup = char(flup);
end
drawnow;
N = size(base, 1);
M = size(flup, 1);
if rem(M, N) ~= 0
    error('make_diffs:numbers_mismatch',...
        'Number of follow-ups (%d) must be multiple of #baselines (%d)',...
        M, N)
end
P = M / N; % number of follow-ups per baseline

% process files
for n = 1:N
    for p = 1:P
        m = n + (p-1)*N; % corresponding follow-up
        bvol = spm_vol(base(n, :));
        fvol = spm_vol(flup(m, :));
        [pth fname ext] = fparts(fvol.fname);
        dfile = fullfile(pth, ['diff_' fname ext]);
        mfile = fullfile(pth, ['mean_' fname ext]);

        bimg = spm_read_vols(bvol);
        fimg = spm_read_vols(fvol);

        % additive mask to place nans where both ims below threshold
        % nanmask = nan(size(bimg));
        % nanmask(bimg >= thr) = 0;
        % nanmask(fimg >= thr) = 0;

        % additive mask to place nans where either im below threshold
        nanmask = zeros(size(bimg));
        nanmask(bimg < thr) = nan;
        nanmask(fimg < thr) = nan;

        dimg = (fimg - bimg) + nanmask;
        dvol = fvol;
        dvol.fname = dfile;
        dvol.descrip = 'Difference image';
        if spm5
            dvol.dt(1) = float32;
        else
            dvol.dim(4) = float32;
        end
        spm_write_vol(dvol, dimg);
    
        if exist('means_too','var') && any(means_too)
            mimg = (fimg + bimg)/2 + nanmask;
            mvol = bvol;
            mvol.fname = mfile;
            mvol.descrip = 'Mean image';
            if spm5
                mvol.dt(1) = float32;
            else
                mvol.dim(4) = float32;
            end
            spm_write_vol(mvol, mimg);
        end
    end
end

if nargout < 1
    clear base; % don't spew to screen if no output args and ';' omitted!
end