function [p, ROI, Cond, G] = roi_percent(varargin)
% FORMAT p = roi_percent(roi,SPMcfg.mat)

% The idea here is to extract from a specified ROI in a
% specified dataset, the percent signal change in the ROI
% for each scan and save that data out in a usable form.
%
% To that end, we first (if the user so chooses)
% preprocess the data temporally from
% the actual uncompressed data images, then extract the basic
% intensity data necessary to compute signal change.
%
% roi: The filename of an .img file specifying an
% ROI to use.
% SPMcfg.mat: a config file containing the design choices for
% the experiment.
%
% p = vector of percent signal change for each scan for this ROI.
%        p also gets saved out as a text and mat file.
%
% ROI = data structure with much info about this ROI, specifically:
%    name        (extension-less name of ROI)
%    img         (filename for ROI .img file)
%    tal         (filename for ROI .tal file)
%    data_imgs (N by 1 vol structs of data images)
%               where N is no of data_imgs
%    XYZ         (3 x n array of ROI voxel locations in voxel coordinates)
%               where n is number of voxels in ROI.
%    XYZmm       (3 x n array of ROI voxel locations in mm coordinates)
%    Yfull_mean      (N by 1 array of means for each data_img)
%    Yfull_var       (N x 1 variances)
%    Yfull_max         (N x 1 max image intensity for each data_img)
%    Yfull_min      (N x 1 min image intensity for each data_img)
%    Ypr_mean       (N x 1 array of pre-processed mean percent signal changes for each data_img)
%    Ypr_mean_collapsed   (M x 1 array of pre-processed mean percent signal changes for each image
%                           averaged across all sessions, 
%                           where M is no. scans/session)
%   Yfull_collapsed    (M x 1 array of mean intensities for each image averaged across all sessions.)
%   Yfull_mean_scaled  (N x 1 array of mean image intensity for each data_img, scaled
%                           by the grand mean intensity session by session.)
%   Yfull_mean_collapsed (same as Ypr_mean_collapsed, but with image intensity)
%   
% 
% Cond - a condition substructure with design information, 
%          of size(Cond) = [1 x #subj]         
% .name   - name of condition
% .iC     - indeces for this conditions in the columns of xX.X
% .iS     - corresponding condition indeces in each Sess
% .iSC    - session specific indeces for this condition
% .onsec  - the onsets for this condition (in sec)
% .onidx  - the onsets for this condition (in scans)
% .num_scans  - (a guess about the) duration of the PST in scans
% .Ypr_avg - the mean percent signal change timeseries of the ROI for that condition
% .Y_avg - the mean intensity timeseries of the ROI for that condition
% .Ypr_std - the standard deviation, at each peristimulus timepoint, for that condition's peristim timecourse
% .Ypr_conint - 1.96 * (Ypr_std / sqrt(n)) - in other words, one-half of the error bar value for this timepoint.
%
%
% G is N x 1 matrix of globals, if required
%
%
% Other major non-obvious data structures used in this include:
%
% Vf   - a V-files substructure with subfields
%     .iS  - V-files indexed by session
%     .G   - global mean activity for the corresponding V-file
% RT   - repetition time (in seconds)



spm2_flag = spm_input('Which SPM version?', 1, 'SPM99 | SPM2', [0 1], 1);
% conveniently, spm decided to change the syntax of
% spm_filter between SPM99 and SPM2.  Backwards
% compatibility is for weaklings!
filter_ver = which('spm_filter');
if ~isempty(strfind(filter_ver, '99'));
    %i.e., we're in SPM99 right now
    run_spm2_flag = 0;             
else
    run_spm2_flag = 1;
end
    

if nargin==0
    files = spm_get(Inf, '*.tal', 'Select ROI .tal files');
else
    files = varargin{1};
end

if nargin < 2
    multisub_flag = spm_input('Single subject or multisubject mode?', 1, 'm', 'Single|Multi', [0;1],1);
    if multisub_flag
        expdir = spm_get(-1, '*', 'Select experiment home directory');

        names = spm_get(-Inf, '*', 'Select subject home directories');
        names = cellstr(names);
        names = deblank(names);
        for i = 1:length(names)
            curr = names{i};
            [dummy_dir, names{i}] = fileparts(curr);
        end

        full_resdir = spm_get(-1, '*', 'Select results directory from first subject');
        [dummy_dir, resdir] = fileparts(full_resdir);
        resdir = deblank(resdir);
        if (exist(fullfile(expdir, names{1}, resdir, 'beta_0001.img')))==0
            % results directory is not at first level of subject dir
            [test_dummy_dir, test_dir] = fileparts(dummy_dir);
            while(~(strcmp(test_dir, names{1})))
                [dummy_dir, newpart] = fileparts(dummy_dir);
                resdir = fullfile(newpart, resdir);
                [test_dummy_dir, test_dir] = fileparts(dummy_dir);
            end
        end
        resdir = deblank(resdir);  %I'm being compulsive; I apologize.
        

        for i = 1:length(names)
            curr_name = names{i};
            curr_res_path = fullfile(expdir, curr_name, resdir);
            if ~spm2_flag  %i.e., SPM99
                if (exist(fullfile(curr_res_path, 'SPMcfg.mat'), 'file'))==0;
                    curr_name
                    curr_cfg_path = spm_get(1, 'SPMcfg.mat', ['Please locate SPMcfg.mat file for ' curr_name '.']);
                else
                    curr_cfg_path = fullfile(curr_res_path, 'SPMcfg.mat');
                end
            else %i.e., SPM2
                if (exist(fullfile(curr_res_path, 'SPM.mat'), 'file'))==0;
                    curr_name
                    curr_cfg_path = spm_get(1, 'SPM.mat', ['Please locate SPM.mat file for ' curr_name '.']);
                else
                    curr_cfg_path = fullfile(curr_res_path, 'SPM.mat');
                end
            end
            cfg_list{i,1} = curr_cfg_path;
        end
        cfg_list = char(cfg_list);
        
    else
        if ~spm2_flag %i.e., SPM99
            cfg_list = spm_get(1, '*SPMcfg.mat', 'Select cfg.mat file');
        else %i.e., SPM2
            cfg_list = spm_get(1, '*SPM.mat', 'Select SPM.mat file');
        end
    end
else
    cfg_list = varargin{2};
end
nsubs = size(cfg_list, 1);



time_baseflag = spm_input('Baseline: steady or time-varying?', '+1', 'b', 'Steady | Vary', [0 1]);

cond_baseflag = 0;

if ~(time_baseflag)
    % Time-varying baseline of condition doesn't make sense... time-varying choice implicitly means varying for all scans
    cond_baseflag = spm_input('Choose baseline: mean of timecourse or condition?', '+1', 'm', 'Mean of timecourse | Mean of condition', [0 1]);
    if cond_baseflag
        if ~spm2_flag %SPM99        
            temp_condlist = kc_get_design('Load', cfg_list(1,:));
            condnames = vertcat(temp_condlist.name);
        else % SPM2
            load(deblank(cfg_list(1,:)));
            condnames = {};
            % accumulate all names first, then pick out
            % unique ones
            for s = 1:length(SPM.Sess)
                condnames = vertcat(condnames, SPM.Sess(s).U.name);
            end
            full_condnames = condnames; 
            [condnames, unique_cond_idx, full_cond_idx] = unique(condnames);
        end
        cond_base_idx = spm_input('Which condition as baseline?', '+1', 'm', condnames, 1:length(condnames));
    end
end

baseregion_list = char({'Same ROI'; 'Grand mean of volume'; 'Mean of user-specified mask'; 'Mean of user-specified inverse mask'});
if time_baseflag
    % Can't have time-varying mean of ROI as baseline - that's what we're extracting!
    baseregion_list(1,:) = [];
end
baseregion_idx = spm_input('Choose region for baseline:', '+1', 'm', baseregion_list);
baseregion = deblank(baseregion_list(baseregion_idx, :));


if strcmp(baseregion, 'Mean of user-specified mask')
    % i.e., we need to choose a mask
    if multisub_flag
        basemask = spm_get(Inf, '*.img', 'Select mask images to extract baseline from, in subject order.');
    else
        basemask = spm_get(1, '*.img', 'Select mask image to extract baseline from:');
    end
elseif strcmp(baseregion, 'Mean of user-specified inverse mask')
    % i.e., we need to choose a negative mask
    % need to choose mask and image to subtract mask from - later on might provide option to supply anatomy and 
    % have program extract itself.
    if multisub_flag
        baseneg = spm_get(length(names), '*.img', 'Select inverse masks to exclude from baseline:');
        basetemplate = spm_get(length(names), '*.img', 'Select original masks to include in baseline:');
    else
        baseneg = spm_get(1, '*.img', 'Select inverse mask to exclude from baseline:');
        basetemplate = spm_get(1, '*.img', 'Select original mask to include in baseline:');
    end
end
    



% -----------------------------------
% Queries for preprocessing
% -----------------------------------

q = spm_input('What preprocessing would you like?',1,'m','All|None|Some',[1;2;3],1);

trim=1;
scale=1;
detr=1;
lowpass=1;
highpass=1;
if q==3
    scale = spm_input('Session-specific scale?','+1','b','Yes|No',[1,0]);
    trim = spm_input('Trim data?','+1','b','Yes|No',[1,0]);
    detr = spm_input('Detrend data by session?','+1','b','Yes|No',[1,0]);
    if ~run_spm2_flag
        lowpass = spm_input('Apply low-pass filter?','+1','b','Yes|No',[1,0]);
    end
    highpass = spm_input('Apply high-pass filter?','+1','b','Yes|No',[1,0]);
else
    if q==2
        trim=0;
        scale=0;
        detr=0;
        lowpass=0;
        highpass=0;
    end
end
% note that we won't even check lowpass down below if we're
% running SPM2.

load(deblank(cfg_list(1,:)));
if ~spm2_flag %SPM99
    nsess = size(Sess,2);
else  %SPM2
    nsess = length(SPM.Sess);
end

if highpass
   	 % use high-pass cutoff from config file
     if ~spm2_flag  %SPM99
        HParam = spm_input('High-pass filter cutoff? (sec)','+1','e',xX.K{1}.HParam,1);
        LFstr = sprintf('[min] Cutoff period %d seconds',min(HParam));
        cLF = 'specify';	
     else
        HParam = spm_input('High-pass filter cutoff? (sec)', '+1', 'e', SPM.xX.K(1).HParam,1);
        LFstr = sprintf('[min Cutoff period %d seconds',min(HParam));
        cLF = 'specify';
     end
else
    HParam = [];
    LFstr = 'none';
    cLF = 'none';
end


if ~run_spm2_flag & lowpass
    % use low-pass cutoff from config file
    low_choice = spm_input('Low-pass filter type?','+1','none|Gaussian|hrf');
    if ~spm2_flag  %SPM99
        def_low_param = xX.K{1}.LParam;
    else % SPM2
        def_low_param = [];
    end
    if strcmp('Gaussian',low_choice)
        PoI_low = spm_input('Gaussian FWHM (secs)?', '+1','r',def_low_param,1);
    	HFstr   = sprintf('Gaussian FWHM %0.1f seconds',PoI_low);
	    LParam  = PoI_low/sqrt(8*log(2));
    else
        LParam = [];
        HFstr = low_choice;
    end
else
    LParam = [];
    low_choice = 'none';
    HFstr = 'none';
end

% ------------------------------------------------
% Query for averaging or one-voxel timeseries
% ------------------------------------------------ 

avg_flag = spm_input('Avg ROI or choose one voxel?',1,'b','Avg|Voxel',[1,0]);

if ~(avg_flag)
    vox_xyz = spm_input('Enter voxel coords (in mm space)','+1','i',[],3);
end

% ------------------------------------------------
% Query for prepending a string to each filename
% ------------------------------------------------ 

if multisub_flag==0
    prependflag = spm_input('Prepend string to output filenames?','+1','b','Yes|No',[1,0]);     
    prepend_list = ' ';
    if prependflag
        prepend_list = spm_input('String to prepend','+1','s');
        prepend_list = [prepend_list '_'];
    end
else
    % if in multi-subject mode, need identifier for each set of ROI files,
    % since they're all going to the same output dir.
    % So we'll use the name of the subject directories
    
    prepend_list = names;
    for i = 1:length(prepend_list)
        prepend_list{i} = [prepend_list{i} '_'];
    end
    prepend_list = char(prepend_list);
end

output_dir = spm_get(-1,'*', 'Select output directory');

%keyboard;
%
% Begin all-encompassing loop
% ------------------------------
%
%keyboard;
for n = 1:nsubs

cfg = deblank(cfg_list(n,:));
prepend = deblank(prepend_list(n,:));
fprintf('\n\nWorking on subject %s', prepend);

% Get design specs for experiment
% -------------------------------
if ~spm2_flag %SPM99
    [Cond, Vf, RT] = kc_get_design('Load',cfg);
    V = Vf.VY;
else
    clear SPM;
    load(cfg);
    [Cond, V, RT] = jc_get_design(SPM);
    % note that Cond will be slightly different from
    % Kalina's here, or slightly clarified - iC and iS are
    % the indices only within the U vector of "actual"
    % conditions, not including things like parametric mods,
    % Volterra kernels, user-spec'ed regressors, etc.
end
    


% If not ROI-based, get timecourse for baseline
% ---------------------------------------------

if strcmp(baseregion, 'Mean of user-specified mask') 
    curr_basemask = deblank(basemask(n, :));
elseif strcmp(baseregion, 'Mean of user-specified inverse mask')
    curr_baseneg = deblank(baseneg(n, :));
    curr_basetemplate = deblank(basetemplate(n, :));
end

switch baseregion
    case 'Grand mean of volume'
        if ~spm2_flag %SPM99
            basic_baseline = Vf.G;
        else %SPM2
            basic_baseline = SPM.xGX.rg;
        end
    case 'Mean of user-specified mask'
        curr_basemask_vol = spm_vol(curr_basemask);
        [base_maskY,base_XYZmm] = spm_read_vols(curr_basemask_vol);
        base_XYZmm = base_XYZmm(:,find(base_maskY==max(max(max(base_maskY)))));

        vinv_base = inv(curr_basemask_vol.mat);  %matrix converting mm->voxels
        base_XYZvox = vinv_base(1:3,:)*[base_XYZmm; ones(1,size(base_XYZmm,2))];
        base_XYZvox = round(base_XYZvox);

        fprintf('\nExtracting %g baseline voxels from %g images\n',length(base_XYZvox),length(V))
        for i = 1:length(V)
            if rem(i,100)==0
                fprintf('\nGot to image %g\n',i)
            end
            basic_baseline(i,1) = mean(spm_sample_vol(V(i), base_XYZvox(1,:), base_XYZvox(2,:), base_XYZvox(3,:),1));   
        end
    case 'Mean of user-specified inverse mask'
        curr_baseneg_vol = spm_vol(curr_baseneg);
        [baseneg_maskY,baseneg_XYZmm] = spm_read_vols(curr_baseneg_vol);
        baseneg_XYZmm = baseneg_XYZmm(:,find(baseneg_maskY==max(max(max(baseneg_maskY)))));

        curr_basetemplate_vol = spm_vol(curr_basetemplate);
        [basetemplate_maskY,basetemplate_XYZmm] = spm_read_vols(curr_basetemplate_vol);
        basetemplate_XYZmm = basetemplate_XYZmm(:,find(basetemplate_maskY==max(max(max(basetemplate_maskY)))));

        base_XYZmm = (setdiff(basetemplate_XYZmm', baseneg_XYZmm', 'rows'))';

        vinv_base = inv(curr_basetemplate_vol.mat);
        base_XYZvox = vinv_base(1:3,:)*[base_XYZmm; ones(1, size(base_XYZmm, 2))];
        base_XYZvox = round(base_XYZvox);

        fprintf('\nExtracting %g baseline voxels from %g images\n',length(base_XYZvox),length(V))
        for i = 1:length(V)
            if rem(i,100)==0
                fprintf('\nGot to image %g\n',i)
            end
            basic_baseline(i,1) = mean(spm_sample_vol(V(i), base_XYZvox(1,:), base_XYZvox(2,:), base_XYZvox(3,:),1));   
        end
    otherwise
        basic_baseline = [];  % will be re-set for each ROI
end


for m = 1:size(files,1)

ROI.file = files(m,:);

[ROI.path,ROI.name,ROI.ext] = fileparts(ROI.file);

ROI.tal = ROI.file;

fprintf('\nWorking on %s\n\n',ROI.file);

    


    % 
    % ------------------------------------------------------------------
    % Assemble ROI structure, other associated vars
    % ------------------------------------------------------------------
    %
%keyboard;    
%v = spm_vol(ROI.space);
%vinv = inv(v.mat); 


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

% This is TERRIBLE style.  I don't even want to talk about it.
% It does the trick, though - allowing one-voxel multi-file mode.

if ~(avg_flag)
    vox_xyz_matrix = repmat(vox_xyz,1,size(ROI.XYZmm,2)); %copies vox_xyz into same size as ROI.XYZ
    if isempty(find(all(vox_xyz_matrix==ROI.XYZmm,1)))
        fprintf('\nCannot find coordinates in current ROI!  Please enter new set.\n');
        ok_flag = 0;
        while ~ok_flag
            vox_xyz = spm_input('Enter voxel coords (in mm space)','+1','i',[],3);
            vox_xyz_matrix = repmat(vox_xyz,1,size(ROI.XYZmm,2)); %copies vox_xyz into same size as ROI.XYZ
            if isempty(find(all(vox_xyz_matrix==ROI.XYZmm,1)))
            fprintf('\nCannot find coordinates in current ROI!  Please enter new set.\n');
            else
                ok_flag = 1;
            end
        end
    end
    %re-set ROI.XYZmm to be, effectively, a one-voxel ROI.
    ROI.XYZmm = vox_xyz;
end  
    

% Convert to ROI.XYZmm to pixel coordinates in ROI.XYZ
% ----------------------------------------------------
if ~spm2_flag %SPM99
    vinv_data = inv(Vf.VY(1).mat);  %matrix converting mm->voxels
else %SPM2
    vinv_data = inv(SPM.xY.VY(1).mat);
end
ROI.XYZ = vinv_data(1:3,:)*[ROI.XYZmm; ones(1,size(ROI.XYZmm,2))];
ROI.XYZ = round(ROI.XYZ);




% Assemble and extract ROI image intensities from images
% -----------------------------------------------------------------

% for each scan image, we need to read in the image and extract the 
% raw image intensity values from the voxels in the ROI

if ~spm2_flag %SPM99
    V = Vf.VY;
else %SPM2
    V = SPM.xY.VY;
end
fprintf('\nRead in vols.\n')

fprintf('\nExtracting %g voxels from %g images\n',size(ROI.XYZ,2),length(V))
for i = 1:length(V)
   if rem(i,100)==0
       fprintf('\nGot to image %g\n',i)
   end
   ROI.Y(i,:) = spm_sample_vol(V(i), ROI.XYZ(1,:), ROI.XYZ(2,:), ROI.XYZ(3,:),1);   
   % leaving trilinear as the default hold, but I don't know why - middle ground, I guess.
   
end


% *****************************


% We now have ROI.Y, a matrix of scaled image intensities in which the columns are co-indexed
% with ROI.XYZ and correspond to the voxel specified at that index in ROI.XYZ.
% The rows of ROI.Y correspond to scans, co-indexed with the 4th dimension of Vol_Y, and the 
% entry at each point in ROI.Y is the image intensity at that voxel for that scan.
% Keep in mind that the Vf.iS structure indexes the scan numbers to sessions, which is nice.

ROI.Yfull_mean = mean(ROI.Y,2);

if avg_flag
    Y = ROI.Yfull_mean;
else
    Y = ROI.Y;
end

scaled_baseline = basic_baseline;

% Y (and ROI_Yfull_mean) is now a column vector in which the entries
% are the global mean intensities of the ROI for each 
% scan in the experiment.
%   We set scaled_baseline equal to basic_baseline here, now that we're
% inside the ROI loop; scaled_baseline will be specific to
% each ROI, as opposed to basic_baseline, which may be
% subject- (but not ROI-specific).  Basic_baseline shouldn't
% be touched inside the ROI loop.


    %
    % -------------------------------------------------------------------
    % Temporal preprocessing
    % -------------------------------------------------------------------
    %
    
 % --------------------------------------------------------------
 % Apply session specific grand-mean scaling, dividing by the 
 % trimmed (%10) grand mean for the entire brain for this session.
 % Necessary for proper interpretation of % signal change, but 
 % we give the option anyways, because, you know, people do strange
 % things sometimes.  We want to scale the baseline as well,
 % so that it doesn't screw up our percent signal
 % calculation later...
 % --------------------------------------------------------------
 
 if scale 
         trimcut = 5; % for the trimmean cut off 10% (5% on each side)
     
    	 for s = 1:nsess
             if ~spm2_flag %SPM99
        		sess_idx = find(Vf.iS==s);
        		trGM = trimmean(Vf.G(sess_idx),trimcut);
            else %SPM2
                sess_idx = SPM.Sess(s).row;
                trGM = trimmean(SPM.xGX.rg(sess_idx), trimcut);
            end
            scaleFactor(s) = 100/trGM;
            Y(sess_idx) = Y(sess_idx)*scaleFactor(s);
            if ~strcmp(baseregion, 'Same ROI')
                scaled_baseline(sess_idx) = scaled_baseline(sess_idx)*scaleFactor(s);
            end
    	 end
      fprintf('Session-specific grand-mean scale complete.\n\n');
  end
      
    


 % ---------------------------------------------
 % Trim outliers by replacing abs(zscore)>2 with 
 % the session specific trimmed mean (10%).  
 % --------------------------------------------
 if trim
 % ---------------------------------------------
	 trimcut = 5; % for the trimmean cut off 10% (5% on each side)
	 prev_idx = 0;
	 for s = 1:nsess
		if ~spm2_flag %SPM99
            sess_idx = find(Vf.iS==s);
        else %SPM2
            sess_idx = SPM.Sess(s).row;
        end
		out_idx = find(abs(zscore(Y(sess_idx))) > 3);
		trM	  = trimmean(Y(sess_idx),trimcut);
		Y(prev_idx+out_idx) = trM;
		prev_idx = sess_idx(end);
	 end
     fprintf('Trim complete.\n\n');
 end
  

% condition flag
% time_baseflag
% Vf.G - global values
% curr_basemask
% curr_baseneg
% curr_basetemplate
% baseregion_list = char({'Same ROI'; 'Grand mean of volume'; 'Mean of user-specified mask'; 'Mean of user-specified inverse mask'});


 % -----------------------------------------------
 % Convert to per cent signal change
 % -----------------------------------------------

% first, get basic_baseline set - should be all set for all non-ROI-based baselines

if strcmp(baseregion, 'Same ROI');
    scaled_baseline = Y;
    % note that this is automatically scaled if we scaled
    % the ROI.
end

% next, find out if we need one number or time-varying
%       - for one number, check if we want mean of whole course or just one condition
%       - for time-varying, need to do array division.


if ~time_baseflag
    % i.e., we have a steady timecourse
    % so check if condition or not
    if cond_baseflag
        cond_base_onsets = Cond(cond_base_idx).onidx;
        if isempty(cond_base_onsets)
            error('Condition has no onsets in this subject.');
        end
        Y_mean_baseline = mean(scaled_baseline(cond_base_onsets));
    else
        Y_mean_baseline = mean(scaled_baseline);
    end
else
    % i.e., we have a time-varying timecourse
    Y_mean_baseline = scaled_baseline;
end

% if Y_mean_baseline is scalar, we can use it in array division, so this covers all our bases. (bases!  ha!)
Y = Y ./ Y_mean_baseline;

Y = Y*100;
Y = Y-100;
% keyboard;

 
        
 % --------------------------------------------------------------
 % Detrend by removing linear trend - separately for each session
 % --------------------------------------------------------------
 if detr
 % --------------------------------------------------------------
	 for s = 1:nsess
        if ~spm2_flag %SPM99
            sess_idx = find(Vf.iS==s);
        else %SPM2
            sess_idx = SPM.Sess(s).row;
        end
		 Y(sess_idx) = detrend(Y(sess_idx));
	 end
     fprintf('Detrending complete.\n\n');
 end
	 


% ----------------------------------------------------------------
% Filtering - SPM-style filtering 
% ----------------------------------------------------------------

if lowpass | highpass
    if ~run_spm2_flag  %running SPM99's version of spm_filter
        for s = 1:nsess
            if ~spm2_flag %SPM99
                K{s} = struct(	'HChoice',	cLF,...
                    'HParam',	HParam,...
                    'LChoice',	low_choice,...
           			'LParam',	LParam,...
           			'row',		find(Vf.iS==s),...
   			        'RT',		RT);
            else %SPM2
                K{s} = struct(	'HChoice',	cLF,...
		            'HParam',	HParam,...
    	            'LChoice',	low_choice,...
                    'LParam',	LParam,...
        	        'row',		SPM.Sess(s).row,...
	                'RT',		RT);
            end
        end
        K = spm_filter('set',K);
    else %running SPM2's version of spm_filter
        for s = 1:nsess
            if ~spm2_flag %SPM99
                K(s) = struct('HChoice', cLF,...
                    'HParam', HParam,...
                    'row', find(Vf.iS==s),...
                    'RT', RT);
            else %SPM2
                K(s) = struct('HChoice', cLF,...
                    'HParam', HParam,...
                    'row', SPM.Sess(s).row,...
                    'RT', RT);
            end
        end
        K = spm_filter(K);
    end 
end

if ~run_spm2_flag & lowpass
    Y = spm_filter('low',K,Y);
    fprintf('Low-pass filtering complete.\n\n');
end
if highpass
    if ~run_spm2_flag %running SPM99 spm_filter
        Y = spm_filter('high',K,Y);
    else % running SPM2 spm_filter
        Y = spm_filter(K, Y);
    end
    fprintf('High-pass filtering complete.\n\n');
end
    
 
            % 
            % ---------------------------------------------
            % Update output structures
            % ---------------------------------------------
            %

% window_length is how long a timecourse to calculate - we
% don't really care how long the condition is, since this is
% just time-locked averaging.  Default is 32 seconds.
window_length = round(32 / RT);

 ROI.Ypr_mean = Y;
  for c = 1:length(Cond)
    if cond_baseflag
        Cond(c).basename = Cond(cond_base_idx).name;
    else
        Cond(c).basename = {'mean'};
    end
	for i=1:window_length
      	tmp_idx = Cond(c).onidx + (i-1);
		tmp_idx = tmp_idx(find(tmp_idx<=length(Y))); % if trials exceed exp end 
        Cond(c).Ypr_avg(i)=mean(ROI.Ypr_mean(tmp_idx));
        Cond(c).Y_avg(i)=mean(ROI.Yfull_mean(tmp_idx));
        Cond(c).Ypr_std(i) = std(ROI.Ypr_mean(tmp_idx));
        Cond(c).Ypr_conint(i) = 1.96*(Cond(c).Ypr_std(i) / sqrt(length(tmp_idx)));
	end
 end 
 

% ROI.Ypr_mean now is a column vector in which the entries are the pre-processed 
% mean percent signal change values for each scan.
% As well, each Condition structure contains its own Ypr_avg vector which 
% consists of the avg ROI % signal change timeseries during that condition across
% all sessions, and a (scaled, but not otherwise processed) Y_avg vector with the mean intensity 
% timeseries.  In order to establish error bars, there are also two vectors the same
% length as Ypr_avg: Ypr_std, which is the standard deviation of the sample of values taken
% at the ith scan after condition c's onset throughout the experiment; and Ypr_conint, which is
% 1.96 * the standard error of that sample, which is provided to aid in plotting error bars.
 
% We can also establish the other elements created by roi_stats, for use in 
% other functions.

% Note that all of these column vectors ARE NOT preprocessed as of now
% EXCEPT for scaling;
% they reports max, min, var, etc. for the scaled intensities.  We might 
% change this in a little while.

%var unfortunately only works along columns...
ROI.Yfull_var = var(ROI.Y');
ROI.Yfull_var = ROI.Yfull_var';

ROI.Yfull_max = max(ROI.Y,[],2);
ROI.Yfull_min = min(ROI.Y,[],2);

% A nice thing to have would be the timeseries for the ROI averaged across
% sessions, rather than just for every individual scan, so we'll compute 
% that as well.  We'll check just to make sure each session has the same 
% number of scans, just in case.  Since each session was scaled individually,
% we won't do this for the image intensities, just the percent signal change.
ROI.Ypr_mean_collapsed = [];
% ROI.Yfull_mean_collapsed = [];
if ~spm2_flag  %SPM99
    ns = size(find(Vf.iS==1),1);
else  %SPM2
    ns = length(SPM.Sess(1).row);
end
skip = 0;
for i = 2:nsess
    if ~spm2_flag  %SPM99
        if (size(find(Vf.iS==i),1)~=ns)
            skip = 1;
            warning('Number of scans per session unequal; skipping between-session averaging.');
        end
    else  %SPM2
        if (length(SPM.Sess(1).row) ~= ns)
            skip = 1;
            warning('Number of scans per session unequal; skipping between-session averaging.');
        end
    end
end
if skip==0
    for i = 1:nsess
        if ~spm2_flag  %SPM99
            sess_vec_pr(:,i) = ROI.Ypr_mean(find(Vf.iS==i));
        else %SPM2
            sess_vec_pr(:,i) = ROI.Ypr_mean(SPM.Sess(i).row);
        end
        %sess_vec(:,i) = ROI.Yfull_mean(find(Vf.iS==i));
    end
    ROI.Ypr_mean_collapsed = mean(sess_vec_pr,2);
    %ROI.Yfull_mean_collapsed = mean(sess_vec,2);
end




%keyboard;
% ---------------------------------
% Save and report these data
% ---------------------------------

p = ROI.Ypr_mean;
p_coll = ROI.Ypr_mean_collapsed;

cdir = pwd;
%cd(ROI.path);


cd(output_dir);
filename = [prepend ROI.name '_percent_signal'];
filename_coll = [filename '_collapsed'];
filename_cond = [filename '_condition'];

     fid = fopen([filename '.txt'],'w');
	     if fid == -1, error(['Unable to create new file - please check permissions', ...
								  ' in current directory.']); 
		  end
   
     fprintf(fid,'%4.4f\n', p);
     fclose(fid);
     
     fid = fopen([filename_coll '.txt'],'w');
	     if fid == -1, error(['Unable to create new file - please check permissions', ...
								  ' in current directory.']); 
		  end
   
     fprintf(fid,'%4.4f\n', p_coll);
     fclose(fid);
     
     fid = fopen([filename_cond '.txt'],'w');
	     if fid == -1, error(['Unable to create new file - please check permissions', ...
								  ' in current directory.']); 
		  end
     fprintf(fid, 'Baseline:  ');
     if cond_baseflag
        fprintf(fid, '%s\n\n', Cond(cond_base_idx).name{1});
     else 
        fprintf(fid, 'mean\n\n');
     end
     for i = 1:length(Cond)
         
         fprintf(fid,'%s\npct:  ',char(Cond(i).name));
         fprintf(fid,'%4.4f  ', Cond(i).Ypr_avg);
         fprintf(fid,'\nstd:  ');
         fprintf(fid,'%4.4f  ', Cond(i).Ypr_std);
         fprintf(fid,'\nconint:  ');
         fprintf(fid,'%4.4f  ', Cond(i).Ypr_conint);
         fprintf(fid,'\nY:  ');
         fprintf(fid,'%4.4f  ', Cond(i).Y_avg);
         fprintf(fid,'\n\n');
     end
     fclose(fid);

save([filename '.mat'], 'p', '-mat');
save([filename_coll '.mat'], 'p_coll', '-mat');
save([filename_cond '.mat'], 'Cond', '-mat');
fprintf('\nThe percent signal change vector was saved as \n%s.mat \nand %s.txt.\n\n',...
    fullfile(pwd,filename), fullfile(pwd,filename));
cd(cdir);
if nargout == 1; varargout{1} = p; end
if nargout == 2; varargout{2} = ROI; end
if nargout == 3; varargout{3} = Cond; end
if nargout > 3
    G(i) = spm_global(V(i)); 
    varargout{4} = G;
end


%
% End all-encompassing loop
% -------------------------
%


clear ROI sess_vec_pr Y V Y_mean_baseline K scaled_baseline
%ROI loop first
end

clear cfg Cond Vf RT 
%then subject loop
end
