

function [ROI, Cond, Vf] = roi_extractYmad(Design,Action,varargin);

% FORMAT roi_extractYmad('Epoch','Plot')
% Plots Bargraph of average intensity value for each condition.
% 
% FORMAT roi_extractYmad('Event','Plot')
% Plots average peristimulus timecourse (PST) for each condition.
%
% Argument values
% Design - 'Epoch','Event'
% Action - '','Plot','Save'
%
% varargin{1} can be 'Ymad','Adj','Fit','Vfiles', but only if
% Design is 'Event'. The variable Extract will be set appropriately.
% Otherwise for Design 'Epoch' data are only extract from Ymad.
%
% varargin{2} can be the full-path results directory, or resdir:
% e.g., '/vol/disk75/kalina/int_ext/No1_Karina_Bryant'
%
% varargin{3} can be the full-path roi file, or roi_path:
% e.g., '/vol/disk75/kalina/int_ext/No1_Karina_Bryant/roi_34_62_0.tal'
%
% varargin{4} can be the preprocessing option: use 'y' for preproc
% 
% Example use could be:
% [ROI, Cond, Vf]=roi_extractYmad('Event','Plot','Fit',resdir,roi_path,'y');
%
% This program uses another routine, kc_get_design.m, in order to
% gather information about the design specifics needed for trial 
% (or block) averaging. For more information, please see 
% http://www-psych.stanford.edu/~kalina/SPM99/Tools/roi.html
% _________________________________________________________________
% @(#)roi_extractYmad.m (roi_plot) 3.2 Kalina Christoff  2001-08-07

if nargin < 2 , error(['Not enough input arguments. Type ''help' ...
						  ' roi_extractYmad'' for more information'])
end

if nargin == 3 & ~strcmpi(Design,'Event'),
  error(['3 input arguments can only be used when Design is ''Event'''])
elseif nargin > 2 & strcmpi(Design,'Event')
  Extract=varargin{1};
elseif nargin == 2 & strcmpi(Design,'Event');   % assign Ymad Extract
  Extract = 'Ymad';                             % as default for 'Event'
end


						  
cwd = pwd;

% Ask for results directory (unless specified) and confirm its name
% -----------------------------------------------------
if nargin == 5 | nargin == 6
  resdir = varargin{2};
else resdir = spm_get(-1,'',...
			'Select directory containing analysis files');
end

cd(resdir);
fprintf('\n\n-----------------------------------------------\n');
fprintf('Working on data in\n%s ...',resdir);
fprintf('\n-----------------------------------------------\n\n');


load SPM; load xCon;

% Get the anatomical .tal ROI file (unless given as input arg)
% --------------------------------------
if nargin == 5 | nargin == 6
  roi_path = varargin{3};
else roi_path=spm_get(1,'*.tal','Select a .tal ROI file');
end


% Assign preprocessing option
if nargin == 6, prepr = varargin{4}; end

 
% Read in coordinates into ROImm
% -------------------
ROI.XYZmm = spm_load(roi_path)';
if size(ROI.XYZmm,1)~=3,error('The .tal should have 3 columns.'),,end


% Convert to ROI.XYZmm to pixel coordinates in ROI.XYZ
% ----------------------------------------------------
iM=inv(M);
ROI.XYZ = iM(1:3,:)*[ROI.XYZmm; ones(1,size(ROI.XYZmm,2))];
ROI.XYZ = round(ROI.XYZ);

% Extract raw data from Y.mad
% ---------------------------
if   exist(fullfile('.','Yidx.mat'),'file') & ...
     exist(fullfile('.','Y.mad'),'file')
     load('Yidx.mat');
else fprintf('\n\nThere is no Y.mad file in \n%s\n\nAborting...\n\n',resdir);
     error('Make sure the fMRI_UFp variable in your spm_defaults is not set to 0.');
end


ROI.Ymad_idx=zeros(1,size(ROI.XYZ,2)); 

   for vox=1:size(ROI.XYZ,2)
     

     % get only the indexes of existing (unmasked) voxels
     %---------------------------------------------------
     xyz_idx            = round(find(...
       (XYZ(1,:)==ROI.XYZ(1,vox)) & ...
       (XYZ(2,:)==ROI.XYZ(2,vox)) & ...
       (XYZ(3,:)==ROI.XYZ(3,vox)) ));

    if ~isempty(round(Yidx)==xyz_idx);
     ymad_idx = round(find(Yidx==xyz_idx));
    end

    if ~isempty(ymad_idx)
     ROI.Ymad_idx(vox)= ymad_idx;
    end

   end
	
	nvox_roi = size(XYZ,2);
	nvox_surv= find(ROI.Ymad_idx);

   % Bail out if there are no suprathreshold voxels
   %---------------------------------------------------------------
   if isempty(find(ROI.Ymad_idx))
    fprintf('\n\nSorry, there are no suprathreshold voxels in this ROI.\n');
    fprintf('It may be best to set fMRI_UFp=1 in your spm_defaults.m\n\n');
    error  ('Aborting...');
   end

   % Eliminate subtreshold and masked voxels from ROI structure
   % ---------------------------------------------------------
	nvox_roi = size(ROI.XYZ,2);

   ROI.XYZmm    = ROI.XYZmm   (:,find(ROI.Ymad_idx));
   ROI.XYZ      = ROI.XYZ     (:,find(ROI.Ymad_idx));
   ROI.Ymad_idx = ROI.Ymad_idx(:,find(ROI.Ymad_idx));
	
	% Announce how many voxels from this roi survived
	nvox_surv= size(ROI.XYZ,2);
	if nvox_surv < nvox_roi
	   fprintf('\n\nThis ROI had %g voxels. Information about %g of them\n',...
			  nvox_roi, nvox_surv);
 	   fprintf(['is available in the Y.mad file. The average intensity \n',...
				'in the plot is an average over the available voxels only.\n',...
				'If this number seems insufficient, consider setting the\n',...
				'fMRI_UFp variable to a higher value and rerunning the', ...
				 ' statistics.\n\n'])
	end
	

ROI.Ymad     = spm_extract('Y.mad',ROI.Ymad_idx); 
ROI.Ymad_avg = mean(ROI.Ymad,2);


% Extract adjusted data or fitted response, if requested
% --------------------------------------------------------
if strcmpi(Extract,'Adj') | strcmpi(Extract,'Fit')
  ROI = roi_fitadj(ROI,xX,xCon,Vbeta,VResMS,1);
end



% ================================================
% Avverage over the PST
% ================================================

   % Get a Cond structure containing design parameters conditions
   % -------------------------------------------------------------
   full_load_path = [resdir,'/SPMcfg.mat'];
   [Cond, Vf, RT] = kc_get_design('Load',full_load_path);

   nVf = length(ROI.Ymad_avg);  % number of V-files

	% For each condition, average across PST for ROI.Ymad_avg
	% ----------------------------------------------------------------
	for c = 1:length(Cond)
     for i=1:Cond(c).num_scans
		  tmp_idx = Cond(c).onidx + (i-1);
		  tmp_idx = tmp_idx(find(tmp_idx<=nVf)); % if trials exceed exp end 
        Cond(c).Ymad_avg(i)=mean(ROI.Ymad_avg(tmp_idx));
     end    
	end
  
	
   % Offer Ymad preprocessing, if Design is 'Event' and Extract
	% ----------------------------------------------------------------
	if strcmpi(Design,'Event') & strcmpi(Extract,'Ymad');
	  
	  if ~exist('prepr','var'),
		 prepr = input('Would you like to preprocess the data [y|n] : ', ...
							's');
	  end

	  if strcmpi(prepr,'y');
		 [ROI, Cond] = roi_preprocYmad({'Trim','Detrend','LP','HP'},ROI,Cond,Vf, ...
												 RT);
	  
		 % Average over the preprocessed data
		 % --------------------------------------------------------
		 for c = 1:length(Cond)
        for i=1:Cond(c).num_scans
      	  tmp_idx = Cond(c).onidx + (i-1);
      	  tmp_idx = tmp_idx(find(tmp_idx<=nVf)); % if trials exceed exp end 
           Cond(c).Ypr_avg(i)=mean(ROI.Ypr(tmp_idx));
        end    
       end
	  
	  end
	  
	end
	
	
	
	% Average the fitted response and adjusted data, if requested
	% ------------------------------------------------------------
	if strcmpi(Extract,'Fit') | strcmpi(Extract,'Adj')
      for c = 1:length(Cond)
        for i=1:Cond(c).num_scans
      	  tmp_idx = Cond(c).onidx + (i-1);
      	  tmp_idx = tmp_idx(find(tmp_idx<=nVf)); % if trials exceed exp end 
           Cond(c).Yfit_avg(i)=mean(ROI.Yfit_avg(tmp_idx));
			  Cond(c).Yadj_avg(i)=mean(ROI.Yadj_avg(tmp_idx));
        end    
      end
	end
	
  
%=======================================================================
switch lower(Action)
%======================================================================

 %-----------------------------------------------------------------------
 case 'save'                                    % save ROI.mat file
 %-----------------------------------------------------------------------
     save ROI ROI Cond resdir roi_path;
	  fprintf('\n\nA file ROI.mat has been saved in \n%s.\n',resdir);
	  fprintf('It can be used later for between subject averaging.\n\n');
  
 %-----------------------------------------------------------------------
 case 'plot'                             
 %-----------------------------------------------------------------------
  h=figure;
  tit1=resdir; if length(tit1)>50, tit1=['...',tit1(end-50:end)]; end
  tit1=['Data in : ',tit1];
  tit2=roi_path; 
  if length(tit2)>50,tit2=['for coordinates in : ...',tit2(end-50:end)];end
  title(char(tit1,tit2),'Interpreter','none');

  %================================================
  switch lower(Design) 
  %================================================
  
  
	  case 'event'       % event-related PST
	  % ------------------------------------- 
	  COL   = ['g','b','r','c','y','m','k']; COL = repmat(COL,1,20);
	  figure(h); hold on
	  
		 for c = 1:length(Cond)  
			  
		     switch lower(Extract)   
			  %-------------------------------
				case 'ymad', 
				   if strcmpi(prepr,'y'), data=Cond(c).Ypr_avg;
				   else data=Cond(c).Ymad_avg;
				   end
				
				case 'fit',  data=Cond(c).Yfit_avg;
				case 'adj',  data=Cond(c).Yadj_avg;
			  end
			  
			 
		     pst = 0:RT:40;
			  pst = pst(1:length(Cond(c).Ymad_avg));
           plot(pst,data,[COL(c) '-'],'LineWidth',1)
			  fprintf('\nThe values for %s are',char(Cond(c).name));
			  data
           leg = [];
           for c = 1:length(Cond)
             tmp=char(Cond(c).name); tmp(findstr(tmp,'_'))='-';
             leg=char(leg, tmp);
           end
           leg = leg(2:end,:);
           legend(leg);
		 end
		 
		 xlabel('Peristimulus time (sec)');

		 
		 
	  case 'epoch'       % blocked design bargraph
	  %------------------------------------------
	  	fprintf('\n============================\n');
		fprintf('\nPlease, define the time span to average over, relative');
		fprintf(' the onset of each block:\n');
		
	      for c = 1:length(Cond)
			  fprintf('\nStart point for %s [sec] : ',char(Cond(c).name));
			  Cond(c).start = input('');
           fprintf('End point for %s [sec] : ',char(Cond(c).name));
			  Cond(c).end = input('');
			  if Cond(c).end > 32, 
				 error ('the maximum end point is 32 seconds');
			  end
			  
			  xlabels{c+1} = char(Cond(c).name);
			  tmp_idx = round(Cond(c).start/RT):round(Cond(c).end/RT);
			  tmp_idx = tmp_idx+1;
			  bar_avg(c) =  mean(Cond(c).Ymad_pst_avg(tmp_idx));	
			  fprintf('The value for %s is \t:  %g\n',...
						 char(Cond(c).name),bar_avg(c));
  		   end
			xlabels{end+1} = '';
         
			if length(bar_avg)==1, vrange=bar_avg;
			elseif length(bar_avg)>1, vrange = max(bar_avg)-min(bar_avg);
			end
			
			axis([0 (length(Cond)+1),...
					min(bar_avg) - 0.1*vrange,...
					max(bar_avg) + 0.1*vrange ])
			figure(h); hold on
			bar(bar_avg,'k');
			set(gca,'XTickLabel',xlabels);
			set(gca,'XTick',0:length(Cond)+1);

			xlabel('Condition');
  
  end

   if strcmpi(Extract,'Ymad') & strcmpi(prepr,'y')
	  %ylabel(['% Signal change (Outliers trimmed, Lin detr, HP.01Hz, LP.125Hz)']); 
      ylabel(['% Signal Change']); 
      
	  else
	  ylabel(['Average fMRI intensity (Scaled to Session Grand Mean; No Global' ...
			  ' Scaling)']);
	end
	

end
			
	  
% Return output arguments, if requested
% -------------------------------------
if nargout >  0; varargout{1} = ROI;  end
if nargout == 2; varargout{2} = Cond; end
if nargout == 3; varargout{2} = Vf; end
	
cd(cwd);










