00001 %script IndiAna
00002 %
00003 %
00004 %
00005 %Purpose:
00006 % Run 3dDeconvolve in an interactive user-friendly way.
00007 %
00008 %Input:
00009 % Major ones are done through keyboard, while rare ones such as -censor , -nforist, etc. are left for user to implement.
00010 %
00011 %Output:
00012 % Wave functions from waver;
00013 % Various output files from 3dDeconvolve
00014 %
00015 %Key Terms:
00016 %
00017 %More Info :
00018 %
00019 % Author : Gang Chen (with help from Ziad Saad)
00020 % Date : Fri Jul 18 13:58:44 EDT 2003
00021 % SSCC/NIMH/ National Institutes of Health, Bethesda MD 20892
00022
00023 clear all;
00024
00025 FuncName = 'Decon';
00026
00027 fprintf(1,'\nOK, let''s roll ...\n\n');
00028
00029 %parameters common to all regressors
00030
00031 %EPI time seris brick
00032
00033 flg1 = 0; %flag for file existence
00034 while flg1 == 0,
00035 InBrik = input('Input dataset (full filename with suffix BRIK): ','s');
00036 fid = fopen (InBrik,'r');
00037 if (fid == -1),
00038 flg1 = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', InBrik);
00039 else flg1 = 1; fclose (fid);
00040 end
00041 end
00042
00043 %Use Ziad's AfniPrefix.m, which uses RemoveExtension.m
00044 [InBrikPrefix, InBrikView, InBrikExt] = AfniPrefix (InBrik);
00045
00046 flg2 = 0;
00047 %Repition time %TR in seconds
00048 while flg2 == 0,
00049 TR = input('TR (second): ');
00050 if (isnumeric(TR) == 0 | isempty(TR)),
00051 flg2 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00052 else flg2 = 1;
00053 end
00054 end
00055
00056 flg1 = 0; %flag for censor file
00057 while flg1 == 0,
00058 Censor = input('\nNeed to censor any time point(s)? (0 - No; 1 - Yes): ');
00059 if (Censor ~= 1 & Censor ~= 0),
00060 flg1 = 0; fprintf(2,'Error: the input is not a valid number. Please try it again.\n');
00061 else flg1 = 1; end
00062 end
00063
00064 if (Censor == 1),
00065 flg1 = 0;
00066 while flg1 == 0,
00067 Censor_fn = input('Censor file name: ','s');
00068 fid2 = fopen (Censor_fn,'r');
00069 if (fid2 == -1),
00070 flg1 = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', Censor_fn);
00071 else flg1 = 1; fclose (fid2); end
00072 end
00073 end
00074
00075 flg3=0;
00076 %Number of tasks/conditions/regressors to be tested
00077 while flg3 == 0,
00078 N_tasks = input('Number of stimuli/tasks/conditions/regressors: ');
00079 if (isnumeric(N_tasks) == 0 | isempty(N_tasks)),
00080 flg3 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00081 else flg3 = 1;
00082 end
00083 end
00084
00085 flg4 = 0;
00086 %Number of runs concatenated
00087 while flg4 == 0,
00088 N_runs = input('Number of concatenated runs: ');
00089 if (isnumeric(N_runs) == 0 | isempty(N_runs)),
00090 flg4 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00091 else flg4 = 1;
00092 end
00093 end
00094
00095 if (N_runs == 1),
00096 run(1).first = 0;
00097 end
00098
00099 %Get the concatenation points and create a concatenation file
00100 if (N_runs > 1),
00101 fprintf('After concatenation, the 1st run usually starts at 0.\n');
00102 Concat_fn = sprintf('%s_concat.1D', InBrikPrefix);
00103 Concat_id = fopen(Concat_fn,'w');
00104 if (Concat_id < 0),
00105 fprintf(2,'Error in creating concatenation file: Failed to write %s to disk.\n', Concat_fn);
00106 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00107 end
00108 % run(1).first = 0;
00109 for (i=1:1:N_runs),
00110 fprintf('Number %i run starts ', i);
00111 flg5 = 0;
00112 while flg5 == 0,
00113 run(i).first = input('at (TR): ');
00114 if (isnumeric(run(i).first) == 0 | isempty(run(i).first)),
00115 flg5 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00116 else flg5 = 1;
00117 end
00118 end
00119 flg5 = 0;
00120 % fprintf('Number %i run ends ', i);
00121 % run(i).last = input('at (TR): ');
00122 % run(i+1).first = run(i).last + 1;
00123 fprintf(Concat_id, '%i\n', run(i).first);
00124 end
00125 end
00126
00127 %Numer of TR's for the whole dataset including concatenated ones
00128 flg6 = 0;
00129 while flg6 == 0,
00130 N_TR = input('Number of TRs for the whole dataset (run 3dinfo on the dataset if you are not sure): ');
00131 if (isnumeric(N_TR) == 0 | isempty(N_TR)),
00132 flg6 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00133 else flg6 = 1;
00134 end
00135 end
00136
00137 for (i=1:1:N_runs),
00138 if (i<N_runs), run(i).last = run(i+1).first - 1;
00139 else run(i).last = N_TR-1; end
00140 end
00141
00142 flg7 = 0;
00143 while flg7 == 0,
00144 fprintf('\nNormalization would automatically covert the regressor coefficients into percent signal change.');
00145 run_norm = input('\nWant to normalize the dataset before running deconvolution/regression? (0 - No; 1 - Yes): ');
00146 if (run_norm ~= 0 & run_norm ~= 1),
00147 flg7 = 0; fprintf(2,'Error: the input has to be 0 or 1. Please try it again.\n');
00148 else flg7 = 1;
00149 end
00150 end
00151
00152 flg7 = 0;
00153 while flg7 == 0,
00154 fprintf('\n3dAutomask uses 3dClipLevel algorithm to find clipping level.');
00155 run_mask.do = input('\nWant to mask the outside of brain? (0 - No; 1 - Yes): ');
00156 if (run_mask.do ~= 0 & run_mask.do ~= 1),
00157 flg7 = 0; fprintf(2,'Error: the input has to be 0 or 1. Please try it again.\n');
00158 else flg7 = 1;
00159 end
00160 end
00161
00162 if (run_mask.do == 1),
00163 flg7 = 0;
00164 while flg7 == 0,
00165 run_mask.N_dilate = input('How many times do you want to dilate the mask outwards? (default - 5): ');
00166 if isempty(run_mask.N_dilate), run_mask.N_dilate = 5; flg7 = 1;
00167 elseif (isnumeric(run_mask.N_dilate) == 0),
00168 flg7 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00169 else flg7 = 1;
00170 end
00171 end
00172 end
00173
00174 flg8 = 0;
00175 %polort = 2; %Read the 3dDeconvolve manual ...
00176 while flg8 == 0,
00177 polort = input('\nOrder of polynormial fitting for baseline (-1 - no basline; 0 - constant; 1 - linear; 2 - quadratic; ...): ');
00178 if (isnumeric(polort) == 0 | isempty(polort)),
00179 flg8 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00180 else flg8 = 1;
00181 end
00182 end
00183
00184 %Analysis type
00185 fprintf('\nYou can run 2 types of analysis: \n\n');
00186 fprintf('(1) Deconvolution + Regression: generating impulse response functions and condition components; \n');
00187 fprintf('(2) Regression: breaking the signal down into various regressor components. \n');
00188
00189 flg9 = 0;
00190 while flg9 == 0,
00191 Run_Type = input('\nWhat type of analysis? 1 - Deconvolution; 2 - Regression: ');
00192 if (Run_Type ~= 1 & Run_Type ~= 2),
00193 flg9 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00194 else flg9 = 1;
00195 end
00196 end
00197
00198 fprintf('\n');
00199 if (Run_Type == 1), %for deconvolution
00200 N_basis = 1;
00201 for (iT=1:1:N_tasks),
00202 flg10 = 0;
00203 while flg10 == 0,
00204 fprintf('Minimum lags for number %i stimulus ',iT);
00205 % Task(iT).MinLags = input('is: ');
00206 Task(iT).BasisOpt_struct(1).minlag = input('is: ');
00207 % if (isnumeric(Task(iT).MinLags) == 0 | isempty(Task(iT).MinLags)),
00208 if (isnumeric(Task(iT).BasisOpt_struct(1).minlag) == 0 | isempty(Task(iT).BasisOpt_struct(1).minlag)),
00209 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again. \n');
00210 else flg10 = 1;
00211 end
00212 end
00213 flg10 = 0;
00214 while flg10 == 0,
00215 fprintf('Maximum lags for number %i stimulus ',iT);
00216 % Task(iT).MaxLags = input('is: ');
00217 Task(iT).BasisOpt_struct(1).maxlag = input('is: ');
00218 % if (isnumeric(Task(iT).MaxLags) == 0 | isempty(Task(iT).MaxLags)),
00219 if (isnumeric(Task(iT).BasisOpt_struct(1).maxlag) == 0 | isempty(Task(iT).BasisOpt_struct(1).maxlag)),
00220 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again. \n');
00221 else flg10 = 1;
00222 end
00223 end
00224 end
00225 end
00226
00227 if (Run_Type == 2), %for regression
00228 fprintf('\nAvailable basis function options are: \n\n');
00229 fprintf('GAM -- Gamma variate \n');
00230 fprintf('tent -- tent function \n');
00231 fprintf('MGH -- used by Massachusetts General Hospital group \n');
00232 fprintf('SPM1 -- double canonical Gamma variates \n\n');
00233 fprintf('SPM2 -- mixture of Gamma variates with time derivative \n\n');
00234
00235 BasisFunc = input('Basis function type (GAM, tent, MGH, SPM1, SPM2, ...): ', 's');
00236
00237 %There is a scale factor 1/1.25^2 with MGH, but it is OK with waver, and -peak 1 should take care of it.
00238 if (strcmp (lower(BasisFunc),'gam') | strcmp (lower(BasisFunc),'mgh')),
00239
00240 flg10 = 0;
00241 while flg10 == 0,
00242 fprintf('\nThe default power paratmeter is: 8.6 - GAM; 2 - MGH.\n');
00243 BasisOpt.gamb = input('Input the Power parameter if different from the default value, otherwise hit Enter: ');
00244 if isempty(BasisOpt.gamb),
00245 if (strcmp (lower(BasisFunc),'gam')), BasisOpt.gamb = 8.6; end
00246 if (strcmp (lower(BasisFunc),'mgh')), BasisOpt.gamb = 2.0; end
00247 flg10 = 1;
00248 elseif isnumeric(BasisOpt.gamb) == 0,
00249 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00250 else flg10 = 1;
00251 end
00252 end
00253
00254 flg10 = 0;
00255 while flg10 == 0,
00256 fprintf('\nThe default scaling paratmeter is: 0.547 - GAM; 1.25 - MGH.\n');
00257 BasisOpt.gamc = input('Input the scaling parameter if different from the default value, otherwise hit Enter:');
00258 if isempty(BasisOpt.gamc),
00259 if (strcmp (lower(BasisFunc),'gam')), BasisOpt.gamc = 0.547; end
00260 if (strcmp (lower(BasisFunc),'mgh')), BasisOpt.gamc = 1.25; end
00261 flg10 = 1;
00262 elseif isnumeric(BasisOpt.gamc) == 0,
00263 flg10 = 0;
00264 fprintf(2,'Error: the input is not a number. Please try it again.\n');
00265 else flg10 = 1;
00266 end
00267 end
00268
00269 flg10 = 0;
00270 while flg10 == 0,
00271 fprintf('\nThe default delay paratmeter is: 0 - GAM; 2.25 - MGH.\n');
00272 BasisOpt.gamd = input('Input the delay parameter if different from the default value, otherwise hit Enter:');
00273 if isempty(BasisOpt.gamd),
00274 if (strcmp (lower(BasisFunc),'gam')), BasisOpt.gamd = 0.0; end
00275 if (strcmp (lower(BasisFunc),'mgh')), BasisOpt.gamd = 2.25; end
00276 flg10 = 1;
00277 elseif isnumeric(BasisOpt.gamd) == 0,
00278 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00279 else
00280 flg10 = 1;
00281 end
00282 end
00283
00284 flg10 = 0;
00285 fprintf('\nCurrently in each condition minlag and maxlag are set to be the same for all basis functions.\n');
00286 fprintf('The default for both is set to 0. Hit Enter if you want to set them 0.\n');
00287 while flg10 == 0,
00288 BasisOpt.minlag = input('Minimum lags (default: 0): '); %Right now same lags for all regressors and basis functions
00289 if isempty(BasisOpt.minlag),
00290 BasisOpt.minlag = 0;
00291 flg10 = 1;
00292 elseif isnumeric(BasisOpt.minlag) == 0,
00293 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00294 else flg10 = 1;
00295 end
00296 end
00297
00298 flg10 = 0;
00299 while flg10 == 0,
00300 BasisOpt.maxlag = input('Maximum lags (default: 0): '); %Right now same lags for all regressors and basis functions
00301 if isempty(BasisOpt.maxlag),
00302 BasisOpt.minlag = 0;
00303 flg10 = 1;
00304 elseif isnumeric(BasisOpt.maxlag) == 0,
00305 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00306 else flg10 = 1;
00307 end
00308 end
00309
00310 end
00311
00312 flg10 = 0;
00313 while flg10 == 0,
00314 if (strcmp (lower(BasisFunc),'gam') | strcmp (lower(BasisFunc),'mgh') | strcmp (lower(BasisFunc),'spm1')), N_basis = 1;
00315 elseif (strcmp (lower(BasisFunc),'spm2')), N_basis = 2;
00316 else N_basis = input('Number of basis functions: ');
00317 end
00318 if (isnumeric(N_basis) == 0 | isempty(N_basis)),
00319 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00320 else flg10 = 1;
00321 end
00322 end
00323
00324 %Total time span of all basis functions (i.e. duration of IRF) in seconds
00325 BasisOpt.tSpan = 0;
00326 if ((N_basis > 1) & (strcmp (lower(BasisFunc),'tent'))),
00327 flg10 = 0;
00328 while flg10 == 0,
00329 BasisOpt.tSpan = input('Total time span of all basis functions in this regressor (sec): ');
00330 if (isnumeric(BasisOpt.tSpan) == 0 | isempty(BasisOpt.tSpan)),
00331 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00332 else flg10 = 1;
00333 end
00334 end
00335 end
00336
00337 if (N_basis>1),
00338 %Recon_dt = 0.1; %Sampling period for reconstruction IRFs
00339 flg10 = 0;
00340 while flg10 == 0,
00341 Recon_dt = input('Sampling period for reconstruction IRFs, i.e., 0.1, (second): ');
00342 if (isnumeric(Recon_dt) == 0 | isempty(Recon_dt)),
00343 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00344 else flg10 = 1;
00345 end
00346 end
00347 end
00348 end
00349
00350 flg10 = 0;
00351 while flg10 == 0,
00352 fprintf('\nThere are 2 kinds of stimulus file:\n');
00353 fprintf('1. Regular timing indexed with 0s and 1s (TR-locked));\n');
00354 fprintf('2. Irregular timing (not synchronized with TR).\n');
00355 Stim_Type = input('\nType of stimulus time files (1 - Regular; 2 - Irregular): ');
00356 if (Stim_Type ~= 1 & Stim_Type ~= 2),
00357 flg10 = 0; fprintf(2,'Error: the input has to be 0 or 1. Please try it again.\n');
00358 else flg10 = 1;
00359 end
00360 end
00361
00362 if (Stim_Type == 2), %Irregular timing
00363 flg10 = 0;
00364 while flg10 == 0,
00365 fprintf('Specify the precision of your timing files. For example, if time points reach the \n');
00366 fprintf('100th decimal, set it as 0.01 seconds.\n');
00367 WAV_dt = input('Irregular timing precision (second): ');
00368 if (isnumeric(WAV_dt) == 0 | isempty(WAV_dt)),
00369 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00370 else flg10 = 1;
00371 end
00372 end
00373 end
00374
00375 %Obtain stimulus files
00376 for (iT=1:1:N_tasks),
00377 flg10 = 0;
00378 while flg10 == 0,
00379 fprintf('Number %i stimulus timing file ',iT);
00380 Task(iT).StimTime = input('is: ','s');
00381 fid = fopen (Task(iT).StimTime,'r');
00382 if (fid == -1),
00383 flg10 = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', Task(iT).StimTime);
00384 else flg10 = 1; fclose (fid);
00385 end
00386 end
00387 flg10 = 0;
00388 while flg10 == 0,
00389 fprintf('Number %i stimulus label ',iT);
00390 Task(iT).Label = input('is: ','s');
00391 if (ischar(Task(iT).Label) == 0 | isempty(Task(iT).Label)),
00392 flg10 = 0; fprintf(2,'Error: the input is not a string. Please try it again.\n');
00393 else flg10 = 1;
00394 end
00395 end
00396 end
00397
00398 fprintf('\nSometimes covariates are considered in the modeling process. For example, if there were head motion');
00399 fprintf('\nduring the scanning, the variation caused by the motion is better accounted for by adding those');
00400 fprintf('\nparameters from 3dvolreg into the analysis. These regressors will be included as part of the baseline');
00401 fprintf('\nso that the finalthe final full F is interpreted as baseline + pure noise versus baseline + real signal + noise.');
00402
00403 flg10 = 0;
00404 while flg10 == 0,
00405 base = input('\n\nDo you want to add some covariates? (1 - yes; 0 - no) ');
00406 if (base ~= 0 & base ~= 1),
00407 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00408 else flg10 = 1;
00409 end
00410 end
00411
00412 N_base = 0;
00413 if (base == 1),
00414 flg10 = 0;
00415 while flg10 == 0,
00416 N_base = input('\nHow many covariates? ');
00417 if (isnumeric(N_base) == 0 | isempty(N_base)),
00418 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00419 else flg10 = 1;
00420 end
00421 end
00422 if (N_base > 0),
00423 for (i = 1:N_base),
00424 flg10 = 0;
00425 while flg10 == 0,
00426 fprintf('Number %i covariate ', i);
00427 stimbase(i).label = input('name is: ', 's');
00428 if (ischar(stimbase(i).label) == 0 | isempty(stimbase(i).label)),
00429 flg10 = 0; fprintf(2,'Error: the input is not a string. Please try it again.\n');
00430 else flg10 = 1;
00431 end
00432 end
00433 flg10 = 0;
00434 while flg10 == 0,
00435 fprintf('Number %i covariate 1D file (or file with column specified) ', i);
00436 stimbase(i).file = input('expression is: ', 's');
00437 if (ischar(stimbase(i).file) == 0 | isempty(stimbase(i).file)),
00438 flg10 = 0; fprintf(2,'Error: the input is not a string. Please try it again.\n');
00439 else flg10 = 1;
00440 end
00441 end
00442 end % for (i = 1:N_base)
00443 end % if (N_baser > 0)
00444 end % if (base == 1)
00445
00446 %Get information for contrast tests from the user
00447
00448 flg10 = 0;
00449 while flg10 == 0,
00450 contrast = input('\nDo you want to run contrast test? (1 - yes; 0 - no) ');
00451 if (contrast ~= 0 & contrast ~= 1),
00452 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00453 else flg10 = 1;
00454 end
00455 end
00456
00457 if (contrast == 1),
00458 flg10 = 0;
00459 while flg10 == 0,
00460 N_contr = input('How many contrast tests do you want to run? ');
00461 if (isnumeric(N_contr) == 0 | isempty(N_contr)),
00462 flg10 = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00463 else flg10 = 1;
00464 end
00465 end
00466
00467 contr = [];
00468 if (N_contr > 0),
00469 for (i = 1:N_contr),
00470 flg10 = 0;
00471 while flg10 == 0,
00472 fprintf('Number %i contrast ', i);
00473 contr(i).label = input('name is: ', 's');
00474 if (ischar(contr(i).label) == 0 | isempty(contr(i).label)),
00475 flg10 = 0; fprintf(2,'Error: the input is not a string. Please try it again.\n');
00476 else flg10 = 1;
00477 end
00478 end
00479 flg10 = 0;
00480 while flg10 == 0,
00481 fprintf('Number %i contrast ', i);
00482 contr(i).expr = input('expression is: ', 's');
00483 if (ischar(contr(i).expr) == 0 | isempty(contr(i).expr)),
00484 flg10 = 0; fprintf(2,'Error: the input is not a string. Please try it again.\n');
00485 else flg10 = 1;
00486 end
00487 end
00488
00489 if ((contr(i).expr(1) ~= '+') | (contr(i).expr(1) ~= '-')),
00490 contr(i).expr =['+' contr(i).expr]; %add + in front of the first label if it does not have one
00491 end
00492 end
00493 end
00494 end
00495 if (contrast == 0), N_contr = 0; end
00496
00497 %Others = input('Any options not listed above (hit Enter if none): ', 's');
00498
00499 %All input from the user is done here.
00500
00501 FoutName_log = sprintf('%s_CommandLog', InBrikPrefix);
00502 foutid_log = fopen(FoutName_log,'w');
00503 if (foutid_log < 0),
00504 fprintf(2,'Error %s: Failed to write %s to disk.\n',FuncName, FoutName_log);
00505 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00506 end
00507
00508 % Start normalization or masking
00509
00510 if (run_norm == 1 | run_mask.do == 1),
00511 [err, preprc_fn] = Norm(InBrik, N_runs, run, run_norm, run_mask, foutid_log);
00512 if (err),
00513 fprintf(2,'Error %s: Failed in normalization.\n', preprc_fn);
00514 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00515 end
00516 FoutName_log = sprintf('%s_CommandLog', preprc_fn);
00517 else FoutName_log = sprintf('%s_CommandLog', InBrikPrefix);
00518 end
00519
00520 if (contrast == 1), %If the user requests for contrast tests
00521 for (i = 1:1:N_contr),
00522 cnt = 0; %number of stimulus types
00523 El(i).labels = contr(i).label;
00524 El(i).sign = [];
00525 El(i).expr = '';
00526 El(i).strln = [];
00527 strlnth = 0; %string length for each label
00528 for (j = 1:1:length(contr(i).expr)),
00529 switch contr(i).expr(j)
00530 case '+'
00531 cnt = cnt + 1;
00532 El(i).strln = [El(i).strln strlnth]; %put the length of previous stimulus label into this array. The 1st one is 0
00533 El(i).sign = [El(i).sign 1]; %put the sign of this stimulus label
00534 strlnth = 0;
00535 case '-'
00536 cnt = cnt + 1;
00537 El(i).strln = [El(i).strln, strlnth]; %put the length of previous stimulus label into this array. The 1st one is 0
00538 El(i).sign = [El(i).sign -1]; %put the sign of this stimulus label
00539 strlnth = 0;
00540 otherwise if (contr(i).expr(j) ~= ' '), %spaces are excluded
00541 strlnth = strlnth+1;
00542 El(i).expr = sprintf('%s%s',El(i).expr, contr(i).expr(j));
00543 end
00544 end
00545 end
00546 El(i).strln = [El(i).strln strlnth]; %add the last which was not done in the loop
00547 El(i).strln = El(i).strln(2:cnt+1); %remove the first number which is 0
00548 El(i).cnt=cnt;
00549 El(i).order = [];
00550 first = 1;
00551 for (k = 1:1:cnt),
00552 next = first+El(i).strln(k)-1;
00553 for (iT=1:1:N_tasks)
00554 if strcmp(El(i).expr(first:next), Task(iT).Label),
00555 El(i).order = [El(i).order iT];
00556 end
00557 end
00558 first = next+1;
00559 end
00560 %Make sure if there is any misspelling of the stimulus labels.
00561 if (size(El(i).order) ~= cnt),
00562 fprintf(2,'Error: Misspelled stimulus label %s.\n', contr(i).expr);
00563 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00564 end
00565 end
00566 end
00567
00568 %some more setup
00569
00570 fprintf (1,'\n');
00571 if (Run_Type == 2), %for regression
00572 %create regressors using waver
00573 if (strcmp (lower(BasisFunc),'spm1') | strcmp(lower(BasisFunc),'spm2')),
00574 for (iT = 1:1:N_tasks),
00575 fprintf (1,'Creating regressor(s) for task %d...\n', iT);
00576 Task(iT).StimReg = zeros(N_basis, N_TR);
00577 [err, Task(iT).BasisOpt_struct] = WaverBasisOption (BasisFunc, N_basis, BasisOpt);
00578 if (err),
00579 fprintf(2,'Error %s: Failed in WaverBasisOption.\n', FuncName);
00580 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00581 end
00582 % p(1) - delay of response (relative to onset) 6 -gamd in waver -GAM
00583 % p(2) - delay of undershoot (relative to onset) 16
00584 % p(3) - dispersion of response 1
00585 % p(4) - dispersion of undershoot 1 X
00586 % p(5) - ratio of response to undershoot 6 X
00587 % p(6) - onset (seconds) 0 -gamd
00588 % p(7) - length of kernel (seconds) 32
00589
00590 % p = [6 16 1 1 6 0 32]; %default parameters
00591 % [bf1 p] = spm_hrf(TR,p); %get gamma hrf from SPM function spm_hrf
00592 % dp = 1;
00593 % p(6) = p(6) + dp; %for calculating time derivative
00594 % bf2 = (bf1 - spm_hrf(TR,p))/dp;
00595 % D = (bf(:,1) - spm_hrf(TR,p))/dp; %time derivative
00596 % bf = [bf D(:)]; %combine the two
00597 % fid = fopen('%s', Task(iT).Name(iB).str, 'w');
00598 vec = load(Task(iT).StimTime);
00599 Noext = RemoveExtension(Task(iT).StimTime, '.1D|.txt');
00600 Task(iT).Name(1).str = sprintf('%s_resp1.1D', Noext);
00601 if (strcmp(lower(BasisFunc),'spm2')), Task(iT).Name(2).str = sprintf('%s_resp2.1D', Noext); end
00602 if (Stim_Type == 1), %regular timing
00603 [err, bf] = spm_bf (TR);
00604 if (err),
00605 fprintf(2,'Error running spm_bf: Failed in running spm_bf. \n');
00606 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00607 end
00608 reg1 = conv(bf(:,1), vec); %use the convolution function in Matlab
00609 reg2 = conv(bf(:,2), vec);
00610 elseif (Stim_Type == 2), %Irregular timing
00611 [err, bf] = spm_bf (WAV_dt);
00612 if (err),
00613 fprintf(2,'Error running spm_bf: Failed in running spm_bf.\n');
00614 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00615 end
00616 irr_vec = zeros(N_TR*TR/WAV_dt, 1); %Make finer stimulus vector
00617 for (i = 1:1:length(vec)),
00618 irr_vec(round(vec(i)/WAV_dt)+1) = 1; %Only at those moments (i.e., 1.5 sec) it is set to 1.
00619 end
00620 reg_f1 = conv(bf(:,1), irr_vec); %use the convolution function in Matlab
00621 reg_f2 = conv(bf(:,2), irr_vec);
00622
00623 reg1 = reg_f1(1:TR/WAV_dt:N_TR*TR/WAV_dt); %Upsample those at TR's
00624 reg2 = reg_f2(1:TR/WAV_dt:N_TR*TR/WAV_dt);
00625
00626 end
00627 [err1, Task(iT).Name(1).str] = wryte3(reg1(1:N_TR), Task(iT).Name(1).str); %Ziad's routine for writing a matrix into a file
00628 if (strcmp(lower(BasisFunc),'spm2')), [err2, Task(iT).Name(2).str] = wryte3(reg2(1:N_TR), Task(iT).Name(2).str);
00629 else err2 = 0; end
00630 if (err1 | err2),
00631 fprintf(2,'Error %s: failed to write into %s\n',Task(iT).Name(1).str);
00632 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00633 end
00634 Task(iT).StimReg(1,:) = load(Task(iT).Name(1).str); %this is for calculating condition number
00635 if (strcmp(lower(BasisFunc),'spm2')), Task(iT).StimReg(2,:) = load(Task(iT).Name(2).str); end
00636 end
00637
00638 else %for other basis function options
00639 for (iT = 1:1:N_tasks),
00640 fprintf (1,'Creating regressor(s) for task %d...\n', iT);
00641 Task(iT).StimReg = zeros(N_basis, N_TR);
00642 [err, Task(iT).BasisOpt_struct] = WaverBasisOption (BasisFunc, N_basis, BasisOpt);
00643 if (err),
00644 fprintf(2,'Error %s: Failed in WaverBasisOption.\n', FuncName);
00645 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00646 end
00647 for (iB = 1:1:N_basis),
00648 if (err),
00649 fprintf(2,'Error %s:Failed in WaverBasisOption\n', FuncName);
00650 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00651 end
00652 Noext = RemoveExtension(Task(iT).StimTime, '.1D|.txt');
00653 Task(iT).Name(iB).str = sprintf('%s_resp%d.1D', Noext, iB);
00654
00655 %form the command for waver, Probably better to put in a modular function later?
00656 if (Stim_Type == 2), %Irregular timing
00657 if (strcmp (lower(BasisFunc),'gam') | strcmp (lower(BasisFunc),'mgh')),
00658 comm_Waver = sprintf('waver -dt %g %s -gamb %g -gamc %g -gamd %g -peak 1 -numout %d -tstim `cat %s` > %s',...
00659 TR, Task(iT).BasisOpt_struct(iB).opt, Task(iT).BasisOpt_struct.power, Task(iT).BasisOpt_struct.scale, ...
00660 Task(iT).BasisOpt_struct.delay, N_TR, Task(iT).StimTime, Task(iT).Name(iB).str);
00661 elseif (strcmp (lower(BasisFunc),'tent')),
00662 comm_Waver = sprintf('waver -dt %g %s -peak 1 -numout %d -tstim `cat %s` > %s',...
00663 TR, Task(iT).BasisOpt_struct(iB).opt, N_TR, Task(iT).StimTime, Task(iT).Name(iB).str);
00664 end
00665 end
00666 if (Stim_Type == 1), %regular timing
00667 if (strcmp (lower(BasisFunc),'gam') | strcmp (lower(BasisFunc),'mgh')),
00668 comm_Waver = sprintf('waver -dt %g %s -gamb %g -gamc %g -gamd %g -peak 1 -numout %d -input %s > %s',...
00669 TR, Task(iT).BasisOpt_struct(iB).opt, Task(iT).BasisOpt_struct.power, Task(iT).BasisOpt_struct.scale, ...
00670 Task(iT).BasisOpt_struct.delay, N_TR, Task(iT).StimTime, Task(iT).Name(iB).str);
00671 elseif (strcmp (lower(BasisFunc),'tent')),
00672 comm_Waver = sprintf('waver -dt %g %s -peak 1 -numout %d -input %s > %s',...
00673 TR, Task(iT).BasisOpt_struct(iB).opt, N_TR, Task(iT).StimTime, Task(iT).Name(iB).str);
00674 end
00675 end
00676 %first log the waver command
00677 fprintf(2, 'Running: %s\n', comm_Waver);
00678 fprintf(foutid_log, '%s\n', comm_Waver);
00679 [s,w] = system(comm_Waver);
00680 if (s),
00681 fprintf(2,'Error %s: calling waver.\n\n\n%s\n',FuncName, w);
00682 while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
00683 end
00684 Task(iT).StimReg(iB,:) = load(Task(iT).Name(iB).str); %this is for calculating condition number
00685 end
00686 end
00687 end
00688 comm_Waver = ''; %useless from this point on
00689
00690 %Now display the results
00691 ShowRegressors(Task);
00692
00693 Butt = questdlg(sprintf('Check the regressors now ...'),...
00694 'Continue ?', 'Yes', 'No', 'No');
00695 switch Butt,
00696 case 'No',
00697 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00698 case 'Yes',
00699 Continue = 1;
00700 end
00701
00702
00703 %fprintf(1,'Check the regressors now, and hit Enter key if you want to continue ...\n');
00704 %pause;
00705
00706 %Perhaps test the quality of their regressors by calculating the condition number
00707
00708 [err,CondNum] = CompCondNum (N_tasks, N_basis, Task, N_TR, polort);
00709 fprintf(2, '\nCondition Number of design matrix is: %g\n. Continue?',CondNum);
00710 Butt = questdlg(sprintf('Condition number of design matrix is %g', CondNum),...
00711 'Continue ?', 'Yes', 'No', 'No');
00712 switch Butt,
00713 case 'No',
00714 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00715 case 'Yes',
00716 Continue = 1;
00717 end
00718
00719 end
00720
00721 %Now create the 3dDeconvolve thing
00722
00723 Mandatory_opt1 = '';
00724
00725 if (run_norm == 1 | run_mask.do == 1),
00726 Mandatory_opt1 = sprintf('%s -input %s%s.BRIK \\\n -num_stimts %d \\\n -polort %d \\\n',...
00727 Mandatory_opt1, preprc_fn, InBrikView, N_tasks*N_basis+N_base, polort);
00728 else Mandatory_opt1 = sprintf('%s -input %s%s.BRIK \\\n -num_stimts %d \\\n -polort %d \\\n',...
00729 Mandatory_opt1, InBrikPrefix, InBrikView, N_tasks*N_basis+N_base, polort);
00730 end
00731 if (N_runs > 1),
00732 Mandatory_opt1 =sprintf('%s -concat %s \\\n', Mandatory_opt1, Concat_fn);
00733 end
00734 if (Censor == 1),
00735 Mandatory_opt1 =sprintf('%s -censor %s \\\n', Mandatory_opt1, Censor_fn);
00736 end
00737
00738 Mandatory_opt2 = '';
00739 Mandatory_opt3 = '';
00740 Mandatory_opt4 = '';
00741 for (iT=1:1:N_tasks),
00742 if (Run_Type == 2), %for regression
00743 for (iB=1:1:N_basis),
00744 Mandatory_opt2 = sprintf('%s -stim_file %d %s \\\n',...
00745 Mandatory_opt2, (iT-1)*N_basis+iB, Task(iT).Name(iB).str);
00746 Mandatory_opt3 = sprintf('%s -stim_label %d %s \\\n',...
00747 Mandatory_opt3, (iT-1)*N_basis+iB, sprintf('%s_b%d', Task(iT).Label, iB));
00748 Mandatory_opt4 = sprintf('%s -stim_minlag %d %d -stim_maxlag %d %d \\\n',...
00749 Mandatory_opt4, (iT-1)*N_basis+iB, Task(iT).BasisOpt_struct(iB).minlag,...
00750 (iT-1)*N_basis+iB, Task(iT).BasisOpt_struct(iB).maxlag);
00751 end
00752 elseif (Run_Type == 1), %for deconvolution
00753 Mandatory_opt2 = sprintf('%s -stim_file %d %s \\\n',...
00754 Mandatory_opt2, iT, Task(iT).StimTime);
00755 Mandatory_opt3 = sprintf('%s -stim_label %d %s \\\n',...
00756 Mandatory_opt3, iT, sprintf('%s', Task(iT).Label));
00757 Mandatory_opt4 = sprintf('%s -stim_minlag %d %d -stim_maxlag %d %d \\\n',...
00758 Mandatory_opt4, iT, Task(iT).BasisOpt_struct(1).minlag, iT, Task(iT).BasisOpt_struct(1).maxlag);
00759 end
00760 end
00761
00762 if (N_base > 0), % for covariates (with stim_base)
00763 for (ii = N_base),
00764 Mandatory_opt2 = sprintf('%s -stim_file %d %s \\\n',...
00765 Mandatory_opt2, N_tasks*N_basis+ii, stimbase(ii).file);
00766 Mandatory_opt3 = sprintf('%s -stim_base %d -stim_label %d %s \\\n',...
00767 Mandatory_opt3, N_tasks*N_basis+ii, N_tasks*N_basis+ii, stimbase(ii).label);
00768 end
00769 end
00770
00771 %create options for task effect test with genereal linear test
00772 %here we assume the user implicitly wants to test all tasks
00773 %first generate glt files by calling function TaskTest.m
00774
00775 % Only for regression(Run_Type == 2); No need for deconvolution (Run_Type == 1)
00776 % since 3dDeconvolve is supposed to automatically take care of that?
00777
00778 gltopt1 = '';
00779 if (Run_Type == 2), %for regression
00780 if (N_basis ~= 1);
00781 [err,taskfile] = TaskTest(N_runs, polort, N_basis, N_tasks, Task, N_base);
00782 gltopt1 = sprintf('%s -num_glt %i\\\n', gltopt1, N_tasks+N_contr);
00783 for (i = 1:1:N_tasks),
00784 % gltopt1 = sprintf('%s -glt 1 %s -glt_label %i %s\\\n', gltopt1, taskfile(i).name, i, Task(i).Label);
00785 gltopt1 = sprintf('%s -glt %i %s -glt_label %i %s\\\n', gltopt1, N_basis, taskfile(i).name, i, Task(i).Label);
00786 % Here a matrix with N_basis rows is set for testing each task
00787 end
00788 end
00789 end
00790
00791 %create glt options for contrast tests
00792
00793 gltopt2 = '';
00794 if (contrast)
00795 [err,contrfile] = ContrastTest (N_runs, polort, N_basis, Task, N_tasks, N_contr, El, N_base);
00796 if (err),
00797 fprintf(2,'Error %s: Failed in ContrastTest.\n', contrfile);
00798 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00799 end
00800 if (N_basis == 1), gltopt2 = sprintf('%s -num_glt %i\\\n', gltopt2, N_contr); end
00801 % The reason for adding the above line is that with N_bais ~= 1 there has been an option of
00802 % -num_glt already implemented in gltopt1 above to get task effects
00803 if (Run_Type == 1), %only for deconvolution
00804 for (i = 1:1:N_contr),
00805 gltopt2 = sprintf('%s -glt 1 %s -glt_label %i %s\\\n', gltopt2, contrfile(i).name, i, contr(i).label);
00806 end
00807 elseif (Run_Type == 2), %only for regression
00808 if (N_basis ~= 1);
00809 for (i = 1:1:N_contr),
00810 gltopt2 = sprintf('%s -glt 1 %s -glt_label %i %s\\\n', gltopt2, contrfile(i).name, i+N_tasks, contr(i).label);
00811 end
00812 else for (i = 1:1:N_contr),
00813 gltopt2 = sprintf('%s -glt 1 %s -glt_label %i %s\\\n', gltopt2, contrfile(i).name, i, contr(i).label);
00814 end
00815 end
00816 end
00817 end
00818
00819 fitsprefix = sprintf('%s_fits', InBrikPrefix); statprefix = sprintf('%s_stat', InBrikPrefix);
00820 FitOpts = sprintf(' -fitts %s\\\n -fout -tout -full_first -bucket %s\\\n',...
00821 fitsprefix, statprefix);
00822 comm_3dDecon = sprintf ('3dDeconvolve %s %s %s %s %s %s %s',...
00823 Mandatory_opt1, Mandatory_opt2, Mandatory_opt3, Mandatory_opt4, gltopt1, gltopt2, FitOpts);
00824
00825 %log the command
00826 fprintf(foutid_log, '%s\n', comm_3dDecon);
00827
00828 %let's run the 3dDeconvolve command
00829 %check for output prefix existence:
00830
00831 Skip = 0;
00832 if (PrefixStatus (sprintf('%s%s', statprefix, InBrikView)) ~= 1 | PrefixStatus (sprintf('%s%s', fitsprefix, InBrikView)) ~= 1),
00833 Butt = questdlg(sprintf('Output dataset %s%s and/or %s%s exists?', statprefix, InBrikView, fitsprefix, InBrikView),...
00834 'Skip ?', 'Yes', 'No', 'No');
00835 switch Butt,
00836 case 'No',
00837 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00838 case 'Yes',
00839 Skip = 1;
00840 end
00841 end
00842
00843 tic,
00844 if (~Skip),
00845 fprintf (1,'Running: %s\n',comm_3dDecon);
00846 [s,w] = system(comm_3dDecon);
00847 if (s),
00848 fprintf(2,'Error %s: calling 3dDeconvolve.\n\n%s\n\n',FuncName,w);
00849 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00850 end
00851 end
00852 t=toc;
00853 fprintf(1, '3dDeconvolve done in %f minutes...\n\n', t/60);
00854
00855 if (N_basis>1),
00856 %Load the header of the output stats brik
00857
00858 fprintf(1,'Recontructing each hemodynamic response function. ');
00859 [err, InfoStat] = BrikInfo(sprintf('%s%s.HEAD',statprefix, InBrikView));
00860 % [err, Task(iT).indx, SubLabel] = WhichSubBricks(InfoStat, 'Coef', Task(iT).Label);
00861 if (err),
00862 fprintf(2,'Error %s: Failed in BrikInfo.\n', FuncName);
00863 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00864 end
00865
00866 %form the 3dcalc command for each task type
00867 for (iT=1:1:N_tasks),
00868 [err, Task(iT).indx, SubLabel] = WhichSubBricks(InfoStat, 'Coef', Task(iT).Label);
00869 if (err),
00870 fprintf(2,'Error %s: Failed in WhichSubBricks.\n', FuncName);
00871 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00872 end
00873 if (isempty(Task(iT).indx)),
00874 fprintf(2,'Error %s: Failed to find subbricks.\n', FuncName);
00875 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00876 end
00877 if (length(Task(iT).indx) ~= N_basis),
00878 fprintf(2,'Error %s: length(Task(iT).indx) ~= N_basis (%d, %d).\n', FuncName, length(Task(iT).indx), N_basis);
00879 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00880 end
00881
00882 input_str='';
00883 expr_str='';
00884
00885 if (strcmp (lower(BasisFunc),'spm1') | strcmp (lower(BasisFunc),'spm2')),
00886 for (iB=1:1:N_basis),
00887 input_str = sprintf('%s -%c ''%s%s[%d]'' \\\n',...
00888 input_str, 96+iB, statprefix, InBrikView, Task(iT).indx(iB) - 1);
00889 end
00890 for (iB=1:1:N_basis),
00891 input_str = sprintf('%s -%c %s \\\n',...
00892 input_str, 96+N_basis+ iB, Task(iT).Name(i).str);
00893 expr_str = sprintf('%s%c*%s+',...
00894 expr_str, 96+iB, 96+N_basis+ iB);
00895 end
00896
00897 %remove trailing +
00898 expr_str = expr_str(1:length(expr_str)-1);
00899
00900 %create the junk.1D file
00901 % fidout = fopen('junk.1D', 'w');
00902 % for (i=1:1:BasisOpt.tSpan./Recon_dt+1),
00903 % fprintf(fidout,'%g\n', i);
00904 % end
00905 % fclose(fidout);
00906 %augment input_str with other options
00907 input_str = sprintf('%s -dt %g -datum float\\\n', input_str, Recon_dt);
00908
00909
00910 else
00911 for (iB=1:1:N_basis),
00912 input_str = sprintf('%s -%c ''%s%s[%d]'' \\\n',...
00913 input_str, 96+iB, statprefix, InBrikView, Task(iT).indx(iB) - 1);
00914 expr_str = sprintf('%s%c*%s+',...
00915 expr_str, 96+iB,Task(iT).BasisOpt_struct(iB).BareOpt);
00916 end
00917
00918 %remove trailing +
00919 expr_str = expr_str(1:length(expr_str)-1);
00920
00921 %create the junk.1D file
00922 fidout = fopen('junk.1D', 'w');
00923 for (i=1:1:BasisOpt.tSpan./Recon_dt+1),
00924 fprintf(fidout,'%g\n', i);
00925 end
00926 fclose(fidout);
00927 %augment input_str with junk.1D
00928 input_str = sprintf('%s -%c junk.1D \\\n -dt %g -datum float\\\n', input_str, 96+N_basis+1, Recon_dt);
00929 end
00930
00931
00932 %prefix of output for this task
00933 Task(iT).IRFprefix = sprintf('%s_irf', Task(iT).Label);
00934
00935 %for the 3dcalc command
00936 comm_3dcalc = sprintf('3dcalc %s -prefix %s \\\n -expr ''%s''\n',...
00937 input_str, Task(iT).IRFprefix, expr_str);
00938
00939 %execute the command
00940 %add command to log file
00941 fprintf(foutid_log, '%s\n', comm_3dcalc);
00942 fprintf (1,'%s: Now running 3dcalc command (see %s)...\n', FuncName, FoutName_log);
00943 [s,w] = system(comm_3dcalc);
00944 if (s),
00945 fprintf(2,'Error %s: calling 3dcalc.\n\n%s\n\n',FuncName,w);
00946 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00947 end
00948
00949 %for the refit
00950 comm_refit = sprintf('3drefit -epan %s%s', Task(iT).IRFprefix, InBrikView);
00951 %add command to log file
00952 fprintf(foutid_log, '%s\n', comm_refit);
00953 fprintf (1,'%s: Now running 3drefit command (see %s)...\n', FuncName, FoutName_log);
00954 [s,w] = system(comm_refit);
00955 if (s),
00956 fprintf(2,'Error %s: calling 3dcalc.\n\n%s\n\n',FuncName,w);
00957 while (1); fprintf(2,'Halted: Ctrl+c to exit!'); pause; end
00958 end
00959 end
00960 end
00961 %done, close log file
00962 fclose(foutid_log);
00963 %fprintf(1,'I am so glad that you have reached this point. Congratulations! Open AFNI and check the results ...\n');
00964
00965 fprintf(1,'%s: Done.\n\n', FuncName);
