1 #include "mrilib.h"
2
3 static int verb = 0 ;
4
5 /*----------------------------------------------------------------------------*/
6 /* Setup to solve a collection of equations of the form
7
8 [z] = [X] [beta] + gamma [b] + delta [c]
9
10 where [z] = data = N vector
11 [X] = fixed N x (M-1) matrix } [X] and [b] together make
12 [b] = fixed N vector } up the N x M matrix [A]
13 [c] = variable N vector
14 [beta] = (M-1) vector = fit coefficients of columns of [X]
15 gamma = scalar fit coefficient of [b]
16 delta = scalar fit coefficient of [c]
17 The LSS goal is to find the value of gamma+delta for a bunch of
18 different [c] vectors. To that end, this function returns an
19 N vector [s] for each input [c], such that gamma+delta = [s] *dot* [z].
20 All the other stuff that COULD be estimated, such as [beta], is ignored
21 for the sake of efficiency.
22
23 If NULL is returned, the inputs are illegal. In particular, the nx element
24 (column length) of the 2 input images must be the same, or you will be
25 chastised and anathematised in public.
26 *//*--------------------------------------------------------------------------*/
27
LSS_setup(MRI_IMAGE * ima,MRI_IMAGE * imc)28 MRI_IMAGE * LSS_setup( MRI_IMAGE *ima , MRI_IMAGE *imc )
29 {
30 int nn, mm, nc, ii, jj, ic , nwarn=0 ;
31 float *cc, *pp, *qq, *qj, *vv, *ss, *pv , cj, cvdot, pc ;
32 MRI_IMAGE *ims , *imp , *imq ;
33 MRI_IMARR *imar ;
34
35 ENTRY("LSS_setup") ;
36
37 if( ima == NULL || imc == NULL ){ /* bad user */
38 ERROR_message("LSS_setup: NULL input image?!") ;
39 RETURN(NULL) ;
40 }
41
42 /* [A] is nn X mm ; [C] is nn x nc */
43
44 nn = ima->nx ; mm = ima->ny ; nc = imc->ny ; cc = MRI_FLOAT_PTR(imc) ;
45
46 if( imc->nx != nn || nn <= mm+2 ){ /* stoopid user */
47 ERROR_message("LSS_setup: ima->nx=%d does not equal imc->nx=%d :-(" ,
48 ima->nx,imc->nx) ;
49 RETURN(NULL) ;
50 }
51
52 /* get imp = [P] = psinv of [A] = mm X nn matrix
53 imq = [Q] = ortproj onto column null space of [A] = nn X nn matrix */
54
55 mri_matrix_psinv_svd(1) ;
56 imar = mri_matrix_psinv_ortproj( ima , 1 ) ;
57
58 if( imar == NULL ){ /* should not happen */
59 ERROR_message("LSS_setup: cannot compute pseudo-inverse :-(") ;
60 RETURN(NULL) ;
61 }
62
63 imp = IMARR_SUBIM(imar,0) ; pp = MRI_FLOAT_PTR(imp) ;
64 imq = IMARR_SUBIM(imar,1) ; qq = MRI_FLOAT_PTR(imq) ;
65
66 /* create output image = [S] = nn X nc
67 Each column of [S] is the vector that we
68 dot into a data vector [z] to get the estimate of
69 gamma+delta for the corresponding column from [C] */
70
71 ims = mri_new(nn,nc,MRI_float) ; ss = MRI_FLOAT_PTR(ims) ;
72
73 /* workspace vectors */
74
75 vv = (float *)malloc(sizeof(float)*nn) ; /* will be [Q] [c] */
76
77 pv = (float *)malloc(sizeof(float)*nn) ; /* last row of [P] */
78 for( ii=0 ; ii < nn ; ii++ ) pv[ii] = pp[ mm-1 + ii*mm ] ;
79
80 /* loop over columns of [C] (and [S]) */
81
82 for( ic=0 ; ic < nc ; ic++,cc+=nn,ss+=nn ){
83
84 /* compute [v] = [Q] [c] */
85
86 for( ii=0 ; ii < nn ; ii++ ) vv[ii] = 0.0f ; /* initialize [v] to 0 */
87 for( jj=0 ; jj < nn ; jj++ ){ /* loop over columns of Q */
88 qj = qq + jj*nn ; /* ptr to jj-th column of Q */
89 cj = cc[jj] ; /* jj-th value of [c] */
90 for( ii=0 ; ii < nn ; ii++ ) vv[ii] += qj[ii] * cj ; /* sum into [v] */
91 }
92
93 /* compute cvdot = [c] *dot* [v]
94 cj = [c] *dot* [c]
95 pc = [c] *dot* {last row of [P] = pv} */
96
97 for( pc=cj=cvdot=ii=0 ; ii < nn ; ii++ ){
98 cvdot += cc[ii]*vv[ii] ; cj += cc[ii]*cc[ii] ; pc += pv[ii]*cc[ii] ;
99 }
100
101 /* initialize [s] = last row of [P] */
102
103 for( ii=0 ; ii < nn ; ii++ ) ss[ii] = pv[ii] ;
104
105 /* If cvdot is zero(ish), this means that the extra column [c]
106 is collinear(ish) with the columns of [A], and we skip the next step.
107 Note that since [Q] is an orthogonal matrix,
108 we are guaranteed that L2norm([Q][c]) == L2norm([c]),
109 which implies that abs(cvdot) <= abs(cj), by the triangle inequality. */
110
111 if( fabsf(cvdot) >= 1.e-5*cj ){
112
113 /* add the proper fraction of [v] into [s] */
114
115 pc = (1.0f - pc) / cvdot ;
116 for( ii=0 ; ii < nn ; ii++ ) ss[ii] += pc * vv[ii] ;
117 } else {
118 nwarn++ ;
119 }
120
121 } /* end of loop over columns of [C] */
122
123 if( verb && nwarn > 0 )
124 INFO_message("%d (out of %d) LSS individual estimators were collinear",
125 nwarn , nc ) ;
126
127 /* toss the trash and return the output set of columns */
128
129 free(pv) ; free(vv) ; DESTROY_IMARR(imar) ; RETURN(ims) ;
130 }
131
132 /*----------------------------------------------------------------------------*/
133 /* Create 2 new matrix images.
134 (0) One where columns jbot..jtop are excised, and the last column is
135 the sum of those columns.
136 (1) The matrix of the excised columns.
137 *//*--------------------------------------------------------------------------*/
138
LSS_mangle_matrix(MRI_IMAGE * ima,int jbot,int jtop)139 MRI_IMARR * LSS_mangle_matrix( MRI_IMAGE *ima, int jbot, int jtop )
140 {
141 int ii , jj , jnew , jn , njj , nn,mm ;
142 MRI_IMAGE *imb , *imc ; MRI_IMARR *imar ;
143 float *aa , *bb , *cc , *acol , *bcol , *ccol ;
144
145 ENTRY("LSS_mangle_matrix") ;
146
147 if( ima == NULL || ima->kind != MRI_float ) RETURN(NULL) ;
148 nn = ima->nx ; mm = ima->ny ;
149 njj = jtop-jbot+1 ; if( njj <= 1 ) RETURN(NULL) ;
150 if( jbot < 0 || jtop >= mm ) RETURN(NULL) ;
151
152 imb = mri_new( nn , mm-njj+1 , MRI_float ) ; /* matrix without jbot..jtop */
153 imc = mri_new( nn , njj , MRI_float ) ; /* matrix with jbot..jtop */
154 aa = MRI_FLOAT_PTR(ima) ;
155 bb = MRI_FLOAT_PTR(imb) ;
156 cc = MRI_FLOAT_PTR(imc) ;
157
158 /* copy non-excised columns into the new imb */
159
160 for( jn=jj=0 ; jj < mm ; jj++ ){
161 if( jj >= jbot && jj <= jtop ) continue ;
162 acol = aa + jj*nn ; /* jj = old column index */
163 bcol = bb + jn*nn ; jn++ ; /* jn = new column index */
164 for( ii=0 ; ii < nn ; ii++ ) bcol[ii] = acol[ii] ;
165 }
166
167 /* copy excised columns into the new imc,
168 and also add them up into the last column of imb */
169
170 bcol = bb + (mm-njj)*nn ; /* last col of imb */
171 for( jn=0,jj=jbot ; jj <= jtop ; jj++ ){
172 acol = aa + jj*nn ;
173 ccol = cc + jn*nn ; jn++ ;
174 for( ii=0 ; ii < nn ; ii++ ){
175 ccol[ii] = acol[ii] ; bcol[ii] += acol[ii] ;
176 }
177 }
178
179 INIT_IMARR(imar) ; ADDTO_IMARR(imar,imb) ; ADDTO_IMARR(imar,imc) ; RETURN(imar) ;
180 }
181
182 /*----------------------------------------------------------------------------*/
183
LSS_help(void)184 void LSS_help(void)
185 {
186 printf(
187 "Usage: 3dLSS [options]\n"
188 "\n"
189 " ** Least-Squares-Sum (LSS) estimation from a -stim_times_IM matrix, as **\n"
190 " * described in the paper: *\n"
191 " * JA Mumford et al. Deconvolving BOLD activation in event-related *\n"
192 " * designs for multivoxel pattern classification analyses. *\n"
193 " * NeuroImage (2011) http://dx.doi.org/10.1016/j.neuroimage.2011.08.076 *\n"
194 " * LSS regression was first mentioned in this poster: *\n"
195 " * B Turner. A comparison of methods for the use of pattern classification *\n"
196 " * on rapid event-related fMRI data. Annual Meeting of the Society for *\n"
197 " ** Neuroscience, San Diego, CA (2010). **\n"
198 "\n"
199 " The method implemented here can be described (by me) as a 'pull one out'\n"
200 " approach. That is, for a single trial in the list of trials, its individual\n"
201 " regressor is pulled out and kept separate, and all the other trials are\n"
202 " combined to give another regressor - so that if there are N trials, only\n"
203 " 2 regressors (instead of N) are used for the response model. This 'pull out'\n"
204 " approach is repeated for each single trial separately (thus doing N separate\n"
205 " regressions), which gives a separate response amplitude (beta coefficient)\n"
206 " for each trial. See the 'Caveats' section below for more information.\n"
207 "\n"
208 "----------------------------------------\n"
209 "Options (the first 'option' is mandatory)\n"
210 "----------------------------------------\n"
211 " -matrix mmm = Read the matrix 'mmm', which should have been\n"
212 " output from 3dDeconvolve via the '-x1D' option, and\n"
213 " should have included exactly one '-stim_times_IM' option.\n"
214 " -->> The 3dLSS algorithm requires that at least 2 different\n"
215 " stimulus times be given in the -stim_times_IM option.\n"
216 " If you have only 1 stim time, this program will not run.\n"
217 " In such a case, the normal '-bucket' output from 3dDeconvolve\n"
218 " (or '-Rbuck' output from 3dREMLfit) will have the single\n"
219 " beta for the single stim time.\n"
220 "\n"
221 " -input ddd = Read time series dataset 'ddd'\n"
222 " ** OR **\n"
223 " -nodata = Just compute the estimator matrix -- to be saved with '-save1D'.\n"
224 " * The number of time points is taken from the matrix header.\n"
225 " * If neither '-input' nor '-nodata' is given, '-nodata' is used.\n"
226 " * If '-input' is used, the number of time points in the dataset\n"
227 " must match the number of time points in the matrix.\n"
228 "\n"
229 " -mask MMM = Read dataset 'MMM' as a mask for the input; voxels outside\n"
230 " the mask will not be fit by the regression model.\n"
231 " -automask = If you don't know what this does by now, please don't use\n"
232 " this program.\n"
233 " * Neither of these options has any meaning for '-nodata'.\n"
234 " * If '-input' is used and neither of these options is given,\n"
235 " then all voxels will be processed.\n"
236 "\n"
237 " -prefix ppp = Prefix name for the output dataset;\n"
238 " this dataset will contain ONLY the LSS estimates of the\n"
239 " beta weights for the '-stim_times_IM' stimuli.\n"
240 " * If you don't use '-prefix', then the prefix is 'LSSout'.\n"
241 "\n"
242 " -save1D qqq = Save the estimator vectors (cf. infra) to a 1D formatted\n"
243 " file named 'qqq'. Each column of this file will be\n"
244 " one estimator vector, the same length as the input\n"
245 " dataset timeseries (after censoring, if any).\n"
246 " * The j-th LSS beta estimate is the dot product of the j-th\n"
247 " column of this file with the data time series (duly censored).\n"
248 " * If you don't use '-save1D', then this file is not saved.\n"
249 "\n"
250 " -verb = Write out progress reports, for fun fun fun in the sun sun sun.\n"
251 "\n"
252 "-------------------\n"
253 "Method == EQUATIONS\n"
254 "-------------------\n"
255 " 3dLSS is fast, since it uses a rank-1 bordering technique to pre-compute\n"
256 " the estimator for each separate stimulus regressor from the fixed part of\n"
257 " the matrix, then applies these estimators to each time series in the input\n"
258 " dataset by a simple dot product. If you wish to peruse the equations, see\n"
259 " https://afni.nimh.nih.gov/pub/dist/doc/misc/3dLSS/3dLSS_mathnotes.pdf \n"
260 " The estimator for each separate beta (as described at '-save1D') is the\n"
261 " N-vector which, when dotted into the N-vector of a voxel's time series,\n"
262 " gives the LSS beta estimate for that voxel.\n"
263 "\n"
264 "---------------------\n"
265 "Caveats == READ THIS!\n"
266 "---------------------\n"
267 " The LSS method produces estimates that tend to have smaller variance than the\n"
268 " LSA method that 3dDeconvolve would produce, but the LSS estimates have greater\n"
269 " bias -- in principle, the LSA method is unbiased if the noise is symmetrically\n"
270 " distributed. For the purpose of using the beta estimates for MVPA (e.g., 3dsvm),\n"
271 " the bias may not matter much and the variance reduction may help improve the\n"
272 " classification, as illustrated in the Mumford paper. For other purposes, the\n"
273 " trade-off might well go the other way -- for ANY application of LSS vs. LSA,\n"
274 " you need to assess the situation before deciding -- probably by the judicious\n"
275 " use of simulation (as in the Mumford paper).\n"
276 "\n"
277 " The bias in the estimate of any given beta is essentially due to the fact\n"
278 " that for any given beta, LSS doesn't use an estimator vector that is orthogonal\n"
279 " to the regressors for other coefficients -- that is what LSA does, using the\n"
280 " pseudo-inverse. Typically, any given LSS-estimated beta will include a mixture\n"
281 " of the betas from neighboring stimuli -- for example,\n"
282 " beta8{LSS} = beta8{LSA} + 0.3*beta7{LSA} - 0.1*beta9{LSA} + smaller stuff\n"
283 " where the weights of the neighbors are larger if the corresponding stimuli\n"
284 " are closer (so the regressors overlap more).\n"
285 "\n"
286 " The LSS betas are NOT biased by including any betas that aren't from the\n"
287 " -stim_times_IM regressors -- the LSS estimator vectors (what '-save1D' gives)\n"
288 " are orthogonal to those 'nuisance' regression matrix columns.\n"
289 "\n"
290 " To investigate these weighting and orthogonality issues yourself, you can\n"
291 " multiply the LSS estimator vectors into the 3dDeconvolve regression matrix\n"
292 " and examine the result -- in the ideal world, the matrix would be all 0\n"
293 " except for 1s on diagonal corresponding to the -stim_times_IM betas. This\n"
294 " calculation can be done in AFNI with commands something like the 'toy' example\n"
295 " below, which has only 6 stimulus times:\n"
296 "\n"
297 " 3dDeconvolve -nodata 50 1.0 -polort 1 -x1D R.xmat.1D -x1D_stop -num_stimts 1 \\\n"
298 " -stim_times_IM 1 '1D: 12.7 16.6 20.1 26.9 30.5 36.5' 'BLOCK(0.5,1)'\n"
299 " 3dLSS -verb -nodata -matrix R.xmat.1D -save1D R.LSS.1D\n"
300 " 1dmatcalc '&read(R.xmat.1D) &transp &read(R.LSS.1D) &mult &write(R.mult.1D)'\n"
301 " 1dplot R.mult.1D &\n"
302 " 1dgrayplot R.mult.1D &\n"
303 "\n"
304 " * 3dDeconvolve is used to setup the matrix into file R.xmat.1D\n"
305 " * 3dLSS is used to compute the LSS estimator vectors into file R.LSS.1D\n"
306 " * 1dmatcalc is used to multiply the '-save1D' matrix into the regression matrix:\n"
307 " [R.mult.1D] = [R.xmat.1D]' [R.LSS.1D]\n"
308 " where [x] = matrix made from columns of numbers in file x, and ' = transpose.\n"
309 " * 1dplot and 1dgrayplot are used to display the results.\n"
310 " * The j-th column in the R.mult.1D file is the set of weights of the true betas\n"
311 " that influence the estimated j-th LSS beta.\n"
312 " * e.g., Note that the 4th and 5th stimuli are close in time (3.6 s), and that\n"
313 " the result is that the LSS estimator for the 4th and 5th beta weights mix up\n"
314 " the 'true' 4th, 5th, and 6th betas. For example, looking at the 4th column\n"
315 " of R.mult.1D, we see that\n"
316 " beta4{LSS} = beta4{LSA} + 0.33*beta5{LSA} - 0.27*beta6{LSA} + small stuff\n"
317 " * The sum of each column of R.mult.1D is 1 (e.g., run '1dsum R.mult.1D'),\n"
318 " and the diagonal elements are also 1, showing that the j-th LSS beta is\n"
319 " equal to the j-th LSA beta plus a weighted sum of the other LSA betas, where\n"
320 " those other weights add up to zero.\n"
321 "\n"
322 "--------------------------------------------------------------------------\n"
323 "-- RWCox - Dec 2011 - Compute fast, abend early, leave a pretty dataset --\n"
324 "--------------------------------------------------------------------------\n"
325 ) ;
326 PRINT_COMPILE_DATE ; exit(0) ;
327 }
328
329 /*----------------------------------------------------------------------------*/
330
main(int argc,char * argv[])331 int main( int argc , char *argv[] )
332 {
333 int iarg , nerr=0 , nvals,nvox , nx,ny,nz , ii,jj,kk ;
334 char *prefix="LSSout" , *save1D=NULL , nbuf[256] ;
335 THD_3dim_dataset *inset=NULL , *outset ;
336 MRI_vectim *inset_mrv=NULL ;
337 byte *mask=NULL ; int mask_nx=0,mask_ny=0,mask_nz=0, automask=0, nmask=0 ;
338 NI_element *nelmat=NULL ; char *matname=NULL ;
339 char *cgl ;
340 int Ngoodlist,*goodlist=NULL , nfull , ncmat,ntime ;
341 NI_int_array *giar ; NI_str_array *gsar ; NI_float_array *gfar ;
342 MRI_IMAGE *imX, *imA, *imC, *imS ; float *Xar, *Sar ; MRI_IMARR *imar ;
343 int nS ; float *ss , *oo , *fv , sum ; int nvec , iv ;
344 int nbstim , nst=0 , jst_bot,jst_top ; char *stlab="LSS" ;
345 int nodata=1 ;
346
347 /*--- help me if you can ---*/
348
349 if( argc < 2 || strcasecmp(argv[1],"-HELP") == 0 ) LSS_help() ;
350
351 /*--- bureaucratic startup ---*/
352
353 PRINT_VERSION("3dLSS"); mainENTRY("3dLSS main"); machdep();
354 AFNI_logger("3dLSS",argc,argv); AUTHOR("RWCox");
355 (void)COX_clock_time() ;
356
357 /**------- scan command line --------**/
358
359 iarg = 1 ;
360 while( iarg < argc ){
361
362 if( strcmp(argv[iarg],"-verb") == 0 ){ verb++ ; iarg++ ; continue ; }
363 if( strcmp(argv[iarg],"-VERB") == 0 ){ verb+=2 ; iarg++ ; continue ; }
364
365 /**========== -mask ==========**/
366
367 if( strcasecmp(argv[iarg],"-mask") == 0 ){
368 THD_3dim_dataset *mset ;
369 if( ++iarg >= argc ) ERROR_exit("Need argument after '-mask'") ;
370 if( mask != NULL || automask ) ERROR_exit("Can't have two mask inputs") ;
371 mset = THD_open_dataset( argv[iarg] ) ;
372 CHECK_OPEN_ERROR(mset,argv[iarg]) ;
373 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
374 mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);
375 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
376 if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%.33s'",argv[iarg]) ;
377 nmask = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;
378 if( verb || nmask < 1 ) INFO_message("Number of voxels in mask = %d",nmask) ;
379 if( nmask < 1 ) ERROR_exit("Mask is too small to process") ;
380 iarg++ ; continue ;
381 }
382
383 if( strcasecmp(argv[iarg],"-automask") == 0 ){
384 if( mask != NULL ) ERROR_exit("Can't have -automask and -mask") ;
385 automask = 1 ; iarg++ ; continue ;
386 }
387
388 /**========== -matrix ==========**/
389
390 if( strcasecmp(argv[iarg],"-matrix") == 0 ){
391 if( ++iarg >= argc ) ERROR_exit("Need argument after '%s'",argv[iarg-1]) ;
392 if( nelmat != NULL ) ERROR_exit("More than 1 -matrix option!?");
393 nelmat = NI_read_element_fromfile( argv[iarg] ) ; /* read NIML file */
394 matname = argv[iarg];
395 if( nelmat == NULL || nelmat->type != NI_ELEMENT_TYPE )
396 ERROR_exit("Can't process -matrix file '%s'!?",matname) ;
397 iarg++ ; continue ;
398 }
399
400 /**========== -nodata ===========**/
401
402 if( strcasecmp(argv[iarg],"-nodata") == 0 ){
403 nodata = 1 ; iarg++ ; continue ;
404 }
405
406 /**========== -input ==========**/
407
408 if( strcasecmp(argv[iarg],"-input") == 0 ){
409 if( inset != NULL ) ERROR_exit("Can't have two -input options!?") ;
410 if( ++iarg >= argc ) ERROR_exit("Need argument after '%s'",argv[iarg-1]) ;
411 inset = THD_open_dataset( argv[iarg] ) ;
412 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
413 nodata = 0 ; iarg++ ; continue ;
414 }
415
416 /**========== -prefix =========**/
417
418 if( strcasecmp(argv[iarg],"-prefix") == 0 ){
419 if( ++iarg >= argc ) ERROR_exit("Need argument after '%s'",argv[iarg-1]) ;
420 prefix = strdup(argv[iarg]) ;
421 if( !THD_filename_ok(prefix) ) ERROR_exit("Illegal string after %s",argv[iarg-1]) ;
422 if( verb && strcmp(prefix,"NULL") == 0 )
423 INFO_message("-prefix NULL ==> no dataset will be written") ;
424 iarg++ ; continue ;
425 }
426
427 /**========== -save1D =========**/
428
429 if( strcasecmp(argv[iarg],"-save1D") == 0 ){
430 if( ++iarg >= argc ) ERROR_exit("Need argument after '%s'",argv[iarg-1]) ;
431 save1D = strdup(argv[iarg]) ;
432 if( !THD_filename_ok(save1D) ) ERROR_exit("Illegal string after %s",argv[iarg-1]) ;
433 iarg++ ; continue ;
434 }
435
436 /***** Loser User *****/
437
438 ERROR_message("Unknown option: %s",argv[iarg]) ;
439 suggest_best_prog_option(argv[0], argv[iarg]);
440 exit(1);
441
442 } /* end of loop over options */
443
444 /*----- check for errors -----*/
445
446 if( nelmat == NULL ){ ERROR_message("No -matrix option!?") ; nerr++ ; }
447 if( nerr > 0 ) ERROR_exit("Can't continue without these inputs!") ;
448
449 if( inset != NULL ){
450 nvals = DSET_NVALS(inset) ; nvox = DSET_NVOX(inset) ;
451 nx = DSET_NX(inset) ; ny = DSET_NY(inset) ; nz = DSET_NZ(inset) ;
452 } else {
453 automask = nvals = 0 ;
454 nvox = nx = ny = nz = nodata = 1 ; /* nodata */
455 mask = NULL ;
456 }
457
458 /*----- masque -----*/
459
460 if( mask != NULL ){ /* check -mask option for compatibility */
461 if( mask_nx != nx || mask_ny != ny || mask_nz != nz )
462 ERROR_exit("-mask dataset grid doesn't match input dataset :-(") ;
463
464 } else if( automask ){ /* create a mask from input dataset */
465 mask = THD_automask( inset ) ;
466 if( mask == NULL )
467 ERROR_message("Can't create -automask from input dataset :-(") ;
468 nmask = THD_countmask( nvox , mask ) ;
469 if( verb || nmask < 1 )
470 INFO_message("Number of voxels in automask = %d (out of %d = %.1f%%)",
471 nmask, nvox, (100.0f*nmask)/nvox ) ;
472 if( nmask < 1 ) ERROR_exit("Automask is too small to process") ;
473
474 } else if( !nodata ) { /* create a 'mask' for all voxels */
475 if( verb )
476 INFO_message("No mask ==> computing for all %d voxels",nvox) ;
477 mask = (byte *)malloc(sizeof(byte)*nvox) ; nmask = nvox ;
478 memset( mask , 1 , sizeof(byte)*nvox ) ;
479
480 }
481
482 /*----- get matrix info from the NIML element -----*/
483
484 if( verb ) INFO_message("extracting matrix info") ;
485
486 ncmat = nelmat->vec_num ; /* number of columns */
487 ntime = nelmat->vec_len ; /* number of rows */
488 if( ntime < ncmat+2 )
489 ERROR_exit("Matrix has too many columns (%d) for number of rows (%d)",ncmat,ntime) ;
490
491 /*--- number of rows in the full matrix (without censoring) ---*/
492
493 cgl = NI_get_attribute( nelmat , "NRowFull" ) ;
494 if( cgl == NULL ) ERROR_exit("Matrix is missing 'NRowFull' attribute!") ;
495 nfull = (int)strtod(cgl,NULL) ;
496 if( nodata ){
497 nvals = nfull ;
498 } else if( nvals != nfull ){
499 ERROR_exit("-input dataset has %d time points, but matrix indicates %d",
500 nvals , nfull ) ;
501 }
502
503 /*--- the goodlist = mapping from matrix row index to time index
504 (which allows for possible time point censoring) ---*/
505
506 cgl = NI_get_attribute( nelmat , "GoodList" ) ;
507 if( cgl == NULL ) ERROR_exit("Matrix is missing 'GoodList' attribute!") ;
508 giar = NI_decode_int_list( cgl , ";," ) ;
509 if( giar == NULL || giar->num < ntime )
510 ERROR_exit("Matrix 'GoodList' badly formatted?!") ;
511 Ngoodlist = giar->num ; goodlist = giar->ar ;
512 if( Ngoodlist != ntime )
513 ERROR_exit("Matrix 'GoodList' incorrect length?!") ;
514 else if( verb > 1 && Ngoodlist < nfull )
515 ININFO_message("censoring reduces time series length from %d to %d",nfull,Ngoodlist) ;
516
517 /*--- extract the matrix from the NIML element ---*/
518
519 imX = mri_new( ntime , ncmat , MRI_float ) ;
520 Xar = MRI_FLOAT_PTR(imX) ;
521
522 if( nelmat->vec_typ[0] == NI_FLOAT ){ /* from 3dDeconvolve_f */
523 float *cd ;
524 for( jj=0 ; jj < ncmat ; jj++ ){
525 cd = (float *)nelmat->vec[jj] ;
526 for( ii=0 ; ii < ntime ; ii++ ) Xar[ii+jj*ntime] = cd[ii] ;
527 }
528 } else if( nelmat->vec_typ[0] == NI_DOUBLE ){ /* from 3dDeconvolve */
529 double *cd ;
530 for( jj=0 ; jj < ncmat ; jj++ ){
531 cd = (double *)nelmat->vec[jj] ;
532 for( ii=0 ; ii < ntime ; ii++ ) Xar[ii+jj*ntime] = cd[ii] ;
533 }
534 } else {
535 ERROR_exit("Matrix file stored with illegal data type!?") ;
536 }
537
538 /*--- find the stim_times_IM option ---*/
539
540 cgl = NI_get_attribute( nelmat , "BasisNstim") ;
541 if( cgl == NULL ) ERROR_exit("Matrix doesn't have 'BasisNstim' attribute!") ;
542 nbstim = (int)strtod(cgl,NULL) ;
543 if( nbstim <= 0 ) ERROR_exit("Matrix 'BasisNstim' attribute is %d",nbstim) ;
544 for( jj=1 ; jj <= nbstim ; jj++ ){
545 sprintf(nbuf,"BasisOption_%06d",jj) ;
546 cgl = NI_get_attribute( nelmat , nbuf ) ;
547 if( cgl == NULL || strcmp(cgl,"-stim_times_IM") != 0 ) continue ;
548 if( nst > 0 )
549 ERROR_exit("More than one -stim_times_IM option was found in the matrix") ;
550 nst = jj ;
551 sprintf(nbuf,"BasisColumns_%06d",jj) ;
552 cgl = NI_get_attribute( nelmat , nbuf ) ;
553 if( cgl == NULL )
554 ERROR_exit("Matrix doesn't have %s attribute!",nbuf) ;
555 jst_bot = jst_top = -1 ;
556 sscanf(cgl,"%d:%d",&jst_bot,&jst_top) ;
557 if( jst_bot < 0 || jst_top < 0 )
558 ERROR_exit("Can't decode matrix attribute %s",nbuf) ;
559 if( jst_bot == jst_top )
560 ERROR_exit("Matrix attribute %s shows only 1 column for -stim_time_IM:\n"
561 " -->> 3dLSS is meant to be used when more than one stimulus\n"
562 " time was given, and then it computes the response beta\n"
563 " for each stim time separately. If you have only one\n"
564 " stim time with -stim_times_IM, you can use the output\n"
565 " dataset from 3dDeconvolve (or 3dREMLfit) to get that\n"
566 " single beta directly.\n" , nbuf ) ;
567 if( jst_bot >= jst_top || jst_top >= ncmat )
568 ERROR_exit("Matrix attribute %s has illegal value: %d:%d (ncmat=%d)",nbuf,jst_bot,jst_top,ncmat) ;
569 sprintf(nbuf,"BasisName_%06d",jj) ;
570 cgl = NI_get_attribute( nelmat , nbuf ) ;
571 if( cgl != NULL ) stlab = strdup(cgl) ;
572 if( verb > 1 )
573 ININFO_message("-stim_times_IM at stim #%d; cols %d..%d",jj,jst_bot,jst_top) ;
574 }
575 if( nst == 0 )
576 ERROR_exit("Matrix doesn't have any -stim_times_IM options inside :-(") ;
577
578 /*--- mangle matrix to segregate IM regressors from the rest ---*/
579
580 if( verb ) INFO_message("setting up LSS vectors") ;
581
582 imar = LSS_mangle_matrix( imX , jst_bot , jst_top ) ;
583 if( imar == NULL )
584 ERROR_exit("Can't compute LSS 'mangled' matrix :-(") ;
585
586 /*--- setup for LSS computations ---*/
587
588 imA = IMARR_SUBIM(imar,0) ;
589 imC = IMARR_SUBIM(imar,1) ;
590 imS = LSS_setup( imA , imC ) ; DESTROY_IMARR(imar) ;
591 if( imS == NULL )
592 ERROR_exit("Can't complete LSS setup :-((") ;
593 nS = imS->ny ; Sar = MRI_FLOAT_PTR(imS) ;
594
595 if( save1D != NULL ){
596 mri_write_1D( save1D , imS ) ;
597 if( verb ) ININFO_message("saved LSS vectors into file %s",save1D) ;
598 } else if( nodata ){
599 WARNING_message("-nodata used but -save1D not used ==> you get no output!") ;
600 }
601
602 if( nodata || strcmp(prefix,"NULL") == 0 ){
603 INFO_message("3dLSS ends since prefix is 'NULL' or -nodata was used") ;
604 exit(0) ;
605 }
606
607 /*----- create output dataset -----*/
608
609 if( verb ) INFO_message("creating output datset in memory") ;
610
611 outset = EDIT_empty_copy(inset) ;
612 EDIT_dset_items( outset ,
613 ADN_prefix , prefix ,
614 ADN_datum_all , MRI_float ,
615 ADN_brick_fac , NULL ,
616 ADN_nvals , nS ,
617 ADN_ntt , nS ,
618 ADN_none ) ;
619 tross_Copy_History( inset , outset ) ;
620 tross_Make_History( "3dLSS" , argc,argv , outset ) ;
621 for( kk=0 ; kk < nS ; kk++ ){
622 EDIT_substitute_brick( outset , kk , MRI_float , NULL ) ;
623 sprintf(nbuf,"%s#%03d",stlab,kk) ;
624 EDIT_BRICK_LABEL( outset , kk , nbuf ) ;
625 }
626
627 /*----- convert input dataset to vectim -----*/
628
629 if( verb ) INFO_message("loading input dataset into memory") ;
630
631 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
632 inset_mrv = THD_dset_to_vectim( inset , mask , 0 ) ;
633 DSET_unload(inset) ;
634
635 /*----- compute dot products, store results -----*/
636
637 if( verb ) INFO_message("computing away, me buckos!") ;
638
639 nvec = inset_mrv->nvec ;
640 for( kk=0 ; kk < nS ; kk++ ){
641 ss = Sar + kk*ntime ;
642 oo = DSET_ARRAY(outset,kk) ;
643 for( iv=0 ; iv < nvec ; iv++ ){
644 fv = VECTIM_PTR(inset_mrv,iv) ;
645 for( sum=0.0f,ii=0 ; ii < ntime ; ii++ )
646 sum += ss[ii] * fv[goodlist[ii]] ;
647 oo[inset_mrv->ivec[iv]] = sum ;
648 }
649 }
650
651 DSET_write(outset) ; WROTE_DSET(outset) ;
652
653 /*-------- Hasta la vista, baby --------*/
654
655 if( verb )
656 INFO_message("3dLSS finished: total CPU=%.2f Elapsed=%.2f",
657 COX_cpu_time() , COX_clock_time() ) ;
658 exit(0) ;
659 }
660