function [roi, G] = xbeta_roi_stats(roi_imgs, data_imgs, flags)
% FORMAT [roi G] = xbeta_roi_stats(roi_imgs, data_imgs, flags)
%
% returns data from ROI
% ROI defined by non-zero voxels in imgs roi_imgs
% data_imgs are images to extract data from
% images can be string or vol struct
% flags contains zero or more of chars
%              c - give 0 25 50 75 100 centile
% default flags is empty
%
% roi is cell array of structures containing fields
%    data_imgs (N by 1 vol structs of data images)
%               where N is no of data_imgs
%    roi_imgs   (vol struct of ROI images)
%    mean      (N by 1 array of means for each data_img)
%    var       (N x 1 variances)
%    n         (N x 1 n of non NaN voxels in ROI)
%   
% G is N x 1 matrix of globals, if required
%
% xbeta_roi_stats adapted from roi_stats by Jeff Cooper, May 2003
  
if nargin < 1
  roi_imgs = '';
end
if nargin < 2
  data_imgs = '';
end
if nargin < 3
  flags = '';
end

if isempty(roi_imgs)
  roi_imgs = spm_get(Inf, 'img', 'Select images defining ROIs');
end
if isempty(data_imgs)
  data_imgs = spm_get(Inf, 'img', 'Select data images');
end
if isempty(flags)
  flags = ' ';
end

if ischar(roi_imgs)
  roi_imgs = spm_vol(roi_imgs);
end
if ischar(data_imgs)
  data_imgs = spm_vol(data_imgs);
end
if ~isstruct(roi_imgs) | ~isstruct(data_imgs)
  error('Input args must be strings or structs')
end

% resampling = trilinear
hold = 1;

% centiles to return
%cents = [0 25 50 75 100];

%%%%%%%%%%%%%%%%%%%%%%%
% begin grand ROI loop
%%%%%%%%%%%%%%%%%%%%%%%

for r = 1:length(roi_imgs)

    curr_roi = roi_imgs(r);
    fprintf('\nWorking on roi %s\n', curr_roi.fname);

    % xyz mm of non-zero points in curr_roi
    [Y,XYZ] = spm_read_vols(curr_roi);
    tmp = find(Y);
    XYZ = [XYZ(:, tmp); ones(1, length(tmp))];
    clear Y

    nimgs = length(data_imgs);

    [roi{r}.max roi{r}.min roi{r}.n roi{r}.mean roi{r}.var] = deal(ones(nimgs,1)*NaN);

    fprintf('\nExtracting %g voxels from %g images\n',size(XYZ,2),nimgs)
    for i = 1:nimgs
       if rem(i,100)==0
           fprintf('\nGot to image %g\n',i)
       end
      % resample data at voxel centres of ROI
      ixyz = data_imgs(i).mat \ XYZ;
      d    = spm_sample_vol(data_imgs(i),...
    			ixyz(1,:),ixyz(2,:),ixyz(3,:),hold);
      
      
%keyboard;
      d = d(find(~(isnan(d))));

      roi{r}.name = spm_str_manip(curr_roi.fname, 'rtdu'); %trim pathname and suffix, remove illegal
      roi{r}.n(i) = length(d);
      if roi{r}.n(i)
        roi{r}.mean(i) = mean(d);
        roi{r}.var(i) = var(d);
        roi{r}.max(i) = max(d);
        roi{r}.min(i) = min(d);
        if nargout > 1
           G(i) = spm_global(data_imgs(i));
        end 
      end
    end

% end roi grand loop
end

