function p = calc_roi_mean(roi_imgs, data_imgs)
%--------------------------------------------------------------------------
% roi_stat_extraction
%
% Extracts mean, median, min, max, and intervals of ROIs
%
% USAGE:
% p = calc_roi_mean(roi_imgs, data_imgs, interval)
%
% INPUT:
% roi_imgs:     char array of ROI image filenames. Means and medians are
%               extracted where ROI~=0
% data_imgs:    char array of data image filenames for which the means and
%               medians are extracted
%
% OUTPUT:
% p.mean_act:   mean timeseries 

% load the default spm settings
spm('Defaults','FMRI')

hold = 0; % 0=nearest neighbor

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use SPM functions to read image information
p.v_roi  = spm_vol(roi_imgs);
p.v_data = spm_vol(data_imgs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% loop through the ROIs
for r = 1:size(roi_imgs,1)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % xyz mm of non-zero points in roi img
    [Y,XYZ] = spm_read_vols(p.v_roi(r));
    tmp     = find(Y);
    XYZ     = [XYZ(:, tmp); ones(1, length(tmp))];
    clear Y
    
    fprintf('Extracting %g voxels from %g images\n', size(XYZ,2), length(p.v_data));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % loop through the data images        
    for i = 1:length(p.v_data)
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % find voxel elements for which to extract stats
        ixyz = p.v_data(i).mat \ XYZ;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % then sample the data from the volume at those locations
        x = spm_sample_vol(p.v_data(i),ixyz(1,:),ixyz(2,:),ixyz(3,:),hold);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % find NaN means
        nans = isnan(x);
        x(nans) = 0;
        n = sum(~nans);
        n(n==0) = NaN;
        m = sum(x) ./ n;
        p.mean_act(i,r) = m;
                
    end
end


