function [roigraph]= spm_roi_compute(roistruct)
%  USAGE:
%  [roigraph]=spm_roi_compute(roistruct)
%
%  ROI selective averaging tool - part of the SPM ROI Toolbox
%  - see http://spm-toolbox.sourceforge.net/ for more details
%
% this version of the ROI graphing tool:
% - encompasses both blocked and event-related designs
% - allows flexible adjustment for confounding conditions
% - allows flexible application of high-pass and/or low-pass filtering
% - computes average using unbiased (finite impulse response) model
%
% NOTE: data are no longer plotted as they were in previous versions
%
% all arguments are now passed through the roistruct structure,
% which has the following fields:
%  REQUIRED FIELDS:
%  roifile: name(s) of image files defining clusters
%  spmfile: location of SPM.mat file from which to compute results
% 
%  OPTIONAL FIELDS:
%  outfilename: name of output file for results (defaults to
%      roiX_avg.txt where X is the index of the ROI)  
%  TER: effective temporal resolution of selective averaging
%      - when TER<TR, data are oversampled (defaults to TR)
%  adjust_flag: controls adjustment for other conditions, taking one
%    of the following values (default='cond'):
%      'raw': perform averaging on raw data with no adjustment
%      'filt': perform averaging on data after high-pass and low-pass
%         filtering as specifified in the statistical model
%      'cond': perform averaging after high/lo pass filtering as well as
%         adjustment for other modeled effects (this is the standard SPM
%         approach)
%  block_summ: if block_summ > 0, integrated signal change will be 
%    recorded to the output file in addition to the timcourse
%  post_time_window: the window after stimulus presentation over which
%    selective averaging is performed
%    - for event-related designs, default = 24 seconds
%    - for blocked designs, default = length of longest block + 10
%    seconds
%  neg_time_window: length of pre-stimulus time window (default = 4 s)
%
%  This software has been developed by Dr. Russell Poldrack and the
%  spm-toolbox development team, based on previous code from SPM99, 
%  and is made available under the terms of the GNU Public License.   
%  This software is meant for research only and should not be used for 
%  clinical purposes.
%

% Input checking
%-----------------------------------------------------------------------
  
%keyboard;
% check for proper roistruct argument and required fields
die_now=0;
try,
  if isempty(roistruct.roifile), die_now=1; end;
  if isempty(roistruct.spmfile), die_now=1; end;
catch,
  die_now=1;
end;
  
if die_now,
  help spm_roi_compute;
  return;
end;

% create the optional fields in roistruct if they don't exist
if ~isfield(roistruct,'outfilename'),roistruct.outfilename=[];end;
if ~isfield(roistruct,'TER'),roistruct.TER=[];end;
if ~isfield(roistruct,'adjust_flag'),roistruct.adjust_flag=[];end;
if ~isfield(roistruct,'block_summ'),roistruct.block_summ=[];end;
if ~isfield(roistruct,'PosTimeWindow'),roistruct.PosTimeWindow=[];end;
if ~isfield(roistruct,'NegTimeWindow'),roistruct.NegTimeWindow=[];end;


% set up default settings for various aspects of analysis
%-----------------------------------------------------------------------
% specify the time window for selective averaging
if isempty(roistruct.NegTimeWindow),roistruct.NegTimeWindow=4;end;

if isempty(roistruct.adjust_flag),
  roistruct.adjust_flag='cond'; % no adjustment
else,
  roistruct.adjust_flag=lower(roistruct.adjust_flag); 
end;



% load stats data to be analyzed
%-----------------------------------------------------------------------
% always use contrast 1, the omnibus F stat
contrastnum=1;

% load data using special version of getSPM that bypasses GUI
[refSPM,refVOL,refxX,refxCon,refxSDM]=spm_roi_getSPM(roistruct.spmfile, ...
						  contrastnum);
% check to make sure raw data are saved in Y.mad
if exist(fullfile(refSPM.swd,'Y.mad')) ~= 2
  spm('alert"',{'No raw data (Y.mad) saved with this analysis:',...
		'Cannot perform ROI analysis...'},...
		mfilename,0,1);
  return;
end;

% check to make sure UFp is high enough
if refxSDM.UFp < 0.99, 
  fprintf('\nWARNING: UFp = %0.3f\ndata for some voxels may be missing\n\n',refxSDM.UFp);
end;

if ~isfield(refxSDM,'Sess')
    fprintf(['Inappropriate dataset for ROI analysis\nSPM.mat must be a' ...
	     ' single-subject (fixed effects) analysis']);
    return;
end

if ~isempty(roistruct.TER),
  if mod(refxX.RT,roistruct.TER)~=0,
    spm('alert',{'TR must be evenly divisible by TER'},mfilename,0,1);
    return;
  end;
