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(). |