computeROI

 computeROI -- compute an ROI average
 
 ------  Background  ------

 computeROI is the main workhorse for computing ROI averages.  It works for
 both blocked and event-related data (but has not been tested for retinotopic
 or abblocked data) as well as gammafit and FIR analyses.  There are a number
 of options that must be set before computeROI can be run.

 To begin, run "ROIs = computeROI('-init')".  This will initialize a data
 structure that is used in the computation of the ROI average.  To get an idea
 of the possible options, type "ROIs" after running the previous command
 to see the setup of the structure (hence the convention "ROIs", which stands
 for "ROI structure").  These paramters are of course going to be different
 for different experimental designs and analyses.

 ------  ROIs structure  ------

 For all designs, the following elements of the ROIs structure must be
 set.  The majority are self-explanatory:

 funcPath, analysis, maskDir (path to your ROI directory), roi, runlist
 (given as a space delimited list of functional directories), numSlices,
 TER, designtype (either 'event-related' or 'blocked'), nconditions,
 parname (name of your paradigm file), nrows, ncols, timewindow, prestim

 For gammafit designs, you additionally have to set the gammafit element.

 ------  computeROI output  ------

 The output of computeROI depends on your type of analysis and experiment
 design.  For blocked designs only, there will exist a non-empty element
 "ROIs.ROImean" which will give you the raw time-course of your signal over
 your ROI.  Note that this is the _raw_ timecourse, in scanner units, as
 FS-FAST does not save the detrended data as an intermediate step.

 For gammafit analyses the elements "ROIhavg", "ROIhstd", and "ROIhstderr"
 will be set; these refer to the average deviation from baseline over the
 ROI, standard devition, and standard error, respectively.  The values
 used in these computations come from the results of FS-FASTs "selxavg"
 program.

 For FIR analyses the elements "ROImean", "ROIstd", and "ROIstderr" will
 be set; these refer to the mean deviation from baseline, standard
 deviation, and standard error, respectively.  The elements are of
 dimension [numberConditions numberTimepoints].

 For both gammafit and FIR analyses the elements "ROIpct", "ROIpctstd",
 and "ROIpctstderr" will be set, containing the percent signal change
 over the ROI, standard deviation, and standard error.  These values are
 calculated using the results of the FS-FAST "selxavg" program and scaling
 by the value in the "h-offset" volume.  For gammafit analyses, the 
 results will be dim(numberConditions,1), while for FIR analyses they are
 dim(numberConditions, numberTimepoints).

 ------  Usage  ------

 ROIs = computeROI('-init') -- initializes a data structure

 ROIs = computeROI('-process', '-blocked', ROIs) -- performs ROI computations
 using the paramters in the ROIs structure for blocked designs.

 ROIs = computeROI('-process', '-event', ROIs) -- performs ROI computations
 using the paramters in the ROIs structure for event-related designs.

 See the scripts within the FSFAST::ROI module, available from 
 http://froi.sourceforge.net, for examples.

 $Id: computeROI.m,v 1.13 2003/09/25 18:28:00 nknouf Exp $

