00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #define PROGRAM_NAME "3dclust"
00015 #define PROGRAM_AUTHOR "RW Cox et al"
00016 #define PROGRAM_DATE "21 Jul 2005"
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 #include <stdio.h>
00064 #include <stdlib.h>
00065
00066 #include "mrilib.h"
00067
00068 static EDIT_options CL_edopt ;
00069 static int CL_ivfim=-1 , CL_ivthr=-1 ;
00070
00071 static int CL_nopt ;
00072
00073 static int CL_summarize = 0 ;
00074
00075 static int CL_noabs = 0;
00076
00077 static int CL_verbose = 0 ;
00078
00079 static int CL_quiet = 0;
00080
00081 static char * CL_prefix = NULL ;
00082
00083 static int CL_do_mni = 0 ;
00084 static int CL_isomode = 0 ;
00085
00086
00087
00088
00089 static THD_coorder CL_cord ;
00090
00091 int compare_cluster( void * , void * ) ;
00092 void CL_read_opts( int , char ** ) ;
00093 #define CL_syntax(str) \
00094 do{ fprintf(stderr,"\n*** %s\a\n",str) ; exit(1) ; } while(0)
00095
00096 void MCW_fc7( float qval , char * buf ) ;
00097
00098
00099
00100 int main( int argc , char * argv[] )
00101 {
00102 float rmm , vmul ;
00103 int iarg ;
00104 THD_3dim_dataset *dset=NULL ;
00105 void * vfim ;
00106 int nx,ny,nz , nxy,nxyz , ivfim ,
00107 iclu , ptmin , ipt , ii,jj,kk , ndet , nopt ;
00108 float dx,dy,dz , xx,yy,zz,mm , xxsum,yysum,zzsum,mmsum , volsum , fimfac ,
00109 xxmax,yymax,zzmax,mmmax , ms, mssum , msmax ,
00110 RLmax, RLmin, APmax, APmin, ISmax, ISmin;
00111 double mean, sem, sqsum, glmmsum, glsqsum, glmssum, glmean, glxxsum, glyysum, glzzsum;
00112 MCW_cluster_array * clar , * clbig ;
00113 MCW_cluster * cl ;
00114 THD_fvec3 fv ;
00115 int nvox_total ;
00116 float vol_total ;
00117 char buf1[16],buf2[16],buf3[16] ;
00118 float dxf,dyf,dzf ;
00119 int do_mni ;
00120
00121 if( argc < 4 || strncmp(argv[1],"-help",4) == 0 ){
00122 printf ("\n\n");
00123 printf ("Program: %s \n", PROGRAM_NAME);
00124 printf ("Author: %s \n", PROGRAM_AUTHOR);
00125 printf ("Date: %s \n", PROGRAM_DATE);
00126 printf ("\n");
00127 printf(
00128 "3dclust - performs simple-minded cluster detection in 3D datasets \n"
00129 " \n"
00130 " This program can be used to find clusters of 'active' voxels and \n"
00131 " print out a report about them. \n"
00132 " * 'Active' refers to nonzero voxels that survive the threshold \n"
00133 " that you (the user) have specified \n"
00134 " * Clusters are defined by a connectivity radius parameter 'rmm' \n"
00135 " \n"
00136 " Note: by default, this program clusters on the absolute values \n"
00137 " of the voxels \n"
00138 "----------------------------------------------------------------------- \n"
00139 "Usage: 3dclust [editing options] [other options] rmm vmul dset ... \n"
00140 "----- \n"
00141 " \n"
00142 "Examples: \n"
00143 "-------- \n"
00144 " \n"
00145 " 3dclust -1clip 0.3 5 2000 func+orig'[1]' \n"
00146 " 3dclust -1noneg -1thresh 0.3 5 2000 func+orig'[1]' \n"
00147 " 3dclust -1noneg -1thresh 0.3 5 2000 func+orig'[1]' func+orig'[3] \n"
00148 " \n"
00149 " 3dclust -noabs -1clip 0.5 -dxyz=1 1 10 func+orig'[1]' \n"
00150 " 3dclust -noabs -1clip 0.5 5 700 func+orig'[1]' \n"
00151 " \n"
00152 " 3dclust -noabs -2clip 0 999 -dxyz=1 1 10 func+orig'[1]' \n"
00153 " \n"
00154 " 3dclust -1clip 0.3 5 3000 func+orig'[1]' \n"
00155 " 3dclust -quiet -1clip 0.3 5 3000 func+orig'[1]' \n"
00156 " 3dclust -summarize -quiet -1clip 0.3 5 3000 func+orig'[1]' \n"
00157 "----------------------------------------------------------------------- \n"
00158 " \n"
00159 "Arguments (must be included on command line): \n"
00160 "--------- \n"
00161 " \n"
00162 " rmm : cluster connection radius (in millimeters). \n"
00163 " All nonzero voxels closer than rmm millimeters \n"
00164 " (center-to-center distance) to the given voxel are \n"
00165 " included in the cluster. \n"
00166 " * If rmm = 0, then clusters are defined by nearest-\n"
00167 " neighbor connectivity \n"
00168 " \n"
00169 " vmul : minimum cluster volume (micro-liters) \n"
00170 " i.e., determines the size of the volume cluster. \n"
00171 " * If vmul = 0, then all clusters are kept. \n"
00172 " * If vmul < 0, then the absolute vmul is the minimum\n"
00173 " number of voxels allowed in a cluster. \n"
00174 " \n"
00175 " dset : input dataset (more than one allowed, but only the \n"
00176 " first sub-brick of the dataset) \n"
00177 " \n"
00178 " The results are sent to standard output (i.e., the screen) \n"
00179 " \n"
00180 "----------------------------------------------------------------------- \n"
00181 " \n"
00182 "Options: \n"
00183 "------- \n"
00184 " \n"
00185 "* Editing options are as in 3dmerge (see 3dmerge -help) \n"
00186 " (including -1thresh, -1dindex, -1tindex, -dxyz=1 options) \n"
00187 " \n"
00188 "* -noabs => Use the signed voxel intensities (not the absolute \n"
00189 " value) for calculation of the mean and Standard \n"
00190 " Error of the Mean (SEM) \n"
00191 " \n"
00192 "* -summarize => Write out only the total nonzero voxel \n"
00193 " count and volume for each dataset \n"
00194 " \n"
00195 "* -nosum => Suppress printout of the totals \n"
00196 " \n"
00197 "* -verb => Print out a progress report (to stderr) \n"
00198 " as the computations proceed \n"
00199 " \n"
00200 "* -quiet => Suppress all non-essential output \n"
00201 " \n"
00202 "* -mni => If the input dataset is in +tlrc coordinates, this \n"
00203 " option will stretch the output xyz-coordinates to the \n"
00204 " MNI template brain. \n"
00205 " \n"
00206 " N.B.1: The MNI template brain is about 5 mm higher (in S), \n"
00207 " 10 mm lower (in I), 5 mm longer (in PA), and tilted \n"
00208 " about 3 degrees backwards, relative to the Talairach- \n"
00209 " Tournoux Atlas brain. For more details, see \n"
00210 " http://www.mrc-cbu.cam.ac.uk/Imaging/mnispace.html \n"
00211 " N.B.2: If the input dataset is not in +tlrc coordinates, \n"
00212 " then the only effect is to flip the output coordinates\n"
00213 " to the 'LPI' (neuroscience) orientation, as if you \n"
00214 " gave the '-orient LPI' option.) \n"
00215 " \n"
00216 "* -isovalue => Clusters will be formed only from contiguous (in the \n"
00217 " rmm sense) voxels that also have the same value. \n"
00218 " \n"
00219 " N.B.: The normal method is to cluster all contiguous \n"
00220 " nonzero voxels together. \n"
00221 " \n"
00222 "* -isomerge => Clusters will be formed from each distinct value \n"
00223 " in the dataset; spatial contiguity will not be \n"
00224 " used (but you still have to supply rmm and vmul \n"
00225 " on the command line). \n"
00226 " \n"
00227 " N.B.: 'Clusters' formed this way may well have components \n"
00228 " that are widely separated! \n"
00229 " \n"
00230 "* -prefix ppp => Write a new dataset that is a copy of the \n"
00231 " input, but with all voxels not in a cluster \n"
00232 " set to zero; the new dataset's prefix is 'ppp' \n"
00233 " \n"
00234 " N.B.: Use of the -prefix option only affects the \n"
00235 " first input dataset \n"
00236 "----------------------------------------------------------------------- \n"
00237 " \n"
00238 "E.g., 3dclust -1clip 0.3 5 3000 func+orig'[1]' \n"
00239 " \n"
00240 " The above command tells 3dclust to find potential cluster volumes for \n"
00241 " dataset func+orig, sub-brick #1, where the threshold has been set \n"
00242 " to 0.3 (i.e., ignore voxels with an activation threshold of >0.3 or \n"
00243 " <-0.3. Voxels must be no more than 5 mm apart, and the cluster volume\n"
00244 " must be at least 3000 micro-liters in size. \n"
00245 " \n"
00246 "Explanation of 3dclust Output: \n"
00247 "----------------------------- \n"
00248 " \n"
00249 " Volume : Number of voxels that make up the volume cluster \n"
00250 " \n"
00251 " CM RL : Center of mass (CM) for the cluster in the Right-Left \n"
00252 " direction (i.e., the coordinates for the CM) \n"
00253 " \n"
00254 " CM AP : Center of mass for the cluster in the \n"
00255 " Anterior-Posterior direction \n"
00256 " \n"
00257 " CM IS : Center of mass for the cluster in the \n"
00258 " Inferior-Superior direction \n"
00259 " \n"
00260 " minRL, maxRL : Bounding box for the cluster, min and max \n"
00261 " coordinates in the Right-Left direction \n"
00262 " \n"
00263 " minAP, maxAP : Min and max coordinates in the Anterior-Posterior \n"
00264 " direction of the volume cluster \n"
00265 " \n"
00266 " minIS, max IS: Min and max coordinates in the Inferior-Superior \n"
00267 " direction of the volume cluster \n"
00268 " \n"
00269 " Mean : Mean value for the volume cluster \n"
00270 " \n"
00271 " SEM : Standard Error of the Mean for the volume cluster \n"
00272 " \n"
00273 " Max Int : Maximum Intensity value for the volume cluster \n"
00274 " \n"
00275 " MI RL : Maximum Intensity value in the Right-Left \n"
00276 " direction of the volume cluster \n"
00277 " \n"
00278 " MI AP : Maximum Intensity value in the Anterior-Posterior \n"
00279 " direction of the volume cluster \n"
00280 " \n"
00281 " MI IS : Maximum Intensity value in the Inferior-Superior \n"
00282 " direction of the volume cluster \n"
00283 "----------------------------------------------------------------------- \n"
00284 " \n"
00285 "Nota Bene: \n"
00286 " \n"
00287 " * The program does not work on complex- or rgb-valued datasets! \n"
00288 " \n"
00289 " * Using the -1noneg option is strongly recommended! \n"
00290 " \n"
00291 " * 3D+time datasets are allowed, but only if you use the \n"
00292 " -1tindex and -1dindex options. \n"
00293 " \n"
00294 " * Bucket datasets are allowed, but you will almost certainly \n"
00295 " want to use the -1tindex and -1dindex options with these. \n"
00296 " \n"
00297 " * SEM values are not realistic for interpolated data sets! \n"
00298 " A ROUGH correction is to multiply the SEM of the interpolated \n"
00299 " data set by the square root of the number of interpolated \n"
00300 " voxels per original voxel. \n"
00301 " \n"
00302 " * If you use -dxyz=1, then rmm should be given in terms of \n"
00303 " voxel edges (not mm) and vmul should be given in terms of \n"
00304 " voxel counts (not microliters). Thus, to connect to only \n"
00305 " 3D nearest neighbors and keep clusters of 10 voxels or more, \n"
00306 " use something like '3dclust -dxyz=1 1.01 10 dset+orig'. \n"
00307 " In the report, 'Volume' will be voxel count, but the rest of \n"
00308 " the coordinate dependent information will be in actual xyz \n"
00309 " millimeters. \n"
00310 " \n"
00311 " * The default coordinate output order is DICOM. If you prefer \n"
00312 " the SPM coordinate order, use the option '-orient LPI' or \n"
00313 " set the environment variable AFNI_ORIENT to 'LPI'. For more \n"
00314 " information, see file README.environment. \n"
00315 ) ;
00316 exit(0) ;
00317 }
00318
00319 mainENTRY("3dclust main"); machdep(); AFNI_logger("3dclust",argc,argv);
00320 PRINT_VERSION("3dclust") ;
00321
00322 THD_coorder_fill( my_getenv("AFNI_ORIENT") , &CL_cord ) ;
00323 CL_read_opts( argc , argv ) ;
00324 nopt = CL_nopt ;
00325
00326 if( CL_do_mni )
00327 THD_coorder_fill( "LPI" , &CL_cord ) ;
00328
00329
00330 if( !CL_quiet ){
00331 printf ("\n\n");
00332 printf ("Program: %s \n", PROGRAM_NAME);
00333 printf ("Author: %s \n", PROGRAM_AUTHOR);
00334 printf ("Date: %s \n", PROGRAM_DATE);
00335 printf ("\n");
00336 }
00337
00338 if( nopt+3 > argc ){
00339 fprintf(stderr,"\n*** No rmm or vmul arguments?\a\n") ;
00340 exit(1) ;
00341 }
00342
00343 rmm = strtod( argv[nopt++] , NULL ) ;
00344 vmul = strtod( argv[nopt++] , NULL ) ;
00345 if( rmm < 0.0 ){
00346 fprintf(stderr,"\n*** Illegal rmm=%f \a\n",rmm) ;
00347 exit(1) ;
00348 }
00349
00350
00351
00352 if( CL_edopt.clust_rmm >= 0.0 ){
00353 fprintf(stderr,"** Warning: -1clust can't be used in 3dclust!\n") ;
00354 CL_edopt.clust_rmm = -1.0 ;
00355 }
00356
00357
00358
00359 dset = NULL ;
00360 for( iarg=nopt ; iarg < argc ; iarg++ ){
00361 if( dset != NULL ) THD_delete_3dim_dataset( dset , False ) ;
00362 dset = THD_open_dataset( argv[iarg] ) ;
00363 if( dset == NULL ){
00364 fprintf(stderr,"** Warning: skipping dataset %s\n",argv[iarg]) ;
00365 continue ;
00366 }
00367 if( DSET_NUM_TIMES(dset) > 1 &&
00368 ( CL_edopt.iv_fim < 0 || CL_edopt.iv_thr < 0 ) ){
00369
00370 fprintf(stderr,
00371 "** Cannot use time-dependent dataset %s\n",argv[iarg]) ;
00372 continue ;
00373 }
00374 THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
00375 if( CL_verbose )
00376 fprintf(stderr,"+++ Loading dataset %s\n",argv[iarg]) ;
00377 THD_load_datablock( dset->dblk ) ;
00378 EDIT_one_dataset( dset , &CL_edopt ) ;
00379
00380
00381
00382 do_mni = (CL_do_mni && dset->view_type == VIEW_TALAIRACH_TYPE) ;
00383
00384 if( CL_ivfim < 0 )
00385 ivfim = DSET_PRINCIPAL_VALUE(dset) ;
00386 else
00387 ivfim = CL_ivfim ;
00388
00389 vfim = DSET_ARRAY(dset,ivfim) ;
00390 fimfac = DSET_BRICK_FACTOR(dset,ivfim) ;
00391 if( vfim == NULL ){
00392 fprintf(stderr,"** Cannot access data in dataset %s\a\n",argv[iarg]) ;
00393 continue ;
00394 }
00395
00396 if( !AFNI_GOOD_FUNC_DTYPE( DSET_BRICK_TYPE(dset,ivfim) ) ||
00397 DSET_BRICK_TYPE(dset,ivfim) == MRI_rgb ){
00398
00399 fprintf(stderr,"** Illegal datum type in dataset %s\a\n",argv[iarg]) ;
00400 continue ;
00401 }
00402
00403 nx = dset->daxes->nxx ; dx = fabs(dset->daxes->xxdel) ;
00404 ny = dset->daxes->nyy ; dy = fabs(dset->daxes->yydel) ;
00405 nz = dset->daxes->nzz ; dz = fabs(dset->daxes->zzdel) ;
00406 nxy = nx * ny ; nxyz = nxy * nz ;
00407
00408 if( CL_edopt.fake_dxyz ){ dxf = dyf = dzf = 1.0 ; }
00409 else { dxf = dx ; dyf = dy ; dzf = dz ; }
00410
00411 if( vmul >= 0.0 )
00412 ptmin = (int) (vmul / (dxf*dyf*dzf) + 0.99) ;
00413 else
00414 ptmin = (int) fabs(vmul) ;
00415
00416
00417 if( !CL_quiet ){
00418
00419 if( CL_summarize != 1 ){
00420 printf( "\n"
00421 "Cluster report for file %s %s\n"
00422 #if 0
00423 "[3D Dataset Name: %s ]\n"
00424 "[ Short Label: %s ]\n"
00425 #endif
00426 "[Connectivity radius = %.2f mm Volume threshold = %.2f ]\n"
00427 "[Single voxel volume = %.1f (microliters) ]\n"
00428 "[Voxel datum type = %s ]\n"
00429 "[Voxel dimensions = %.3f mm X %.3f mm X %.3f mm ]\n",
00430 argv[iarg] ,
00431 do_mni ? "[MNI coords]" : "" ,
00432 #if 0
00433 dset->self_name , dset->label1 ,
00434 #endif
00435 rmm , ptmin*dx*dy*dz , dx*dy*dz ,
00436 MRI_TYPE_name[ DSET_BRICK_TYPE(dset,ivfim) ] ,
00437 dx,dy,dz );
00438
00439 if( CL_edopt.fake_dxyz )
00440 printf("[Fake voxel dimen = %.3f mm X %.3f mm X %.3f mm ]\n",
00441 dxf,dyf,dzf) ;
00442
00443 if (CL_noabs)
00444 printf ("Mean and SEM based on Signed voxel intensities: \n\n");
00445 else
00446 printf ("Mean and SEM based on Absolute Value "
00447 "of voxel intensities: \n\n");
00448
00449 printf (
00450 "Volume CM %s CM %s CM %s min%s max%s min%s max%s min%s max%s Mean SEM Max Int MI %s MI %s MI %s\n"
00451 "------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ------- ------- ------- ----- ----- -----\n",
00452
00453 ORIENT_tinystr[ CL_cord.xxor ] ,
00454 ORIENT_tinystr[ CL_cord.yyor ] ,
00455 ORIENT_tinystr[ CL_cord.zzor ] ,
00456 ORIENT_tinystr[ CL_cord.xxor ] , ORIENT_tinystr[ CL_cord.xxor ] ,
00457 ORIENT_tinystr[ CL_cord.yyor ] , ORIENT_tinystr[ CL_cord.yyor ] ,
00458 ORIENT_tinystr[ CL_cord.zzor ] , ORIENT_tinystr[ CL_cord.zzor ] ,
00459 ORIENT_tinystr[ CL_cord.xxor ] ,
00460 ORIENT_tinystr[ CL_cord.yyor ] ,
00461 ORIENT_tinystr[ CL_cord.zzor ]
00462 ) ;
00463
00464 } else {
00465 if (CL_noabs)
00466 printf ("Mean and SEM based on Signed voxel intensities: \n");
00467 else
00468 printf ("Mean and SEM based on Absolute Value "
00469 "of voxel intensities: \n");
00470 printf("Cluster summary for file %s %s\n" ,
00471 argv[iarg] , do_mni ? "[MNI coords]" : "");
00472 printf("Volume CM %s CM %s CM %s Mean SEM \n",ORIENT_tinystr[ CL_cord.xxor ],ORIENT_tinystr[ CL_cord.yyor ] ,ORIENT_tinystr[ CL_cord.zzor ]);
00473 }
00474 }
00475
00476
00477
00478 clar = NIH_find_clusters( nx,ny,nz , dxf,dyf,dzf ,
00479 DSET_BRICK_TYPE(dset,ivfim) , vfim , rmm , CL_isomode ) ;
00480
00481
00482
00483 PURGE_DSET( dset ) ;
00484
00485 if( clar == NULL || clar->num_clu == 0 ){
00486 printf("** NO CLUSTERS FOUND ***\n") ;
00487 if( clar != NULL ) DESTROY_CLARR(clar) ;
00488 continue ;
00489 }
00490
00491
00492
00493 INIT_CLARR(clbig) ;
00494 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
00495 cl = clar->clar[iclu] ;
00496 if( cl != NULL && cl->num_pt >= ptmin ){
00497 ADDTO_CLARR(clbig,cl) ;
00498 clar->clar[iclu] = NULL ;
00499 }
00500 }
00501 DESTROY_CLARR(clar) ;
00502 clar = clbig ;
00503 if( clar == NULL || clar->num_clu == 0 ){
00504 printf("** NO CLUSTERS FOUND ***\n") ;
00505 if( clar != NULL ) DESTROY_CLARR(clar) ;
00506 continue ;
00507 }
00508
00509
00510
00511
00512
00513 if( iarg == nopt && CL_prefix != NULL ){
00514 int qv ; byte *mmm ;
00515
00516
00517
00518 mmm = (byte *) calloc(sizeof(byte),nxyz) ;
00519 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
00520 cl = clar->clar[iclu] ; if( cl == NULL ) continue ;
00521 for( ipt=0 ; ipt < cl->num_pt ; ipt++ ){
00522 ii = cl->i[ipt] ; jj = cl->j[ipt] ; kk = cl->k[ipt] ;
00523 mmm[ii+jj*nx+kk*nxy] = 1 ;
00524 }
00525 }
00526
00527 DSET_load( dset ) ;
00528
00529 EDIT_dset_items( dset ,
00530 ADN_prefix , CL_prefix ,
00531 ADN_none ) ;
00532
00533 tross_Make_History( "3dclust" , argc , argv , dset ) ;
00534
00535
00536
00537 for( qv=0 ; qv < DSET_NVALS(dset) ; qv++ ){
00538
00539 switch( DSET_BRICK_TYPE(dset,qv) ){
00540
00541 case MRI_short:{
00542 short *bar = (short *) DSET_ARRAY(dset,qv) ;
00543 for( ii=0 ; ii < nxyz ; ii++ )
00544 if( mmm[ii] == 0 ) bar[ii] = 0 ;
00545 }
00546 break ;
00547
00548 case MRI_byte:{
00549 byte *bar = (byte *) DSET_ARRAY(dset,qv) ;
00550 for( ii=0 ; ii < nxyz ; ii++ )
00551 if( mmm[ii] == 0 ) bar[ii] = 0 ;
00552 }
00553 break ;
00554
00555 case MRI_int:{
00556 int *bar = (int *) DSET_ARRAY(dset,qv) ;
00557 for( ii=0 ; ii < nxyz ; ii++ )
00558 if( mmm[ii] == 0 ) bar[ii] = 0 ;
00559 }
00560 break ;
00561
00562 case MRI_float:{
00563 float *bar = (float *) DSET_ARRAY(dset,qv) ;
00564 for( ii=0 ; ii < nxyz ; ii++ )
00565 if( mmm[ii] == 0 ) bar[ii] = 0.0 ;
00566 }
00567 break ;
00568
00569 case MRI_double:{
00570 double *bar = (double *) DSET_ARRAY(dset,qv) ;
00571 for( ii=0 ; ii < nxyz ; ii++ )
00572 if( mmm[ii] == 0 ) bar[ii] = 0.0 ;
00573 }
00574 break ;
00575
00576 case MRI_complex:{
00577 complex *bar = (complex *) DSET_ARRAY(dset,qv) ;
00578 for( ii=0 ; ii < nxyz ; ii++ )
00579 if( mmm[ii] == 0 ) bar[ii].r = bar[ii].i = 0.0 ;
00580 }
00581 break ;
00582
00583 case MRI_rgb:{
00584 byte *bar = (byte *) DSET_ARRAY(dset,qv) ;
00585 for( ii=0 ; ii < nxyz ; ii++ )
00586 if( mmm[ii] == 0 ) bar[3*ii] = bar[3*ii+1] = bar[3*ii+2] = 0 ;
00587 }
00588 break ;
00589 }
00590 }
00591
00592
00593
00594 DSET_write(dset) ;
00595 fprintf(stderr,"++ Wrote dataset %s\n",DSET_BRIKNAME(dset)) ;
00596 PURGE_DSET(dset) ; free(mmm) ;
00597 }
00598
00599
00600
00601 if( clar->num_clu < 1000 ){
00602 SORT_CLARR(clar) ;
00603 } else if( CL_summarize != 1 ){
00604 printf("** TOO MANY CLUSTERS TO SORT BY VOLUME ***\n") ;
00605 }
00606 ndet = 0 ;
00607
00608 vol_total = nvox_total = 0 ;
00609 glmmsum = glmssum = glsqsum = glxxsum = glyysum = glzzsum = 0;
00610
00611 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
00612 cl = clar->clar[iclu] ;
00613 if( cl == NULL || cl->num_pt < ptmin ) continue ;
00614
00615 volsum = cl->num_pt * dxf*dyf*dzf ;
00616 xxsum = yysum = zzsum = mmsum = mssum = 0.0 ;
00617 xxmax = yymax = zzmax = mmmax = msmax = 0.0 ;
00618 sqsum = sem = 0;
00619
00620
00621 RLmax = APmax = ISmax = -1000;
00622 RLmin = APmin = ISmin = 1000;
00623
00624 for( ipt=0 ; ipt < cl->num_pt ; ipt++ ){
00625
00626 #if 0
00627
00628 IJK_TO_THREE( cl->ijk[ipt] , ii,jj,kk , nx,nxy ) ;
00629 #endif
00630 ii = cl->i[ipt] ; jj = cl->j[ipt] ; kk = cl->k[ipt] ;
00631
00632 fv = THD_3dind_to_3dmm( dset , TEMP_IVEC3(ii,jj,kk) ) ;
00633 fv = THD_3dmm_to_dicomm( dset , fv ) ;
00634 xx = fv.xyz[0] ; yy = fv.xyz[1] ; zz = fv.xyz[2] ;
00635 if( !do_mni )
00636 THD_dicom_to_coorder( &CL_cord , &xx,&yy,&zz ) ;
00637 else
00638 THD_3tta_to_3mni( &xx , &yy , &zz ) ;
00639
00640 ms = cl->mag[ipt];
00641 mm = fabs(ms);
00642
00643 mssum += ms;
00644 mmsum += mm;
00645 sqsum += mm * mm;
00646 xxsum += mm * xx ; yysum += mm * yy ; zzsum += mm * zz ;
00647 if( mm > mmmax ){
00648 xxmax = xx ; yymax = yy ; zzmax = zz ;
00649 mmmax = mm ; msmax = ms ;
00650 }
00651
00652
00653 if ( xx > RLmax )
00654 RLmax = xx;
00655 if ( xx < RLmin )
00656 RLmin = xx;
00657 if ( yy > APmax )
00658 APmax = yy;
00659 if ( yy < APmin )
00660 APmin = yy;
00661 if ( zz > ISmax )
00662 ISmax = zz;
00663 if ( zz < ISmin )
00664 ISmin = zz;
00665
00666 }
00667 if( mmsum == 0.0 ) continue ;
00668
00669 glmssum += mssum;
00670 glmmsum += mmsum;
00671 glsqsum += sqsum ;
00672 glxxsum += xxsum;
00673 glyysum += yysum;
00674 glzzsum += zzsum;
00675
00676 ndet++ ;
00677 xxsum /= mmsum ; yysum /= mmsum ; zzsum /= mmsum ;
00678
00679 if (CL_noabs) mean = mssum / cl->num_pt;
00680 else mean = mmsum / cl->num_pt;
00681
00682 if( fimfac != 0.0 )
00683 { mean *= fimfac; msmax *= fimfac;
00684 sqsum *= fimfac*fimfac; }
00685
00686
00687
00688
00689
00690 if (cl->num_pt > 1)
00691 {
00692 sem = (sqsum - (cl->num_pt * mean * mean)) / (cl->num_pt - 1);
00693 if (sem > 0.0) sem = sqrt( sem / cl->num_pt ); else sem = 0.0;
00694 }
00695 else
00696 sem = 0.0;
00697
00698 if( CL_summarize != 1 ){
00699 MCW_fc7(mean, buf1) ;
00700 MCW_fc7(msmax,buf2) ;
00701 MCW_fc7(sem ,buf3) ;
00702
00703 printf("%6.0f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %7s %7s %7s %5.1f %5.1f %5.1f \n",
00704 volsum, xxsum, yysum, zzsum, RLmin, RLmax, APmin, APmax, ISmin, ISmax, buf1, buf3, buf2, xxmax, yymax, zzmax ) ;
00705 }
00706
00707 nvox_total += cl->num_pt ;
00708 vol_total += volsum ;
00709
00710 }
00711
00712 DESTROY_CLARR(clar) ;
00713 if( ndet == 0 )
00714 printf("** NO CLUSTERS FOUND ABOVE THRESHOLD VOLUME ***\n") ;
00715
00716
00717
00718
00719 if (CL_noabs) glmean = glmssum / nvox_total;
00720 else glmean = glmmsum / nvox_total;
00721
00722
00723 if( fimfac != 0.0 ){ glsqsum *= fimfac*fimfac ; glmean *= fimfac ; }
00724 if (nvox_total > 1)
00725 {
00726 sem = (glsqsum - (nvox_total*glmean*glmean)) / (nvox_total - 1);
00727 if (sem > 0.0) sem = sqrt( sem / nvox_total ); else sem = 0.0;
00728 }
00729 else
00730 sem = 0.0;
00731
00732 glxxsum /= glmmsum ; glyysum /= glmmsum ; glzzsum /= glmmsum ;
00733
00734
00735 if( CL_summarize == 1 )
00736 { if( !CL_quiet )
00737 printf( "------ ----- ----- ----- -------- -------- \n");
00738 printf("%6.0f %5.1f %5.1f %5.1f %8.1f %6.3f\n" , vol_total, glxxsum, glyysum, glzzsum, glmean, sem ) ;
00739 }
00740 else if( ndet > 1 && CL_summarize != -1 )
00741 {
00742 MCW_fc7(glmean ,buf1) ;
00743 MCW_fc7(sem ,buf3) ;
00744 if( !CL_quiet )
00745 printf ("------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ------- ------- ------- ----- ----- -----\n");
00746 printf ("%6.0f %5.1f %5.1f %5.1f %7s %7s \n",
00747 vol_total, glxxsum, glyysum, glzzsum, buf1, buf3 ) ;
00748 }
00749
00750 }
00751
00752 exit(0) ;
00753 }
00754
00755
00756
00757
00758
00759
00760
00761 #ifdef CLDEBUG
00762 # define DUMP1 fprintf(stderr,"ARG: %s\n",argv[nopt])
00763 # define DUMP2 fprintf(stderr,"ARG: %s %s\n",argv[nopt],argv[nopt+1])
00764 # define DUMP3 fprintf(stderr,"ARG: %s %s %s\n",argv[nopt],argv[nopt+1],argv[nopt+2])
00765 #else
00766 # define DUMP1
00767 # define DUMP2
00768 # define DUMP3
00769 #endif
00770
00771 void CL_read_opts( int argc , char * argv[] )
00772 {
00773 int nopt = 1 ;
00774 int ival ;
00775
00776 INIT_EDOPT( &CL_edopt ) ;
00777
00778 while( nopt < argc && argv[nopt][0] == '-' ){
00779
00780
00781
00782 ival = EDIT_check_argv( argc , argv , nopt , &CL_edopt ) ;
00783 if( ival > 0 ){
00784 nopt += ival ;
00785 continue ;
00786 }
00787
00788
00789
00790 if( strcmp(argv[nopt],"-isovalue") == 0 ){
00791 CL_isomode = ISOVALUE_MODE ;
00792 nopt++ ; continue ;
00793 }
00794
00795 if( strcmp(argv[nopt],"-isomerge") == 0 ){
00796 CL_isomode = ISOMERGE_MODE ;
00797 nopt++ ; continue ;
00798 }
00799
00800
00801
00802 if( strcmp(argv[nopt],"-mni") == 0 ){
00803 CL_do_mni = 1 ;
00804 nopt++ ; continue ;
00805 }
00806
00807
00808
00809 if( strcmp(argv[nopt],"-prefix") == 0 ){
00810 if( ++nopt >= argc ){
00811 fprintf(stderr,"need an argument after -prefix!\n") ; exit(1) ;
00812 }
00813 CL_prefix = argv[nopt] ;
00814 if( !THD_filename_ok(CL_prefix) ){
00815 fprintf(stderr,"-prefix string is illegal: %s\n",CL_prefix); exit(1);
00816 }
00817 nopt++ ; continue ;
00818 }
00819
00820
00821
00822 if( strncmp(argv[nopt],"-1dindex",5) == 0 ){
00823 if( ++nopt >= argc ){
00824 fprintf(stderr,"need an argument after -1dindex!\n") ; exit(1) ;
00825 }
00826 CL_ivfim = CL_edopt.iv_fim = (int) strtod( argv[nopt++] , NULL ) ;
00827 continue ;
00828 }
00829
00830 if( strncmp(argv[nopt],"-1tindex",5) == 0 ){
00831 if( ++nopt >= argc ){
00832 fprintf(stderr,"need an argument after -1tindex!\n") ; exit(1) ;
00833 }
00834 CL_ivthr = CL_edopt.iv_thr = (int) strtod( argv[nopt++] , NULL ) ;
00835 continue ;
00836 }
00837
00838
00839
00840 if( strncmp(argv[nopt],"-summarize",5) == 0 ){
00841 CL_summarize = 1 ;
00842 nopt++ ; continue ;
00843 }
00844
00845
00846
00847 if( strncmp(argv[nopt],"-nosum",6) == 0 ){
00848 CL_summarize = -1 ;
00849 nopt++ ; continue ;
00850 }
00851
00852
00853
00854 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00855 CL_verbose = CL_edopt.verbose = 1 ;
00856 nopt++ ; continue ;
00857 }
00858
00859
00860
00861 if( strncmp(argv[nopt],"-quiet",5) == 0 ){
00862 CL_quiet = 1 ;
00863 nopt++ ; continue ;
00864 }
00865
00866
00867
00868 if( strncmp(argv[nopt],"-orient",5) == 0 ){
00869 THD_coorder_fill( argv[++nopt] , &CL_cord ) ;
00870 nopt++ ; continue ;
00871 }
00872
00873
00874
00875 if( strncmp(argv[nopt],"-noabs",6) == 0 ){
00876 CL_noabs = 1 ;
00877 nopt++ ; continue ;
00878 }
00879
00880
00881
00882 fprintf(stderr,"** Unrecognized option %s\a\n",argv[nopt]) ;
00883 exit(1) ;
00884
00885 }
00886
00887 #ifdef CLDEBUG
00888 printf("** Finished with options\n") ;
00889 #endif
00890
00891 CL_nopt = nopt ;
00892 return ;
00893 }
00894
00895
00896
00897 void MCW_fc7( float qval , char * buf )
00898 {
00899 float aval = fabs(qval) ;
00900 int lv , il ;
00901 char lbuf[16] ;
00902
00903
00904
00905 lv = (int) qval ;
00906
00907 if( qval == lv && abs(lv) < 100000 ){
00908 if( lv >= 0 ) sprintf( buf , " %d" , lv ) ;
00909 else sprintf( buf , "%d" , lv ) ;
00910 return ;
00911 }
00912
00913
00914
00915 #define BSTRIP \
00916 for( il=6 ; il>1 && lbuf[il]=='0' ; il-- ) lbuf[il] = '\0'
00917
00918
00919
00920 lv = (int) (10.0001 + log10(aval)) ;
00921
00922 switch( lv ){
00923
00924 default:
00925 if( qval > 0.0 ) sprintf( lbuf , "%7.1e" , qval ) ;
00926 else sprintf( lbuf , "%7.0e" , qval ) ;
00927 break ;
00928
00929 case 7:
00930 case 8:
00931 case 9:
00932 case 10:
00933 sprintf( lbuf , "%7.4f" , qval ) ; BSTRIP ; break ;
00934
00935 case 11:
00936 sprintf( lbuf , "%7.3f" , qval ) ; BSTRIP ; break ;
00937
00938 case 12:
00939 sprintf( lbuf , "%7.2f" , qval ) ; BSTRIP ; break ;
00940
00941 case 13:
00942 sprintf( lbuf , "%7.1f" , qval ) ; BSTRIP ; break ;
00943
00944 case 14:
00945 sprintf( lbuf , "%7.0f" , qval ) ; break ;
00946
00947 }
00948 strcpy(buf,lbuf) ;
00949 return ;
00950 }
00951
