1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 /*----------------------------- macros and prototypes ------------------------*/
8
9 #undef DR_workspace
10 #define DR_workspace(na,nb) (double *)malloc(sizeof(double)*(2*(na)+(nb)+3)*(nb))
11
12 static void vec_dimred( int nt , int nv , int nfixed , float *inar , int rdim ,
13 int polort , float dt , float fbot , float ftop ,
14 int despike , int vnorm ,
15 double *ws , float *umm , float *svv ) ;
16
17 static void preproc_dimred( int nt , int nvv , float *inar ,
18 int polort ,
19 float dt , float fbot , float ftop ,
20 int despike , int vnorm ) ;
21
22 static void DR_set_svd_sort( int ss ) ;
23 static void DR_svd_double( int m, int n, double *a, double *s, double *u, double *v ) ;
24
25 /*----------------------------------------------------------------------------*/
26
27 static char *chunk_fnam = NULL ; /* for -chunk temp storage */
28 static FILE *chunk_file = NULL ;
29
DR_atexit(void)30 static void DR_atexit(void) /*-- called by exit(): delete chunk_file --*/
31 {
32 if( chunk_file != NULL ){ fclose(chunk_file) ; chunk_file = NULL ; }
33 if( chunk_fnam != NULL ){
34 INFO_message("Deleting -chunk file %s",chunk_fnam) ;
35 remove(chunk_fnam) ; chunk_fnam = NULL ;
36 }
37 return ;
38 }
39
40 /*----------------------------------------------------------------------------*/
41
append_to_prefix(char * pre,char * app)42 static char * append_to_prefix( char *pre , char *app )
43 {
44 char *prA , *prB , *oot ; int lll ;
45
46 prA = strdup(pre) ; lll = strlen(prA) ;
47
48 if( STRING_HAS_SUFFIX(prA,".nii") && lll > 4 ){
49 prA[lll-4] = '\0' ;
50 prB = strdup(".nii") ;
51 } else if( STRING_HAS_SUFFIX(prA,".nii.gz") && lll > 7 ){
52 prA[lll-7] = '\0' ;
53 prB = strdup(".nii.gz") ;
54 } else {
55 prB = strdup("\0") ;
56 }
57
58 oot = (char *)malloc(sizeof(char)*(lll+strlen(app)+9)) ;
59 strcpy(oot,prA) ;
60 strcat(oot,"_") ;
61 strcat(oot,app) ;
62 strcat(oot,prB) ;
63 free(prB) ; free(prA) ; return oot ;
64 }
65
66 /*----------------------------------------------------------------------------*/
67
DR_syntax(void)68 void DR_syntax(void)
69 {
70 printf(
71 "Program 3dDimRed does 'dimensional reduction' on a collection of input 3D+time\n"
72 "datasets, and produces as output a smaller collection of 3D+time datasets.\n"
73 "The process of dimensional reduction is via the Singular Value Decomposition\n"
74 "(SVD), and the goal is to produce (at each voxel) a collection of time series\n"
75 "vectors whose linear span captures most of the variability present in the\n"
76 "input. The number of output datasets is specified by the user, via the\n"
77 "'-rdim' option.\n"
78 "\n"
79 "Usage: 3dDimRed [options] inset1 inset2 ...\n"
80 "\n"
81 "where 'inset1' is the first input 3D+time dataset, 'inset2' is the second\n"
82 "input 3D+time dataset, and so forth. All input datasets must be on the\n"
83 "same 3D+time grid, and there must be at least 2 input datasets.\n"
84 "\n"
85 "--------\n"
86 "OPTIONS:\n"
87 "--------\n"
88 " -rdim rr = The number 'rr' is the number of dimensions to use for the\n"
89 " output = the number of output 3D+time datasets.\n"
90 " * 'rr' must be more than 0 and less than the number of\n"
91 " input datasets.\n"
92 " * If 'rr' is set to 0, or is not given, then no output\n"
93 " 3D+time datasets will be produced!\n"
94 "\n"
95 " -prefix pp = The string 'pp' is the root prefix for the output datasets.\n"
96 " The first one will get the prefix 'pp_001', etc.\n"
97 " * If '-prefix' is not given, then the default root prefix\n"
98 " is 'dimred'.\n"
99 "\n"
100 " -sing = Save the singular values at each voxel into a dataset with\n"
101 " prefix 'pp_sing', where 'pp' is the root prefix.\n"
102 " * If this option is not given, the singular values are not\n"
103 " saved.\n"
104 "\n"
105 " -1Dcols a b ... = This option lets you specify some extra time series to be\n"
106 " included in the dimension reduction, which are fixed for each\n"
107 " input voxel. Each column of each 1D file following '-1Dcols'\n"
108 " is added to the collection of vectors to be processed at each\n"
109 " voxel.\n"
110 "\n"
111 " -input x y ... = Alternate way to input 3D+time datasets.\n"
112 "\n"
113 " -polort qq = Detrend the time series with polynomial basis of order 'qq'\n"
114 " prior to further processing.\n"
115 " * You cannot use this option with '-band'!\n"
116 " * If neither '-polort' nor '-band' is given, then the mean\n"
117 " of each time series is removed (e.g., '-polort 0').\n"
118 "\n"
119 " -band fbot ftop = This option specifies to bandpass the time series prior\n"
120 " to further processing (as in program 3dBandpass).\n"
121 " * You cannot use this option with '-polort'!\n"
122 "\n"
123 " -dt dd = Set time step to 'dd' seconds (for use with '-band').\n"
124 " *OR* * Usually the time step is taken from the 3D+time dataset\n"
125 " -TR dd header, but '-dt' (or '-TR') lets you over-ride that.\n"
126 " * If all inputs are 1D files, then this option is needed\n"
127 " to set the right time step; otherwise, TR is taken as 1.\n"
128 "\n"
129 " -despike = Despike each input time series before other processing.\n"
130 " * Hopefully, you don't need to do this, which is why it\n"
131 " is optional.\n"
132 "\n"
133 " -novnorm = By default, just before the SVD, each time series vector is\n"
134 " normalized to L2 magnitude 1. Use '-novnorm' to turn this\n"
135 " step off.\n"
136 "\n"
137 " -mask mset = Mask dataset\n"
138 " -automask = Create mask from first input dataset\n"
139 "\n"
140 " -chunk = 3dDimRed can use up a LOT of memory. This option lets\n"
141 " you cut down on that by having it operate on the voxels\n"
142 " in chunks. The downside to doing this is more disk I/O.\n"
143 " * '-chunk' alone means do 10,000 voxels at one time.\n"
144 " * '-chunk 999' means to do 999 voxels at one time (etc).\n"
145 "\n"
146 " -quiet = Suppress some of the fun informative progress messages.\n"
147 "\n"
148 "------\n"
149 "NOTES:\n"
150 "------\n"
151 " * You can input ONLY 1D files, via the -1Dcols option. In that case,\n"
152 " there is only 1 'voxel', and instead of 'rdim' datasets being output\n"
153 " output to hold the results, a single 1D dataset with 'rdim' columns\n"
154 " will be produced.\n"
155 " * The output vectors are the first 'rdim' left singular vectors, as\n"
156 " would be produced by '1dsvd -1Dleft' (for example). They are intended\n"
157 " to be used in further processing scripts, and are probably not useful\n"
158 " by themselves.\n"
159 " * Written by Zhark the Singular -- April 2012\n"
160 ) ;
161 PRINT_AFNI_OMP_USAGE("3dDimRed",NULL) ;
162 PRINT_COMPILE_DATE ; exit(0) ;
163 }
164
165 /*----------------------------------------------------------------------------*/
166
main(int argc,char * argv[])167 int main( int argc , char *argv[] )
168 {
169 int nopt=1 , verb=1 ;
170 int do_despike=0 , do_vnorm=1 , do_automask=0 , do_sing=0 , polort=-1 ;
171 float fbot=-666.0f, ftop=-999.9f , dt=0.0f ; int have_freq=0 ;
172 byte *mask=NULL ; int mask_nx=0,mask_ny=0,mask_nz=0,nmask=0 , *imask=NULL ;
173 int nx,ny,nz,nt,nvox , nfft=0 , rdim=0 , kk,ii,qq,pp , nbad=0 ;
174 int ndset=0,nim=0,nvim=0 , ndim=0 , chunk=0,im,imbot,imtop ;
175 char *prefix="dimred" ;
176 MRI_IMARR *imarc ;
177 RwcPointer_array *dsar ;
178 THD_3dim_dataset *iset=NULL,*kset ; MRI_IMAGE *tim ;
179 float *cfixv=NULL , *cvect=NULL , *csave=NULL ;
180 size_t msize , ntsiz , ndsiz ;
181 char *prefout , prefapp[32] ;
182
183 /*-- startup --*/
184
185 mainENTRY("3dDimRed"); machdep();
186 AFNI_logger("3dDimRed",argc,argv);
187 PRINT_VERSION("3dDimRed"); AUTHOR("Uncle John's Band");
188 enable_mcw_malloc() ;
189 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
190
191 INIT_XTARR(dsar) ; INIT_IMARR(imarc) ;
192 while( nopt < argc && argv[nopt][0] == '-' ){
193
194 if( strcmp(argv[nopt],"-") == 0 ){ nopt++ ; continue ; }
195
196 if( strcasecmp(argv[1],"-help") == 0 ) DR_syntax() ;
197
198 if( strcasecmp(argv[nopt],"-1Dcols") == 0 ){
199 MRI_IMAGE *qim ;
200 if( ++nopt >= argc ) ERROR_exit("need an argument after %s",argv[nopt-1]) ;
201 if( argv[nopt][0] == '-' ) ERROR_exit("Illegal argument after %s",argv[nopt-1]) ;
202 for( ; nopt < argc && argv[nopt][0] != '-' ; nopt++ ){
203 qim = mri_read_1D(argv[nopt]) ;
204 if( qim == NULL ) ERROR_exit("Can't read from -1Dcols file '%s'",argv[nopt]) ;
205 mri_add_name(argv[nopt],qim) ;
206 ADDTO_IMARR(imarc,qim) ;
207 }
208 continue ;
209 }
210
211 if( strcasecmp(argv[nopt],"-rdim") == 0 ){
212 if( ++nopt >= argc ) ERROR_exit("need an argument after -rdim!") ;
213 rdim = (int)strtod(argv[nopt],NULL) ;
214 if( rdim <= 0 )
215 WARNING_message("-rdim value means no output time series datasets") ;
216 nopt++ ; continue ;
217 }
218
219 if( strcasecmp(argv[nopt],"-chunk") == 0 ){
220 if( ++nopt < argc && isdigit(argv[nopt][0]) ){
221 chunk = (int)strtod(argv[nopt++],NULL) ;
222 if( chunk < 666 ) chunk = 666 ; /* Lame user */
223 } else {
224 chunk = 10000 ;
225 }
226 nopt++ ; continue ;
227 }
228
229 if( strcasecmp(argv[nopt],"-despike") == 0 ){
230 do_despike++ ; nopt++ ; continue ;
231 }
232
233 if( strcasecmp(argv[nopt],"-sing") == 0 ){
234 do_sing++ ; nopt++ ; continue ;
235 }
236
237 if( strcasecmp(argv[nopt],"-prefix") == 0 ){
238 if( ++nopt >= argc ) ERROR_exit("need an argument after -prefix!") ;
239 prefix = strdup(argv[nopt]) ;
240 if( !THD_filename_ok(prefix) ) ERROR_exit("bad -prefix option!") ;
241 nopt++ ; continue ;
242 }
243
244 if( strcasecmp(argv[nopt],"-automask") == 0 ){
245 if( mask != NULL ) ERROR_exit("Can't use -mask AND -automask!") ;
246 do_automask = 1 ; nopt++ ; continue ;
247 }
248
249 if( strcasecmp(argv[nopt],"-mask") == 0 ){
250 THD_3dim_dataset *mset ;
251 if( ++nopt >= argc ) ERROR_exit("Need argument after '-mask'") ;
252 if( mask != NULL || do_automask ) ERROR_exit("Can't have two mask inputs") ;
253 mset = THD_open_dataset( argv[nopt] ) ;
254 CHECK_OPEN_ERROR(mset,argv[nopt]) ;
255 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
256 mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);
257 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
258 if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[nopt]) ;
259 nmask = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;
260 if( verb ) INFO_message("Number of voxels in mask = %d",nmask) ;
261 if( nmask < 1 ) ERROR_exit("Mask is empty: nothing to process :-(") ;
262 nopt++ ; continue ;
263 }
264
265 if( strcasecmp(argv[nopt],"-novnorm") == 0 ){
266 do_vnorm = 0 ; nopt++ ; continue ;
267 }
268
269 if( strcasecmp(argv[nopt],"-quiet") == 0 ){
270 verb = 0 ; nopt++ ; continue ;
271 }
272
273 if( strcasecmp(argv[nopt],"-dt") == 0 || strcasecmp(argv[nopt],"-TR") == 0 ){
274 if( ++nopt >= argc ) ERROR_exit("need an argument after %s",argv[nopt-1]) ;
275 dt = (float)strtod(argv[nopt],NULL) ;
276 if( dt <= 0.0f ) WARNING_message("value after %s illegal!",argv[nopt-1]) ;
277 nopt++ ; continue ;
278 }
279
280 if( strcasecmp(argv[nopt],"-polort") == 0 ){
281 if( ++nopt >= argc ) ERROR_exit("need an argument after -polort!") ;
282 if( have_freq ) ERROR_exit("Can't use -polort and -band together!") ;
283 polort = (int)strtod(argv[nopt],NULL) ;
284 nopt++ ; continue ;
285 }
286
287 if( strcasecmp(argv[nopt],"-band") == 0 ){
288 if( ++nopt >= argc-1 ) ERROR_exit("need 2 arguments after -band!") ;
289 if( polort >= 0 ) ERROR_exit("Can't use -band and -polort together!") ;
290 fbot = strtod(argv[nopt++],NULL) ;
291 ftop = strtod(argv[nopt++],NULL) ;
292 if( fbot < 0.0f || ftop <= fbot ) WARNING_message("values after -band are illegal!") ;
293 else have_freq = 1 ;
294 continue ;
295 }
296
297 if( strcasecmp(argv[nopt],"-input") == 0 ){
298 if( ++nopt >= argc ) ERROR_exit("Need argument after '%s'",argv[nopt-1]);
299 if( argv[nopt][0] == '-' ) ERROR_exit("Illegal argument after '%s'",argv[nopt-1]) ;
300 for( ; nopt < argc && argv[nopt][0] != '-' ; nopt++ ){
301 kset = THD_open_dataset(argv[nopt]) ;
302 if( kset == NULL ) ERROR_exit("Can't read dataset '%s'",argv[nopt]) ;
303 ADDTO_XTARR(dsar,kset) ;
304 ii = XTARR_NUM(dsar)-1 ; XTARR_IC(dsar,ii) = IC_DSET ;
305 }
306 continue ;
307 }
308
309 ERROR_message("Unknown option: %s\n",argv[nopt]) ;
310 suggest_best_prog_option(argv[0], argv[nopt]);
311 exit(1);
312
313 } /*-- end of loop over options --*/
314
315 if( argc < 2 ) DR_syntax() ;
316
317 /*--- read rest of args as input datasets ---*/
318
319 for( ; nopt < argc ; nopt++ ){
320 kset = THD_open_dataset(argv[nopt]) ;
321 if( kset == NULL ) ERROR_exit("Can't read dataset '%s'",argv[nopt]) ;
322 ADDTO_XTARR(dsar,kset) ;
323 ii = XTARR_NUM(dsar)-1 ; XTARR_IC(dsar,ii) = IC_DSET ;
324 }
325
326 /*---- check stuff ----*/
327
328 ndset = XTARR_NUM(dsar) ; nim = IMARR_COUNT(imarc) ;
329 if( ndset < 2 && nim <= 0 )
330 ERROR_exit("Inadequate number of input datasets! What am I supposed to do?!") ;
331 if( ndset == 0 && nim == 1 && IMARR_SUBIM(imarc,0)->ny == 1 )
332 ERROR_exit("Inadequate number of input 1D columns! What am I supposed to do?!") ;
333
334 if( rdim <= 0 && !do_sing )
335 ERROR_exit("No output specified by either -rdim or -sing :-(") ;
336
337 /* check input datasets and/or 1D files for matches of various kinds */
338
339 if( ndset > 0 ){
340 iset = (THD_3dim_dataset *)XTARR_XT(dsar,0) ;
341 nx = DSET_NX(iset); ny = DSET_NY(iset); nz = DSET_NZ(iset); nvox = nx*ny*nz;
342 nt = DSET_NVALS(iset) ;
343 if( dt <= 0.0f ) dt = DSET_TR(iset) ;
344 if( nt < 2 ){
345 ERROR_message("Dataset %s has only 1 point in the time direction!",DSET_HEADNAME(iset)) ;
346 nbad++ ;
347 }
348 if( dt <= 0.0f && have_freq ){
349 WARNING_message("TR not present in first input dataset: using 1 s") ;
350 dt = 1.0f ;
351 }
352
353 for( kk=1 ; kk < ndset ; kk++ ){
354 kset = (THD_3dim_dataset *)XTARR_XT(dsar,kk) ;
355 if( DSET_NVOX(kset) != nvox ){
356 ERROR_message("Dataset %s has %d voxels, but first dataset has %d",
357 DSET_HEADNAME(kset) , DSET_NVOX(kset) , nvox ) ;
358 nbad++ ;
359 }
360 if( DSET_NVALS(kset) != nt ){
361 ERROR_message("Dataset %s has %d time points, but first dataset has %d",
362 DSET_HEADNAME(kset) , DSET_NVALS(kset) , nt ) ;
363 nbad++ ;
364 }
365 }
366 } else { /* no actual input datasets! */
367 tim = IMARR_SUBIM(imarc,0) ;
368 nvox = 1 ; nt = tim->nx ;
369 if( nt < 2 ){
370 ERROR_message("1D file %s has only 1 time point!",tim->name) ;
371 nbad++ ;
372 }
373 if( mask != NULL ){
374 WARNING_message("skipping -mask since no 3D+time datasets were input") ;
375 free(mask) ; mask = NULL ;
376 }
377 }
378
379 for( kk=0 ; kk < nim ; kk++ ){
380 tim = IMARR_SUBIM(imarc,kk) ;
381 if( tim->nx != nt ){
382 ERROR_message("1D file %s has %d time points, but should have %d",
383 tim->name , tim->nx , nt ) ;
384 nbad++ ;
385 }
386 nvim += tim->ny ; /* count of number of 1D time series vectors */
387 }
388
389 if( mask != NULL ){
390 if( mask_nx != nx || mask_ny != ny || mask_nz != nz ){
391 ERROR_message("-mask dataset grid doesn't match input dataset") ;
392 nbad++ ;
393 }
394 } else if( do_automask ){
395 if( ndset == 0 ){
396 WARNING_message("ignoring -automask since no 3D+time datasets were input") ;
397 } else {
398 mask = THD_automask(iset) ;
399 if( mask == NULL ){
400 ERROR_message("Can't create -automask from first input dataset?") ; nbad++ ;
401 }
402 nmask = THD_countmask( nvox , mask ) ;
403 if( verb ) INFO_message("Number of voxels in automask = %d",nmask);
404 if( nmask < 1 ){ ERROR_message("Automask is too small to process"); nbad++; }
405 }
406 }
407 if( mask == NULL ){
408 mask = (byte *)malloc(sizeof(byte)*nvox) ; nmask = nvox ;
409 memset(mask,1,sizeof(byte)*nvox) ;
410 if( verb && nvox > 1 ) INFO_message("No mask ==> processing all %d voxels",nvox);
411 }
412
413 /* make list of voxels to process, from mask */
414
415 imask = (int *)malloc(sizeof(int)*nmask) ;
416 for( qq=im=0 ; im < nvox ; im++ ){
417 if( mask[im] ) imask[qq++] = im ;
418 }
419 free(mask) ; mask = NULL ;
420
421 ndim = ndset + nvim ; /* total number of input dimensions */
422 if( ndim <= rdim ){
423 ERROR_message("Input dimension = %d is not bigger than output = %d",ndim,rdim) ;
424 nbad++ ;
425 }
426
427 ndsiz = sizeof(float)*ndim ; /* size of singular values vector */
428 ntsiz = sizeof(float)*nt ; /* size of time series vectors */
429
430 if( nbad > 0 )
431 ERROR_exit("Can't continue after above ERRORs") ; /*** BAD USER ***/
432
433 /*---- create temp file for usufruct of chunk -----*/
434
435 if( chunk > 0 && ndset == 0 ) chunk = 0 ;
436
437 if( chunk > 0 && chunk < nmask ){
438 chunk_fnam = UNIQ_idcode() ;
439 chunk_fnam[0] = 'D' ; chunk_fnam[1] = 'R' ; chunk_fnam[2] = 'X' ;
440 chunk_file = fopen( chunk_fnam , "w+b" ) ;
441 if( chunk_file == NULL )
442 ERROR_exit("Unable to create -chunk temporary file %s",chunk_fnam) ;
443 INFO_message("Creating -chunk temporary file %s",chunk_fnam) ;
444 atexit(DR_atexit) ;
445 }
446
447 if( chunk <= 0 || chunk > nmask ) chunk = nmask ;
448
449 if( verb )
450 INFO_message("Chunk size = %s voxels",commaized_integer_string(nmask)) ;
451
452 DR_set_svd_sort(-1) ;
453
454 /*--- allocate space to hold results from one chunk ---*/
455 /*--- for each voxel, rdim vectors of length nt, ---*/
456 /*--- plus the singular values, of length ndim. ---*/
457
458 /* macro == vi-th output nt-vector for the ic-th voxel */
459
460 #undef CSUU
461 #define CSUU(ic,vi) (csave+((ic)*(rdim*nt+ndim)+(vi)*nt))
462
463 /* macro == output ndim-vector of singular values for the ic-th voxel */
464
465 #undef CSSV
466 #define CSSV(ic) (csave+((ic)*(rdim*nt+ndim)+rdim*nt))
467
468 /* macro == time series nt-vector from the kd-th dataset, ic-th voxel */
469
470 #undef CVD
471 #define CVD(kd,ic) (cvect+(kd+(ic)*ndset)*nt)
472
473 /* macro == ic-th fixed time series nt-vector */
474
475 #undef CFV
476 #define CFV(ic) (cfixv+(ic)*nt)
477
478 msize = (rdim*ntsiz+ndsiz)*chunk ;
479 csave = (float *)malloc(msize) ;
480 if( verb || csave == NULL )
481 ININFO_message("Output memory chunk = %s bytes (%s)",
482 commaized_integer_string(msize) ,
483 approximate_number_string((double)msize) ) ;
484 if( csave == NULL )
485 ERROR_exit("Can't allocate memory for output chunk") ;
486
487 /* memory to hold dataset inputs from one chunk: ndset vectors of length nt */
488
489 if( ndset > 0 ){
490 msize = ntsiz*chunk*ndset ;
491 cvect = (float *)malloc(msize) ;
492 if( verb || cvect == NULL )
493 ININFO_message("Dataset input memory chunk = %s bytes (%s)",
494 commaized_integer_string(msize) ,
495 approximate_number_string((double)msize) ) ;
496 if( cvect == NULL )
497 ERROR_exit("Can't allocate memory for input chunk") ;
498 }
499
500 /* memory for fixed vectors, and pre-process them now */
501
502 if( nvim > 0 ){
503 MRI_IMAGE *tim ; float *tar ;
504 cfixv = (float *)malloc(ntsiz*nvim) ;
505 if( cfixv == NULL )
506 ERROR_exit("Can't allocate memory for fixed vectors") ;
507 for( qq=kk=0 ; kk < nim ; kk++ ){
508 tim = IMARR_SUBIM(imarc,kk) ; tar = MRI_FLOAT_PTR(tim) ;
509 for( pp=0 ; pp < tim->ny ; pp++,qq++ )
510 AAmemcpy( CFV(qq) , tar+pp*nt , ntsiz ) ;
511 }
512 DESTROY_IMARR(imarc) ;
513 preproc_dimred( nt , nvim , cfixv ,
514 polort , dt , fbot , ftop , do_despike , do_vnorm ) ;
515 }
516
517 #ifdef USE_OMP
518 #pragma omp parallel
519 {
520 if( verb && omp_get_thread_num() == 0 )
521 INFO_message("OpenMP thread count = %d",omp_get_num_threads()) ;
522 }
523 #else
524 if( verb && nmask > 1 ) INFO_message("Start main voxel loop") ;
525 #endif
526
527 /********** loop over chunks **********/
528
529 for( imbot=0 ; imbot < nmask ; imbot+=chunk ){
530
531 imtop = imbot + chunk ; if( imtop > nmask ) imtop = nmask ;
532
533 if( verb && chunk < nmask )
534 ININFO_message("start chunk: voxels %d..%d",imbot,imtop-1) ;
535
536 /* extract dataset data into cvect */
537
538 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
539 for( kk=0 ; kk < ndset ; kk++ ){
540 kset = (THD_3dim_dataset *)XTARR_XT(dsar,kk) ;
541 ININFO_message("Load dataset %s",DSET_HEADNAME(kset)) ;
542 DSET_load(kset) ;
543 if( !DSET_LOADED(kset) )
544 ERROR_exit("Can't load data from dataset %s",DSET_BRIKNAME(kset)) ;
545 for( ii=imbot ; ii < imtop ; ii++ )
546 THD_extract_array( imask[ii] , kset , 0 , CVD(kk,ii-imbot) ) ;
547 DSET_unload(kset) ;
548 }
549 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
550
551 AFNI_OMP_START ;
552 #pragma omp parallel if( imtop-imbot > 66 )
553 { float *tsar , *umat , *sval ; double *ws ; int iv,kd ;
554 tsar = (float *)malloc(ntsiz*ndim) ; /* input vectors */
555 umat = (rdim <= 0) ? NULL
556 : (float *)malloc(sizeof(float)*ntsiz*rdim); /* output vectors */
557 sval = (float *)malloc(ndsiz) ; /* output sing vals */
558 ws = DR_workspace(nt,ndim) ;
559 #pragma omp for
560 for( iv=imbot ; iv < imtop ; iv++ ){
561 /* copy fixed data into tsar */
562 if( cfixv != NULL ){
563 ININFO_message("iv = %d: cfixv -> tsar",iv) ;
564 AAmemcpy( tsar , cfixv , ntsiz*nvim ) ;
565 }
566 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
567 /* copy dataset data into tsar */
568 if( cvect != NULL ){
569 ININFO_message("iv = %d: cvect -> tsar",iv) ;
570 for( kd=0 ; kd < ndset ; kd++ )
571 AAmemcpy( tsar+(nvim+kd)*nt , CVD(kd,iv-imbot) , ntsiz ) ;
572 }
573 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
574 /* process it into umat and sval */
575 ININFO_message("iv = %d: call vec_dimred",iv) ;
576 vec_dimred( nt , ndim , nvim , tsar , rdim ,
577 polort , dt , fbot , ftop , do_despike , do_vnorm ,
578 ws , umat , sval ) ;
579 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
580 /* copy umat and sval into csave */
581 if( umat != NULL ){
582 ININFO_message("iv = %d: umat -> csave",iv) ;
583 AAmemcpy( CSUU(iv-imbot,0) , umat , ntsiz*rdim ) ;
584 }
585 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
586 ININFO_message("iv = %d: sval -> csave",iv) ;
587 AAmemcpy( CSSV(iv-imbot) , sval , ndsiz ) ;
588 } /* end of parallel-ized loop over voxels */
589 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
590 free(ws) ; free(sval) ; free(tsar) ; if( umat != NULL ) free(umat) ;
591 } /* end of parallel-ized code */
592 AFNI_OMP_END ;
593
594 /* if using temporary file, write csave to disk */
595
596 if( chunk_file != NULL ){
597 msize = fwrite( csave, (rdim*ntsiz+ndsiz), imtop-imbot, chunk_file ) ;
598 if( msize < imtop-imbot )
599 ERROR_exit("Failure to write to -chunk file %s -- is disk full?",
600 chunk_fnam ) ;
601 }
602 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
603
604 } /* end of loop over chunks of voxels */
605
606 if( cvect != NULL ) free(cvect) ;
607 if( cfixv != NULL ) free(cfixv) ;
608
609 if( ndset > 0 && nvox > 1 ){ /* create output datasets for vectors */
610
611 for( kk=0 ; kk < rdim ; kk++ ){
612
613 /* create the empty dataset */
614
615 sprintf(prefapp,"%03d",kk+1) ;
616 prefout = append_to_prefix( prefix , prefapp ) ;
617 iset = (THD_3dim_dataset *)XTARR_XT(dsar,0) ;
618 kset = EDIT_empty_copy(iset) ;
619 EDIT_dset_items( kset , ADN_prefix , prefout , NULL ) ;
620 tross_Copy_History( iset , kset ) ;
621 tross_Make_History( "3dDimRed" , argc,argv , kset ) ;
622 free(prefout) ;
623
624 /* put in the sub-bricks */
625 for( pp=0 ; pp < nt ; pp++ )
626 EDIT_substitute_brick( kset , pp , MRI_float , NULL ) ;
627
628 if( chunk_file != NULL ) rewind(chunk_file) ;
629
630 /* loop over chunks, fill dataset time series */
631
632 for( imbot=0 ; imbot < nmask ; imbot += chunk ){
633
634 imtop = imbot + chunk ; if( imtop > nmask ) imtop = nmask ;
635
636 if( chunk_file != NULL ){ /* get data back from disk */
637 msize = fread( csave, (rdim*ntsiz+ndsiz), imtop-imbot, chunk_file ) ;
638 if( msize < imtop-imbot ) /* should never happen */
639 ERROR_exit("Failure to read from -chunk file %s :-(",chunk_fnam) ;
640 }
641
642 for( im=imbot ; im < imtop ; im++ ) /* shove it into the time series */
643 THD_insert_series( imask[im], kset, nt, MRI_float, CSUU(im-imbot,kk),0 ) ;
644 }
645
646 DSET_write(kset) ; WROTE_DSET(kset) ; DSET_delete(kset) ;
647 }
648
649 } else if( rdim > 0 ){ /* create a 1D output file */
650
651 MRI_IMAGE *qim ; float *qar ; char *prefout ;
652
653 qim = mri_new(nt,rdim,MRI_float) ; qar = MRI_FLOAT_PTR(qim) ;
654 AAmemcpy( qar , CSUU(0,0) , ntsiz*rdim ) ;
655 prefout = append_to_prefix( prefix , "vect.1D" ) ;
656 mri_write_1D(prefout,qim) ; mri_free(qim) ;
657 INFO_message("Wrote vector file %s",prefout) ; free(prefout) ;
658
659 } /* end of writing out vectors */
660
661 /* write the singular values dataset? */
662
663 if( do_sing ){
664 if( ndset > 0 && nvox > 1 ){
665 prefout = append_to_prefix( prefix , "sing" ) ;
666 iset = (THD_3dim_dataset *)XTARR_XT(dsar,0) ;
667 kset = EDIT_empty_copy(iset) ;
668 EDIT_dset_items( kset ,
669 ADN_prefix , prefout ,
670 ADN_nvals , ndim ,
671 NULL ) ;
672 tross_Copy_History( iset , kset ) ;
673 tross_Make_History( "3dDimRed" , argc,argv , kset ) ;
674 free(prefout) ;
675
676 for( pp=0 ; pp < ndim ; pp++ )
677 EDIT_substitute_brick( kset , pp , MRI_float , NULL ) ;
678
679 if( chunk_file != NULL ) rewind(chunk_file) ;
680
681 for( imbot=0 ; imbot < nmask ; imbot += chunk ){
682
683 imtop = imbot + chunk ; if( imtop > nmask ) imtop = nmask ;
684
685 if( chunk_file != NULL ){
686 msize = fread( csave, (rdim*ntsiz+ndsiz), imtop-imbot, chunk_file ) ;
687 if( msize < imtop-imbot ) /* should never happen */
688 ERROR_exit("Failure to read from -chunk file %s :-(",chunk_fnam) ;
689 }
690
691 for( im=imbot ; im < imtop ; im++ )
692 THD_insert_series( imask[im], kset, ndim, MRI_float, CSSV(im-imbot),0 ) ;
693
694 }
695
696 DSET_write(kset) ; WROTE_DSET(kset) ; DSET_delete(kset) ;
697
698 } else { /* write a 1D file with the singular values */
699
700 MRI_IMAGE *qim ; float *qar ; char *prefout ;
701
702 qim = mri_new(ndim,1,MRI_float) ; qar = MRI_FLOAT_PTR(qim) ;
703 AAmemcpy( qar , CSSV(0) , ndsiz ) ;
704 prefout = append_to_prefix( prefix , "sing.1D" ) ;
705 mri_write_1D(prefout,qim) ; mri_free(qim) ;
706 INFO_message("Wrote singular values file %s",prefout) ; free(prefout) ;
707 }
708
709 } /* end of dosing */
710
711 free(csave) ;
712 exit(0) ;
713 }
714
715 /*----------------------------------------------------------------------------*/
716
717 #undef INVEC
718 #define INVEC(k) (inar+(k)*nt)
719
preproc_dimred(int nt,int nvv,float * inar,int polort,float dt,float fbot,float ftop,int despike,int vnorm)720 static void preproc_dimred( int nt , int nvv , float *inar ,
721 int polort ,
722 float dt , float fbot , float ftop ,
723 int despike , int vnorm )
724 {
725 static double *pcc=NULL ; /* for polort detrending */
726 static float **pref=NULL ;
727 int ii,jj,kk ;
728
729 /*-- cleanup the trash call? --*/
730
731 #pragma omp critical
732 { if( nt <= 0 ){
733 if( pcc != NULL ){ free(pcc); pcc = NULL; }
734 if( pref != NULL ){
735 for( ii=0 ; ii <= polort ; ii++ ) free(pref[ii]) ;
736 free(pref) ; pref = NULL ;
737 }
738 }
739 }
740 if( nt <= 0 || nvv <= 0 ) return ;
741
742 /*-- first time in: create polort detrending stuff? --*/
743
744 #pragma omp critical
745 { if( pcc == NULL && polort > 0 ){
746 pref = THD_build_polyref(polort+1,nt) ;
747 pcc = startup_lsqfit( nt , NULL , polort+1 , pref ) ;
748 }
749 }
750
751 /*-- despike? --*/
752
753 if( despike ){
754 for( kk=0 ; kk < nvv ; kk++ ) THD_despike9(nt,INVEC(kk)) ;
755 }
756
757 /*-- polort detrend? --*/
758
759 if( polort == 0 ){
760 for( kk=0 ; kk < nvv ; kk++ ) THD_const_detrend(nt,INVEC(kk),NULL) ;
761 } else if( polort > 0 && pcc != NULL ){
762 float *coef , *vec , *pr , cf ;
763 for( kk=0 ; kk < nvv ; kk++ ){
764 vec = INVEC(kk) ;
765 coef = delayed_lsqfit( nt , vec , polort+1 , pref , pcc ) ;
766 if( coef != NULL ){
767 for( jj=0 ; jj <= polort ; jj++ ){
768 cf = coef[jj] ; pr = pref[jj] ;
769 for( ii=0 ; ii < nt ; ii++ ) vec[ii] -= cf * pr[ii] ;
770 }
771 free(coef) ;
772 }
773 }
774 }
775
776 /*-- bandpass? --*/
777
778 if( dt > 0.0f && fbot >= 0.0f && fbot < ftop ){
779 float **vvar = malloc(sizeof(float *)*nvv) ;
780 for( kk=0 ; kk < nvv ; kk++ ) vvar[kk] = INVEC(kk) ;
781 THD_bandpass_vectors( nt , nvv , vvar , dt , fbot , ftop , 2 , 0 , NULL ) ;
782 free(vvar) ;
783 }
784
785 /*-- L2 normalize? --*/
786
787 if( vnorm ){
788 for( kk=0 ; kk < nvv ; kk++ ) THD_normalize(nt,INVEC(kk)) ;
789 }
790
791 return ;
792 }
793
794 /*----------------------------------------------------------------------------*/
795 /* assume the first nfixed vectors are already pre-processed */
796
vec_dimred(int nt,int nv,int nfixed,float * inar,int rdim,int polort,float dt,float fbot,float ftop,int despike,int vnorm,double * ws,float * umm,float * svv)797 static void vec_dimred( int nt , int nv , int nfixed , float *inar , int rdim ,
798 int polort , float dt , float fbot , float ftop ,
799 int despike , int vnorm ,
800 double *ws , float *umm , float *svv )
801 {
802 int ntv=nt*nv , ntr=nt*rdim , ii,jj,pp,nnz ;
803 float *vv ;
804 double *wss , *amat , *sval , *umat , *vmat ;
805 MRI_IMAGE *outim ; float *outar ;
806
807 if( umm == NULL && svv == NULL ) return ; /* WTF? */
808
809 ININFO_message(" preproc") ;
810 preproc_dimred( nt , nv-nfixed , inar+nfixed*nt ,
811 polort , dt,fbot,ftop , despike , vnorm ) ;
812 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
813
814 wss = ws ; if( wss == NULL ) wss = DR_workspace(nt,nv) ;
815
816 amat = wss ; /* nt*nv long */
817 umat = amat + ntv ; /* nt*nv long */
818 vmat = umat + nv*nv ; /* nv*nv long */
819 sval = vmat + ntv ; /* nv long */
820
821 /** copy vectors into amat, removing all zero columns **/
822
823 ININFO_message(" copy data in") ;
824 #undef AM
825 #undef IM
826 #define AM(i,j) amat[(i)+(j)*nt]
827 #define IM(i,j) inar[(i)+(j)*nt]
828 for( jj=pp=0 ; jj < nv ; jj++ ){
829 for( nnz=ii=0 ; ii < nt ; ii++ ){
830 AM(ii,pp) = (double)IM(ii,jj) ; nnz += ( AM(ii,pp) != 0.0 ) ;
831 }
832 if( nnz > 0 ) pp++ ;
833 }
834 #undef AM
835 #undef IM
836 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
837
838 if( pp < nv )
839 ININFO_message("removed %d all zero vectors",nv-pp) ;
840
841 /* initialize outputs to 0? */
842
843 if( pp < nv ){
844 if( umm != NULL ) AAmemset(umm,0,sizeof(float)*ntr) ;
845 if( svv != NULL ) AAmemset(svv,0,sizeof(float)*nv ) ;
846 if( pp == 0 ){ /* if all inputs 0, output 0 */
847 if( wss != ws ) free(wss) ;
848 return ;
849 }
850 }
851
852 /* second dimension is now pp, possibly less than nv */
853
854 ININFO_message("**svd") ;
855 DR_svd_double( nt , pp , amat , sval , umat , vmat ) ;
856 ININFO_message("**svd done") ;
857 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
858
859 if( umm != NULL ){
860 ININFO_message(" umat -> umm") ;
861 nnz = nt*MIN(pp,rdim) ;
862 for( ii=0 ; ii < nnz ; ii++ ) umm[ii] = (float)umat[ii] ;
863 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
864 }
865 if( svv != NULL ){
866 ININFO_message(" sval -> svv") ;
867 for( ii=0 ; ii < pp ; ii++ ) svv[ii] = (float)sval[ii] ;
868 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
869 }
870
871 if( wss != ws ) free(wss) ;
872 return ;
873 }
874
875 /*--------------------------------------------------------------------*/
876
877 #include "eispack.h"
878
879 #ifdef isfinite
880 # define IS_GOOD_FLOAT(x) isfinite(x)
881 #else
882 # define IS_GOOD_FLOAT(x) finite(x)
883 # define isfinite finite
884 #endif
885
886 /*--------------------------------------------------------------------*/
887
888 #define CHECK_SVD
889
890 #undef CHK
891 #ifdef CHECK_SVD
892 # define CHK 1
893 # define A(i,j) aa[(i)+(j)*m]
894 # define U(i,j) uu[(i)+(j)*m]
895 # define V(i,j) vv[(i)+(j)*n]
896 #else
897 # define CHK 0
898 #endif
899
900 /** setup for sorting SVD values:
901 0 = no sort (whatever the function returns)
902 +1 = sort in increasing order of singular values
903 -1 = sort in descending order of singular values **/
904
905 static int svd_sort = 0 ;
DR_set_svd_sort(int ss)906 static void DR_set_svd_sort( int ss ){ svd_sort = ss; }
907
908 /*----------------------------------------------------------------------------*/
909 /*! Compute SVD of double precision matrix: T
910 [a] = [u] diag[s] [v]
911 - m = # of rows in a = length of each column
912 - n = # of columns in a = length of each row
913 - a = pointer to input matrix; a[i+j*m] has the (i,j) element
914 (m X n matrix, stored in column-first order)
915 - s = pointer to output singular values; length = n (cannot be NULL)
916 - u = pointer to output matrix, if desired; length = m*n (m X n matrix)
917 - v = pointer to output matrix, if desired; length = n*n (n x n matrix)
918
919 Modified 10 Jan 2007 to add sorting of s and corresponding columns of u & v.
920 ------------------------------------------------------------------------------*/
921
DR_svd_double(int m,int n,double * a,double * s,double * u,double * v)922 static void DR_svd_double( int m, int n, double *a, double *s, double *u, double *v )
923 {
924 integer mm,nn , lda,ldu,ldv , ierr ;
925 doublereal *aa, *ww , *uu , *vv , *rv1 ;
926 logical matu , matv ;
927
928 if( a == NULL || s == NULL || m < 1 || n < 1 ) return ;
929
930 mm = m ;
931 nn = n ;
932 aa = a ;
933 lda = m ;
934 ww = s ;
935
936 /* make space for u matrix, if not supplied */
937
938 if( u == NULL ){
939 matu = (logical) CHK ;
940 uu = (doublereal *)calloc(sizeof(double),m*n) ;
941 } else {
942 matu = (logical) 1 ;
943 uu = u ;
944 }
945 ldu = m ;
946
947 /* make space for v matrix if not supplied */
948
949 if( v == NULL ){
950 matv = (logical) CHK ;
951 vv = (CHK) ? (doublereal *)calloc(sizeof(double),n*n) : NULL ;
952 } else {
953 matv = (logical) 1 ;
954 vv = v ;
955 }
956 ldv = n ;
957
958 rv1 = (double *)calloc(sizeof(double),n) ; /* workspace */
959
960 /** the actual SVD **/
961
962 (void) svd_( &mm , &nn , &lda , aa , ww ,
963 &matu , &ldu , uu , &matv , &ldv , vv , &ierr , rv1 ) ;
964
965 #ifdef CHECK_SVD
966 /** back-compute [A] from [U] diag[ww] [V]'
967 and see if it is close to the input matrix;
968 if not, compute the results in another function;
969 this is needed because the svd() function compiles with
970 rare computational errors on some compilers' optimizers **/
971 { register int i,j,k ; register doublereal aij ; double err=0.0,amag=1.e-11 ;
972 for( j=0 ; j < n ; j++ ){
973 for( i=0 ; i < m ; i++ ){
974 aij = A(i,j) ; amag += fabs(aij) ;
975 for( k=0 ; k < n ; k++ ) aij -= U(i,k)*V(j,k)*ww[k] ;
976 err += fabs(aij) ;
977 }}
978 amag /= (m*n) ; /* average absolute value of matrix elements */
979 err /= (m*n) ; /* average absolute error per matrix element */
980 if( err >= 1.e-5*amag || !IS_GOOD_FLOAT(err) ){
981 fprintf(stderr,"SVD avg err=%g; recomputing ...",err) ;
982
983 #if 1 /* mangle all zero columns */
984 { double arep=1.e-11*amag , *aj ;
985 for( j=0 ; j < nn ; j++ ){
986 aj = aa + j*mm ;
987 for( i=0 ; i < mm ; i++ ) if( aj[i] != 0.0 ) break ;
988 if( i == mm ){
989 ININFO_message("mangling all zero column %d",j) ;
990 for( i=0 ; i < mm ; i++ ) aj[i] = (drand48()-0.5)*arep ;
991 }
992 }
993 }
994 #endif
995
996 /* svd_slow is compiled without optimization */
997
998 (void) svd_slow_( &mm , &nn , &lda , aa , ww ,
999 &matu , &ldu , uu , &matv , &ldv , vv , &ierr , rv1 ) ;
1000 err = 0.0 ;
1001 for( j=0 ; j < n ; j++ ){
1002 for( i=0 ; i < m ; i++ ){
1003 aij = A(i,j) ;
1004 for( k=0 ; k < n ; k++ ) aij -= U(i,k)*V(j,k)*ww[k] ;
1005 err += fabs(aij) ;
1006 }}
1007 err /= (m*n) ;
1008 fprintf(stderr," new avg err=%g %s\n",
1009 err , (err >= 1.e-5*amag || !IS_GOOD_FLOAT(err)) ? "**BAD**" : "**OK**" ) ;
1010 }
1011 }
1012 #endif
1013
1014 free((void *)rv1) ;
1015
1016 /* discard [u] and [v] spaces if not needed for output */
1017
1018 if( u == NULL && uu != NULL ) free((void *)uu) ;
1019 if( v == NULL && vv != NULL ) free((void *)vv) ;
1020
1021 /*--- 10 Jan 2007: sort the singular values and columns of U and V ---*/
1022
1023 if( n > 1 && svd_sort != 0 ){
1024 double *sv ; int *iv , jj,kk ;
1025 sv = (double *)malloc(sizeof(double)*n) ;
1026 iv = (int *) malloc(sizeof(int) *n) ;
1027 for( kk=0 ; kk < n ; kk++ ){
1028 iv[kk] = kk ; sv[kk] = (svd_sort > 0) ? s[kk] : -s[kk] ;
1029 }
1030 ININFO_message("sorting sv") ;
1031 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
1032 qsort_doubleint( n , sv , iv ) ;
1033 ININFO_message(" -- sorted") ;
1034 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
1035 if( u != NULL ){
1036 double *cc ;
1037 ININFO_message("malloc-ing cc: n=%d",n) ;
1038 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
1039 cc = (double *)calloc(sizeof(double),m*n) ;
1040 ININFO_message("copying u") ;
1041 AAmemcpy( cc , u , sizeof(double)*m*n ) ;
1042 for( jj=0 ; jj < n ; jj++ ){
1043 ININFO_message(" u[%d] <- cc[%d]",jj,kk) ;
1044 kk = iv[jj] ; /* where the new jj-th col came from */
1045 AAmemcpy( u+jj*m , cc+kk*m , sizeof(double)*m ) ;
1046 }
1047 ININFO_message("freeing cc") ;
1048 free((void *)cc) ;
1049 ININFO_message("cc is freed") ;
1050 }
1051 if( v != NULL ){
1052 double *cc ;
1053 ININFO_message("malloc-ing cc: n=%d",n) ;
1054 ININFO_message("MCW_MALLOC_status = %s",MCW_MALLOC_status) ;
1055 cc = (double *)calloc(sizeof(double),n*n) ;
1056 ININFO_message("copying cc <- v") ;
1057 AAmemcpy( cc , v , sizeof(double)*n*n ) ;
1058 for( jj=0 ; jj < n ; jj++ ){
1059 kk = iv[jj] ;
1060 AAmemcpy( v+jj*n , cc+kk*n , sizeof(double)*n ) ;
1061 }
1062 free((void *)cc) ;
1063 }
1064 ININFO_message("getting s back") ;
1065 for( kk=0 ; kk < n ; kk++ )
1066 s[kk] = (svd_sort > 0) ? sv[kk] : -sv[kk] ;
1067 free((void *)iv) ; free((void *)sv) ;
1068 }
1069
1070 return ;
1071 }
1072