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 <math.h>
00008 #include <stdio.h>
00009 #include <stdlib.h>
00010 #include <string.h>
00011 
00012 #include "mrilib.h"
00013 #include "pcor.h"
00014 #include "ts.h"
00015 
00016 #ifndef MAX
00017 #  define MAX(x,y)  (((x)>(y)) ? (x) : (y))
00018 #endif
00019 
00020 #ifndef MIN
00021 #  define MIN(x,y)  (((x)<(y)) ? (x) : (y))
00022 #endif
00023 
00024 #ifndef MALLOC_ERR
00025 #  define MALLOC_ERR(str) {fprintf(stderr,"malloc error for %s\a\n",str);exit(1);}
00026 #endif
00027 
00028 #define SCALE 10000
00029 #define BLAST 33333.0
00030 
00031 #define DFILT_NONE   0           /* no derivative filtering */
00032 
00033 #define DFILT_SPACE  1           /* spatial mode */
00034 #define DFILT_SPACE0 11
00035 
00036 #define DFILT_TIME   2           /* temporal mode */
00037 #define DFILT_TIME0  21
00038 
00039 #define DFILT_BOTH   3           /* both, separately */
00040 #define DFILT_BOTH0  31
00041 
00042 #define DFILT_SPACETIME  4       /* both, together */
00043 #define DFILT_SPACETIME0 41
00044 
00045 #define CLIP_FRAC 0.10           /* for clipping nonbrain (low intensity) pixels */
00046 
00047 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00048 
00049 /** structure type to hold results from scanning the command line options **/
00050 
00051 typedef struct {
00052     char *fname_corr ,          /* filename for correlation image */
00053          *fname_fim ,           /* filename for activation image */
00054          *fname_cnr ,           /* filename for contrast-to-noise image */
00055          *fname_sig ,           /* filename for standard deviation image */
00056          *fname_fit  ;          /* filename root for fit coefficients */
00057     MRI_IMARR * ims ;           /* array of images themselves */
00058     MRI_IMARR * rims ;          /* array of registered images */
00059     float scale_fim ,           /* scale factor for alpha image */
00060           thresh_pcorr ,        /* correlation coeff. threshold */
00061           thresh_report ;       /* reporting threshold for output */
00062     int nxim , nyim ,           /* image dimensions */
00063         ntime ;                 /* number of time points */
00064     time_series *weight ;       /* pointer to weighting vector */
00065     time_series_array * refts,  /* array of reference (detrending) vectors */
00066                       * idts  ; /* array of ideal vectors */
00067     int flim ,                  /* if TRUE, write floating point images */
00068         norm_fim ,              /* if TRUE, normalize alpha image */
00069         scale_output ;          /* if TRUE, force scale data into outputs' corners */
00070     int dfilt_code ;            /* one of the DFILT_ consts above */
00071     int reg_bilinear ;          /* 1 for bilinear registration, 0 for bicubic */
00072     MRI_IMAGE * first_flim ;    /* first image to use, in flim format */
00073     int do_clip , debug , quiet ;
00074     char * fname_subort ;
00075 } line_opt ;
00076 
00077 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00078 
00079 /*** internal prototypes ***/
00080 
00081 void get_line_opt( int argc , char *argv[] , line_opt *opt ) ;
00082 
00083 void Syntax( char * ) ;
00084 
00085 void write_images( line_opt * , char * , MRI_IMAGE * ,
00086                    float      , char * , MRI_IMAGE *  ) ;
00087 
00088 time_series * edit_weight( time_series * , time_series * ) ;
00089 
00090 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00091 
00092 int main( int argc , char *argv[] )
00093 {
00094    line_opt   opt ;        /* holds data constructed from command line */
00095    references *ref ;       /* holds reference covariance computations */
00096    voxel_corr *voxcor ;    /* holds voxel-specific covariances */
00097    float      *pcorr ;     /* holds partial correlation coefficients */
00098    float      *alpha ;     /* holds activation levels */
00099    float      *refvec ,    /* values of reference vectors at given time */
00100               *voxvec ;    /* values of voxel data at given time */
00101    MRI_IMAGE  *voxim ;     /* image data (contains voxvec) */
00102    MRI_IMAGE  *imcor ,     /* correlation image (pcorr is part of this) */
00103               *imalp  ;    /* activation level image (alpha is in this ) */
00104    int nref , nideal ;
00105 
00106    float      wt ;         /* weight factor at given time */
00107 
00108    /* duplicates for rims */
00109 
00110    references *ref_reg ;
00111    voxel_corr *voxcor_reg ;
00112    float      *pcorr_reg , *alpha_reg ;
00113    MRI_IMAGE  *imcor_reg , *imalp_reg  ;
00114    int nref_reg ;
00115 
00116    int   itime , numvox , good , vox , ii , ii_idts ;
00117    int   do_cnr , do_sig , do_fit , do_subort ;
00118 
00119    time_series * wtnew ;
00120 
00121    MRI_IMARR * cor_ims , * alp_ims , * cnr_ims , * sig_ims ;
00122    MRI_IMARR ** fit_ims ;  /* a whole bunch of arrays of images! */
00123    MRI_IMAGE * best_im ;
00124    int       * best ;
00125    float     * cornew ;
00126    float       cnew , cold ;
00127 
00128    int   npass_pos , npass_neg ;   /* number of voxels that pass the test */
00129    float cormin , cormax ,         /* min and max correlations seen */
00130          alpmin , alpmax ;         /* ditto for activations */
00131 
00132 /*** read command line, and set up as it bids ***/
00133 
00134    machdep() ;
00135 
00136    get_line_opt( argc , argv , &opt ) ;
00137 
00138    numvox = opt.nxim * opt.nyim ;
00139    nref   = opt.refts->num ;
00140    nideal = opt.idts->num ;
00141 
00142    do_cnr    = opt.fname_cnr    != NULL ;
00143    do_sig    = opt.fname_sig    != NULL ;
00144    do_fit    = opt.fname_fit    != NULL ;
00145    do_subort = opt.fname_subort != NULL ;
00146 
00147    /** create place to put output images for each ideal **/
00148 
00149    INIT_IMARR(cor_ims) ;
00150    INIT_IMARR(alp_ims) ;
00151    if( do_cnr ) INIT_IMARR(cnr_ims) ;
00152    if( do_sig ) INIT_IMARR(sig_ims) ;
00153    if( do_fit || do_subort )
00154       fit_ims = (MRI_IMARR **) malloc( sizeof(MRI_IMARR *) * nideal ) ;
00155 
00156    /*** June 1995: set up to clip on intensities of first_flim ***/
00157 
00158    if( CLIP_FRAC > 0.0 && opt.do_clip ){
00159       float ftop , fbot , clipper , val ;
00160       float * far ;
00161       int nclipper ;
00162 
00163       ftop = mri_max( opt.first_flim ) ;
00164       fbot = mri_min( opt.first_flim ) ;
00165       ftop = (fabs(ftop) > fabs(fbot)) ? fabs(ftop) : fabs(fbot) ;
00166       ftop = CLIP_FRAC * ftop ;
00167       far  = MRI_FLOAT_PTR( opt.first_flim ) ;
00168 
00169       clipper  = 0.0 ;
00170       nclipper = 0 ;
00171       for( vox=0 ; vox < numvox ; vox++ ){
00172          val = fabs(far[vox]) ;
00173          if( val >= ftop ){ clipper += val ; nclipper++ ; }
00174       }
00175       clipper  = CLIP_FRAC * clipper / nclipper ;
00176       nclipper = 0 ;
00177       for( vox=0 ; vox < numvox ; vox++ ){
00178          val = fabs(far[vox]) ;
00179          if( val < clipper ){ far[vox] = 0.0 ; nclipper++ ; }
00180       }
00181       if( nclipper > 0 && !opt.quiet )
00182          printf("CLIPPING %d voxels to zero for low intensity in base image!\n",nclipper) ;
00183    }
00184 
00185    refvec = (float *) malloc( sizeof(float) * nref ) ;
00186    if( refvec == NULL ) MALLOC_ERR("refvec") ;
00187 
00188    /** July 1995: master loop over multiple ideals **/
00189 
00190    for( ii_idts=0 ; ii_idts < nideal ; ii_idts++ ){
00191 
00192       opt.refts->tsarr[nref-1] = opt.idts->tsarr[ii_idts] ;          /* last ref = new ideal */
00193       wtnew = edit_weight( opt.weight , opt.idts->tsarr[ii_idts] ) ; /* make weight vector */
00194 
00195       if( wtnew == NULL ){
00196          fprintf(stderr,"** bad weight at ideal # %d -- end of run!\a\n",ii_idts+1) ;
00197          exit(1) ;
00198       }
00199 
00200       ref    = new_references( nref ) ;
00201       voxcor = new_voxel_corr( numvox , nref ) ;
00202 
00203       if( opt.rims != NULL ){
00204          nref_reg   = nref - 3 ;                    /* don't use the 1st 3 refs */
00205          ref_reg    = new_references( nref_reg ) ;
00206          voxcor_reg = new_voxel_corr( numvox , nref_reg ) ;
00207       }
00208 
00209       /*** loop through time,
00210            getting new reference vector data, and
00211            new images to correlate with the references vectors ***/
00212 
00213       itime = 0 ;
00214       do{
00215          int wt_not_one ;
00216 
00217          /*** find weighting factor for this time ***/
00218 
00219          wt = wtnew->ts[itime] ;
00220 
00221          /*** process new image for this time (if any weight is given to it) ***/
00222 
00223          if( wt != 0.0 ){
00224             wt_not_one = (wt != 1.0) ;
00225             voxim      = (wt_not_one) ? mri_to_float( opt.ims->imarr[itime] )  /* copy */
00226                                       : opt.ims->imarr[itime] ;                /* pointer */
00227             voxvec = MRI_FLOAT_PTR( voxim ) ;
00228             for( ii=0 ; ii < nref ; ii++ )
00229                refvec[ii] = opt.refts->tsarr[ii]->ts[itime] ;
00230             if( wt_not_one ){
00231                for( vox=0 ; vox < nref   ; vox++ ) refvec[vox] *= wt ;
00232                for( vox=0 ; vox < numvox ; vox++ ) voxvec[vox] *= wt ;
00233             }
00234             update_references( refvec , ref ) ;
00235             update_voxel_corr( voxvec , ref , voxcor ) ;
00236             if( wt_not_one ) mri_free( voxim ) ;
00237 
00238             /*** same for registered images, but don't use refs 0, 1, or 2 ***/
00239 
00240             if( opt.rims != NULL ){
00241                voxim      = (wt_not_one) ? mri_to_float( opt.rims->imarr[itime] )  /* copy */
00242                                          : opt.rims->imarr[itime] ;                /* pointer */
00243                voxvec = MRI_FLOAT_PTR( voxim ) ;
00244                for( ii=0 ; ii < nref_reg ; ii++ )
00245                   refvec[ii] = opt.refts->tsarr[ii+3]->ts[itime] ;
00246                if( wt_not_one ){
00247                   for( vox=0 ; vox < nref_reg ; vox++ ) refvec[vox] *= wt ;
00248                   for( vox=0 ; vox < numvox   ; vox++ ) voxvec[vox] *= wt ;
00249                }
00250                update_references( refvec , ref_reg ) ;
00251                update_voxel_corr( voxvec , ref_reg , voxcor_reg ) ;
00252                if( wt_not_one ) mri_free( voxim ) ;
00253             }
00254          }
00255 
00256       } while( ++itime < opt.ntime ) ;
00257 
00258       /*** compute outputs ***/
00259 
00260       imcor = mri_new( opt.nxim , opt.nyim , MRI_float ) ;   /* output images */
00261       imalp = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
00262       pcorr = MRI_FLOAT_PTR( imcor ) ;                       /* output arrays */
00263       alpha = MRI_FLOAT_PTR( imalp ) ;
00264 
00265       get_pcor( ref , voxcor , pcorr ) ;   /* compute partial correlation */
00266       get_coef( ref , voxcor , alpha ) ;  /* compute activation levels */
00267 
00268       /*** if have images AND registered images (DFBOTH), merge them ***/
00269 
00270       if( opt.rims != NULL ){
00271          float pc , pcreg ;
00272 
00273          imcor_reg  = mri_new( opt.nxim , opt.nyim , MRI_float ) ;   /* output images */
00274          imalp_reg  = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
00275          pcorr_reg  = MRI_FLOAT_PTR( imcor ) ;                       /* output arrays */
00276          alpha_reg  = MRI_FLOAT_PTR( imalp ) ;
00277 
00278          get_pcor( ref_reg , voxcor_reg , pcorr_reg ) ;
00279          get_coef( ref_reg , voxcor_reg , alpha_reg ) ;
00280 
00281          for( ii=0 ; ii < numvox ; ii++ ){
00282             pc = pcorr[ii] ; pcreg = pcorr_reg[ii] ;
00283             if( fabs(pc) > fabs(pcreg) ){
00284                pcorr[ii] = pcreg ; alpha[ii] = alpha_reg[ii] ;
00285             }
00286          }
00287 
00288          mri_free( imalp_reg ) ; mri_free( imcor_reg ) ;
00289          free_references( ref_reg ) ; free_voxel_corr( voxcor_reg ) ;
00290       }
00291 
00292       /*** June 1995: clip on intensities of first_flim ***/
00293 
00294       if( CLIP_FRAC > 0.0 && opt.do_clip ){
00295          float * far = MRI_FLOAT_PTR( opt.first_flim ) ;
00296          for( vox=0 ; vox < numvox ; vox++ )
00297             if( far[vox] == 0.0 ){ pcorr[vox] = alpha[vox] = 0.0 ; }
00298       }
00299 
00300       /*** store images for later processing ***/
00301 
00302       ADDTO_IMARR( cor_ims , imcor ) ;
00303       ADDTO_IMARR( alp_ims , imalp ) ;
00304 
00305       /*** compute CNR and SIG, if desired ***/
00306 
00307       if( do_cnr || do_sig ){
00308          MRI_IMAGE * imcnr , * imsig ;
00309          float     * cnr   , * sig ;
00310          float rbot,rtop , scale , sbot ;
00311          int ii , first , its ;
00312 
00313          first = 1 ;
00314          rbot  = rtop = 0 ;
00315          its   = nref - 1 ;         /* index of ideal */
00316 
00317          for( ii=0 ; ii < opt.ntime ; ii++ ){
00318             if( wtnew->ts[ii] > 0.0 ){
00319                if( first ){
00320                   rtop  = rbot = opt.refts->tsarr[its]->ts[ii] ;
00321                   first = 0 ;
00322                } else {
00323                   rbot = MIN( opt.refts->tsarr[its]->ts[ii] , rbot ) ;
00324                   rtop = MAX( opt.refts->tsarr[its]->ts[ii] , rtop ) ;
00325                }
00326             }
00327          }
00328          scale = rtop-rbot ;
00329 
00330          imsig = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
00331          sig   = MRI_FLOAT_PTR(imsig) ;
00332          get_variance( voxcor , sig ) ;
00333          sbot = 0.0 ;
00334          for( vox=0 ; vox < numvox ; vox++ )
00335             if( sig[vox] > 0.0 ){ sig[vox] = sqrt(sig[vox]) ; sbot += sig[vox] ; }
00336             else                  sig[vox] = 0.0 ;
00337          sbot = 0.001 * sbot / numvox ;   /* for clipping cnr */
00338 
00339          if( do_cnr ){
00340             imcnr = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
00341             cnr   = MRI_FLOAT_PTR(imcnr) ;
00342             for( vox=0 ; vox < numvox ; vox++ )
00343                if( sig[vox] > sbot ) cnr[vox] = scale * alpha[vox] / sig[vox] ;
00344                else                  cnr[vox] = 0.0 ;
00345             ADDTO_IMARR( cnr_ims , imcnr ) ;
00346 
00347 #ifdef DEBUG
00348    { char buf[64] ;
00349      printf("ideal %d: sbot = %g\n",ii_idts,sbot) ;
00350      sprintf(buf,"cnr.%02d",ii_idts) ; mri_write(buf,imcnr) ;
00351      sprintf(buf,"sig.%02d",ii_idts) ; mri_write(buf,imsig) ;
00352      sprintf(buf,"alp.%02d",ii_idts) ; mri_write(buf,imalp) ;
00353    }
00354 #endif
00355          }
00356 
00357          if( do_sig ) ADDTO_IMARR( sig_ims , imsig ) ;
00358          else         mri_free(imsig) ;
00359       }
00360 
00361       /** save fit coefficients for the -fitfile option **/
00362 
00363       if( do_fit || do_subort ){
00364          MRI_IMARR * fitim ;
00365          MRI_IMAGE * tim ;
00366          float ** fitar ;
00367          int ir ;
00368 
00369          INIT_IMARR(fitim) ;  /* create array of fit coefficients */
00370                               /* (one for each ref vector) */
00371 
00372          fitar = (float **) malloc( sizeof(float *) * nref ) ;
00373          for( ir=0 ; ir < nref ; ir++ ){
00374             tim = mri_new( opt.nxim , opt.nyim , MRI_float ) ; /* ir-th fit image */
00375             fitar[ir] = MRI_FLOAT_PTR(tim) ;                   /* ir-th array */
00376             ADDTO_IMARR(fitim,tim) ;
00377          }
00378 
00379          get_lsqfit( ref , voxcor , fitar ) ; /* compute all fit arrays at once */
00380          fit_ims[ii_idts] = fitim ;           /* add a whole new image array */
00381 
00382          free( fitar ) ;  /* don't need these pointers, are stored in fit_ims */
00383       }
00384 
00385       /*** cleanup for next ideal ***/
00386 
00387       free_references( ref ) ; free_voxel_corr( voxcor ) ;
00388       RWC_free_time_series( wtnew ) ;
00389 
00390    }  /*** end of loop over ideals (ii_idts) ***/
00391 
00392    /** release images and other data (unless they are needed below) **/
00393 
00394    if( !do_subort )       DESTROY_IMARR(opt.ims) ;
00395    if( opt.rims != NULL ) DESTROY_IMARR(opt.rims) ;
00396 
00397    opt.refts->tsarr[nref-1] = NULL ;           /* this one is in opt.idts */
00398    if( !do_subort ) DESTROY_TSARR(opt.refts) ;
00399    DESTROY_TSARR(opt.idts) ;
00400 
00401    mri_free( opt.first_flim ) ;
00402    free( refvec ) ;
00403    RWC_free_time_series( opt.weight ) ;
00404 
00405    /*************** scan through all ideals and pick the "best" ***************/
00406 
00407    if( nideal == 1 ){
00408       imcor = cor_ims->imarr[0] ;
00409       imalp = alp_ims->imarr[0] ;
00410    } else {
00411       imcor   = mri_to_float( cor_ims->imarr[0] ) ;          /* to be best correlation */
00412       pcorr   = MRI_FLOAT_PTR( imcor ) ;
00413 
00414       best_im = mri_new( opt.nxim , opt.nxim , MRI_int ) ;   /* to be index of best */
00415       best    = MRI_INT_PTR(best_im) ;
00416       for( vox=0 ; vox < numvox ; vox++ ) best[vox] = 0 ;    /* start at first ideal */
00417 
00418       /** find best correlation image **/
00419 
00420       for( ii_idts=1 ; ii_idts < nideal ; ii_idts++ ){
00421          cornew = MRI_FLOAT_PTR( cor_ims->imarr[ii_idts] ) ;  /* ii_idts'th correlation */
00422          for( vox=0 ; vox < numvox ; vox++ ){
00423             cnew = cornew[vox] ; cold = pcorr[vox] ;
00424             if( fabs(cnew) > fabs(cold) ){
00425                best[vox]  = ii_idts ; pcorr[vox] = cnew ;
00426             }
00427          }
00428       }
00429 
00430       /** load best alpha image **/
00431 
00432       imalp = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
00433       alpha = MRI_FLOAT_PTR( imalp ) ;
00434       for( vox=0 ; vox < numvox ; vox++ )
00435          alpha[vox] = MRI_FLOAT_PTR( alp_ims->imarr[best[vox]] )[vox] ;
00436    }
00437 
00438    /*** write correlation and alpha out ***/
00439 
00440    write_images( &opt ,            opt.fname_corr , imcor ,
00441                  opt.thresh_pcorr, opt.fname_fim  , imalp  ) ;
00442 
00443    /*** write a report to the screen ***/
00444 
00445    npass_pos = npass_neg = 0 ;
00446 
00447    for( vox=0 ; vox < numvox ; vox++ )
00448       if( fabs(pcorr[vox]) >= opt.thresh_pcorr ){
00449          if( pcorr[vox] > 0 ) npass_pos++ ;
00450          if( pcorr[vox] < 0 ) npass_neg++ ;
00451       }
00452 
00453    cormin = mri_min( imcor ) ; cormax = mri_max( imcor ) ;
00454    alpmin = mri_min( imalp ) ; alpmax = mri_max( imalp ) ;
00455 
00456    if( !opt.quiet ){
00457       printf( "minimum activation  = %11.4g  maximum = %11.4g\n", alpmin,alpmax);
00458       printf( "minimum correlation = %11.4g  maximum = %11.4g\n", cormin,cormax);
00459       printf( "number of voxels with correlation >=  %5.3f is %d\n",
00460               opt.thresh_pcorr , npass_pos ) ;
00461       printf( "number of voxels with correlation <= -%5.3f is %d\n",
00462               opt.thresh_pcorr , npass_neg ) ;
00463    }
00464 
00465    DESTROY_IMARR( cor_ims ) ;
00466    DESTROY_IMARR( alp_ims ) ;
00467    if( nideal > 1 ){ mri_free(imcor) ; mri_free(imalp) ; }
00468 
00469    /*** write CNR out ***/
00470 
00471    if( do_cnr ){
00472       MRI_IMAGE * imcnr ;
00473       float * cnr ;
00474 
00475       if( nideal == 1 ){
00476          imcnr = cnr_ims->imarr[0] ;
00477       } else {
00478          imcnr = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
00479          cnr   = MRI_FLOAT_PTR( imcnr ) ;
00480          for( vox=0 ; vox < numvox ; vox++ )
00481             cnr[vox] = MRI_FLOAT_PTR( cnr_ims->imarr[best[vox]] )[vox] ;
00482       }
00483 
00484       if( opt.flim ){
00485          mri_write( opt.fname_cnr , imcnr ) ;
00486       } else {
00487          MRI_IMAGE * shim ;
00488          shim = mri_to_short( 100.0 , imcnr ) ;
00489          mri_write( opt.fname_cnr , shim ) ;
00490          mri_free( shim ) ;
00491       }
00492 
00493       DESTROY_IMARR( cnr_ims ) ;
00494       if( nideal > 1 ) mri_free(imcnr) ;
00495    }
00496 
00497    /*** write SIG out ***/
00498 
00499    if( do_sig ){
00500       MRI_IMAGE * imsig ;
00501       float * sig ;
00502 
00503       if( nideal == 1 ){
00504          imsig = sig_ims->imarr[0] ;
00505       } else {
00506          imsig = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
00507          sig   = MRI_FLOAT_PTR( imsig ) ;
00508          for( vox=0 ; vox < numvox ; vox++ )
00509             sig[vox] = MRI_FLOAT_PTR( sig_ims->imarr[best[vox]] )[vox] ;
00510       }
00511 
00512       mri_write( opt.fname_sig , imsig ) ;  /* always flim */
00513       DESTROY_IMARR( sig_ims ) ;
00514       if( nideal > 1 ) mri_free(imsig) ;
00515    }
00516 
00517    /*** write FIT out ***/
00518 
00519    if( do_fit || do_subort ){
00520       char root[128] , fname[128] ;
00521       int ir , ib ;
00522       MRI_IMAGE * tim ;
00523       float * tar , * qar ;
00524       float ortval ;
00525 
00526       if( do_fit ){
00527          strcpy(root,opt.fname_fit) ; ir = strlen(root) ;
00528          if( root[ir-1] != '.' ){ root[ir]   = '.' ; root[ir+1] = '\0' ; }
00529       }
00530 
00531       /* if have more than 1 ideal, must pick best fit and put into new image */
00532 
00533       if( nideal > 1 ){
00534          tim = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
00535          tar = MRI_FLOAT_PTR( tim ) ;
00536       }
00537 
00538       for( ir=0 ; ir < nref ; ir++ ){
00539          if( nideal == 1 ){
00540             tim = fit_ims[0]->imarr[ir] ;  /* 1 ideal --> output the 1 fit */
00541             tar = MRI_FLOAT_PTR( tim ) ;   /* ptr to coefficients */
00542          } else {
00543             for( vox=0 ; vox < numvox ; vox++ )
00544                tar[vox] = MRI_FLOAT_PTR( fit_ims[best[vox]]->imarr[ir] )[vox] ;
00545          }
00546 
00547          if( do_fit ){
00548             sprintf(fname,"%s%02d",root,ir+1) ;
00549             mri_write( fname , tim ) ;          /* always flim */
00550          }
00551 
00552          if( do_subort && ir < nref-1 ){  /* subtract ort # ir from images */
00553 
00554             for( itime=0 ; itime < opt.ntime ; itime++ ){   /* loop over time */
00555                ortval = opt.refts->tsarr[ir]->ts[itime] ;
00556                if( fabs(ortval) >= BLAST ) continue ;       /* skip this one */
00557                qar = MRI_FLOAT_PTR(opt.ims->imarr[itime]) ; /* current image */
00558                for( vox=0 ; vox < numvox ; vox++ )          /* loop over space */
00559                   qar[vox] -= ortval * tar[vox] ;           /* subtract ort */
00560 
00561             } /* end of loop over time */
00562          }
00563       } /* end of loop over ref vectors */
00564 
00565       for( ii_idts=0 ; ii_idts < nideal ; ii_idts++ )  /* don't need fits no more */
00566          DESTROY_IMARR(fit_ims[ii_idts]) ;
00567       free(fit_ims) ;
00568 
00569       if( nideal > 1 ) mri_free(tim) ;  /* don't need best fit image no more */
00570 
00571       if( do_subort ){
00572          strcpy(root,opt.fname_subort) ; ir = strlen(root) ;
00573          if( root[ir-1] != '.' ){ root[ir]   = '.' ; root[ir+1] = '\0' ; }
00574 
00575          for( itime=0 ; itime < opt.ntime ; itime++ ){
00576             sprintf(fname,"%s%04d",root,itime+1) ;
00577             mri_write( fname , opt.ims->imarr[itime] ) ; /* always flim */
00578          }
00579 
00580          DESTROY_IMARR(opt.ims) ;    /* no longer used */
00581          DESTROY_TSARR(opt.refts) ;
00582       }
00583    }
00584 
00585    /******** all done ********/
00586 
00587    if( nideal > 1 ) mri_free( best_im ) ;
00588    exit(0) ;
00589 }
00590 
00591 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00592 
00593 /*** write images out
00594      (inputs must be floating mrilib format) ***/
00595 
00596 void write_images( line_opt *opt ,    char *fname_imcor , MRI_IMAGE *imcor ,
00597                    float thresh_cor , char *fname_imalp , MRI_IMAGE *imalp  )
00598 {
00599    int vox ,
00600        numvox = imcor->nx * imcor->ny ;
00601 
00602    float *pcorr = imcor->im.float_data ,
00603          *alpha = imalp->im.float_data  ;
00604 
00605    MRI_IMAGE *shim ;  /* image of shorts to actually write */
00606    short     *shar ;  /* array inside shim */
00607 
00608    float sfac , alpmin , alpmax ;
00609 
00610 /*** threshold alpha on pcorr ***/
00611 
00612    if( thresh_cor > 0.0 ){
00613       for( vox=0 ; vox < numvox ; vox++ )
00614          if( fabs(pcorr[vox]) < thresh_cor ) alpha[vox] = 0.0 ;
00615    }
00616 
00617 /*** write output images ***/
00618 
00619    if( opt->flim ){  /* write floating point images [an advanced user :-)] */
00620 
00621       if( fname_imcor != NULL ){            /* write correlation image */
00622          mri_write( fname_imcor , imcor ) ;
00623       }
00624 
00625       if( fname_imalp != NULL ){            /* write activation image */
00626          mri_write( fname_imalp , imalp ) ;
00627       }
00628 
00629    } else {  /* write short images [a primitive user :-(] */
00630 
00631    /*** scale and write correlation image ***/
00632 
00633       if( fname_imcor != NULL ){
00634 
00635          sfac = SCALE ;
00636          shim = mri_to_short( sfac , imcor ) ;   /* scale to shorts */
00637 
00638          if( opt->scale_output ){
00639             shar    = MRI_SHORT_PTR( shim ) ;    /* access array in shim */
00640             shar[0] = 0 ;                           /* for display purposes */
00641             shar[1] = -SCALE ;
00642             shar[2] =  SCALE ;
00643          }
00644 
00645          mri_write( fname_imcor , shim ) ;   /* write short image */
00646          mri_free( shim ) ;                      /* toss it away */
00647       }
00648 
00649    /*** scale and write activation image ***/
00650 
00651       if( fname_imalp != NULL ){
00652 
00653          alpmin = mri_min( imalp ) ; alpmin = fabs(alpmin) ;  /* find */
00654          alpmax = mri_max( imalp ) ; alpmax = fabs(alpmax) ;  /* biggest */
00655          alpmax = (alpmin < alpmax) ? (alpmax) : (alpmin) ;   /* alpha */
00656 
00657          /*** default scaling (normalization) ***/
00658 
00659          if( opt->norm_fim ){  /* normalize alpha to +/-SCALE range */
00660             if( alpmax==0.0 ){   /* no data range! */
00661                sfac = SCALE ;
00662             } else {
00663                sfac = SCALE / alpmax ;
00664             }
00665          } else {
00666 
00667          /*** user input scale factor ***/
00668 
00669             sfac = opt->scale_fim ;
00670          }
00671 
00672          shim = mri_to_short( sfac , imalp ) ;    /* do the scaling */
00673          shar = MRI_SHORT_PTR( shim ) ;        /* get array of shorts */
00674 
00675          for( vox=0 ; vox < numvox ; vox++ )
00676                  if( shar[vox] >  SCALE ) shar[vox] =  SCALE ;
00677             else if( shar[vox] < -SCALE ) shar[vox] = -SCALE ;
00678 
00679          if( opt->scale_output ){
00680             shar    = MRI_SHORT_PTR( shim ) ;    /* access array in shim */
00681             shar[0] = 0 ;                           /* for display purposes */
00682             shar[1] = -SCALE ;
00683             shar[2] =  SCALE ;
00684          }
00685 
00686          mri_write( fname_imalp , shim ) ;      /* write short image */
00687          mri_free( shim ) ;                     /* and toss it */
00688       }
00689 
00690    }  /* endif (opt->flim) */
00691 
00692    return ;
00693 }
00694 
00695 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00696 
00697 /*** read and interpret command line arguments ***/
00698 
00699 #ifdef DEBUG
00700 #  define DBOPT(str) fprintf(stderr,"at option %s: %s\n",argv[nopt],str)
00701 #else
00702 #  define DBOPT(str) /* nothing */
00703 #endif
00704 
00705 void get_line_opt( int argc , char *argv[] , line_opt *opt )
00706 {
00707    int nopt , ii , pp , rr , bad ;
00708    int first_im , do_corr , num_im , num_time , polref ;
00709    float ftemp , fscal ;
00710    MRI_IMAGE *im ;
00711    time_series *ideal , *vec ;
00712    char err_note[128] , * regbase ;
00713 
00714 /*** give help if needed ***/
00715 
00716    if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) Syntax(NULL) ;
00717 
00718 /*** setup opt with defaults ***/
00719 
00720    opt->fname_corr    = NULL ;   /* no output names yet */
00721    opt->fname_fim     = NULL ;
00722    opt->fname_cnr     = NULL ;
00723    opt->fname_fit     = NULL ;
00724    opt->fname_subort  = NULL ;
00725    opt->fname_sig     = NULL ;
00726    opt->thresh_pcorr  = 0.70 ;   /* default 70% threshold */
00727    opt->thresh_report = 10.0 ;   /* not currently used */
00728    opt->nxim          = 0 ;
00729    opt->nyim          = 0 ;
00730    opt->ntime         = 0 ;
00731    opt->weight        = NULL ;   /* no weight vector yet */
00732    opt->ims           = NULL ;
00733    opt->rims          = NULL ;
00734    opt->flim          = FALSE ;  /* default: don't write float images */
00735    opt->scale_fim     = 1.0 ;    /* scale factor for writing short images */
00736                                  /* (used only if norm_fim is FALSE) */
00737    opt->norm_fim      = TRUE ;   /* default: normalize alpha image */
00738    opt->scale_output  = TRUE ;   /* default: force scales into output corners */
00739 
00740    opt->dfilt_code    = DFILT_NONE ;    /* default: do neither of these things */
00741    opt->reg_bilinear  = 0 ;             /* default: bicubic */
00742    opt->do_clip       = 0 ;
00743    opt->debug         = 0 ;
00744    opt->quiet         = 0 ;
00745 
00746    /*** initialize array of time series data ***/
00747 
00748    INIT_TSARR(opt->refts) ;
00749    INIT_TSARR(opt->idts) ;
00750 
00751 /*** set defaults in local variables ***/
00752 
00753    polref   = 0 ;         /* maximum order of polynomial reference vectors */
00754    first_im = 1 ;         /* first image to actually use */
00755    num_im   = 999999 ;    /* maximum number of images to use */
00756    do_corr  = FALSE ;     /* write correlation image out? */
00757    ideal    = NULL ;      /* time_series of ideal response vector */
00758    regbase  = NULL ;      /* pointer to name of registration image */
00759 
00760 /*** read command line arguments and process them:
00761       coding technique is to examine each argv, and if it matches
00762       something we expect (e.g., -flim), process it, then skip to
00763       the next loop through ('continue' statements in each strncmp if).
00764       Note that the 'while' will increment the argument index (nopt) ***/
00765 
00766    nopt = 1 ;
00767    do{  /* a while loop, continuing until we reach the MRI filenames */
00768 
00769       /** clip option **/
00770 
00771       if( strncmp(argv[nopt],"-clip",5) == 0 ){
00772          DBOPT("-clip") ;
00773          opt->do_clip = 1 ;
00774          continue ;
00775       }
00776 
00777       /** quiet option **/
00778 
00779       if( strncmp(argv[nopt],"-q",2) == 0 ){
00780          DBOPT("-q") ;
00781          opt->quiet = 1 ;
00782          continue ;
00783       }
00784 
00785       /** debug option **/
00786 
00787       if( strncmp(argv[nopt],"-debug",5) == 0 ){
00788          DBOPT("-debug") ;
00789          opt->debug = 1 ;
00790          continue ;
00791       }
00792 
00793       /** anti-motion filtering options **/
00794 
00795       if( strncmp(argv[nopt],"-bilinear",6) == 0 ){
00796          DBOPT("-bilinear") ;
00797          opt->reg_bilinear = 1 ;
00798          continue ;
00799       }
00800 
00801       if( strncmp(argv[nopt],"-regbase",6) == 0 ){
00802          DBOPT("-regbase") ;
00803          if( ++nopt >= argc ) Syntax("-regbase needs a filename!") ;
00804          regbase = argv[nopt] ;
00805          continue ;
00806       }
00807 
00808 #ifdef ALLOW_DFTIME
00809       if( strstr(argv[nopt],"-dftime") != NULL ){
00810          DBOPT("-dftime") ;
00811          opt->dfilt_code = DFILT_TIME ;
00812          if( strstr(argv[nopt],":0") != NULL ) opt->dfilt_code = DFILT_TIME0 ;
00813          continue ;
00814       }
00815 #endif
00816 
00817       if( strstr(argv[nopt],"-dfspace") != NULL && strstr(argv[nopt],"-dfspacetime") == NULL ){
00818          DBOPT("-dfspace") ;
00819          opt->dfilt_code = DFILT_SPACE ;
00820          if( strstr(argv[nopt],":0") != NULL ) opt->dfilt_code = DFILT_SPACE0 ;
00821          continue ;
00822       }
00823 
00824 #ifdef ALLOW_DFTIME
00825       if( strstr(argv[nopt],"-dfboth") != NULL ){
00826          DBOPT("-dfboth") ;
00827          opt->dfilt_code = DFILT_BOTH ;
00828          if( strstr(argv[nopt],":0") != NULL ) opt->dfilt_code = DFILT_BOTH0 ;
00829          continue ;
00830       }
00831 #endif
00832 
00833 #ifdef ALLOW_DFTIME
00834       if( strstr(argv[nopt],"-dfspacetime") != NULL ){
00835          DBOPT("-dfspacetime") ;
00836          opt->dfilt_code = DFILT_SPACETIME ;
00837          if( strstr(argv[nopt],":0") != NULL ) opt->dfilt_code = DFILT_SPACETIME0 ;
00838          continue ;
00839       }
00840 #endif
00841 
00842       /** -clean **/
00843 
00844       if( strncmp(argv[nopt],"-clean",5) == 0 ){
00845          DBOPT("-clean") ;
00846          opt->scale_output = FALSE;
00847          continue ;
00848       }
00849 
00850       /*** -flim  ==> output floating point images (default is shorts) ***/
00851 
00852       if( strncmp(argv[nopt],"-flim",5) == 0 ){
00853          DBOPT("-flim") ;
00854          opt->flim = TRUE ;
00855          continue ;
00856       }
00857 
00858       /*** -non  ==> don't normalize alpha image ***/
00859 
00860       if( strncmp(argv[nopt],"-non",4) == 0 ){
00861          DBOPT("-non") ;
00862          opt->norm_fim = FALSE ;
00863          continue ;
00864       }
00865 
00866       /*** -coef #  ==>  coefficient for normalizing alpha image ***/
00867 
00868       if( strncmp(argv[nopt],"-coe",4) == 0 ){
00869          DBOPT("-coef") ;
00870          if( ++nopt >= argc ) Syntax("-coef needs an argument") ;
00871          ftemp = strtod( argv[nopt] , NULL ) ;
00872          if( ftemp <= 0 ){
00873             sprintf( err_note , "illegal -coef value: %f" , ftemp ) ;
00874             Syntax(err_note) ;
00875          }
00876          opt->scale_fim = ftemp ;
00877          continue ;
00878       }
00879 
00880       /*** -list #  ==>  threshold for report listing (not used now) ***/
00881 
00882       if( strncmp(argv[nopt],"-list",4) == 0 ){
00883          DBOPT("-list") ;
00884          if( ++nopt >= argc ) Syntax("-list needs an argument") ;
00885          ftemp = strtod( argv[nopt] , NULL ) ;
00886 #if 0
00887          if( ftemp <= 0 ){
00888             sprintf( err_note , "illegal -list value: %f" , ftemp ) ;
00889             Syntax(err_note) ;
00890          }
00891 #else
00892          fprintf(stderr,"WARNING: -list option is ignored in fim2!\n") ;
00893 #endif
00894          opt->thresh_report = ftemp ;
00895          continue ;
00896       }
00897 
00898       /*** -polref #  ==>  order of polynomial reference functions ***/
00899 
00900       if( strncmp(argv[nopt],"-polref",5) == 0 || strncmp(argv[nopt],"-polort",5) == 0 ){
00901          DBOPT("-polref") ;
00902          if( ++nopt >= argc ) Syntax("-polref needs an argument") ;
00903          ftemp = strtod( argv[nopt] , NULL ) ;
00904          if( ftemp > 3 ){
00905             fprintf( stderr , "WARNING: large -polref value %f\n" , ftemp ) ;
00906          }
00907          polref = (int) ftemp ;
00908          continue ;
00909       }
00910 
00911       /*** -pcnt #  ==>  percent of deviation from perfect correlation ***/
00912 
00913       if( strncmp(argv[nopt],"-pcnt",4) == 0 ){
00914          DBOPT("-pcnt") ;
00915          if( ++nopt >= argc ) Syntax("-pcnt needs an argument") ;
00916          ftemp = strtod( argv[nopt] , NULL ) ;
00917          if( ftemp < 0.0 || ftemp > 100.0 ){
00918             sprintf( err_note , "illegal -pcnt value: %f" , ftemp ) ;
00919             Syntax(err_note) ;
00920          }
00921          opt->thresh_pcorr = 1.0 - 0.01 * ftemp ;
00922          continue ;
00923       }
00924 
00925       /*** -pcthresh #  ==>  actual threshold on \rho ***/
00926 
00927       if( strncmp(argv[nopt],"-pcthresh",5) == 0 ){
00928          DBOPT("-pcthresh") ;
00929          if( ++nopt >= argc ) Syntax("-pcthresh needs an argument") ;
00930          ftemp = strtod( argv[nopt] , NULL ) ;
00931          if( ftemp < 0.0 || ftemp > 1.0 ){
00932             sprintf( err_note , "illegal -pcthresh value: %f" , ftemp ) ;
00933             Syntax(err_note) ;
00934          }
00935          opt->thresh_pcorr = ftemp ;
00936          continue ;
00937       }
00938 
00939       /*** -im1 #  ==>  index (1...) of first image to actually use ***/
00940 
00941       if( strncmp(argv[nopt],"-im1",4) == 0 ){
00942          DBOPT("-im1") ;
00943          if( ++nopt >= argc ) Syntax("-im1 needs an argument") ;
00944          ftemp = strtod( argv[nopt] , NULL ) ;
00945          if( ftemp < 1.0 ){
00946             sprintf( err_note , "illegal -im1 value: %f" , ftemp ) ;
00947             Syntax(err_note) ;
00948          }
00949          first_im = (int)(ftemp+0.499) ;
00950          continue ;
00951       }
00952 
00953       /*** -num #  ==>  maximum number of images to use from command line ***/
00954 
00955       if( strncmp(argv[nopt],"-num",4) == 0 ){
00956          DBOPT("-num") ;
00957          if( ++nopt >= argc ) Syntax("-num needs an argument") ;
00958          ftemp = strtod( argv[nopt] , NULL ) ;
00959          if( ftemp < 2 ){
00960             sprintf( err_note , "illegal -num value: %f" , ftemp ) ;
00961             Syntax(err_note) ;
00962          }
00963          num_im = (int)(ftemp+0.499) ;
00964          continue ;
00965       }
00966 
00967 #define OPENERS "[{/%"
00968 #define CLOSERS "]}/%"
00969 
00970 #define IS_OPENER(sss) (strlen((sss))==1 && strstr(OPENERS,(sss))!=NULL)
00971 #define IS_CLOSER(sss) (strlen((sss))==1 && strstr(CLOSERS,(sss))!=NULL)
00972 
00973       /*** -ort file  ==>  reference time series ***/
00974 
00975       if( strncmp(argv[nopt],"-ort",4) == 0 ){
00976          DBOPT("-ort") ;
00977          if( ++nopt >= argc ) Syntax("-ort needs a filename") ;
00978 
00979          /** July 1995: read a bunch of orts using [ a b c ... ] format **/
00980 
00981          if( ! IS_OPENER(argv[nopt]) ){                                 /* one file */
00982             vec = RWC_read_time_series( argv[nopt] ) ;
00983             if( vec == NULL ){ fprintf( stderr , "cannot continue\a\n" ) ; exit(1) ; }
00984             ADDTO_TSARR( opt->refts , vec ) ;
00985          } else {                                                           /* many */
00986             for( nopt++ ; !IS_CLOSER(argv[nopt]) && nopt<argc ; nopt++ ){
00987                vec = RWC_read_time_series( argv[nopt] ) ;
00988                if( vec == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
00989                ADDTO_TSARR( opt->refts , vec ) ;
00990             }
00991             if( nopt == argc ) Syntax("-ort never finishes!") ;
00992          }
00993          continue ;
00994       }
00995 
00996       /*** -ideal file  ==>  ideal response vector time series ***/
00997 
00998       if( strncmp(argv[nopt],"-ideal",5) == 0 ){
00999          DBOPT("-ideal") ;
01000          if( ++nopt >= argc ) Syntax("-ideal needs a filename") ;
01001 
01002          /** July 1995: read a bunch of ideals using [ a b c ... ] format **/
01003 
01004          if( ! IS_OPENER(argv[nopt]) ){                                 /* one file */
01005             ideal = RWC_read_time_series( argv[nopt] ) ;
01006             if( ideal == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
01007             ADDTO_TSARR( opt->idts , ideal ) ;
01008          } else {                                                           /* many */
01009             for( nopt++ ; !IS_CLOSER(argv[nopt]) && nopt<argc ; nopt++ ){
01010                ideal = RWC_read_time_series( argv[nopt] ) ;
01011                if( ideal == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
01012                ADDTO_TSARR( opt->idts , ideal ) ;
01013             }
01014             if( nopt == argc ) Syntax("-ideal never finishes!") ;
01015          }
01016          continue ;
01017       }
01018 
01019       /*** -weight file  ==>  weight vector time series ***/
01020 
01021       if( strncmp(argv[nopt],"-weight",5) == 0 ){
01022          DBOPT("-weight") ;
01023          if( ++nopt >= argc ) Syntax("-weight needs a filename") ;
01024          if( opt->weight != NULL ){
01025             fprintf( stderr , "cannot have two weight vectors!\a\n" ) ;
01026             exit(1) ;
01027          }
01028          opt->weight = RWC_read_time_series( argv[nopt] ) ;
01029          if( opt->weight == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
01030          DBOPT("read in OK") ;
01031          continue ;
01032       }
01033 
01034       /*** -fimfile file  ==>  name of output file ***/
01035 
01036       if( strncmp(argv[nopt],"-fimfile",5) == 0 ){
01037          DBOPT("-fimfile") ;
01038          if( ++nopt >= argc ) Syntax("-fimfile needs a filename") ;
01039          if( opt->fname_fim != NULL ){
01040             fprintf( stderr , "cannot have two fim output files!\a\n" ) ;
01041             exit(1) ;
01042          }
01043          opt->fname_fim = argv[nopt] ;
01044          DBOPT("stored as fimfile") ;
01045          continue ;
01046       }
01047 
01048       /*** -corr  ==>  write correlation file as fimfile.CORR ***/
01049 
01050       if( strncmp(argv[nopt],"-corr",5) == 0 ){
01051          DBOPT("-corr") ;
01052          do_corr = TRUE ;
01053          continue ;
01054       }
01055 
01056       /*** -corfile file  ==>  write correlation file as this name ***/
01057 
01058       if( strncmp(argv[nopt],"-corfile",5)  == 0 ||
01059           strncmp(argv[nopt],"-corrfile",6) == 0   ){
01060 
01061          DBOPT("-corfile") ;
01062          if( ++nopt >= argc ) Syntax("-corfile needs a filename") ;
01063          if( opt->fname_corr != NULL ){
01064             fprintf( stderr , "cannot have two corr output files!\a\n" ) ;
01065             exit(1) ;
01066          }
01067          opt->fname_corr = argv[nopt] ;
01068          do_corr         = TRUE ;
01069          DBOPT("stored as corfile") ;
01070          continue ;
01071       }
01072 
01073       /*** -cnrfile file  ==>  write cnr file as this name ***/
01074 
01075       if( strncmp(argv[nopt],"-cnrfile",5)  == 0 ){
01076 
01077          DBOPT("-cnrfile") ;
01078          if( ++nopt >= argc ) Syntax("-cnrfile needs a filename") ;
01079          if( opt->fname_cnr != NULL ){
01080             fprintf( stderr , "cannot have two cnr output files!\a\n" ) ;
01081             exit(1) ;
01082          }
01083          opt->fname_cnr = argv[nopt] ;
01084          DBOPT("stored as cnrfile") ;
01085          continue ;
01086       }
01087 
01088       /*** -sigfile file  ==>  write sig file as this name ***/
01089 
01090       if( strncmp(argv[nopt],"-sigfile",5)  == 0 ){
01091 
01092          DBOPT("-sigfile") ;
01093          if( ++nopt >= argc ) Syntax("-sigfile needs a filename") ;
01094          if( opt->fname_sig != NULL ){
01095             fprintf( stderr , "cannot have two sig output files!\a\n" ) ;
01096             exit(1) ;
01097          }
01098          opt->fname_sig = argv[nopt] ;
01099          DBOPT("stored as sigfile") ;
01100          continue ;
01101       }
01102 
01103       /*** -fitfile file  ==>  write fit files as this name ***/
01104 
01105       if( strncmp(argv[nopt],"-fitfile",5)  == 0 ){
01106 
01107          DBOPT("-fitfile") ;
01108          if( ++nopt >= argc ) Syntax("-fitfile needs a filename") ;
01109          if( opt->fname_fit != NULL ){
01110             fprintf( stderr , "cannot have two fit output filenames!\a\n" ) ;
01111             exit(1) ;
01112          }
01113          opt->fname_fit = argv[nopt] ;
01114          DBOPT("stored as fitfile") ;
01115          continue ;
01116       }
01117 
01118       /*** -subort file ==> compute time series with orts removed ***/
01119 
01120       if( strncmp(argv[nopt],"-subort",5) == 0 ){
01121 
01122          DBOPT("-subort") ;
01123          if( ++nopt >= argc ) Syntax("-subort needs a filename") ;
01124          if( opt->fname_subort != NULL ){
01125             fprintf( stderr , "cannot have two subort output filenames!\a\n" ) ;
01126             exit(1) ;
01127          }
01128          opt->fname_subort = argv[nopt] ;
01129          DBOPT("stored as subortfile") ;
01130          continue ;
01131       }
01132 
01133       /*** get to here, and start with a '-'  ==>  bad news city, Arizona ***/
01134 
01135       if( strncmp(argv[nopt],"-",1) == 0 ){
01136          sprintf( err_note , "unknown option %s" , argv[nopt] ) ;
01137          Syntax(err_note) ;
01138       }
01139 
01140    /*** finished with switches.  Anything that makes it here is a filename ***/
01141 
01142       if( opt->idts->num == 0 ){   /* if ideal timeseries not given yet, is first */
01143          DBOPT("will be ideal") ;
01144          ideal = RWC_read_time_series( argv[nopt] ) ;
01145          if( ideal == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
01146          ADDTO_TSARR( opt->idts , ideal ) ;
01147          continue ;
01148       }
01149 
01150       /*** this point is the start of image filenames;
01151            break out of the loop and process them       ***/
01152 
01153       DBOPT("1st image file") ;
01154 
01155       break ;   /* time to leave */
01156 
01157    } while( ++nopt < argc ) ;  /* loop until all args are read, or we break */
01158 
01159 /***** when this loop ends, nopt = index of first image filename in argv[] *****/
01160 
01161    /*** check the situation for reasonabilityositiness ***/
01162 
01163    if( opt->idts->num == 0 ) Syntax("no ideal response vector is defined!") ;
01164 
01165    /*** compute the number of image files ***/
01166 
01167    num_time = argc - nopt ;                   /* # of arguments left */
01168    if( opt->fname_fim == NULL ){
01169       num_time -- ;                           /* one less, if need be */
01170       opt->fname_fim = argv[argc-1] ;         /* last name = output file */
01171    }
01172 
01173    /*** July 1995: read all images in now! ***/
01174 
01175    if( num_time < 1 ) Syntax("No input files given on command line?!") ;
01176 
01177    opt->ims = mri_read_many_files( num_time , argv+nopt ) ;
01178    if( opt->ims == NULL ) Syntax("Cannot read all images!") ;
01179    num_time = MIN( opt->ims->num , num_im ) ;
01180    if( num_time < 2 ) Syntax("must have at least 2 images to correlate!") ;
01181    opt->ntime = num_time ;
01182    opt->nxim  = opt->ims->imarr[0]->nx ;
01183    opt->nyim  = opt->ims->imarr[0]->ny ;
01184 
01185 #ifdef DEBUG
01186    fprintf(stderr,"num_time = %d  nx = %d  ny = %d\n",num_time,opt->nxim,opt->nyim) ;
01187 #endif
01188 
01189    /** convert images to float format **/
01190 
01191    for( ii=0 ; ii < num_time ; ii++ ){
01192       if( opt->ims->imarr[ii]->kind != MRI_float ){
01193          im = mri_to_float( opt->ims->imarr[ii] ) ;
01194          mri_free( opt->ims->imarr[ii] ) ;
01195          opt->ims->imarr[ii] = im ;
01196       }
01197       if( opt->ims->imarr[ii]->nx != opt->nxim ||
01198           opt->ims->imarr[ii]->ny != opt->nyim   ){
01199 
01200          fprintf(stderr,"** Image size mismatch at image # %d -- end of run!\a\n",ii+1) ;
01201          exit(1) ;
01202       }
01203    }
01204 
01205    if( first_im < 1 || first_im >= num_time ){
01206       sprintf( err_note ,
01207                "-im1 was %d, but only have %d images" , first_im , num_time ) ;
01208       Syntax(err_note) ;
01209    }
01210 
01211    /*** replace earliest images with copies of "first_im", if needed ***/
01212 
01213    for( ii=0 ; ii < first_im-1 ; ii++ ){
01214       im = mri_to_float( opt->ims->imarr[first_im-1] ) ;  /* copy data, not just pointer */
01215       mri_free( opt->ims->imarr[ii] ) ;
01216       opt->ims->imarr[ii] = im ;
01217    }
01218 
01219    if( do_corr && opt->fname_corr == NULL ){
01220 #ifdef DEBUG
01221       fprintf(stderr,"creating .CORR filename\n");
01222 #endif
01223       ii = strlen( opt->fname_fim ) + 7 ;
01224       opt->fname_corr = (char *) malloc( sizeof(char) * ii ) ;
01225       if( opt->fname_corr == NULL ) MALLOC_ERR(".CORR filename") ;
01226       strcpy( opt->fname_corr , opt->fname_fim ) ;
01227       strcat( opt->fname_corr , ".CORR" ) ;
01228    }
01229 
01230    /*** put the polynomial responses into opt->refts ***/
01231 
01232    if( polref >= 0 ){
01233 
01234 #ifdef DEBUG
01235       fprintf(stderr,"creating polynomial references; polref=%d\n",polref);
01236 #endif
01237 
01238       for( pp=0 ; pp <= polref ; pp++ ){
01239          vec = RWC_blank_time_series( num_time ) ;
01240 #if 1
01241          if( pp == 0 ){
01242             for( ii=0 ; ii < num_time ; ii++ ) vec->ts[ii] = 1.0 ;
01243          } else {
01244             fscal = 1.0 / num_time ;
01245             for( ii=0 ; ii < num_time ; ii++ ) vec->ts[ii] = pow(fscal*ii,pp) ;
01246          }
01247 #else
01248          ftemp = 0.5 * num_time - 0.4321 ;  /* 0.4321 ensures (ii-ftemp)!=0 */
01249          fscal = 2.0 / num_time ;           /* range of arg to pow is -1..1 */
01250          for( ii=0 ; ii < num_time ; ii++ )
01251             vec->ts[ii] = pow( fscal*(ii-ftemp) , pp ) ;
01252 #endif
01253 
01254          ADDTO_TSARR( opt->refts , vec ) ;
01255       }
01256    }
01257 
01258 /*** if no weight vector input (via -weight filename), make one up ***/
01259 
01260    if( opt->weight == NULL ){
01261 
01262 #ifdef DEBUG
01263       fprintf(stderr,"creating default weight\n");
01264 #endif
01265 
01266       vec = RWC_blank_time_series( num_time ) ;
01267 
01268       for( ii=0 ; ii < num_time   ; ii++ ) vec->ts[ii] = 1.0 ;  /* weight = all ones */
01269 
01270       opt->weight = vec ;
01271    }
01272    for( ii=0 ; ii < first_im-1 ; ii++ ) opt->weight->ts[ii] = 0.0 ;  /* except skipped images */
01273 
01274    /*** make space for the ideal vector in refts (but don't but it there yet) ***/
01275 
01276    ADDTO_TSARR( opt->refts , NULL ) ;
01277 
01278 /*** check all time series for being long enough ***/
01279 
01280    bad = FALSE ;
01281 
01282    for( ii=0 ; ii < opt->refts->num-1 ; ii++ ){  /* note not checking last one */
01283 #ifdef DEBUG
01284       fprintf(stderr,"checking ref %d for size\n",ii) ;
01285 #endif
01286       if( opt->refts->tsarr[ii]->len < num_time ){
01287          fprintf( stderr ,
01288                   "reference vector %s has %d elements: too short!\a\n" ,
01289                   opt->refts->tsarr[ii]->fname , opt->refts->tsarr[ii]->len ) ;
01290          bad = TRUE ;
01291       }
01292    }
01293 
01294    for( ii=0 ; ii < opt->idts->num-1 ; ii++ ){
01295 #ifdef DEBUG
01296       fprintf(stderr,"checking ideal %d for size\n",ii) ;
01297 #endif
01298       if( opt->idts->tsarr[ii]->len < num_time ){
01299          fprintf( stderr ,
01300                   "ideal vector %s has %d elements: too short!\a\n" ,
01301                   opt->idts->tsarr[ii]->fname , opt->idts->tsarr[ii]->len ) ;
01302          bad = TRUE ;
01303       }
01304    }
01305 
01306    if( opt->weight->len < num_time ){
01307       fprintf( stderr ,
01308                "weight vector %s has %d elements: too short!\a\n" ,
01309                opt->weight->fname , opt->weight->len ) ;
01310       bad = TRUE ;
01311    }
01312 
01313    if( bad ) exit(1) ;
01314 
01315 /*** zero out weight at any time where a refts time_series is blasted ***/
01316 
01317 #ifdef DEBUG
01318    fprintf(stderr,"blasting away ... ") ;
01319 #endif
01320 
01321    for( ii=0 ; ii < num_time ; ii++ ){
01322 
01323       bad = (opt->weight->ts[ii] >= BLAST) || (opt->weight->ts[ii] < 0.0) ;
01324 
01325       for( rr=0 ; !bad && rr < opt->refts->num-1 ; rr++ )
01326          bad = (opt->refts->tsarr[rr]->ts[ii] >= BLAST) ;
01327 
01328       if( bad ){
01329          opt->weight->ts[ii] = 0 ;
01330          for( rr=0 ; rr < opt->refts->num-1 ; rr++ )  /* note not checking last one */
01331             opt->refts->tsarr[rr]->ts[ii] = 0 ;
01332       }
01333    }
01334 
01335 #ifdef DEBUG
01336    fprintf(stderr,"\n") ;
01337 #endif
01338 
01339 /*** check to see if the edited vectors are still nonzero enough ***/
01340 
01341    bad = FALSE ;
01342 
01343 #ifdef DEBUG
01344    fprintf(stderr,"checking weight for nonnegativity\n") ;
01345 #endif
01346 
01347    ftemp = RWC_norm_ts( num_time , opt->weight ) ;
01348    if( ftemp < 1.e-9 ){
01349       fprintf( stderr , "there is no time with nonzero weighting!\n" ) ;
01350       bad = TRUE ;
01351    }
01352 
01353    for( rr=0 ; rr < opt->refts->num-1 ; rr++ ){  /* note not checking last one */
01354 #ifdef DEBUG
01355       fprintf(stderr,"checking ref %d for nonzeroness\n",rr) ;
01356 #endif
01357       ftemp = RWC_norm_ts( num_time , opt->refts->tsarr[rr] ) ;
01358       if( ftemp < 1.e-9 ){
01359          fprintf( stderr , "reference vector %d is all zeroes\n" , rr+1 ) ;
01360          bad = TRUE ;
01361       }
01362    }
01363    if( bad ) exit(1) ;
01364 
01365    /*** June 1995: get first image with nonzero weighting ***/
01366 
01367    if( regbase == NULL ){
01368       for( ii=0 ; ii < num_time ; ii++ )
01369          if( opt->weight->ts[ii] != 0.0 && ideal->ts[ii] < BLAST ) break ;
01370 
01371       if( ii == num_time ){ fprintf(stderr,"FIRST_IM: scan error!\a\n");exit(1); }
01372 
01373       opt->first_flim = mri_to_float( opt->ims->imarr[ii] ) ;  /* copy it */
01374    } else {
01375       MRI_IMAGE * qim ;
01376       qim = mri_read_just_one( regbase ) ;
01377       if( qim == NULL ) Syntax("Couldn't read -regbase image!") ;
01378       if( qim->kind == MRI_float ) opt->first_flim = qim ;
01379       else {
01380          opt->first_flim = mri_to_float( qim ) ;
01381          mri_free( qim ) ;
01382       }
01383       if( opt->first_flim->nx != opt->nxim || opt->first_flim->ny != opt->nyim ){
01384          fprintf(stderr,"-regbase: image size mismatch!\a\n") ; exit(1) ;
01385       }
01386    }
01387 
01388    /*** Setup for differential filtering (registration) ***/
01389 
01390    if( opt->dfilt_code != DFILT_NONE ){
01391       int alcode ;
01392       MRI_IMARR * tims ;
01393       time_series * dxts , * dyts , * phits ;
01394 
01395       switch( opt->dfilt_code ){
01396          case DFILT_TIME:   alcode = 0                                       ; break ;
01397 
01398          case DFILT_TIME0:  alcode = ALIGN_NOITER_CODE                       ; break ;
01399 
01400          case DFILT_SPACETIME:
01401          case DFILT_BOTH:
01402          case DFILT_SPACE:  alcode = ALIGN_REGISTER_CODE                     ; break ;
01403 
01404          case DFILT_SPACETIME0:
01405          case DFILT_BOTH0:
01406          case DFILT_SPACE0: alcode = ALIGN_REGISTER_CODE | ALIGN_NOITER_CODE ; break ;
01407 
01408          default:
01409             Syntax("Internal error: opt->dfilt_code illegal!") ;
01410       }
01411 
01412       dxts  = RWC_blank_time_series( num_time ) ;  /* space for motion params */
01413       dyts  = RWC_blank_time_series( num_time ) ;
01414       phits = RWC_blank_time_series( num_time ) ;
01415 
01416       if( opt->reg_bilinear ) alcode |= ALIGN_BILINEAR_CODE ;
01417       if( opt->debug        ) alcode |= ALIGN_DEBUG_CODE ;
01418 
01419       tims = mri_align_dfspace( opt->first_flim , NULL , opt->ims ,
01420                                 alcode , dxts->ts , dyts->ts , phits->ts ) ;
01421 
01422       switch( opt->dfilt_code ){
01423 
01424          case DFILT_TIME:
01425          case DFILT_TIME0:
01426             ADDTO_TSARR( opt->refts , NULL ) ;  /* make 3 blank spots at end */
01427             ADDTO_TSARR( opt->refts , NULL ) ;
01428             ADDTO_TSARR( opt->refts , NULL ) ;
01429 
01430             /* move previously existing time series up by 3 */
01431 
01432             for( ii=opt->refts->num-4 ; ii >=0 ; ii-- )
01433                opt->refts->tsarr[ii+3] = opt->refts->tsarr[ii] ;
01434 
01435             /* put dx,dy,phi at beginning of list */
01436 
01437             opt->refts->tsarr[0] = dxts ;
01438             opt->refts->tsarr[1] = dyts ;
01439             opt->refts->tsarr[2] = phits ;
01440          break ;
01441 
01442          case DFILT_SPACE:
01443          case DFILT_SPACE0:
01444             DESTROY_IMARR( opt->ims ) ;   /* put registered images in */
01445             opt->ims = tims ;             /* place of the input images */
01446 
01447 #if 0
01448 if( opt->debug ){
01449    char fff[64] ;
01450    for( ii=0 ; ii < IMARR_COUNT(opt->ims) ; ii++ ){
01451       sprintf(fff,"rrr.%04d",ii+1) ;
01452       mri_write( fff , IMARR_SUBIMAGE(opt->ims,ii) ) ;
01453    }
01454 }
01455 #endif
01456 
01457             RWC_free_time_series( dxts ) ;
01458             RWC_free_time_series( dyts ) ;
01459             RWC_free_time_series( phits ) ;
01460          break ;
01461 
01462          case DFILT_BOTH:
01463          case DFILT_BOTH0:
01464             ADDTO_TSARR( opt->refts , NULL ) ;
01465             ADDTO_TSARR( opt->refts , NULL ) ;
01466             ADDTO_TSARR( opt->refts , NULL ) ;
01467 
01468             for( ii=opt->refts->num-4 ; ii >=0 ; ii-- )
01469                opt->refts->tsarr[ii+3] = opt->refts->tsarr[ii] ;
01470 
01471             opt->refts->tsarr[0] = dxts ;
01472             opt->refts->tsarr[1] = dyts ;
01473             opt->refts->tsarr[2] = phits ;
01474 
01475             opt->rims = tims ;   /* keep registered images as a separate data set */
01476          break ;
01477 
01478          case DFILT_SPACETIME:
01479          case DFILT_SPACETIME0:
01480             ADDTO_TSARR( opt->refts , NULL ) ;
01481             ADDTO_TSARR( opt->refts , NULL ) ;
01482             ADDTO_TSARR( opt->refts , NULL ) ;
01483 
01484             for( ii=opt->refts->num-4 ; ii >=0 ; ii-- )
01485                opt->refts->tsarr[ii+3] = opt->refts->tsarr[ii] ;
01486 
01487             opt->refts->tsarr[0] = dxts ;
01488             opt->refts->tsarr[1] = dyts ;
01489             opt->refts->tsarr[2] = phits ;
01490 
01491             DESTROY_IMARR( opt->ims ) ;   /* put registered images in */
01492             opt->ims = tims ;             /* place of the input images */
01493          break ;
01494       }
01495    }
01496 
01497    /*** Done! ***/
01498 
01499    return ;
01500 }
01501 
01502 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
01503 
01504 /*** give some help and exit ***/
01505 
01506 void Syntax( char *note )
01507 {
01508    static char *mesg[] = {
01509    "Usage: fim2 [options] image_files ..." ,
01510    "where 'image_files ...' is a sequence of MRI filenames," ,
01511    " " ,
01512    "options are:",
01513    "-pcnt #         correlation coeff. threshold will be 1 - 0.01 * #",
01514    "-pcthresh #     correlation coeff. threshold will be #",
01515    "-im1 #          index of image file to use as first in time series;",
01516    "                  default is 1; previous images are filled with this",
01517    "                  image to synchronize with the reference time series",
01518    "-num #          number of images to actually use, if more than this",
01519    "                  many are specified on the command line;  default is",
01520    "                  to use all images",
01521    "-non            this option turns off the default normalization of",
01522    "                  the output activation image;  the user should provide",
01523    "                  a scaling factor via '-coef #', or '1' will be used",
01524    "-coef #         the scaling factor used to convert the activation output",
01525    "                  from floats to short ints (if -non is also present)",
01526    " ",
01527    "-ort fname      fname = filename of a time series to which the image data",
01528    "                  will be orthogonalized before correlations are computed;",
01529    "                  any number of -ort options (from 0 on up) may be used",
01530    "-ideal fname    fname = filename of a time series to which the image data",
01531    "                  is to be correlated; exactly one such time series is",
01532    "                  required;  if the -ideal option is not used, then the",
01533    "                  first filename after all the options will be used",
01534    "      N.B.: This version of fim2 allows the specification of more than",
01535    "            one ideal time series file.  Each one is separately correlated",
01536    "            with the image time series and the one most highly correlated",
01537    "            is selected for each pixel.  Multiple ideals are specified",
01538    "            using more than one '-ideal fname' option, or by using the",
01539    "            form '-ideal [ fname1 fname2 ... ]' -- this latter method",
01540    "            allows the use of wildcarded ideal filenames.",
01541    "            The '[' character that indicates the start of a group of",
01542    "            ideals can actually be any ONE of these: " OPENERS ,
01543    "            and the ']' that ends the group can be:  " CLOSERS ,
01544    " ",
01545    "      [Format of ort and ideal time series files:",
01546    "       ASCII; one number per line;",
01547    "       Same number of lines as images in the time series;",
01548    "       Value over 33333 --> don't use this image in the analysis]",
01549    " ",
01550    "-polref #       use polynomials of order 0..# as extra 'orts';",
01551    "[or -polort #]    default is 0 (yielding a constant vector).",
01552    "                  Use # = -1 to suppress this feature.",
01553 #if 0
01554    "-weight fname   fname = filename of a times series to be used as weights",
01555    "                  in the correlation calculation;  all time series",
01556    "                  (orts, ideal, and image) will be weighted at time i",
01557    "                  by the weight at that time;  if not given, defaults to",
01558    "                  1 at all times (but any ort or ideal >= 33333 gives 0)",
01559 #endif
01560    " ",
01561    "-fimfile fname  fname = filename to save activation magnitudes in;",
01562    "                  if not given, the last name on the command line will",
01563    "                  be used",
01564    "-corr           if present, indicates to write correlation output to",
01565    "                  image file 'fimfile.CORR' (next option is better)",
01566    "-corfile fname  fname = filename to save correlation image in;",
01567    "                  if not present, and -corr is not present, correlation",
01568    "                  image is not saved.",
01569    "-cnrfile fname  fname = filename to save contrast-to-noise image in;" ,
01570    "                  if not present, will not be computed or saved;" ,
01571    "                  CNR is scaled by 100 if images are output as shorts" ,
01572    "                  and is written 'as-is' if output as floats (see -flim)." ,
01573    "                  [CNR is defined here to be alpha/sigma, where" ,
01574    "                   alpha = amplitude of normalized ideal in a pixel" ,
01575    "                   sigma = standard deviation of pixel after removal" ,
01576    "                           of orts and ideal" ,
01577    "                   normalized ideal = ideal scaled so that trough-to-peak" ,
01578    "                     height is one.]" ,
01579    "-sigfile fname  fname = filename to save standard deviation image in;" ,
01580    "                  the standard deviation is of what is left after the" ,
01581    "                  least squares removal of the -orts, -polrefs, and -ideal." ,
01582    "                 N.B.: This is always output in the -flim format!" ,
01583    "-fitfile fname  Image files of the least squares fit coefficients of" ,
01584    "                  all the -ort and -polref time series that" ,
01585    "                  are projected out of the data time series before" ,
01586    "                  the -ideal is fit.  The actual filenames will" ,
01587    "                  be fname.01 fname.02 ...." ,
01588    "                  Their order is -orts, then -polrefs, and last -ideal." ,
01589    "                 N.B.: These are always output in the -flim format!" ,
01590    "-subort fname   A new timeseries of images is written to disk, with",
01591    "                  names of the form 'fname.0001', etc.  These images",
01592    "                  have the orts and polrefs (but not ideals) subtracted out.",
01593    "                 N.B.: These are always output in the -flim format!" ,
01594    "-flim           If present, write outputs in mrilib 'float' format,",
01595    "                  rather than scale+convert to integers",
01596    "                  [The 'ftosh' program can later convert to short integers]",
01597    "-clean          if present, then output images won't have the +/- 10000",
01598    "                  values forced into their corners for scaling purposes.",
01599    "-clip           if present, output correlations, etc., will be set to",
01600    "                  zero in regions of low intensity.",
01601    "-q              if present, indicates 'quiet' operation.",
01602    "-dfspace[:0]    Indicates to use the 'dfspace' filter (a la imreg) to",
01603    "                  register the images spatially before filtering." ,
01604 #ifdef ALLOW_DFTIME
01605    "-dftime[:0]     Indicates to use the 'dftime' filter (a la imreg) to",
01606    "                  produce 3 additional orts, in an attempt to reduce",
01607    "                  motion artifacts.",
01608    "-dfboth[:0]     Indicates to use both -dftime and -dfspace (separately)," ,
01609    "                  then take as the resulting correlation the smaller of the",
01610    "                  two results in each pixel (the conservative approach: " ,
01611    "                  to be above -pcthresh, both calculations must 'hit')." ,
01612    "                  The resulting fim is the one corresponding to the chosen" ,
01613    "                  correlation in each pixel." ,
01614    "-dfspacetime:[0] Indicates to use -dfspace and then -dftime together",
01615    "                  (not separately) on the time series of images." ,
01616 #endif
01617    "-regbase fname   Indicates to read image in file 'fname' as the base",
01618    "                  image for registration.  If not given, the first image",
01619    "                  in the time series that is used in the correlation",
01620    "                  computations will be used.  This is also the image",
01621    "                  that is used to define 'low intensity' for the -clip option.",
01622    " "
01623    } ;
01624 
01625    int nsize , ii ;
01626 
01627    if( note != NULL && (nsize=strlen(note)) > 0 ){
01628       fprintf(stderr,"\n") ;
01629       for(ii=0;ii<nsize;ii++) fprintf(stderr,"*") ;
01630       fprintf(stderr,"\a\n%s\n", note ) ;
01631       for(ii=0;ii<nsize;ii++) fprintf(stderr,"*") ;
01632       fprintf(stderr,"\ntry fim2 -help for instructions\a\n") ;
01633       exit(-1) ;
01634    }
01635 
01636    for( ii=0 ; ii < sizeof(mesg)/sizeof(char *) ; ii++ ){
01637       printf( " %s\n" , mesg[ii] );
01638    }
01639    exit(0) ;
01640 }
01641 
01642 /*----------------------------------------------------------------------------*/
01643 
01644 time_series * edit_weight( time_series * wt , time_series * ideal )
01645 {
01646    time_series * wtnew ;
01647    int ii , ntime ;
01648    float ftemp ;
01649 
01650    ntime = MIN(wt->len,ideal->len) ;
01651    wtnew = RWC_blank_time_series( ntime ) ;
01652 
01653    for( ii=0 ; ii < ntime ; ii++ )
01654       wtnew->ts[ii] = (ideal->ts[ii] >= BLAST) ? (0.0) : (wt->ts[ii]) ;
01655 
01656    ftemp = RWC_norm_ts( ntime , wtnew ) ;
01657    if( ftemp < 1.e-9 ){
01658       fprintf(stderr,"** bad ideal: no times with nonzero weight!\n") ;
01659       RWC_free_time_series( wtnew ) ;
01660       return NULL ;
01661    }
01662 
01663    return wtnew ;
01664 }