Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
imstat.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include <string.h>
00008 #include "mrilib.h"
00009
00010 int main( int argc , char *argv[] )
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 ;
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 ;
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 }