00011 {
00012    MRI_IMAGE *im , *flim ;
00013    MRI_IMARR *imar ;
00014    int kim ;
00015 
00016    int ii , npix , imk , nzero ;
00017 
00018    float im_min,im_min2 , im_max,im_max2 , im_ave , im_std ;
00019    float *flar ;
00020 
00021    int nopt=1 , dolabel=TRUE , doquiet = FALSE ;
00022    char * pix_prefix = NULL ;
00023 
00024    float xcm , ycm , mmm ;  /* 08 Nov 2001 */
00025 
00026    /*-------*/
00027 
00028    if( argc < 2 || strncmp(argv[1],"-h",2) == 0 ){
00029       printf(
00030         "Calculation of statistics of one or more images.\n"
00031         "Usage: imstat [-nolabel] [-pixstat prefix] [-quiet] image_file ...\n"
00032         "  -nolabel        = don't write labels on each file's summary line\n"
00033         "  -quiet          = don't print statistics for each file\n"
00034         "  -pixstat prefix = if more than one image file is given, then\n"
00035         "                     'prefix.mean' and 'prefix.sdev' will be written\n"
00036         "                     as the pixel-wise statistics images of the whole\n"
00037         "                     collection.  These images will be in the 'flim'\n"
00038         "                     floating point format.  [This option only works\n"
00039         "                     on 2D images!]\n"
00040       ) ;
00041       exit(0) ;
00042    }
00043 
00044    machdep() ;
00045 
00046    while( nopt < argc && argv[nopt][0] == '-' ){
00047       if( strncmp(argv[nopt],"-nolabel",4) == 0 ){
00048          dolabel = FALSE ;
00049          nopt++ ; continue ;
00050       }
00051       if( strncmp(argv[nopt],"-quiet",4) == 0 ){
00052          doquiet = TRUE ;
00053          nopt++ ; continue ;
00054       }
00055       if( strncmp(argv[nopt],"-pixstat",4) == 0 ){
00056          pix_prefix = argv[++nopt] ;
00057          nopt++ ; continue ;
00058       }
00059       fprintf(stderr,"Unknown option %s\n",argv[nopt]) ; exit(1) ;
00060    }
00061    if( nopt >= argc ){
00062       fprintf(stderr,"No filename arguments?!\n") ; exit(1) ;
00063    }
00064 
00065    /*--------*/
00066 
00067    imar = mri_read_many_files( argc-nopt , argv+nopt ) ;
00068    if( imar == NULL || IMARR_COUNT(imar) < 1 ){
00069       fprintf(stderr,"Can't read any images!?\n") ; exit(1) ;
00070    }
00071 
00072    for( kim=0 ; kim < IMARR_COUNT(imar) ; kim++ ){
00073 
00074       im  = IMARR_SUBIMAGE(imar,kim) ;
00075       imk = im->kind ;
00076 
00077       if( im->kind == MRI_float ){
00078          flim = im ;
00079       } else {
00080          flim = mri_to_float( im ) ; mri_free( im ) ;
00081       }
00082 
00083       npix = flim->nvox ;
00084       flar = MRI_FLOAT_PTR( flim ) ;
00085 
00086       xcm = ycm = mmm = 0.0 ;    /* 08 Nov 2001 */
00087 
00088       im_max = im_min = flar[0] ;
00089       im_ave = 0.0 ;
00090       for( ii=0 ; ii < npix ; ii++ ){
00091               if( flar[ii] > im_max ) im_max  = flar[ii] ;
00092          else if( flar[ii] < im_min ) im_min  = flar[ii] ;
00093          im_ave += flar[ii] ;
00094 
00095          xcm += fabs(flar[ii]) * (ii/flim->nx) ;
00096          ycm += fabs(flar[ii]) * (ii%flim->nx) ;
00097          mmm += fabs(flar[ii]) ;
00098       }
00099       if( mmm > 0.0 ){ xcm /= mmm; ycm /= mmm; }
00100 
00101       im_ave /= npix ;
00102 
00103       im_max2 = 2 * im_min - im_max ;
00104       im_min2 = 2 * im_max - im_min ;
00105       for( ii=0 ; ii < npix ; ii++ ){
00106          if( flar[ii] > im_max2 && flar[ii] < im_max ) im_max2 = flar[ii] ;
00107          if( flar[ii] < im_min2 && flar[ii] > im_min ) im_min2 = flar[ii] ;
00108       }
00109 
00110       im_std = 0.0 ;
00111       nzero  = 0 ;
00112       for( ii=0 ; ii < npix ; ii++ ){
00113          im_std += (flar[ii]-im_ave) * (flar[ii]-im_ave) ;
00114          if( flar[ii] == 0.0 ) nzero++ ;
00115       }
00116       im_std = sqrt(im_std/npix) ;
00117 
00118       if( ! doquiet ){
00119          if( dolabel){
00120             printf( "\nfile = %s  nx = %d  ny = %d  data type = %s\n" ,
00121                     flim->name , flim->nx , flim->ny , MRI_TYPE_name[imk] ) ;
00122             printf( "min =%11.4g  next min=%11.4g  max=%11.4g  next max=%11.4g\n" ,
00123                     im_min,im_min2 , im_max,im_max2 ) ;
00124             printf( "mean=%11.4g  std.dev.=%11.4g  number of zero pixels = %d\n" ,
00125                     im_ave,im_std,nzero ) ;
00126             if( mmm > 0.0 ) printf("x-CM=%11.4g     y-CM=%11.4g\n",xcm,ycm) ;
00127          } else {
00128             printf( "%d %d %d " , flim->nx , flim->ny , imk ) ;
00129             printf( "%11.4g %11.4g %11.4g %11.4g %11.4g %11.4g %d\n" ,
00130                     im_min,im_min2 , im_max,im_max2 , im_ave,im_std,nzero ) ;
00131          }
00132       }
00133 
00134       if( pix_prefix != NULL ) (void)mri_stat_seq( flim ) ;
00135 
00136       mri_free(flim) ;
00137    }
00138 
00139    if( pix_prefix != NULL ){
00140       char fname[128] ;
00141       MRI_IMAGE ** sim ;
00142       float * mar , * sar ;
00143       float cmax ;
00144 
00145       sim = mri_stat_seq( NULL ) ;
00146       ii = strlen(pix_prefix) ;
00147 
00148       if( pix_prefix[ii-1] == '.' ) sprintf( fname , "%smean" , pix_prefix ) ;
00149       else                          sprintf( fname , "%s.mean" , pix_prefix ) ;
00150       mri_write( fname , sim[0] ) ;
00151       printf("-- Wrote mean image to %s\n",fname) ;
00152 
00153       if( pix_prefix[ii-1] == '.' ) sprintf( fname , "%ssdev" , pix_prefix ) ;
00154       else                          sprintf( fname , "%s.sdev" , pix_prefix ) ;
00155       mri_write( fname , sim[1] ) ;
00156       printf("-- Wrote standard deviation image to %s\n",fname) ;
00157 
00158       if( pix_prefix[ii-1] == '.' ) sprintf( fname , "%scvar" , pix_prefix ) ;
00159       else                          sprintf( fname , "%s.cvar" , pix_prefix ) ;
00160       npix = sim[0]->nx * sim[0]->ny ;
00161       mar  = MRI_FLOAT_PTR(sim[0]) ;
00162       sar  = MRI_FLOAT_PTR(sim[1]) ;
00163       for( ii=0 ; ii < npix ; ii++ )
00164          if( mar[ii] != 0.0 ) mar[ii] = sar[ii] / fabs(mar[ii]) ;
00165       cmax = mri_max( sim[0] ) ;
00166       for( ii=0 ; ii < npix ; ii++ )
00167          if( mar[ii] == 0.0 ) mar[ii] = cmax ;
00168       mri_write( fname , sim[0] ) ;
00169       printf("-- Wrote coefficient of variation image to %s\n",fname) ;
00170    }
00171 
00172    exit(0) ;
00173 }