else,
  roistruct.TER=refxX.RT;
end;

if roistruct.TER ~= refxX.RT & ~exist('resample'),
  fprintf(['WARNING: Oversampling requires the resample function\n from' ...
	   ' the MATLAB Signal Processing Toolbox\n Setting TER=TR\n\n']);
  
  roistruct.TER=refxX.RT;
end;

%-Find nearest voxel [Euclidean distance] in point list & update GUI
%-----------------------------------------------------------------------
if ~length(refSPM.XYZ)
   spm('alert!','No suprathreshold voxels!',mfilename,0);
   return;
end


maxroi=size(roistruct.roifile,1);

for roicount=1:maxroi
     fprintf('Processing roi #%d of %d: %s\n',roicount,maxroi,roistruct.roifile(roicount,:));
     empty_log=fopen('empty_voxels.log','a');
     fprintf(empty_log, 'ROI #%d of %d: %s\n',roicount,maxroi,roistruct.roifile(roicount,:));
     fclose(empty_log);

    % set these using the reference each time around, because they can be
    % changed by intervening processing inside the roicount loop
    SPM=refSPM;
    VOL=refVOL;
    xX=refxX;
    xCon=refxCon;
    xSDM=refxSDM;
    
    % load roi image

    roiV=spm_vol(roistruct.roifile(roicount,:));
    [roiY,roiXYZ]=spm_read_vols(roiV);
    roiXYZ=roiXYZ(:,find(roiY>0));

    numroivox=size(roiXYZ,2);
    cluster_SPM=[];
    counter=0;
    empty_count=0;
    XLAB=[];
    YLAB=[];
 
    for i=1:numroivox, % loop through all voxels in the roi image
      point=tal2vox(roiXYZ(:,i),VOL);
      v=find(SPM.XYZ(1,:)==point(1) & SPM.XYZ(2,:)==point(2) & SPM.XYZ(3,:)==point(3));

      % check for empty voxels
      if ~isempty(v) & SPM.QQ(v)>0
	counter=counter+1;
	cluster_SPM(counter)=v;
      else,
        empty_log=fopen('empty_voxels.log','a');
        fprintf(empty_log,'voxel %d in ROI %d (%d,%d,%d): empty\n',i,...
		roicount,roiXYZ(1,i),roiXYZ(2,i),roiXYZ(3,i));
        fprintf('voxel %d in ROI %d (%d,%d,%d): empty\n',i,...
		roicount,roiXYZ(1,i),roiXYZ(2,i),roiXYZ(3,i));
        fclose(empty_log);
        empty_count=empty_count+1;
      end;
    end

    %-Extract required data from results files
    %=======================================================================
    %-Get (approximate) raw data y from Y.mad file
    %-NB: Data in Y.mad file is compressed, and therefore not fully accurate
    %     Therefore, parameters & ResMS should be read from the image files,
    %     rather than recomputing them on the basis of the Y.mad data.
    %-----------------------------------------------------------------------

    sizeinfo=spm_extract(fullfile(SPM.swd,'Y.mad'));
    nscans=sizeinfo(1);
    npoints=sizeinfo(2);
    
    
    ss    = length(xSDM.Sess); % # of sessions
    tr=[];  %                    # of conditions per session
    for s = ss
	tr    = [tr length(xSDM.Sess{s}.pst)]; % # conds per sess
    end
    tr    = max(tr);
    
    if ~isempty(roistruct.PosTimeWindow),
       NPosEst = round(roistruct.PosTimeWindow/roistruct.TER);
    else,
      % Determine time window based on basis functions
      fprintf('determining window length from basis functions...\n');
      bf_window=zeros(ss,tr);
      for s=ss,
        for t=length(xSDM.Sess{s}.pst),
	  tmp_bf=xSDM.Sess{s}.bf{tr};
	  if size(tmp_bf,2)>1,tmp_bf=sum(tmp_bf,2);end; 
	  pos_bf_timepoints=find(tmp_bf>0);
	  bf_window(s,t)=max(pos_bf_timepoints)*refxX.dt;
        end;
      end;
      bf_window=max(max(bf_window));
      % add 8 seconds to end of window
      NPosEst = round(bf_window/roistruct.TER) + round(8/roistruct.TER); ...
		

    end;

    NNegEst = round(roistruct.NegTimeWindow/roistruct.TER); 
    ntp_est=NPosEst + NNegEst;  % total # of timepoints estimated
    
    % set up variables for saving results of each voxel
    time_course=zeros(size(cluster_SPM,2),ntp_est,ss,tr);
    avg_time_course=zeros(ntp_est,ss,tr);


    fprintf('Calculating graphs for %d voxels (of %d voxels in ROI)\n',...
	    counter,numroivox);

    for roi_voxel=1:counter
        % load data from Y.mad - use try/catch to deal with errors
        % reading file
        try
	   y=spm_extract(fullfile(SPM.swd,'Y.mad'),SPM.QQ(cluster_SPM(roi_voxel)));
        catch,
           fprintf('caught error reading Y.mad! aborting!\n');
           return;
        end;
	if isempty(y), 
	  fprintf('no data for voxel %d\n',roi_voxel);
	end;
	
	% apply the pre-specified filters (from the analysis) to the raw
        % data - this may include HPF and LPF depending upon analysis
	if ~strcmp(roistruct.adjust_flag,'raw'),
	  refy=spm_filter('apply',xX.K,y);
	else,
	  refy=y;
	end;
	
	% find the voxel in the data that is nearest to the voxel in roi
        % space
	
	talcoord=vox2tal(SPM.XYZ(:,cluster_SPM(roi_voxel)),VOL);
	[xyz,i] = spm_XYZreg('NearestXYZ',talcoord,SPM.XYZmm);    

	fprintf('y for voxel %d (%d %d %d): min %3.1f, max %3.1f\n',...
		roi_voxel,xyz,min(y),max(y));

	%-Reset pointer, compute voxel indices, compute location string
	%-----------------------------------------------------------------------
	rcp     = VOL.iM(1:3,:)*[xyz;1];
	if roi_voxel==1, XYZstr  = sprintf(' at [%g, %g, %g]',xyz);end;


	%-Get parameter estimates, ResMS, (compute) fitted data & residuals
	%=======================================================================
	%     Fitted data & residuals are for temporally smoothed model

	%-Parameter estimates: beta = xX.pKX*xX.K*y;
	%-----------------------------------------------------------------------
	beta  = ones(length(xSDM.Vbeta),1);
	for i = 1:length(beta)
		beta(i) = spm_sample_vol(xSDM.Vbeta(i),rcp(1),rcp(2),rcp(3),0);
	end

	%-Compute residuals
	%-----------------------------------------------------------------------
	if isempty(y)

		% make R = NaN so it will not be plotted
		%---------------------------------------------------------------
		R   = NaN*ones(size(xX.xKXs.X,1),1);

	else
		% residuals
		%---------------------------------------------------------------
		R   = spm_filter('apply',xX.K,y) - xX.xKXs.X*beta;

	end

	%-Residual mean square: ResMS = sum(R.^2)/xX.trRV;
	%-----------------------------------------------------------------------

	ResMS = spm_sample_vol(xSDM.VResMS,rcp(1),rcp(2),rcp(3),0);
	SE  = sqrt(ResMS*diag(xX.Bcov));

	Cplot='Event/epoch-related responses';

	% modeling evoked responses based on Sess
	%----------------------------------------------------------------------

	% get sessions
	%--------------------------------------------------------------
	Sess  = xSDM.Sess;
	rep   = Sess{1}.rep;  % replicated?
	if ss > 1
                        ss=[1:ss];
	end

	% get index for conditions
	%--------------------------------------------------------------
	  tr=[];
	  for s = ss
		tr    = [tr length(Sess{s}.pst)]; % # conds per sess
	  end
	  tr    = max(tr);
	  if tr > 1
                tr=[1:tr];
	  end
       


	% cycle over selected events
	%--------------------------------------------------------------
	dx      = xX.dt; % time bin

	% reconstruct response with filtering
	%--------------------------------------------------------------

	for s = ss % sessions
	    for t = tr % conditions

              % basis functions, filter and parameters
	      %------------------------------------------------------
	      i      = Sess{s}.row(:); % index vector for images
	      cond_exists=1;
              try,
		% this is just a trick to see if there is data
		% only take the canonical bfs, not the derivs
		  j      = 1:size(Sess{s}.sf{t},2):length(Sess{s}.ind{t});
              catch,
		  fprintf('No data for Session %d, condition %d\n',s,t);
		  cond_exists=0;
	      end;
	      
	      if (cond_exists>0)
		% this takes a slightly different approach from the
                % original code which was based on spm_graph.m
		% in this version I get the regressors from the
                % design matrix and use that to recreate the signal
		% rather than recreating the basis functions
		% rather than using the basis functions from Sess.bf
		
		if strcmp(roistruct.adjust_flag,'cond'),
		  j      = Sess{s}.col(Sess{s}.ind{t}(j)); % use only
                                                           % columns for
                                                           % this
                                                           % specific
                                                           % effect
		else 
		  j      = Sess{s}.col(:); % use entire
                                                        % desmtx for this
                                                        % session
							
		end;
		B      = beta(j); % parameter estimate for column j

		X = xX.X(Sess{s}.row,j);
		% fitted responses, adjusted data and standard error
		%------------------------------------------------------

		Y      = X*B; % param est * basis fns
		
		% this stuff is necessary for later plotting but is no
                % longer used to create adjusted data
		
		X      = Sess{s}.bf{t}; % basis function for cond t
		q      = 1:size(X,1);   % index vector for basis fn
		x      = (q - 1)*dx;    % basis fn in steps of dx
		se = 0;
		pst    = Sess{s}.pst{t};
		bin    = round(pst/dx);
		q      = find((bin >= 0) & (bin < size(X,1)));

		
		if ~strcmp(roistruct.adjust_flag,'raw'),
		  y      = Y + R(i);   % reconstruct adjusted data by
				       % adding residuals back to estimate
		else,         
		 y=refy(i');           % no adjustment - use raw data
		end;

		% oversample the reconstructed data to TER
		if (roistruct.TER~=xX.RT),
		  y=resample(y,xX.RT/roistruct.TER,1);
		end;
	
		% create FIR design matrix
		FIRconf=struct(...
			'NPosEst', NPosEst,...
			'NNegEst', NNegEst,...
			'condition', t,...
		        'TER',roistruct.TER);
	   
		FIR_X=spm_roi_fir_desmtx2(xX,Sess{s},FIRconf);
		FIR_X=[FIR_X ones(size(FIR_X,1),1)];
		
		% note: since the FIR design matrix will always have full
                % rank, I use the simple (and faster) form of the
                % pseudoinverse rather than the SVD-based version using
                % pinv
		
		FIRbeta=inv(FIR_X'*FIR_X)*FIR_X'*y;
		time_course(roi_voxel,1:size(FIRbeta,1)-1,s,t)=FIRbeta(1:(size(FIRbeta,1)-1));
		increment=roistruct.TER;
		anchor=[-1*FIRconf.NNegEst*roistruct.TER:roistruct.TER:(FIRconf.NPosEst-1)*roistruct.TER];
		time_point=ntp_est+1;
		trial_time=(NPosEst + 1)*roistruct.TER;
		   
	      else, % if cond does not exist - put zeros into timecourse
		     time_course(roi_voxel,:,s,t)=0;
	   end; % if cond_exists

           try,
              cond_name=Sess{s}.name{t};
           catch,
              cond_name='missing';
           end;
	   

	  end;
	end;

    end % looping through voxels

      % average timecourses over voxels
      for s = ss
	for t = tr
	    for pt = 1:time_point-1
	       avg_time_course(pt,s,t)=mean(time_course(:,pt,s,t));
	    end
	end
      end

    keyboard;
    foo=roistruct.roifile(roicount);
    if isempty(roistruct.outfilename),
      save_file=sprintf('roi%d_avg.txt',roicount);	
    else,,
       save_file=sprintf('%s',roistruct.outfilename);
    end;

    try,
     fid = fopen(save_file,'w');
     fprintf(fid,'%%%s\n%%%s\n',roistruct.roifile(roicount,:),roistruct.spmfile);
     fprintf(fid,'Sess\tCond\t');
     for s=1:time_point-1,fprintf(fid,'%0.3f\t',anchor(s));end;
     fprintf(fid,'\n');
   
     for s = ss
	for t = tr
           try,
              cond_name=Sess{s}.name{t};
           catch,
              cond_name='missing';
           end;
          
          fprintf(fid,'%d\t%s\t',s,cond_name);
          for tp=1:time_point-1,
	   fprintf(fid,'%0.3f\t',avg_time_course(tp,s,t));
          end
          fprintf(fid,'\n');
	end
     end
     fclose(fid);
    catch,
      fprintf('cannot open output file %s!\n',save_file);
    end;
    
    fprintf('done with ROI #%d\n',roicount);
    if empty_count > 0,
       fprintf('WARNING: %d voxels were not included in the ROI analysis\nCheck empty_voxels.log for more information\n',empty_count);
       try,
        empty_log=fopen('empty_voxels.log','a');
        fprintf(empty_log, 'Found total of %d empty voxels\n\n',empty_count);
        fclose(empty_log);
       catch,
	 fprintf('cannot open empty voxels log\n');
       end;
    else
	 empty_log=fopen('empty_voxels.log','a');
	 fprintf(empty_log, 'no empty voxels found\n\n');
	 fclose(empty_log);
    end
    
  end % for roifiles

fprintf('done with last ROI #%d of %d\n',roicount,maxroi);


roigraph=struct(...
'avg_time_course',    avg_time_course,...
'time_course',        time_course, ...
'SPM',                SPM,...
'VOL',                VOL,...
'xX',                 xX,...
'xSDM',               xSDM,...
'xCon',               xCon);









