00772 {
00773    int nopt = 1 ;
00774    int ival ;
00775 
00776    INIT_EDOPT( &CL_edopt ) ;
00777 
00778    while( nopt < argc && argv[nopt][0] == '-' ){
00779 
00780       /**** check for editing options ****/
00781 
00782       ival = EDIT_check_argv( argc , argv , nopt , &CL_edopt ) ;
00783       if( ival > 0 ){
00784          nopt += ival ;
00785          continue ;
00786       }
00787 
00788       /**** 30 Apr 2002: -isovalue and -isomerge ****/
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       /**** 30 Apr 2002: -mni ****/
00801 
00802       if( strcmp(argv[nopt],"-mni") == 0 ){
00803          CL_do_mni = 1 ;
00804          nopt++ ; continue ;
00805       }
00806 
00807       /**** 29 Nov 2001: -prefix ****/
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       /**** Sep 16 1999: -1tindex and -1dindex ****/
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       /**** -summarize ****/
00839 
00840       if( strncmp(argv[nopt],"-summarize",5) == 0 ){
00841          CL_summarize = 1 ;
00842          nopt++ ; continue ;
00843       }
00844 
00845       /**** -nosum [05 Jul 2001] ****/
00846 
00847       if( strncmp(argv[nopt],"-nosum",6) == 0 ){
00848          CL_summarize = -1 ;
00849          nopt++ ; continue ;
00850       }
00851 
00852       /**** -verbose ****/
00853 
00854       if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00855          CL_verbose = CL_edopt.verbose = 1 ;
00856          nopt++ ; continue ;
00857       }
00858 
00859       /**** -quiet ****/
00860 
00861       if( strncmp(argv[nopt],"-quiet",5) == 0 ){
00862          CL_quiet = 1 ;
00863          nopt++ ; continue ;
00864       }
00865 
00866       /**** -orient code ****/
00867 
00868       if( strncmp(argv[nopt],"-orient",5) == 0 ){
00869          THD_coorder_fill( argv[++nopt] , &CL_cord ) ; /* July 1997 */
00870          nopt++ ; continue ;
00871       }
00872 
00873       /**** -noabs ****/                                /* BDW  19 Jan 1999 */
00874 
00875       if( strncmp(argv[nopt],"-noabs",6) == 0 ){
00876          CL_noabs = 1 ;
00877          nopt++ ; continue ;
00878       }
00879 
00880       /**** unknown switch ****/
00881 
00882       fprintf(stderr,"** Unrecognized option %s\a\n",argv[nopt]) ;
00883       exit(1) ;
00884 
00885    }  /* end of loop over options */
00886 
00887 #ifdef CLDEBUG
00888 printf("** Finished with options\n") ;
00889 #endif
00890 
00891    CL_nopt = nopt ;
00892    return ;
00893 }

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 ;                  /* 24 Jan 2001: for -dxyz=1 option */
00119    int do_mni ;                         /* 30 Apr 2002 */
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 ) ; /* July 1997 */
00323    CL_read_opts( argc , argv ) ;
00324    nopt = CL_nopt ;
00325 
00326    if( CL_do_mni )
00327      THD_coorder_fill( "LPI" , &CL_cord ) ;  /* 30 Apr 2002 */
00328 
00329  /*----- Identify software -----*/
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    /* BDW  26 March 1999  */
00351 
00352    if( CL_edopt.clust_rmm >= 0.0 ){  /* 01 Nov 1999 */
00353       fprintf(stderr,"** Warning: -1clust can't be used in 3dclust!\n") ;
00354       CL_edopt.clust_rmm  = -1.0 ;
00355    }
00356 
00357    /**-- loop over datasets --**/
00358 
00359    dset = NULL ;
00360    for( iarg=nopt ; iarg < argc ; iarg++ ){
00361       if( dset != NULL ) THD_delete_3dim_dataset( dset , False ) ; /* flush old   */
00362       dset = THD_open_dataset( argv[iarg] ) ;                      /* open new    */
00363       if( dset == NULL ){                                          /* failed?     */
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 ) ){      /* no time     */
00369 
00370          fprintf(stderr,                                           /* dependence! */
00371                  "** Cannot use time-dependent dataset %s\n",argv[iarg]) ;
00372          continue ;
00373       }
00374       THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ; /* don't mmap  */
00375       if( CL_verbose )
00376          fprintf(stderr,"+++ Loading dataset %s\n",argv[iarg]) ;
00377       THD_load_datablock( dset->dblk  ) ;                          /* read in     */
00378       EDIT_one_dataset( dset , &CL_edopt ) ;                       /* editing?    */
00379 
00380       /* 30 Apr 2002: check if -mni should be used here */
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) ;                     /* useful data */
00386       else
00387          ivfim  = CL_ivfim ;                                       /* 16 Sep 1999 */
00388 
00389       vfim   = DSET_ARRAY(dset,ivfim) ;                            /* ptr to data */
00390       fimfac = DSET_BRICK_FACTOR(dset,ivfim) ;                     /* scl factor  */
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 ; }         /* 24 Jan 2001 */
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) ;  /* 30 Apr 2002 */
00415 
00416       /*-- print report header --*/
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]" : "" ,  /* 30 Apr 2002 */
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 )  /* 24 Jan 2001 */
00440                printf("[Fake voxel dimen    = %.3f mm X %.3f mm X %.3f mm ]\n",
00441                       dxf,dyf,dzf) ;
00442 
00443             if (CL_noabs)                                   /* BDW  19 Jan 1999 */
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)                                   /* BDW  19 Jan 1999 */
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       } /* end of report header */
00475 
00476       /*-- actually find the clusters in the dataset */
00477 
00478       clar = NIH_find_clusters( nx,ny,nz , dxf,dyf,dzf ,
00479                                 DSET_BRICK_TYPE(dset,ivfim) , vfim , rmm , CL_isomode ) ;
00480 
00481       /*-- don't need dataset data any more --*/
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 ;                               /* next dataset */
00489       }
00490 
00491       /** edit for volume (June 1995) **/
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 ){ /* big enough */
00497             ADDTO_CLARR(clbig,cl) ;               /* copy pointer */
00498             clar->clar[iclu] = NULL ;             /* null out original */
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       /** end of June 1995 addition **/
00510 
00511       /*-- 29 Nov 2001: write out an edited dataset? --*/
00512 
00513       if( iarg == nopt && CL_prefix != NULL ){
00514         int qv ; byte *mmm ;
00515 
00516         /* make a mask of voxels to keep */
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 ) ;             /* reload data from disk */
00528 
00529         EDIT_dset_items( dset ,         /* rename dataset internally */
00530                            ADN_prefix , CL_prefix ,
00531                          ADN_none ) ;
00532 
00533         tross_Make_History( "3dclust" , argc , argv , dset ) ;
00534 
00535         /* mask out each sub-brick */
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           } /* end of switch over sub-brick type */
00590         } /* end of loop over sub-bricks */
00591 
00592         /* write dataset out */
00593 
00594         DSET_write(dset) ;
00595         fprintf(stderr,"++ Wrote dataset %s\n",DSET_BRIKNAME(dset)) ;
00596         PURGE_DSET(dset) ; free(mmm) ;
00597       }
00598 
00599       /** sort clusters by size, to make a nice report **/
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 ;  /* no good */
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          /* These should be pegged at whatever actual max/min values are */
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 /** this is obsolete and nonfunctional code **/
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 ) ;  /* July 1997 */
00637             else
00638               THD_3tta_to_3mni( &xx , &yy , &zz ) ;           /* 30 Apr 2002 */
00639 
00640             ms = cl->mag[ipt];                           /* BDW  18 Jan 1999 */
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             /* Dimensions: */
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;     /* BDW  19 Jan 1999 */
00680          else            mean = mmsum / cl->num_pt;
00681 
00682          if( fimfac != 0.0 )
00683            { mean  *= fimfac;  msmax *= fimfac;
00684              sqsum *= fimfac*fimfac; }                      /* BDW 12/27/96 */
00685 
00686          /* MSB 11/1/96  Calculate SEM using SEM^2=s^2/N,
00687             where s^2 = (SUM Y^2)/N - (Ymean)^2
00688             where sqsum = (SUM Y^2 ) */
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       /* MSB 11/1/96  Calculate global SEM */
00718 
00719       if (CL_noabs)   glmean = glmssum / nvox_total;    /* BDW  19 Jan 1999 */
00720       else            glmean = glmmsum / nvox_total;
00721 
00722       /* BDW 12/27/96 */
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       /* MSB 11/1/96 Modified so that mean and SEM would print in correct column */
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 }

