function [roi, G] = roi_stats(roi_img, data_imgs, flags)
% FORMAT [roi G] = roi_stats(roi_img, data_imgs, flags)
%
% returns data from ROI
% ROI defined by non-zero voxels in img roi_img
% 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 structure containing fields
%    data_imgs (N by 1 vol structs of data images)
%               where N is no of data_imgs
%    roi_img   (vol struct of ROI image)
%    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
%
  
if nargin < 1
  roi_img = '';
end
if nargin < 2
  data_imgs = '';
end
if nargin < 3
  flags = '';
end

if isempty(roi_img)
  roi_img = spm_get(1, 'img', 'Select image defining ROI');
end
if isempty(data_imgs)
  data_imgs = spm_get(Inf, 'img', 'Select data images');
end
if isempty(flags)
  flags = ' ';
end

if ischar(roi_img)
  roi_img = spm_vol(roi_img);
end
if ischar(data_imgs)
  data_imgs = spm_vol(data_imgs);
end
if ~isstruct(roi_img) | ~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];

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

nimgs = length(data_imgs);
%[roi.mean roi.var roi.max roi.min roi.n G] = deal(ones(nimgs,1)*NaN); 
[roi.max roi.min roi.n ] = deal(ones(nimgs,1)*NaN);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%roi.cent = ones(nimgs, length(cents))*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);
  
  % mask with NaNrep
  NaNrep = spm_type(data_imgs(i).dim(4),'nanrep');
  data = d(d~=NaNrep);
  d = data;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  d = d(find(~isnan(d)));
 %%%%%%%% 

  roi.n(i) = length(d);
  if roi.n(i)
    roi.mean(i) = mean(d);
    roi.var(i) = var(d);
    roi.max(i) = max(d);
    roi.min(i) = min(d);

    
       
    if nargout > 1
      G(i) = spm_global(data_imgs(i));
    end
    
  end
  
end

%roi.data_imgs = data_imgs;
%roi.roi_img = roi_img;
