00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include <stdlib.h>
00008 #include <stdio.h>
00009 #include <string.h>
00010 #include <math.h>
00011 #include <ctype.h>
00012 
00013 #include "mrilib.h"
00014 
00015 #define MAX_NAME 64
00016 
00017 typedef struct {
00018       int         count ;
00019       char        name[MAX_NAME] ;
00020       MRI_IMAGE * avim , * sdim ;
00021 } SF_interval ;
00022 
00023 /*** global data!!! ***/
00024 
00025 static int           SF_numint = 0 ;
00026 static SF_interval * SF_int = NULL ;
00027 static char          SF_bname[MAX_NAME] = "rest"  ;
00028 static char          SF_pname[MAX_NAME] = "sfim." ;
00029 static char          SF_iname[MAX_NAME] = "sfint" ;
00030 static MRI_IMARR *   SF_imts = NULL ;
00031 static int           SF_localbase = 0 ;
00032 
00033 /*** prototypes ***/
00034 
00035 void SFIM_getopts( int , char * argv[] ) ;
00036 SF_interval * SFIM_load_intervals( char * ) ;
00037 void SFIM_syntax( char * ) ;
00038 void SFIM_write_avs() ;
00039 
00040 /**********************************************************************/
00041 
00042 int main( int argc , char * argv[] )
00043 {
00044    int lin , kim , kbot,ktop , nx,ny , npix , ii ,
00045        lbase , lup,ldown ;
00046    MRI_IMAGE ** stat_ret ;
00047    MRI_IMAGE * imb ;
00048    float     * bar , * bav ;
00049 
00050    printf(
00051     "MCW SFIM: Stepwise Functional IMages, by RW Cox\n") ;
00052 
00053    if( argc < 2 ) SFIM_syntax("type sfim -help for usage details") ;
00054    else if( strcmp(argv[1],"-help") == 0 ) SFIM_syntax(NULL) ;
00055 
00056    machdep() ;
00057 
00058    SFIM_getopts( argc , argv ) ;
00059 
00060    /*----- average over each interval -----*/
00061 
00062    nx = SF_imts->imarr[0]->nx ;
00063    ny = SF_imts->imarr[0]->ny ; npix = nx * ny ;
00064 
00065    lin = 0 ; kbot = 0 ;
00066    do {
00067       ktop = kbot + SF_int[lin].count ;
00068       if( ktop > SF_imts->num ) ktop = SF_imts->num ;
00069 
00070       if( isalpha(SF_int[lin].name[0]) ){     /* average if a good name */
00071          for( kim=kbot ; kim < ktop ; kim++ )
00072            (void) mri_stat_seq( SF_imts->imarr[kim] ) ;
00073 
00074          stat_ret         = mri_stat_seq( NULL ) ;
00075          SF_int[lin].avim = stat_ret[0] ;
00076          SF_int[lin].sdim = stat_ret[1] ;
00077       }
00078 
00079       kbot = ktop ; lin ++ ;
00080    } while( SF_int[lin].count > 0 && kbot < SF_imts->num ) ;
00081    SF_numint = lin ;
00082 
00083    /*----- find the number of base intervals -----*/
00084 
00085    lbase = 0 ;
00086    for( lin=0 ; lin < SF_numint ; lin++ )
00087       if( strcmp(SF_int[lin].name,SF_bname) == 0 ) lbase++ ;
00088 
00089    /* no bases --> write averages out now and quit */
00090 
00091    if( lbase <= 0 ){
00092       printf("** no 'base' intervals --> task means not adjusted\n") ;
00093       SFIM_write_avs() ;
00094       exit(0) ;
00095    }
00096 
00097    /* bases yes, but not localbase --> compute global average of bases */
00098 
00099    if( ! SF_localbase ){
00100       int knum = 0 ;
00101 
00102       kbot = 0 ;
00103       for( lin=0 ; lin < SF_numint ; lin++ ){
00104          ktop = kbot + SF_int[lin].count ;
00105          if( ktop > SF_imts->num ) ktop = SF_imts->num ;
00106          if( strcmp(SF_int[lin].name,SF_bname) == 0 ){  /* if a base */
00107             for( kim=kbot ; kim < ktop ; kim++ ){
00108                (void) mri_stat_seq( SF_imts->imarr[kim] ) ; /* average in */
00109                knum ++ ;
00110             }
00111          }
00112          kbot = ktop ;
00113       }
00114       stat_ret = mri_stat_seq( NULL ) ;
00115       imb      = stat_ret[0] ;           /* average of all bases */
00116       mri_free( stat_ret[1] ) ;          /* don't keep st. dev.  */
00117 
00118       printf("** global base = average of %d images\n",knum) ;
00119    }
00120 
00121    /*----- for each non-base interval,
00122            subtract the relevant base average -----*/
00123 
00124    for( lin=0 ; lin < SF_numint ; lin++ ){
00125       int free_imb = 0 ;
00126 
00127       if( !isalpha(SF_int[lin].name[0]) ||
00128           strcmp(SF_int[lin].name,SF_bname) == 0 ) continue ;  /* skip this */
00129 
00130       if( SF_localbase ){
00131          for( lup=lin+1 ; lup < SF_numint ; lup++ )  /* look for a base above */
00132             if( strcmp(SF_int[lup].name,SF_bname) == 0 ) break ;
00133 
00134          for( ldown=lin-1 ; ldown >=0 ; ldown-- )    /* look for a base below */
00135             if( strcmp(SF_int[ldown].name,SF_bname) == 0 ) break ;
00136 
00137          if( ldown < 0 && lup >= SF_numint ){  /* no base?  an error! */
00138             fprintf(stderr,"*** can't find base above or below at lin=%d\n",lin) ;
00139             SFIM_syntax("INTERNAL ERROR -- should not occur!") ;
00140          }
00141 
00142          /* if only have one neighbor, use it, otherwise make average */
00143 
00144          if( ldown <  0         ){
00145             imb = SF_int[lup].avim ; free_imb = 0 ;
00146             printf("** local base for %s = average of %d images above\n",
00147                    SF_int[lin].name , SF_int[lup].count ) ;
00148          }
00149          else if( lup   >= SF_numint ){
00150             imb = SF_int[ldown].avim ; free_imb = 0 ;
00151             printf("** local base for %s = average of %d images below\n",
00152                    SF_int[lin].name , SF_int[ldown].count ) ;
00153          }
00154          else {
00155             float * bup , * bdown ;
00156             bup   = mri_data_pointer( SF_int[lup].avim ) ;
00157             bdown = mri_data_pointer( SF_int[ldown].avim ) ;
00158             imb   = mri_new( nx , ny , MRI_float ) ; free_imb = 1 ;
00159             bar   = mri_data_pointer( imb ) ;
00160             for( ii=0 ; ii < npix ; ii++ )
00161                bar[ii] = 0.5 * ( bup[ii] + bdown[ii] ) ;
00162 
00163             printf("** local base for %s = average of %d below, %d above\n",
00164                    SF_int[lin].name, SF_int[ldown].count, SF_int[lup].count ) ;
00165          }
00166       }
00167 
00168       /* subtract imb (base average) from current interval average */
00169 
00170       bar = mri_data_pointer( imb ) ;
00171       bav = mri_data_pointer( SF_int[lin].avim ) ;
00172       for( ii=0 ; ii < npix ; ii++ ) bav[ii] -= bar[ii] ;
00173 
00174       if( SF_localbase && free_imb ) mri_free( imb ) ;
00175    }
00176 
00177    /*----- now write the averages out -----*/
00178 
00179    SFIM_write_avs() ;
00180    exit(0) ;
00181 }
00182 
00183 /**********************************************************************/
00184 
00185 void SFIM_write_avs()
00186 {
00187    char name[256] ;
00188    int lin , ldown , tnum ;
00189 
00190    for( lin=0 ; lin < SF_numint ; lin++ ){
00191       if( !isalpha(SF_int[lin].name[0]) ||
00192           SF_int[lin].avim == NULL         ) continue ;  /* skip */
00193 
00194       tnum = 1 ;
00195       for( ldown=lin-1 ; ldown >=0 ; ldown-- )
00196          if( strcmp(SF_int[lin].name,SF_int[ldown].name) == 0 ) tnum++ ;
00197 
00198       sprintf( name , "%s%s.%04d" , SF_pname , SF_int[lin].name , tnum ) ;
00199 
00200       printf("-- writing file %s\n",name) ;
00201       mri_write( name , SF_int[lin].avim ) ;
00202    }
00203 }
00204 
00205 /**********************************************************************/
00206 
00207 void SFIM_getopts( int argc , char * argv[] )
00208 {
00209    int nopt = 1 ;
00210    int kk , nx , ny , kount ;
00211 
00212    /*--- do options ---*/
00213 
00214    while( nopt < argc && argv[nopt][0] == '-' ){
00215 
00216       /** -sfint iname **/
00217 
00218       if( strncmp(argv[nopt],"-sfint",5) == 0 ){
00219          if( ++nopt >= argc ) SFIM_syntax("-sfint needs a name after it!") ;
00220          strcpy(SF_iname,argv[nopt]) ;
00221          nopt++ ; continue ;
00222       }
00223 
00224       /** -base bname **/
00225 
00226       if( strncmp(argv[nopt],"-base",5) == 0 ){
00227          if( ++nopt >= argc ) SFIM_syntax("-base needs a name after it!") ;
00228          strcpy(SF_bname,argv[nopt]) ;
00229          nopt++ ; continue ;
00230       }
00231 
00232       /** -prefix pname **/
00233 
00234       if( strncmp(argv[nopt],"-prefix",5) == 0 ){
00235          if( ++nopt >= argc ) SFIM_syntax("-prefix needs a name after it!") ;
00236          strcpy(SF_pname,argv[nopt]) ;
00237          kk = strlen(SF_pname) ;
00238          if( SF_pname[kk-1] != '.' ){
00239             SF_pname[kk]   = '.' ;
00240             SF_pname[kk+1] = '\0' ;
00241          }
00242          nopt++ ; continue ;
00243       }
00244 
00245       /** -localbase **/
00246 
00247       if( strncmp(argv[nopt],"-localbase",5) == 0 ){
00248          SF_localbase = 1 ;
00249          nopt++ ; continue ;
00250       }
00251 
00252       /** illegal option **/
00253 
00254       fprintf(stderr,"*** illegal command line option: %s\n",argv[nopt]) ;
00255       SFIM_syntax("type sfim -help for more details") ;
00256 
00257    }
00258 
00259    /*--- do images ---*/
00260 
00261    SF_imts = mri_read_many_files( argc-nopt , argv+nopt ) ;
00262    if( SF_imts == NULL ) SFIM_syntax("cannot continue without input images!") ;
00263 
00264    /* check images for consistency */
00265 
00266    nx = SF_imts->imarr[0]->nx ;
00267    ny = SF_imts->imarr[0]->ny ;
00268 
00269    for( kk=1 ; kk < SF_imts->num ; kk++ ){
00270       if( nx != SF_imts->imarr[kk]->nx || ny != SF_imts->imarr[kk]->ny ){
00271          fprintf(stderr,"*** image %d not conformant to image 0\n",kk) ;
00272          SFIM_syntax("cannot continue with images whose sizes differ!") ;
00273       }
00274    }
00275 
00276    /*--- read intervals ---*/
00277 
00278    SF_int = SFIM_load_intervals( SF_iname ) ;
00279    if( SF_int[0].count <= 0 ) SFIM_syntax("cannot read interval definitions") ;
00280 
00281    kount = kk = 0 ;
00282    do {
00283       kount += SF_int[kk++].count ;
00284    } while( SF_int[kk].count > 0 ) ;
00285 
00286    if( kount < SF_imts->num ){
00287       fprintf(stderr,"*** # images = %d but only %d images in intervals!\n",
00288               SF_imts->num , kount ) ;
00289       SFIM_syntax("must have at least as many in intervals as actual images!") ;
00290    }
00291 
00292    return ;
00293 }
00294 
00295 /**********************************************************************/
00296 
00297 void SFIM_syntax( char * str )
00298 {
00299    if( str != NULL ){
00300       fprintf(stderr,"*** %s\n",str) ;
00301       exit(-1) ;
00302    }
00303 
00304    printf(
00305     "\n"
00306     "Usage: sfim [options] image_files ...\n"
00307     "\n"
00308     "  + image_files are in the same format AFNI accepts\n"
00309     "  + options are from the following:\n"
00310     "\n"
00311     "  -sfint iname:   'iname' is the name of a file which has\n"
00312     "                  the interval definitions; an example is\n"
00313     "                    3*# 5*rest 4*A 5*rest 4*B 5*rest 4*A 5*rest\n"
00314     "                  which says:\n"
00315     "                    - ignore the 1st 3 images\n"
00316     "                    - take the next 5 as being in task state 'rest'\n"
00317     "                    - take the next 4 as being in task state 'A'\n"
00318     "                    and so on;\n"
00319     "                  task names that start with a nonalphabetic character\n"
00320     "                  are like the '#' above and mean 'ignore'.\n"
00321     "              *** the default 'iname' is 'sfint'\n"
00322     "\n"
00323     "  -base bname:    'bname' is the task state name to use as the\n"
00324     "                  baseline; other task states will have the mean\n"
00325     "                  baseline state subtracted; if there are no task\n"
00326     "                  states from 'iname' that match 'bname', this\n"
00327     "                  subtraction will not occur.\n"
00328     "              *** the default 'bname' is 'rest'\n"
00329     "\n"
00330     "  -localbase:     if this option is present, then each non-base\n"
00331     "                  task state interval has the mean of the two\n"
00332     "                  nearest base intervals subtracted instead of the\n"
00333     "                  grand mean of all the base task intervals.\n"
00334     "\n"
00335     "  -prefix pname:  'pname' is the prefix for output image filenames for\n"
00336     "                  all states:  the i'th interval with task state name\n"
00337     "                  'fred' will be writen to file 'pname.fred.i'.\n"
00338     "              *** the default 'pname' is 'sfim'\n"
00339     "\n"
00340     "  Output files are the base-mean-removed averages for each non-base\n"
00341     "  task interval, and simply the mean for each base task interval.\n"
00342     "  Output images are in the 'flim' (floating pt. image) format, and\n"
00343     "  may be converted to 16 bit shorts using the program 'ftosh'.\n"
00344    ) ;
00345 
00346    exit(0) ;
00347 }
00348 
00349 /**********************************************************************/
00350 
00351 SF_interval * SFIM_load_intervals( char * fname )
00352 {
00353    SF_interval * sfint ;
00354    FILE * fd ;
00355    int num_int , num_all , ic , cc ;
00356    char nn[MAX_NAME] ;
00357 #define INC_SFINT 8
00358 
00359    if( fname == NULL || strlen(fname) == 0 ){
00360       fprintf(stderr,"SFIM_load_intervals illegal filename\n") ;
00361       exit(-1) ;
00362    }
00363    fd = fopen( fname , "r" ) ;
00364    if( fd == NULL ){
00365       fprintf(stderr,"SFIM_load_intervals cannot open filename %s\n",fname) ;
00366       exit(-1) ;
00367    }
00368 
00369    sfint = (SF_interval *) malloc( sizeof(SF_interval) * INC_SFINT ) ;
00370    if( sfint == NULL ){
00371      fprintf(stderr,"SFIM_load_intervals fails to malloc!\n") ;
00372      exit(-1) ;
00373    }
00374    num_int = 0 ;
00375    num_all = INC_SFINT ;
00376 
00377    do {
00378       ic = fscanf( fd , "%d*%s" , &cc , nn ) ;
00379       if( num_int == num_all ){
00380          num_all += INC_SFINT ;
00381          sfint    = (SF_interval *) realloc( (void *) sfint ,
00382                                              sizeof(SF_interval) * num_all ) ;
00383          if( sfint == NULL ){
00384             fprintf(stderr,"SFIM_load_intervals fails to realloc!\n") ;
00385             exit(-1) ;
00386          }
00387       }
00388       if( ic != 2 ){ /* bad scanning */
00389          sfint[num_int].count = -1 ;
00390          break ;
00391       }
00392       sfint[num_int].count = cc ; /* put results into interval */
00393       sfint[num_int].avim  = NULL ;
00394       sfint[num_int].sdim  = NULL ;
00395       strcpy( sfint[num_int].name , nn ) ;
00396 
00397       num_int ++ ;
00398    } while( cc > 0 ) ;
00399 
00400    fclose( fd ) ;
00401    return sfint ;
00402 }