0001 function [ROIdata] = computeROI(varargin)
0002 % computeROI -- compute an ROI average
0003 %
0004 % ------  Background  ------
0005 %
0006 % computeROI is the main workhorse for computing ROI averages.  It works for
0007 % both blocked and event-related data (but has not been tested for retinotopic
0008 % or abblocked data) as well as gammafit and FIR analyses.  There are a number
0009 % of options that must be set before computeROI can be run.
0010 %
0011 % To begin, run "ROIs = computeROI('-init')".  This will initialize a data
0012 % structure that is used in the computation of the ROI average.  To get an idea
0013 % of the possible options, type "ROIs" after running the previous command
0014 % to see the setup of the structure (hence the convention "ROIs", which stands
0015 % for "ROI structure").  These paramters are of course going to be different
0016 % for different experimental designs and analyses.
0017 %
0018 % ------  ROIs structure  ------
0019 %
0020 % For all designs, the following elements of the ROIs structure must be
0021 % set.  The majority are self-explanatory:
0022 %
0023 % funcPath, analysis, maskDir (path to your ROI directory), roi, runlist
0024 % (given as a space delimited list of functional directories), numSlices,
0025 % TER, designtype (either 'event-related' or 'blocked'), nconditions,
0026 % parname (name of your paradigm file), nrows, ncols, timewindow, prestim
0027 %
0028 % For gammafit designs, you additionally have to set the gammafit element.
0029 %
0030 % ------  computeROI output  ------
0031 %
0032 % The output of computeROI depends on your type of analysis and experiment
0033 % design.  For blocked designs only, there will exist a non-empty element
0034 % "ROIs.ROImean" which will give you the raw time-course of your signal over
0035 % your ROI.  Note that this is the _raw_ timecourse, in scanner units, as
0036 % FS-FAST does not save the detrended data as an intermediate step.
0037 %
0038 % For gammafit analyses the elements "ROIhavg", "ROIhstd", and "ROIhstderr"
0039 % will be set; these refer to the average deviation from baseline over the
0040 % ROI, standard devition, and standard error, respectively.  The values
0041 % used in these computations come from the results of FS-FASTs "selxavg"
0042 % program.
0043 %
0044 % For FIR analyses the elements "ROImean", "ROIstd", and "ROIstderr" will
0045 % be set; these refer to the mean deviation from baseline, standard
0046 % deviation, and standard error, respectively.  The elements are of
0047 % dimension [numberConditions numberTimepoints].
0048 %
0049 % For both gammafit and FIR analyses the elements "ROIpct", "ROIpctstd",
0050 % and "ROIpctstderr" will be set, containing the percent signal change
0051 % over the ROI, standard deviation, and standard error.  These values are
0052 % calculated using the results of the FS-FAST "selxavg" program and scaling
0053 % by the value in the "h-offset" volume.  For gammafit analyses, the
0054 % results will be dim(numberConditions,1), while for FIR analyses they are
0055 % dim(numberConditions, numberTimepoints).
0056 %
0057 % ------  Usage  ------
0058 %
0059 % ROIs = computeROI('-init') -- initializes a data structure
0060 %
0061 % ROIs = computeROI('-process', '-blocked', ROIs) -- performs ROI computations
0062 % using the paramters in the ROIs structure for blocked designs.
0063 %
0064 % ROIs = computeROI('-process', '-event', ROIs) -- performs ROI computations
0065 % using the paramters in the ROIs structure for event-related designs.
0066 %
0067 % See the scripts within the FSFAST::ROI module, available from
0068 % http://froi.sourceforge.net, for examples.
0069 %
0070 % $Id: computeROI.m,v 1.13 2003/09/25 18:28:00 nknouf Exp $
0071 %
0072 
0073 % License:    Perl Artistic License
0074 % History:    9/10/03, Initial Release
0075 
0076 %% if there are no arguments, print usage
0077 if (nargin == 0)
0078     printUsage;
0079     return;
0080 end
0081 
0082 %% if there is only one argument, it means we should initialize the structure
0083 %% or show the help information
0084 if (nargin == 1)
0085     flag = deblank(varargin{1});
0086     if (strncmp(flag, '-init', 5))
0087         ROIdata = computeROI_struct;
0088         return;
0089     elseif (strncmp(flag, '-help', 5))
0090         showHelp;
0091         return;
0092     end
0093 %% if there are three arguments, it means that we should
0094 %% actually compute the ROI average
0095 elseif (nargin == 3)
0096     flag = deblank(varargin{1});
0097     if (strncmp(flag, '-process', 8))
0098         %% if the first second argument is '-event', we do the event-realted
0099         %% analysis
0100         flag = deblank(varargin{2});
0101         if (strncmp(flag, '-event', 6))
0102             ROIdata = computeROI_event(varargin{3});
0103             return;
0104         %% else, we do the blocked analysis
0105         elseif (strncmp(flag, '-blocked', 8))
0106             ROIdata = computeROI_blocked(varargin{3});
0107             return;
0108         end
0109     end
0110 end
0111 
0112 %% this function computes the ROI average per run for a blocked design
0113 function ROIs = computeROI_blocked(ROIs)
0114 
0115     %% setup path for output
0116     outputPath = sprintf('%s/%s_%s', ROIs.maskDir, ROIs.analysis, ROIs.roi);
0117 
0118     %% get number of runs
0119     [junk, numRuns] = size(ROIs.runlist);
0120 
0121     %% get the appropriate time scale
0122     ROIs.time = [-ROIs.prestim:ROIs.TER:(ROIs.timewindow-2*ROIs.prestim)]';
0123 
0124     %% read in the mask
0125     fprintf('compute_roi: reading in the mask named %s...\n', ROIs.roi);
0126     maskDir = sprintf('%s/%s.roi', ROIs.maskDir, ROIs.roi);
0127     mask = textread(maskDir);
0128     %%mask = expandROI(mask, [ROIs.numSlices ROIs.nrows ROIs.ncols]);
0129 
0130     %% read in the hemodynamic volume
0131     fprintf(1, '\ncompute_roi: reading in the hemodynamic volume...\n');
0132     hemoDir = sprintf('%s%s/h', ROIs.funcPath, ROIs.analysis);
0133     %%hvol = bfloat2mat(hemoDir, ROIs.numSlices);
0134     hvol = fmri_ldbvolume(hemoDir);
0135 
0136     %% read in the hemodynamic offset volume
0137     fprintf(1, '\ncompute_roi: reading in the hemodynamic offset volume...\n');
0138     hemoOffsetDir = sprintf('%s%s/h-offset', ROIs.funcPath, ROIs.analysis);
0139     %%hoffsetvol = bfloat2mat(hemoOffsetDir, ROIs.numSlices);
0140     hoffsetvol = fmri_ldbvolume(hemoOffsetDir);
0141 
0142     %% read in h.dat file using fsfast toolbox code
0143     fprintf(1, '\ncompute_roi: reading in hemodynamic data file...\n');
0144     hdatName = [hemoDir '.dat'];
0145     hdat = fmri_lddat2(hdatName);
0146 
0147     %% taken directly from yakview.m
0148     %% reading in the number of degrees of freedom (from design matrix) for use
0149     %% in stderr calculations later
0150     d = diag(hdat.SumXtX);
0151     nndof = d(1:hdat.Nh:length(d))';
0152     dof = [hdat.DOF-sum(nndof) nndof];
0153 
0154     %% because the hemodynamic data is stored as one continuous string of data
0155     %% in the h volumes, we need to split it apart by condition, splitting
0156     %% the means and standard deviations.  This is done with the function
0157     %% convert_hvolume, which was taken almost verbatim from an fsfast program
0158     %% written by doug greve (yak.m?)
0159     fprintf(1, '\ncompute_roi: converting hemodynamic data...\n');
0160     [havg, hstd, numSamples] = convert_hvolume(hvol, ROIs.nconditions);
0161 
0162     %% get the number of voxels in the ROI
0163     numVoxels = size(mask,1);
0164 
0165     %% since this is blocked data, we want to get the average ROI time course for each
0166     %% run, so here we loop through the number of runs
0167     for i=1:numRuns
0168         tempROI = [];
0169         run = ROIs.runlist(1,i);
0170 
0171         %% we need to format the string correctly for reading in the
0172         %% funcitonal data
0173         fprintf('\ncompute_roi: reading in raw data for run %d...\n', run);
0174         if (run < 10)
0175             funcDir = sprintf('%s00%s/%s', ROIs.funcPath, num2str(run), ROIs.funcStem);
0176             output = sprintf('%s_00%s', ROIs.outputStem, num2str(run));
0177         else
0178             funcDir = sprintf('%s0%s/%s', ROIs.funcPath, num2str(run), ROIs.funcStem);
0179             output = sprintf('%s_0%s', ROIs.outputStem, num2str(run));
0180         end
0181         
0182         data = fmri_ldbvolume(funcDir);
0183         fprintf('compute_roi: done reading in raw data for run %d\n', run);
0184         fprintf('compute_roi: calculating intersection of ROI and raw data...\n');
0185 
0186         %% get the data from each voxel in the ROI
0187         for j=1:numVoxels
0188             slice = mask(j,1);
0189             x = mask(j,2);
0190             y = mask(j,3);
0191             value = max((100/abs(hoffsetvol(slice,x,y)))*[squeeze(havg(slice, x, y, :, :))]);
0192             
0193             if ((value <= ROIs.upper_thresh) | (value >= ROIs.lower_thresh))
0194                 tempROI = [tempROI squeeze(data(slice, x, y, :))];
0195             end
0196         end
0197 
0198         %% make sure the ROI data has the correct dimensions
0199             tempROI = squeeze(tempROI);
0200             tempROI = tempROI';
0201 
0202         %% do as the information says
0203             fprintf('compute_roi: calculating mean and standard deviation for run %d...\n', run)
0204             ROIs.ROImean = [ROIs.ROImean; mean(tempROI)];
0205             ROIs.ROIstd = [ROIs.ROIstd; std(tempROI)];
0206 
0207         %% write the data to file
0208         if ROIs.text
0209             filename = sprintf('%s_%s.txt', outputPath, output); 
0210             fid = fopen(filename, 'w', 'b');
0211             fprintf(fid, '%f\n', ROIs.ROImean(i,:)');
0212             fclose(fid);
0213         end
0214     end
0215 
0216     tempROIhavg = [];
0217     tempROIpct = [];
0218     tempROIhavgstd = [];
0219     tempROIpctstd = [];
0220       
0221     %% loop through each voxel in each slices to check if
0222     %% the mask is positive for that particular voxel
0223     %% if it is, add that time-course to the ROI
0224     %% there is probably a more efficient way of doing this,
0225     %% but at the moment, this works fast enough on our server
0226     fprintf('\ncompute_roi: calculating ROI mean and percent signal change from estimated responses...\n');
0227     included = 0;
0228     for i=1:numVoxels
0229         slice = mask(i,1);
0230         x = mask(i,2);
0231         y = mask(i,3);
0232             
0233         value = max((100/abs(hoffsetvol(slice,x,y)))*[squeeze(havg(slice, x, y, :, :))]);
0234 
0235         if ((value > ROIs.upper_thresh) | (value < ROIs.lower_thresh))
0236             ROIs.excludedVoxels = ROIs.excludedVoxels + 1;
0237         elseif ((value <= ROIs.upper_thresh) | (value >= ROIs.lower_thresh))
0238             included = included + 1;
0239             tempROIhavg(included, :, :) = [squeeze(havg(slice, x, y, :, :))];
0240             tempROIpct(included, :, :) = (100/abs(hoffsetvol(slice,x,y)))*[squeeze(havg(slice, x, y, :, :))];  %% divide by the offset to get the percent signal change
0241             tempROIhstd(included, :, :) = [squeeze(hstd(slice, x, y, :, :))];
0242                   tempROIpctstd(included, :, :) = (100/abs(hoffsetvol(slice,x,y)))*[squeeze(hstd(slice, x, y, :, :))];
0243         end
0244     end
0245         
0246         %% record the number of voxels
0247         ROIs.numVoxels = numVoxels;
0248 
0249         %% now compute the data for the h volumes
0250             ROIs.ROIhavg = squeeze(mean(tempROIhavg,1))';
0251             ROIs.ROIpct = squeeze(mean(tempROIpct,1))';
0252 
0253         %% compute the standard deviation
0254         variance = tempROIhstd.^2;
0255         variance = squeeze(mean(variance,1));
0256         variancepct = tempROIpctstd.^2;
0257         variancepct = squeeze(mean(variancepct,1));
0258         
0259         %% let's assume the noise is perfectly uncorrelated, then we divide by the number of voxels in the ROI
0260         %% if we don't want to make this assumption, we would divide by 1 (assuming perfectly correlated noise)
0261         variance = variance/sqrt(numVoxels);
0262         ROIs.ROIhstd = sqrt(variance)';
0263         variancepct = variance/sqrt(numVoxels);
0264         ROIs.ROIpctstd = sqrt(variancepct)';
0265         clear variance, variancepct;
0266 
0267         %% now let's compute the standard error, using a formula from doug greve and fsfast
0268         repeated = repmat(dof, [numSamples 1])';    
0269         ROIs.ROIhstderr = ROIs.ROIhstd./sqrt(abs(ROIs.ROIhstd.*repeated));
0270         ROIs.ROIpctstderr = ROIs.ROIpctstd./sqrt(abs(ROIs.ROIpctstd.*repeated));
0271 
0272     %% finally, lets write some legends
0273     for i = 1:numRuns
0274         runName = sprintf('Run %s', num2str(ROIs.runlist(1,i)));
0275         ROIs.legendRun = {ROIs.legendRun{:} runName};
0276     end
0277 
0278     for i = 1:ROIs.nconditions
0279         condName = sprintf('Condition %s', num2str(i));
0280         ROIs.legendCond = {ROIs.legendCond{:} condName};
0281     end
0282 
0283     %% finally, write everything to file
0284         %% write the data to file
0285         if ROIs.writeROI
0286         ROIhavg = ROIs.ROIhavg;
0287         ROIhstd = ROIs.ROIhstd;
0288         ROIhstderr = ROIs.ROIhstderr;
0289         ROIpct = ROIs.ROIpct;
0290         ROIpctstd = ROIs.ROIpctstd;
0291         ROIpctstderr = ROIs.ROIpctstderr;
0292         if ROIs.text
0293                     filename = sprintf('%s_havg.txt', outputPath); 
0294             command = sprintf('save %s ROIhavg -ASCII', filename);
0295             eval(command);
0296                        filename = sprintf('%s_hstd.txt', outputPath);
0297             command = sprintf('save %s ROIhstd -ASCII', filename);
0298             eval(command);
0299                    filename = sprintf('%s_hstderr.txt', outputPath);
0300             command = sprintf('save %s ROIhstderr -ASCII', filename);
0301             eval(command);
0302                    filename = sprintf('%s_pct.txt', outputPath);
0303             command = sprintf('save %s ROIpct -ASCII', filename);
0304             eval(command);
0305                    filename = sprintf('%s_pctstd.txt', outputPath);
0306             command = sprintf('save %s ROIpctstd -ASCII', filename);
0307             eval(command);
0308                    filename = sprintf('%s_pctstderr.txt', outputPath);
0309             command = sprintf('save %s ROIpctstderr -ASCII', filename);
0310             eval(command);
0311                     filename = sprintf('%s_all.txt', outputPath);
0312             ROImean = ROIs.ROImean;
0313             command = sprintf('save %s ROImean -ASCII', filename);
0314             eval(command);
0315         end
0316                 filename = sprintf('%s_struct', outputPath);
0317         command = sprintf('save %s ROIs', filename);
0318         eval(command);
0319         end
0320 
0321     fprintf('compute_roi: done \n\n');
0322 return;
0323 
0324 %% the following function computes ROIs for event related data
0325 %% its functionality differs whether or not the gammafit vector is non-empty
0326 %% if it is non-empty, then computeROI_event will load the h files in the analysis
0327 %% directory and get only the estimated max responses.  If it is empty, then
0328 %% computeROI_event will be able to get the complete hemodynamic curve over
0329 %% the timewindow
0330 function ROIs = computeROI_event(ROIs);
0331     
0332     %% setup path for output
0333     outputPath = sprintf('%s/%s_%s', ROIs.maskDir, ROIs.analysis, ROIs.roi);
0334 
0335     %% get the number of runs
0336     [junk, numRuns] = size(ROIs.runlist);
0337 
0338     %% get the appropriate time scale
0339     ROIs.time = [-ROIs.prestim:ROIs.TER:(ROIs.timewindow-2*ROIs.prestim)]';
0340 
0341     %% read in the mask
0342     fprintf('compute_roi: reading in the ROI named %s...\n', ROIs.roi);
0343     maskDir = sprintf('%s/%s.roi', ROIs.maskDir, ROIs.roi);
0344     mask = textread(maskDir);
0345     %%mask = expandROI(mask, [ROIs.numSlices ROIs.nrows ROIs.ncols]);
0346 
0347     %% read in the hemodynamic volume
0348     fprintf('\ncompute_roi: reading in the hemodynamic volume...\n');
0349     hemoDir = sprintf('%s%s/h', ROIs.funcPath, ROIs.analysis);
0350     hvol = fmri_ldbvolume(hemoDir);
0351 
0352     %% read in the hemodynamic offset volume
0353     fprintf('\ncompute_roi: reading in the hemodynamic offset volume...\n');
0354     hemoOffsetDir = sprintf('%s%s/h-offset', ROIs.funcPath, ROIs.analysis);
0355     hoffsetvol = fmri_ldbvolume(hemoOffsetDir);
0356 
0357     %% read in h.dat file using fsfast toolbox code
0358     fprintf('\ncompute_roi: reading in hemodynamic data file...\n');
0359     hdatName = [hemoDir '.dat'];
0360     hdat = fmri_lddat2(hdatName);
0361 
0362     %% taken directly from yakview.m
0363     %% this gets us a DOF matrix for use in computing the standard
0364     %% error values later on
0365     d = diag(hdat.SumXtX);
0366     nndof = d(1:hdat.Nh:length(d))';
0367     dof = [hdat.DOF-sum(nndof) nndof];
0368 
0369     %% because the hemodynamic data is stored as one continuous string of data
0370     %% in the h volumes, we need to split it apart by condition, splitting
0371     %% the means and standard deviations.  This is done with the function
0372     %% convert_hvolume, which was taken almost verbatim from an fsfast program
0373     %% written by doug greve (yak.m?)
0374     fprintf(1, '\ncompute_roi: converting hemodynamic data...\n');
0375     [havg, hstd, numSamples] = convert_hvolume(hvol, ROIs.nconditions);
0376 
0377     [numSlices, xsize, ysize, numCond, numSamples] = size(havg);
0378 
0379     %% setup blank variables
0380     tempROIhavg = [];
0381     tempROIpct = [];
0382     tempROIhstd = [];
0383     tempROIhstderr = [];
0384 
0385     fprintf(1, '\ncompute_roi: find intersection of mask and data...\n');
0386     %% loop through each voxel in each slices to check if
0387     %% the mask is positive for that particular voxel
0388     %% if it is, add that time-course to the ROI
0389     %% there is probably a more efficient way of doing this,
0390     %% but at the moment, this works fast enough on our server
0391 
0392     %% we shouldn't need to loop through an entire volume here, simply
0393     %% take the voxels that are below our upper threshold and include them
0394     %% in the calculations
0395     numVoxels = size(mask, 1);
0396     included = 0;
0397     for i=1:numVoxels
0398         slice = mask(i,1);
0399         x = mask(i,2);
0400         y = mask(i,3);
0401             
0402         value = max(max((100/abs(hoffsetvol(slice,x,y)))*[squeeze(havg(slice, x, y, :, :))]));
0403 
0404         if ((value > ROIs.upper_thresh) | (value < ROIs.lower_thresh))
0405             ROIs.excludedVoxels = ROIs.excludedVoxels + 1;
0406         elseif ((value <= ROIs.upper_thresh) | (value >= ROIs.lower_thresh))
0407             included = included + 1;        
0408                     tempROIhavg(included, :, :) = [squeeze(havg(slice, x, y, :, :))];
0409                            tempROIpct(included, :, :) = (100/abs(hoffsetvol(slice,x,y)))*[squeeze(havg(slice, x, y, :, :))];
0410                            tempROIhstd(included, :, :) = [squeeze(hstd(slice, x, y, :, :))];
0411                            tempROIpctstd(included, :, :) = (100/abs(hoffsetvol(slice,x,y)))*[squeeze(hstd(slice, x, y, :, :))];
0412         end
0413     end
0414 
0415     %% record the number of voxels
0416     ROIs.numVoxels = numVoxels;
0417 
0418         fprintf('\ncompute_roi: calculating mean and standard deviation...\n')
0419         ROIs.ROImean = squeeze(mean(tempROIhavg,1))';
0420         ROIs.ROIpct = squeeze(mean(tempROIpct,1))';
0421 
0422     %% compute the standard deviation
0423     variance = tempROIhstd.^2;
0424     variance = squeeze(mean(variance,1));
0425     variancepct = tempROIpctstd.^2;
0426     variancepct = squeeze(mean(variancepct,1));
0427         
0428     %% let's assume the noise is perfectly uncorrelated, then we divide by the number of voxels in the ROI
0429     %% if we don't want to make this assumption, we would divide by 1 (assuming perfectly correlated noise)
0430     variance = variance/sqrt(numVoxels);
0431     ROIs.ROIstd = sqrt(variance)';
0432     variancepct = variance/sqrt(numVoxels);
0433     ROIs.ROIpctstd = sqrt(variancepct)';
0434     clear variance, variancepct;
0435 
0436     %% now let's compute the standard error, using a formula from doug greve and fsfast
0437     repeated = repmat(dof, [numSamples 1]);
0438     ROIs.ROIstderr = ROIs.ROIstd./sqrt(ROIs.ROIstd.*repeated);
0439     ROIs.ROIpctstderr = ROIs.ROIpctstd./sqrt(ROIs.ROIpctstd.*repeated);
0440        
0441     %% finally, lets write some legends
0442     for i = 1:numRuns
0443         runName = sprintf('Run %s', num2str(ROIs.runlist(1,i)));
0444         ROIs.legendRun = {ROIs.legendRun{:} runName};
0445     end
0446 
0447     for i = 1:ROIs.nconditions
0448         condName = sprintf('Condition %s', num2str(i));
0449         ROIs.legendCond = {ROIs.legendCond{:} condName};
0450     end
0451  
0452     %% write the data to file
0453         if ROIs.writeROI
0454         ROImean = ROIs.ROImean;
0455         ROImeanstd = ROIs.ROIstd;
0456         ROImeanstderr = ROIs.ROIstderr;
0457         ROIpct = ROIs.ROIpct;
0458         ROIpctstd = ROIs.ROIpctstd;
0459         ROIpctstderr = ROIs.ROIpctstderr;
0460         if ROIs.text
0461                     filename = sprintf('%s_mean.txt', outputPath);
0462             command = sprintf('save %s ROImean -ASCII', filename);
0463             eval(command);
0464                   filename = sprintf('%s_meanstd.txt', outputPath);
0465             command = sprintf('save %s ROImeanstd -ASCII', filename);
0466             eval(command);
0467                   filename = sprintf('%s_meanstderr.txt', outputPath);
0468             command = sprintf('save %s ROImeanstderr -ASCII', filename);
0469             eval(command);
0470                   filename = sprintf('%s_pct.txt', outputPath);
0471             command = sprintf('save %s ROIpct -ASCII', filename);
0472             eval(command);
0473                    filename = sprintf('%s_pctstd.txt', outputPath);
0474             command = sprintf('save %s ROIpctstd -ASCII', filename);
0475             eval(command);
0476                    filename = sprintf('%s_pctstderr.txt', outputPath);
0477             command = sprintf('save %s ROIpctstderr -ASCII', filename);
0478             eval(command);
0479         end
0480         filename = sprintf('%s_struct', outputPath);
0481         command = sprintf('save %s ROIs', filename);
0482         eval(command);
0483         end
0484 return;
0485 
0486 %% setup a default data structure (called with the '-init' parameter)
0487 %% this data structure is what is modified when setting up the
0488 %% ROI computation and when actually computing the ROI
0489 function ROIs = computeROI_struct
0490     ROIs.funcPath = '';
0491     ROIs.analysis = '';
0492     ROIs.maskDir = '';
0493     ROIs.roi = '';
0494     ROIs.runlist = [];
0495     ROIs.numSlices = 0;
0496     ROIs.segThresh = 0.5;
0497     ROIs.upper_thresh = 1000000;
0498     ROIs.lower_thresh= -1000000;
0499     ROIs.numVoxels = 0;
0500     ROIs.excludedVoxels = 0;
0501     ROIs.timewindow = 16.0;
0502     ROIs.prestim = 4.0;
0503     ROIs.TER = 2.0000;
0504     ROIs.time = [];
0505     ROIs.gammafit = [];
0506     ROIs.rescale = 1000;
0507     ROIs.designtype = 'choose';
0508     ROIs.nconditions = 0;
0509     ROIs.parname = 'para.para';
0510     ROIs.runlistfile = 'runlist';
0511     ROIs.funcStem = 'fmc';
0512     ROIs.writeROI = 1;
0513     ROIs.outputStem = 'fmc';
0514     ROIs.nrows = 64;
0515     ROIs.ncols = 64;
0516     ROIs.legendCond = {'fixation'};
0517     ROIs.legendRun = {};
0518     ROIs.ROImean = [];
0519     ROIs.ROIpct = [];
0520     ROIs.ROIstd = [];
0521     ROIs.ROIstderr = [];
0522     ROIs.ROIpctstd = [];
0523     ROIs.ROIpctstderr = [];
0524     ROIs.ROIhavg = [];
0525     ROIs.ROIhstd = [];
0526     ROIs.ROIhstderr = [];
0527 return; 
0528 
0529 %% print basic usage information when no parameters are given
0530 function printUsage()
0531     fprintf('\nType \"help computeROI\" for detailed information\n\n');
0532     fprintf(1, 'computeROI(''-help'') -- gives detailed help information\n');
0533     fprintf(1, 'ROIs = computeROI(''-init'') -- gives basic data structure needed for further analysis\n');
0534     fprintf(1, 'ROIs = computeROI(''-process'', ''-blocked'', ROIs) -- computes ROI for blocked data\n');
0535     fprintf(1, 'ROIs = computeROI(''-process'', ''-event'', ROIs) -- computes ROI for event-related data\n');
0536 
0537     return;
0538 
0539 function showHelp()
0540     fprintf(1, '$Id: computeROI.m,v 1.13 2003/09/25 18:28:00 nknouf Exp $\n\n');
0541     fprintf(1, 'this will show some help information someday\n');
0542 
0543     return;