Doxygen Source Code Documentation
3dpc.c File Reference
#include <string.h>#include "mrilib.h"#include <stdlib.h>#include <ctype.h>#include "cs.h"Go to the source code of this file.
Defines | |
| #define | PC_lprin_calc PC_brnum |
| #define | XX(i, j) PC_brickdata[(j)][(i)] |
| #define | AA(i, j) (aa[(i)+(j)*adim]) |
| #define | VV(i, j) (zout[(i)+(j)*adim]) |
Functions | |
| void | PC_read_opts (int, char **) |
| void | PC_syntax (char *) |
| void | PC_read_opts (int argc, char *argv[]) |
| int | main (int argc, char *argv[]) |
Variables | |
| int | PC_dmean = 0 |
| int | PC_vmean = 0 |
| int | PC_vnorm = 0 |
| int | PC_normalize = 0 |
| int | PC_lprin_save = 0 |
| int | PC_be_quiet = 1 |
| int | PC_do_float = 0 |
| char ** | PC_dsname = NULL |
| int | PC_dsnum = 0 |
| int | PC_brnum = 0 |
| int | PC_1ddum = 0 |
| THD_3dim_dataset ** | PC_dset = NULL |
| char | PC_prefix [THD_MAX_PREFIX] = "pc" |
| float ** | PC_brickdata |
| byte * | PC_mask = NULL |
| int | PC_mask_nvox = 0 |
| int | PC_mask_hits = 0 |
Define Documentation
|
|
i,j element of covariance matrix * Definition at line 49 of file 3dpc.c. Referenced by main(). |
|
|
Definition at line 29 of file 3dpc.c. Referenced by main(). |
|
|
i'th element of j'th eigenvector * Definition at line 53 of file 3dpc.c. Referenced by apply_yshear(), apply_zshear(), flip_xy(), flip_xz(), flip_yz(), and main(). |
|
|
i'th element of j'th input brick * Definition at line 45 of file 3dpc.c. Referenced by huber_func(), linear_filter_extend(), linear_filter_trend(), linear_filter_zero(), and main(). |
Function Documentation
|
||||||||||||
|
---------- Adapted from 3dZeropad.c by RWCox - 08 Aug 2001 ----------* Definition at line 343 of file 3dpc.c. References AA, addto_args(), ADN_func_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, AFNI_logger(), ANAT_BUCK_TYPE, argc, DSET_BRIKNAME, DSET_delete, DSET_HEADNAME, DSET_NX, DSET_NY, DSET_NZ, DSET_write, EDIT_BRICK_LABEL, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), fout, free, FUNC_BUCK_TYPE, ISANAT, machdep(), mainENTRY, malloc, MRI_FLOAT_PTR, mri_free(), mri_new(), mri_write_ascii(), nz, PC_1ddum, PC_be_quiet, PC_brickdata, PC_brnum, PC_dsnum, PC_lprin_calc, PC_lprin_save, PC_mask, PC_mask_hits, PC_prefix, PC_read_opts(), PC_syntax(), PRINT_VERSION, symeig_double(), THD_is_file(), THD_MAX_NAME, tross_Copy_History(), tross_Make_History(), VV, and XX.
00344 {
00345 int nx,ny,nz , nxyz , ii,jj,ll , nn,mm,mmmm , npos,nneg ;
00346 float fmax , ftem ;
00347 THD_3dim_dataset * dset , * new_dset ;
00348 double * dsdev = NULL ;
00349 float * fout , * perc ;
00350 short * bout ;
00351 register float sum ;
00352 register double dsum ;
00353 register int kk ;
00354 int ifirst,ilast,idel ;
00355 int ierror; /* number of errors in editing data */
00356 MRI_IMAGE * vecim ;
00357 char vname[THD_MAX_NAME] ;
00358
00359 /** data for eigenvalue routine (some only used with LAPACK) **/
00360
00361 double * aa , * wout , * zout , * work ;
00362 double abstol , atrace ;
00363 int adim , il , iu , mout , lwork , info ;
00364 int * iwork , * ifail ;
00365
00366 /*-- read command line arguments --*/
00367
00368 if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) PC_syntax(NULL) ;
00369
00370 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
00371
00372 mainENTRY("3dpc main"); machdep(); PRINT_VERSION("3dpc") ;
00373
00374 { int new_argc ; char ** new_argv ;
00375 addto_args( argc , argv , &new_argc , &new_argv ) ;
00376 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00377 }
00378
00379 AFNI_logger("3dpc",argc,argv) ;
00380
00381 PC_read_opts( argc , argv ) ;
00382
00383 /*-- get dimensions --*/
00384
00385 dset = PC_dset[0] ;
00386 nx = DSET_NX(dset) ;
00387 ny = DSET_NY(dset) ;
00388 nz = DSET_NZ(dset) ;
00389 nxyz = nx * ny * nz ;
00390
00391 nn = nxyz ; /* vector length */
00392 mm = PC_brnum ; /* number of vectors */
00393
00394 /*-- space for eigenvalue computations --*/
00395
00396 adim = mm ;
00397 aa = (double *) malloc( sizeof(double) * adim * adim ) ; /* matrix */
00398 wout = (double *) malloc( sizeof(double) * adim ) ; /* evals */
00399
00400 if( aa == NULL || wout == NULL )
00401 PC_syntax("can't malloc space for covariance matrix") ;
00402
00403 #ifdef USE_LAPACK
00404 il = adim + 1 - PC_lprin_calc ; /* lowest index */
00405 iu = adim ; /* upper index */
00406 abstol = 0.0 ;
00407 zout = (double *) malloc( sizeof(double) * adim * PC_lprin_calc ) ;
00408 lwork = 32 * adim ;
00409 work = (double *) malloc( sizeof(double) * lwork ) ;
00410 iwork = (int *) malloc( sizeof(int) * 6 * adim ) ;
00411 ifail = (int *) malloc( sizeof(int) * adim ) ;
00412
00413 if( zout == NULL || work == NULL || iwork == NULL || ifail == NULL )
00414 PC_syntax("can't malloc eigen workspace") ;
00415 #endif
00416
00417 /*-- remove means, if ordered --*/
00418
00419 if( PC_vmean ){
00420 float * vxmean = (float *) malloc( sizeof(float) * nn ) ; /* voxel means */
00421
00422 if( !PC_be_quiet ){ printf("--- remove timeseries means"); fflush(stdout); }
00423
00424 for( kk=0 ; kk < nn ; kk++ ) vxmean[kk] = 0.0 ;
00425
00426 for( jj=0 ; jj < mm ; jj++ ){
00427 for( kk=0 ; kk < nn ; kk++ ) vxmean[kk] += XX(kk,jj) ;
00428 if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00429 }
00430
00431 sum = 1.0 / mm ;
00432 for( kk=0 ; kk < nn ; kk++ ) vxmean[kk] *= sum ;
00433
00434 for( jj=0 ; jj < mm ; jj++ ){
00435 for( kk=0 ; kk < nn ; kk++ ) XX(kk,jj) -= vxmean[kk] ;
00436 if( !PC_be_quiet ){ printf("-"); fflush(stdout); }
00437 }
00438
00439 free(vxmean) ;
00440 if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00441
00442 } else if( PC_dmean ){
00443 if( !PC_be_quiet ){ printf("--- remove brick means"); fflush(stdout); }
00444
00445 if( PC_mask == NULL ){
00446 for( jj=0 ; jj < mm ; jj++ ){
00447 sum = 0.0 ;
00448 for( kk=0 ; kk < nn ; kk++ ) sum += XX(kk,jj) ;
00449 if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00450 sum /= nn ;
00451 for( kk=0 ; kk < nn ; kk++ ) XX(kk,jj) -= sum ;
00452 if( !PC_be_quiet ){ printf("-"); fflush(stdout); }
00453 }
00454 } else { /* 15 Sep 1999 */
00455 for( jj=0 ; jj < mm ; jj++ ){
00456 sum = 0.0 ;
00457 for( kk=0 ; kk < nn ; kk++ ) if( PC_mask[kk] ) sum += XX(kk,jj) ;
00458 if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00459 sum /= PC_mask_hits ;
00460 for( kk=0 ; kk < nn ; kk++ ) if( PC_mask[kk] ) XX(kk,jj) -= sum ;
00461 if( !PC_be_quiet ){ printf("-"); fflush(stdout); }
00462 }
00463 }
00464 if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00465 }
00466
00467 /*-- 07 July 1999: vnorm --*/
00468
00469 if( PC_vnorm ){
00470 float * vxnorm = (float *) malloc( sizeof(float) * nn ) ; /* voxel norms */
00471
00472 if( !PC_be_quiet ){ printf("--- normalize timeseries"); fflush(stdout); }
00473
00474 for( kk=0 ; kk < nn ; kk++ ) vxnorm[kk] = 0.0 ;
00475
00476 for( jj=0 ; jj < mm ; jj++ ){
00477 for( kk=0 ; kk < nn ; kk++ ) vxnorm[kk] += XX(kk,jj) * XX(kk,jj) ;
00478 if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00479 }
00480
00481 for( kk=0 ; kk < nn ; kk++ )
00482 if( vxnorm[kk] > 0.0 ) vxnorm[kk] = 1.0 / sqrt(vxnorm[kk]) ;
00483
00484 for( jj=0 ; jj < mm ; jj++ ){
00485 for( kk=0 ; kk < nn ; kk++ ) XX(kk,jj) *= vxnorm[kk] ;
00486 if( !PC_be_quiet ){ printf("*"); fflush(stdout); }
00487 }
00488
00489 free(vxnorm) ;
00490 if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00491 }
00492
00493 /*-- load covariance matrix
00494 (very short code that takes a long time to run!) --*/
00495
00496 if( !PC_be_quiet ){ printf("--- compute covariance matrix"); fflush(stdout); }
00497
00498 idel = 1 ; /* ii goes forward */
00499 for( jj=0 ; jj < mm ; jj++ ){
00500
00501 ifirst = (idel==1) ? 0 : jj ; /* back and forth in ii to */
00502 ilast = (idel==1) ? jj+1 : -1 ; /* maximize use of cache/RAM */
00503
00504 for( ii=ifirst ; ii != ilast ; ii += idel ){
00505 dsum = 0.0 ;
00506 if( PC_mask == NULL ){
00507 for( kk=0 ; kk < nn ; kk++ ) dsum += XX(kk,ii) * XX(kk,jj) ;
00508 } else {
00509 for( kk=0 ; kk < nn ; kk++ )
00510 if( PC_mask[kk] ) dsum += XX(kk,ii) * XX(kk,jj) ;
00511 }
00512 AA(ii,jj) = AA(jj,ii) = dsum ;
00513 }
00514
00515 if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00516
00517 idel = -idel ; /* reverse direction of ii */
00518 }
00519 if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00520
00521 /*-- check diagonal for OK-ness --**/
00522
00523 atrace = 0.0 ;
00524 ii = 0 ;
00525 for( jj=0 ; jj < mm ; jj++ ){
00526 if( AA(jj,jj) <= 0.0 ){
00527 fprintf(stderr,"*** covariance diagonal (%d,%d) = %g\n",
00528 jj+1,jj+1,AA(jj,jj)) ;
00529 ii++ ;
00530 }
00531 atrace += AA(jj,jj) ;
00532 }
00533 if( ii > 0 ){ printf("*** program exiting right here and now!\n"); exit(1); }
00534
00535 if( !PC_be_quiet ){ printf("--- covariance trace = %g\n",atrace); fflush(stdout); }
00536
00537 /*-- normalize, if desired --*/
00538
00539 if( PC_normalize ){
00540 if( !PC_be_quiet ){ printf("--- normalizing covariance"); fflush(stdout); }
00541
00542 dsdev = (double *) malloc( sizeof(double) * mm ) ; /* brick stdev */
00543
00544 for( jj=0 ; jj < mm ; jj++ ) dsdev[jj] = sqrt( AA(jj,jj) ) ;
00545
00546 for( jj=0 ; jj < mm ; jj++ )
00547 for( ii=0 ; ii < mm ; ii++ ) AA(ii,jj) /= (dsdev[ii]*dsdev[jj]) ;
00548
00549 atrace = mm ;
00550
00551 if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00552 }
00553
00554 if( !PC_be_quiet ){ printf("--- compute eigensolution\n"); fflush(stdout); }
00555
00556 #ifdef USE_LAPACK
00557 (void) dsyevx_( "V" , /* eigenvalues and vectors */
00558 "I" , /* a subrange of eigenvalues */
00559 "U" , /* use upper triangle of A */
00560 &adim , /* dimension of A */
00561 aa , /* the matrix A */
00562 &adim , /* leading dimension of A */
00563 NULL , /* not used */
00564 NULL , /* not used */
00565 &il , /* lowest eigen-index */
00566 &iu , /* highest eigen-index */
00567 &abstol , /* tolerance */
00568 &mout , /* output # of eigenvalues */
00569 wout , /* output eigenvalues */
00570 zout , /* output eigenvectors */
00571 &adim , /* leading dimension of zout */
00572 work , /* double work array */
00573 &lwork , /* size of work array */
00574 iwork , /* another work array */
00575 ifail , /* output failure list */
00576 &info /* results code */
00577 ) ;
00578
00579 free(aa) ; free(work) ; free(iwork) ; free(ifail) ;
00580
00581 if( info != 0 ){
00582 fprintf(stderr,"*** DSYEVX returns error code info=%d\n",info);
00583 if( info < 0 ) exit(1) ;
00584 }
00585 #else
00586 symeig_double( mm , aa , wout ) ; /* eigenvectors go over aa */
00587 zout = aa ; /* eigenvalues go into wout */
00588 #endif
00589
00590 if( !PC_be_quiet ) printf("\n") ;
00591
00592 sum = 0.0 ;
00593
00594 perc = (float *) malloc( sizeof(float) * PC_lprin_calc ) ;
00595
00596 printf("Num. --Eigenvalue-- -Var.Fraction- -Cumul.Fract.-\n") ;
00597 for( jj=0 ; jj < PC_lprin_calc ; jj++ ){
00598 ll = PC_lprin_calc - 1-jj ; /* reversed order of eigensolution! */
00599 perc[jj] = wout[ll]/atrace ;
00600 sum += perc[jj] ;
00601 printf("%4d %14.7g %14.7g %14.7g\n",
00602 jj+1 , wout[ll] , perc[jj] , sum ) ;
00603 }
00604 fflush(stdout) ;
00605
00606 /*--- form and save output dataset ---*/
00607
00608 dset = PC_dset[0] ;
00609 new_dset = EDIT_empty_copy( dset ) ;
00610
00611 if( PC_dsnum == 1 ) tross_Copy_History( dset , new_dset ) ;
00612 tross_Make_History( "3dpc" , argc,argv , new_dset ) ;
00613
00614 EDIT_dset_items( new_dset,
00615 ADN_prefix , PC_prefix ,
00616 ADN_nvals , PC_lprin_save ,
00617 ADN_ntt , 0 ,
00618 ADN_func_type , ISANAT(dset) ? ANAT_BUCK_TYPE
00619 : FUNC_BUCK_TYPE ,
00620 ADN_none ) ;
00621
00622 if( THD_is_file(DSET_HEADNAME(new_dset)) ){
00623 fprintf(stderr,
00624 "\n*** Output dataset %s already exists--will be destroyed!\n",
00625 DSET_HEADNAME(new_dset) ) ;
00626
00627 } else if( !PC_be_quiet ){
00628 printf("--- output dataset %s" , DSET_BRIKNAME(new_dset) ) ;
00629 fflush(stdout) ;
00630 }
00631
00632 fout = (float *) malloc( sizeof(float) * nn ) ; /* output buffer */
00633
00634 for( jj=0 ; jj < PC_lprin_save ; jj++ ){
00635 ll = PC_lprin_calc - 1-jj ;
00636
00637 /** output = weighted sum of input datasets,
00638 with weights from the ll'th eigenvector **/
00639
00640 for( kk=0 ; kk < nn ; kk++ ) fout[kk] = 0.0 ;
00641
00642 for( ii=0 ; ii < mm ; ii++ ){
00643 sum = VV(ii,ll) ; if( PC_normalize ) sum /= dsdev[ii] ;
00644
00645 for( kk=0 ; kk < nn ; kk++ ) fout[kk] += XX(kk,ii) * sum ;
00646 }
00647
00648 fmax = 0.0 ; npos = nneg = 0 ;
00649 for( kk=0 ; kk < nn ; kk++ ){
00650 ftem = fabs(fout[kk]) ; if( fmax < ftem ) fmax = ftem ;
00651 if( fout[kk] > 0 ) npos++ ;
00652 else if( fout[kk] < 0 ) nneg++ ;
00653 }
00654
00655 if( PC_do_float ){
00656 if( nneg > npos )
00657 for( kk=0 ; kk < nn ; kk++ ) fout[kk] = -fout[kk] ;
00658
00659 EDIT_substitute_brick( new_dset , jj , MRI_float , fout ) ;
00660
00661 fout = (float *) malloc( sizeof(float) * nn ) ; /* new buffer */
00662 } else {
00663
00664 bout = (short *) malloc( sizeof(short) * nn ) ; /* output buffer */
00665 if( fmax != 0.0 ){
00666 fmax = 10000.49/fmax ; if( nneg > npos ) fmax = -fmax ;
00667 for( kk=0 ; kk < nn ; kk++ ) bout[kk] = fmax * fout[kk] ;
00668 } else {
00669 for( kk=0 ; kk < nn ; kk++ ) bout[kk] = 0.0 ;
00670 }
00671 EDIT_substitute_brick( new_dset , jj , MRI_short , bout ) ;
00672 }
00673
00674 sprintf(vname,"var=%6.3f%%" , 100.0*perc[jj]+0.499 ) ;
00675 EDIT_BRICK_LABEL( new_dset , jj , vname ) ;
00676
00677 if( !PC_be_quiet ){ printf(".") ; fflush(stdout); }
00678 }
00679 free(fout) ;
00680
00681 DSET_write(new_dset) ;
00682 fprintf(stderr,"++ output dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
00683 DSET_delete(new_dset) ;
00684 if( !PC_be_quiet ){ printf("!\n") ; fflush(stdout); }
00685
00686 /*-- write eigenvectors also --*/
00687
00688 mmmm = mm + PC_1ddum ;
00689 vecim = mri_new( PC_lprin_save , mmmm , MRI_float ) ;
00690 fout = MRI_FLOAT_PTR(vecim) ;
00691 for( jj=0 ; jj < PC_lprin_save ; jj++ ){
00692 ll = PC_lprin_calc - 1-jj ;
00693 for( ii=0 ; ii < PC_1ddum ; ii++ )
00694 fout[jj + ii*PC_lprin_save] = 999999.0 ;
00695 for( ii=0 ; ii < mm ; ii++ )
00696 fout[jj + (ii+PC_1ddum)*PC_lprin_save] = VV(ii,ll) ;
00697 }
00698 sprintf(vname,"%s.1D",PC_prefix) ;
00699 mri_write_ascii( vname, vecim ) ;
00700 mri_free(vecim) ;
00701
00702 for( jj=0 ; jj < PC_lprin_save ; jj++ ){
00703 ll = PC_lprin_calc - 1-jj ;
00704 vecim = mri_new( 1 , mmmm , MRI_float ) ;
00705 fout = MRI_FLOAT_PTR(vecim) ;
00706 for( ii=0 ; ii < PC_1ddum ; ii++ ) fout[ii] = 999999.0 ;
00707 for( ii=0 ; ii < mm ; ii++ ) fout[ii+PC_1ddum] = VV(ii,ll) ;
00708 sprintf(vname,"%s%02d.1D",PC_prefix,jj) ;
00709 mri_write_ascii( vname, vecim ) ;
00710 mri_free(vecim) ;
00711 }
00712
00713 #if 0
00714 free(PC_dset) ;
00715 for( ii=0 ; ii < mm ; ii++ ) free(PC_brickdata[ii]) ;
00716 free(PC_brickdata) ;
00717 free(aa) ; free(wout) ; free(perc) ; if( dsdev!=NULL ) free(dsdev) ;
00718 #endif
00719
00720 exit(0) ;
00721 }
|
|
||||||||||||
|
Use EISPACK instead Definition at line 81 of file 3dpc.c. References argc, DSET_ARRAY, DSET_BRICK_TYPE, DSET_delete, DSET_load, DSET_LOADED, DSET_NVALS, DSET_NVOX, DSET_NX, DSET_NY, DSET_NZ, DSET_unload, DSET_unload_one, EDIT_coerce_type(), free, ISVALID_3DIM_DATASET, malloc, mc, MCW_strncpy, nz, PC_1ddum, PC_be_quiet, PC_brickdata, PC_brnum, PC_dmean, PC_do_float, PC_dsname, PC_dsnum, PC_lprin_save, PC_mask, PC_mask_hits, PC_mask_nvox, PC_normalize, PC_prefix, PC_syntax(), PC_vmean, PC_vnorm, THD_makemask(), THD_MAX_PREFIX, and THD_open_dataset(). Referenced by main().
00082 {
00083 int nopt = 1 ;
00084 float val ;
00085 int kk, nx, ny, nz, nxyz, mm,nn ;
00086 THD_3dim_dataset * dset ;
00087
00088 while( nopt < argc && argv[nopt][0] == '-' ){
00089
00090 /**** -verbose ****/
00091
00092 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00093 PC_be_quiet = 0 ;
00094 nopt++ ; continue ;
00095 }
00096
00097 /**** -float ****/
00098
00099 if( strncmp(argv[nopt],"-float",6) == 0 ){
00100 PC_do_float = 1 ;
00101 nopt++ ; continue ;
00102 }
00103
00104 /**** -prefix prefix ****/
00105
00106 if( strncmp(argv[nopt],"-prefix",6) == 0 ){
00107 nopt++ ;
00108 if( nopt >= argc ) PC_syntax("need argument after -prefix!") ;
00109 MCW_strncpy( PC_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00110 continue ;
00111 }
00112
00113 /**** -dmean ****/
00114
00115 if( strncmp(argv[nopt],"-dmean",6) == 0 ){
00116 PC_dmean = 1 ;
00117 nopt++ ; continue ;
00118 }
00119
00120 /**** -vmean ****/
00121
00122 if( strncmp(argv[nopt],"-vmean",6) == 0 ){
00123 PC_vmean = 1 ;
00124 nopt++ ; continue ;
00125 }
00126
00127 /**** -vnorm ****/
00128
00129 if( strncmp(argv[nopt],"-vnorm",6) == 0 ){
00130 PC_vnorm = 1 ;
00131 nopt++ ; continue ;
00132 }
00133
00134 /**** -normalize ****/
00135
00136 if( strncmp(argv[nopt],"-normalize",6) == 0 ){
00137 PC_normalize = 1 ;
00138 nopt++ ; continue ;
00139 }
00140
00141 /**** -pcsave # ****/
00142
00143 if( strncmp(argv[nopt],"-pcsave",6) == 0 ){
00144 nopt++ ;
00145 if( nopt >= argc ) PC_syntax("need argument after -pcsave!") ;
00146 PC_lprin_save = strtol( argv[nopt] , NULL , 10 ) ;
00147 if( PC_lprin_save < 0 ) PC_syntax("value after -pcsave is illegal!") ;
00148 nopt++ ; continue ;
00149 }
00150
00151 /**** -1ddum # ****/
00152
00153 if( strncmp(argv[nopt],"-1ddum",6) == 0 ){
00154 nopt++ ;
00155 if( nopt >= argc ) PC_syntax("need argument after -1ddum!") ;
00156 PC_1ddum = strtol( argv[nopt] , NULL , 10 ) ;
00157 if( PC_1ddum < 0 ) PC_syntax("value after -1ddum is illegal!") ;
00158 nopt++ ; continue ;
00159 }
00160
00161 /**** -mask mset [15 Sep 1999] ****/
00162
00163 if( strncmp(argv[nopt],"-mask",5) == 0 ){
00164 THD_3dim_dataset * mset ; int ii,mc ;
00165 nopt++ ;
00166 if( nopt >= argc ) PC_syntax("need argument after -mask!") ;
00167 mset = THD_open_dataset( argv[nopt] ) ;
00168 if( mset == NULL ) PC_syntax("can't open -mask dataset!") ;
00169 PC_mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
00170 PC_mask_nvox = DSET_NVOX(mset) ;
00171 DSET_delete(mset) ;
00172 if( PC_mask == NULL ) PC_syntax("can't use -mask dataset!") ;
00173 for( ii=mc=0 ; ii < PC_mask_nvox ; ii++ ) if( PC_mask[ii] ) mc++ ;
00174 if( mc == 0 ) PC_syntax("mask is all zeros!") ;
00175 printf("--- %d voxels in mask\n",mc) ;
00176 PC_mask_hits = mc ;
00177 nopt++ ; continue ;
00178 }
00179
00180 /**** unknown switch ****/
00181
00182 fprintf(stderr,"\n*** unrecognized option %s\n",argv[nopt]) ;
00183 exit(1) ;
00184
00185 } /* end of loop over options */
00186
00187 /*--- a simple consistency check ---*/
00188
00189 if( PC_vmean && PC_dmean )
00190 PC_syntax("can't have both -dmean and -vmean!") ;
00191
00192 /*--- rest of inputs are dataset names ---*/
00193
00194 PC_dsnum = argc - nopt ;
00195 if( PC_dsnum < 1 ) PC_syntax("no input dataset names?") ;
00196
00197 PC_dsname = (char **) malloc( sizeof(char *) * PC_dsnum ) ;
00198 for( kk=0 ; kk < PC_dsnum ; kk++ ) PC_dsname[kk] = argv[kk+nopt] ;
00199
00200 PC_dset = (THD_3dim_dataset **) malloc( sizeof(THD_3dim_dataset *) * PC_dsnum ) ;
00201 for( kk=0 ; kk < PC_dsnum ; kk++ ){
00202 PC_dset[kk] = dset = THD_open_dataset( PC_dsname[kk] ) ; /* allow for selector */
00203 if( !ISVALID_3DIM_DATASET(dset) ){
00204 fprintf(stderr,"\n*** can't open dataset file %s\n",PC_dsname[kk]) ;
00205 exit(1) ;
00206 }
00207 PC_brnum += DSET_NVALS(dset) ;
00208 }
00209 if( PC_brnum < 2 ) PC_syntax("need at least 2 input bricks!") ;
00210
00211 /*--- another consistency check ---*/
00212
00213 if( PC_lprin_save <= 0 )
00214 PC_lprin_save = PC_brnum ;
00215 else if( PC_lprin_save > PC_brnum )
00216 PC_syntax("can't save more components than input bricks!") ;
00217
00218 /*--- load bricks for all input datasets ---*/
00219
00220 nx = DSET_NX(PC_dset[0]) ;
00221 ny = DSET_NY(PC_dset[0]) ;
00222 nz = DSET_NZ(PC_dset[0]) ;
00223 nxyz = nx * ny * nz ; /* Total number of voxels per brick */
00224
00225 /* 15 Sep 1999: check if mask is right size */
00226
00227 if( PC_mask_nvox > 0 && PC_mask_nvox != nxyz )
00228 PC_syntax("mask and input dataset bricks don't match in size!") ;
00229
00230 PC_brickdata = (float **) malloc( sizeof(float *) * PC_brnum ) ;
00231
00232 nn = 0 ; /* current brick index */
00233
00234 if( !PC_be_quiet ){ printf("--- read dataset bricks"); fflush(stdout); }
00235
00236 for( kk=0 ; kk < PC_dsnum ; kk++ ){
00237
00238 if( DSET_NVOX(PC_dset[kk]) != nxyz ) {
00239 fprintf(stderr,
00240 "\n*** dataset in file %s nonconformant with first dataset!\n"
00241 " nx1=%d ny1=%d nz1=%d nx=%d ny=%d nz=%d\n",
00242 PC_dsname[kk], nx, ny, nz ,
00243 DSET_NX(PC_dset[kk]), DSET_NY(PC_dset[kk]), DSET_NZ(PC_dset[kk]) ) ;
00244 exit(1) ;
00245 }
00246
00247 DSET_load( PC_dset[kk] ) ;
00248 if( !DSET_LOADED( PC_dset[kk] ) ){
00249 fprintf(stderr,"\n*** Can't load dataset %s BRIK from disk!\n",PC_dsname[kk]) ;
00250 exit(1) ;
00251 }
00252 if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00253
00254 /* copy brick data into float storage */
00255
00256 for( mm=0 ; mm < DSET_NVALS(PC_dset[kk]) ; mm++,nn++ ){
00257
00258 PC_brickdata[nn] = (float *) malloc( sizeof(float) * nxyz ) ;
00259
00260 if( PC_brickdata[nn] == NULL )
00261 PC_syntax("*** can't malloc intermediate storage") ;
00262
00263 EDIT_coerce_type( nxyz , DSET_BRICK_TYPE(PC_dset[kk],mm) ,
00264 DSET_ARRAY(PC_dset[kk],mm) ,
00265 MRI_float , PC_brickdata[nn] ) ;
00266
00267 DSET_unload_one( PC_dset[kk] , mm ) ;
00268
00269 if( PC_mask != NULL ){ /* 15 Sep 1999 */
00270 int kk ;
00271 for( kk=0 ; kk < nxyz ; kk++ )
00272 if( !PC_mask[kk] ) PC_brickdata[nn][kk] = 0.0 ;
00273 }
00274
00275 if( !PC_be_quiet ){ printf("."); fflush(stdout); }
00276 }
00277
00278 if( kk == 0 ){
00279 DSET_unload( PC_dset[kk] ) ; /* don't need dataset's data anymore */
00280 } else {
00281 DSET_delete( PC_dset[kk] ) ; /* don't need this at all anymore */
00282 PC_dset[kk] = NULL ;
00283 }
00284
00285 }
00286 if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00287
00288 free( PC_dsname ) ; PC_dsname = NULL ;
00289 return ;
00290 }
|
|
||||||||||||
|
|
|
|
Definition at line 294 of file 3dpc.c. References MASTER_SHORTHELP_STRING. Referenced by main(), and PC_read_opts().
00295 {
00296 if( msg != NULL ){ fprintf(stderr,"\n*** %s\n",msg) ; exit(1) ; }
00297
00298 printf(
00299 "Principal Component Analysis of 3D Datasets\n"
00300 "Usage: 3dpc [options] dataset dataset ...\n"
00301 "\n"
00302 "Each input dataset may have a sub-brick selector list.\n"
00303 "Otherwise, all sub-bricks from a dataset will be used.\n"
00304 "\n"
00305 "OPTIONS:\n"
00306 " -dmean = remove the mean from each input brick (across space)\n"
00307 " -vmean = remove the mean from each input voxel (across bricks)\n"
00308 " [N.B.: -dmean and -vmean are mutually exclusive]\n"
00309 " [default: don't remove either mean]\n"
00310 " -vnorm = L2 normalize each input voxel time series\n"
00311 " [occurs after the de-mean operations above,]\n"
00312 " [and before the brick normalization below. ]\n"
00313 " -normalize = L2 normalize each input brick (after mean subtraction)\n"
00314 " [default: don't normalize]\n"
00315 " -pcsave sss = 'sss' is the number of components to save in the output;\n"
00316 " it can't be more than the number of input bricks\n"
00317 " [default = all of them = number of input bricks]\n"
00318 " -prefix pname = Name for output dataset (will be a bucket type);\n"
00319 " also, the eigen-timeseries will be in 'pname'.1D\n"
00320 " (all of them) and in 'pnameNN.1D' for eigenvalue\n"
00321 " #NN individually (NN=00 .. 'sss'-1, corresponding\n"
00322 " to the brick index in the output dataset)\n"
00323 " [default value of pname = 'pc']\n"
00324 " -1ddum ddd = Add 'ddd' dummy lines to the top of each *.1D file.\n"
00325 " These lines will have the value 999999, and can\n"
00326 " be used to align the files appropriately.\n"
00327 " [default value of ddd = 0]\n"
00328 " -verbose = Print progress reports during the computations\n"
00329 " -float = Save eigen-bricks as floats\n"
00330 " [default = shorts, scaled so that |max|=10000]\n"
00331 " -mask mset = Use the 0 sub-brick of dataset 'mset' as a mask\n"
00332 " to indicate which voxels to analyze (a sub-brick\n"
00333 " selector is allowed) [default = use all voxels]\n"
00334 ) ;
00335
00336 printf("\n" MASTER_SHORTHELP_STRING ) ;
00337
00338 exit(0) ;
00339 }
|
Variable Documentation
|
|
Definition at line 27 of file 3dpc.c. Referenced by main(), and PC_read_opts(). |
|
|
Definition at line 21 of file 3dpc.c. Referenced by main(), and PC_read_opts(). |
|
|
Definition at line 35 of file 3dpc.c. Referenced by main(), and PC_read_opts(). |
|
|
Definition at line 26 of file 3dpc.c. Referenced by main(), and PC_read_opts(). |
|
|
inputs * Definition at line 16 of file 3dpc.c. Referenced by PC_read_opts(). |
|
|
Definition at line 22 of file 3dpc.c. Referenced by PC_read_opts(). |
|
|
|
|
|
Definition at line 24 of file 3dpc.c. Referenced by PC_read_opts(). |
|
|
Definition at line 25 of file 3dpc.c. Referenced by main(), and PC_read_opts(). |
|
|
Definition at line 20 of file 3dpc.c. Referenced by main(), and PC_read_opts(). |
|
|
Definition at line 37 of file 3dpc.c. Referenced by main(), and PC_read_opts(). |
|
|
Definition at line 39 of file 3dpc.c. Referenced by main(), and PC_read_opts(). |
|
|
Definition at line 38 of file 3dpc.c. Referenced by PC_read_opts(). |
|
|
Definition at line 19 of file 3dpc.c. Referenced by PC_read_opts(). |
|
|
Definition at line 33 of file 3dpc.c. Referenced by main(), and PC_read_opts(). |
|
|
Definition at line 17 of file 3dpc.c. Referenced by PC_read_opts(). |
|
|
Definition at line 18 of file 3dpc.c. Referenced by PC_read_opts(). |