00898 {
00899    float aval = fabs(qval) ;
00900    int lv , il ;
00901    char lbuf[16] ;
00902 
00903    /* special case if the value is an integer */
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 /* macro to strip trailing zeros from output */
00914 
00915 #define BSTRIP \
00916    for( il=6 ; il>1 && lbuf[il]=='0' ; il-- ) lbuf[il] = '\0'
00917 
00918    /* noninteger: choose floating format based on magnitude */
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:  /* 0.001 -0.01  */
00930       case  8:  /* 0.01  -0.1   */
00931       case  9:  /* 0.1   -1     */
00932       case 10:  /* 1     -9.99  */
00933          sprintf( lbuf , "%7.4f" , qval ) ; BSTRIP ; break ;
00934 
00935       case 11:  /* 10-99.9 */
00936          sprintf( lbuf , "%7.3f" , qval ) ; BSTRIP ; break ;
00937 
00938       case 12:  /* 100-999.9 */
00939          sprintf( lbuf , "%7.2f" , qval ) ; BSTRIP ; break ;
00940 
00941       case 13:  /* 1000-9999.9 */
00942          sprintf( lbuf , "%7.1f" , qval ) ; BSTRIP ; break ;
00943 
00944       case 14:  /* 10000-99999.9 */
00945          sprintf( lbuf , "%7.0f" , qval ) ; break ;
00946 
00947    }
00948    strcpy(buf,lbuf) ;
00949    return ;
00950 }