function output_file=fmri2pcnt(stat_file, input_file, ....
   exclude, mask_thresh, output_file_base)

%FMRI2PCNT converts stat files to percent of fmri data.
%
% OUTPUT_FILE = FMRI2PCNT( STAT_FILE, INPUT_FILE [, EXCLUDE 
%    [, MASK_THRESH [, OUTPUT_FILE_BASE]]] )
% 
% STAT_FILE is the file to be converted, usually an _ef, _sd, or _resid 
% file for magnitudes (t stat files and delays are of course unaffected).
%
% INPUT_FILE is the fMRI image file.
% 
% Both STAT_FILE and INPUT_FILE can be a single 4D file with multiple  
% frames, or a matrix of image file names, each with a single 3D frame,
% either ANALYZE (.img) or MINC (.mnc) format. Extra blanks are ignored. 
% File separator can be / or \ on Windows. Gzipped files are gunzipped on     
% unix. 
%
% MASK_THRESH: all voxels where the first frame of INPUT_FILE > MASK_THRESH
% are converted to %, the rest are set to zero. If MASK_THRESH=[], then it  
% uses fmri_mask_thresh.m to find a suitable value. Default is [].
% 
% OUTPUT_FILE_BASE is a matrix, one row for each row of STAT_FILE; 
% default is STAT_FILE minus extension plus _pcnt.
%
% EXCLUDE is a list of frames that should be excluded from the
% analysis. This must be used with Siemens EPI scans to remove the
% first few frames, which do not represent steady-state images.
% Default is [1].
%
% OUTPUT_FILE is a matrix whose rows are the output file names, 
% OUTPUT_FILE_BASE.mnc or .img,
 
%############################################################################
% COPYRIGHT:   Copyright 2002 K.J. Worsley
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center, 
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada. 
%              worsley@math.mcgill.ca, liao@math.mcgill.ca
%
%              Permission to use, copy, modify, and distribute this
%              software and its documentation for any purpose and without
%              fee is hereby granted, provided that the above copyright
%              notice appear in all copies.  The author and McGill University
%              make no representations about the suitability of this
%              software for any purpose.  It is provided "as is" without
%              express or implied warranty.
%############################################################################

% Defaults:

if nargin < 3
   exclude=[1];
end
if nargin < 4
   mask_thresh=[];
end

numstatfiles=size(stat_file,1);
cid=min(find(stat_file(1,:)=='.'));
ext=stat_file(1,cid+(0:3));

output_file=[];
if nargin < 5
   nid=size(stat_file,2)+9;
   for i=1:numstatfiles
      cid=min(find(stat_file(i,:)=='.'));
      id=[deblank(stat_file(i,1:(cid-1))) '_pcnt' ext];
      output_file=[output_file; [id repmat(' ',1,nid-cid)]];
   end
else
   nid=size(output_file_base,2)+9
   for i=1:size(output_file_base,1)
      id=[deblank(output_file_base(i,:)) '_pcnt' ext]
      output_file=[output_file; [id repmat(' ',1,nid-length(id))]]
   end
end
output_file=deblank(output_file);

d=fmris_read_image(deblank(stat_file(1,:)),0,0);
numstatframes=d.dim(4);

numinputfiles=size(input_file,1);
d=fmris_read_image(input_file(1,:),0,0);
numinputframes=max(d.dim(4),numinputfiles);

allpts = 1:numinputframes;
allpts(exclude) = zeros(1,length(exclude));
keep = allpts( find( allpts ) );
n=length(keep);

if isempty(mask_thresh)
   mask_thresh=fmri_mask_thresh(input_file);
end

for slice=1:d.dim(3)
   slice
   if numinputfiles==1
      d=fmris_read_image(input_file,slice,keep);
      average=squeeze(mean(d.data,4));
   else
      average=zeros(d.dim(1:2));
      for i=1:n
         d=fmris_read_image(deblank(input_file(keep(i),:)),slice,1);
         average=average+d.data;
      end
      average=average/n;
   end
   
   for i=2:numstatfiles
      if numstatfiles==1
         d=fmris_read_image(deblank(stat_file),slice,i);
      else
         d=fmris_read_image(deblank(stat_file(i,:)),slice,1);
      end
      d.data=d.data./(average+(average<=0)).*(average>mask_thresh)*100;
      d.file_path='';
      d.file_name=deblank(output_file(i,:));
      d.parent_file=deblank(stat_file(i,:));
      d.dim(4)=numstatframes*(numstatframes>1);
      fmris_write_image(d,slice,1:numstatframes);
   end
end
         
return


output_file_base=['ab    '; 'cde   ']

output_file_base=['c:/keith/results/ha_100326_hot';
                  'c:/keith/results/ha_100326_wrm';
                  'c:/keith/results/ha_100326_hmw']
stat_file=['c:/keith/results1/ha_093923_hmw_mag_ef.mnc';
                    'c:/keith/results1/ha_100326_hmw_mag_ef.mnc';
                    'c:/keith/results1/ha_101410_hmw_mag_ef.mnc';
                    'c:/keith/results1/ha_102703_hmw_mag_ef.mnc'];
input_file='c:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc';
exclude=[1 2]
output_file=fmri2pcnt(stat_file, input_file, [1 2]);