/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mrilib.h"
#include "pcor.h"
#include "ts.h"
#ifndef MAX
# define MAX(x,y) (((x)>(y)) ? (x) : (y))
#endif
#ifndef MIN
# define MIN(x,y) (((x)<(y)) ? (x) : (y))
#endif
#ifndef MALLOC_ERR
# define MALLOC_ERR(str) {fprintf(stderr,"malloc error for %s\a\n",str);exit(1);}
#endif
#define SCALE 10000
#define BLAST 33333.0
#define DFILT_NONE 0 /* no derivative filtering */
#define DFILT_SPACE 1 /* spatial mode */
#define DFILT_SPACE0 11
#define DFILT_TIME 2 /* temporal mode */
#define DFILT_TIME0 21
#define DFILT_BOTH 3 /* both, separately */
#define DFILT_BOTH0 31
#define DFILT_SPACETIME 4 /* both, together */
#define DFILT_SPACETIME0 41
#define CLIP_FRAC 0.10 /* for clipping nonbrain (low intensity) pixels */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/** structure type to hold results from scanning the command line options **/
typedef struct {
char *fname_corr , /* filename for correlation image */
*fname_fim , /* filename for activation image */
*fname_cnr , /* filename for contrast-to-noise image */
*fname_sig , /* filename for standard deviation image */
*fname_fit ; /* filename root for fit coefficients */
MRI_IMARR * ims ; /* array of images themselves */
MRI_IMARR * rims ; /* array of registered images */
float scale_fim , /* scale factor for alpha image */
thresh_pcorr , /* correlation coeff. threshold */
thresh_report ; /* reporting threshold for output */
int nxim , nyim , /* image dimensions */
ntime ; /* number of time points */
time_series *weight ; /* pointer to weighting vector */
time_series_array * refts, /* array of reference (detrending) vectors */
* idts ; /* array of ideal vectors */
int flim , /* if TRUE, write floating point images */
norm_fim , /* if TRUE, normalize alpha image */
scale_output ; /* if TRUE, force scale data into outputs' corners */
int dfilt_code ; /* one of the DFILT_ consts above */
int reg_bilinear ; /* 1 for bilinear registration, 0 for bicubic */
MRI_IMAGE * first_flim ; /* first image to use, in flim format */
int do_clip , debug , quiet ;
char * fname_subort ;
} line_opt ;
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*** internal prototypes ***/
void get_line_opt( int argc , char *argv[] , line_opt *opt ) ;
void Syntax( char * ) ;
void write_images( line_opt * , char * , MRI_IMAGE * ,
float , char * , MRI_IMAGE * ) ;
time_series * edit_weight( time_series * , time_series * ) ;
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
int main( int argc , char *argv[] )
{
line_opt opt ; /* holds data constructed from command line */
references *ref ; /* holds reference covariance computations */
voxel_corr *voxcor ; /* holds voxel-specific covariances */
float *pcorr ; /* holds partial correlation coefficients */
float *alpha ; /* holds activation levels */
float *refvec , /* values of reference vectors at given time */
*voxvec ; /* values of voxel data at given time */
MRI_IMAGE *voxim ; /* image data (contains voxvec) */
MRI_IMAGE *imcor , /* correlation image (pcorr is part of this) */
*imalp ; /* activation level image (alpha is in this ) */
int nref , nideal ;
float wt ; /* weight factor at given time */
/* duplicates for rims */
references *ref_reg ;
voxel_corr *voxcor_reg ;
float *pcorr_reg , *alpha_reg ;
MRI_IMAGE *imcor_reg , *imalp_reg ;
int nref_reg ;
int itime , numvox , good , vox , ii , ii_idts ;
int do_cnr , do_sig , do_fit , do_subort ;
time_series * wtnew ;
MRI_IMARR * cor_ims , * alp_ims , * cnr_ims , * sig_ims ;
MRI_IMARR ** fit_ims ; /* a whole bunch of arrays of images! */
MRI_IMAGE * best_im ;
int * best ;
float * cornew ;
float cnew , cold ;
int npass_pos , npass_neg ; /* number of voxels that pass the test */
float cormin , cormax , /* min and max correlations seen */
alpmin , alpmax ; /* ditto for activations */
/*** read command line, and set up as it bids ***/
machdep() ;
get_line_opt( argc , argv , &opt ) ;
numvox = opt.nxim * opt.nyim ;
nref = opt.refts->num ;
nideal = opt.idts->num ;
do_cnr = opt.fname_cnr != NULL ;
do_sig = opt.fname_sig != NULL ;
do_fit = opt.fname_fit != NULL ;
do_subort = opt.fname_subort != NULL ;
/** create place to put output images for each ideal **/
INIT_IMARR(cor_ims) ;
INIT_IMARR(alp_ims) ;
if( do_cnr ) INIT_IMARR(cnr_ims) ;
if( do_sig ) INIT_IMARR(sig_ims) ;
if( do_fit || do_subort )
fit_ims = (MRI_IMARR **) malloc( sizeof(MRI_IMARR *) * nideal ) ;
/*** June 1995: set up to clip on intensities of first_flim ***/
if( CLIP_FRAC > 0.0 && opt.do_clip ){
float ftop , fbot , clipper , val ;
float * far ;
int nclipper ;
ftop = mri_max( opt.first_flim ) ;
fbot = mri_min( opt.first_flim ) ;
ftop = (fabs(ftop) > fabs(fbot)) ? fabs(ftop) : fabs(fbot) ;
ftop = CLIP_FRAC * ftop ;
far = MRI_FLOAT_PTR( opt.first_flim ) ;
clipper = 0.0 ;
nclipper = 0 ;
for( vox=0 ; vox < numvox ; vox++ ){
val = fabs(far[vox]) ;
if( val >= ftop ){ clipper += val ; nclipper++ ; }
}
clipper = CLIP_FRAC * clipper / nclipper ;
nclipper = 0 ;
for( vox=0 ; vox < numvox ; vox++ ){
val = fabs(far[vox]) ;
if( val < clipper ){ far[vox] = 0.0 ; nclipper++ ; }
}
if( nclipper > 0 && !opt.quiet )
printf("CLIPPING %d voxels to zero for low intensity in base image!\n",nclipper) ;
}
refvec = (float *) malloc( sizeof(float) * nref ) ;
if( refvec == NULL ) MALLOC_ERR("refvec") ;
/** July 1995: master loop over multiple ideals **/
for( ii_idts=0 ; ii_idts < nideal ; ii_idts++ ){
opt.refts->tsarr[nref-1] = opt.idts->tsarr[ii_idts] ; /* last ref = new ideal */
wtnew = edit_weight( opt.weight , opt.idts->tsarr[ii_idts] ) ; /* make weight vector */
if( wtnew == NULL ){
fprintf(stderr,"** bad weight at ideal # %d -- end of run!\a\n",ii_idts+1) ;
exit(1) ;
}
ref = new_references( nref ) ;
voxcor = new_voxel_corr( numvox , nref ) ;
if( opt.rims != NULL ){
nref_reg = nref - 3 ; /* don't use the 1st 3 refs */
ref_reg = new_references( nref_reg ) ;
voxcor_reg = new_voxel_corr( numvox , nref_reg ) ;
}
/*** loop through time,
getting new reference vector data, and
new images to correlate with the references vectors ***/
itime = 0 ;
do{
int wt_not_one ;
/*** find weighting factor for this time ***/
wt = wtnew->ts[itime] ;
/*** process new image for this time (if any weight is given to it) ***/
if( wt != 0.0 ){
wt_not_one = (wt != 1.0) ;
voxim = (wt_not_one) ? mri_to_float( opt.ims->imarr[itime] ) /* copy */
: opt.ims->imarr[itime] ; /* pointer */
voxvec = MRI_FLOAT_PTR( voxim ) ;
for( ii=0 ; ii < nref ; ii++ )
refvec[ii] = opt.refts->tsarr[ii]->ts[itime] ;
if( wt_not_one ){
for( vox=0 ; vox < nref ; vox++ ) refvec[vox] *= wt ;
for( vox=0 ; vox < numvox ; vox++ ) voxvec[vox] *= wt ;
}
update_references( refvec , ref ) ;
update_voxel_corr( voxvec , ref , voxcor ) ;
if( wt_not_one ) mri_free( voxim ) ;
/*** same for registered images, but don't use refs 0, 1, or 2 ***/
if( opt.rims != NULL ){
voxim = (wt_not_one) ? mri_to_float( opt.rims->imarr[itime] ) /* copy */
: opt.rims->imarr[itime] ; /* pointer */
voxvec = MRI_FLOAT_PTR( voxim ) ;
for( ii=0 ; ii < nref_reg ; ii++ )
refvec[ii] = opt.refts->tsarr[ii+3]->ts[itime] ;
if( wt_not_one ){
for( vox=0 ; vox < nref_reg ; vox++ ) refvec[vox] *= wt ;
for( vox=0 ; vox < numvox ; vox++ ) voxvec[vox] *= wt ;
}
update_references( refvec , ref_reg ) ;
update_voxel_corr( voxvec , ref_reg , voxcor_reg ) ;
if( wt_not_one ) mri_free( voxim ) ;
}
}
} while( ++itime < opt.ntime ) ;
/*** compute outputs ***/
imcor = mri_new( opt.nxim , opt.nyim , MRI_float ) ; /* output images */
imalp = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
pcorr = MRI_FLOAT_PTR( imcor ) ; /* output arrays */
alpha = MRI_FLOAT_PTR( imalp ) ;
get_pcor( ref , voxcor , pcorr ) ; /* compute partial correlation */
get_coef( ref , voxcor , alpha ) ; /* compute activation levels */
/*** if have images AND registered images (DFBOTH), merge them ***/
if( opt.rims != NULL ){
float pc , pcreg ;
imcor_reg = mri_new( opt.nxim , opt.nyim , MRI_float ) ; /* output images */
imalp_reg = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
pcorr_reg = MRI_FLOAT_PTR( imcor ) ; /* output arrays */
alpha_reg = MRI_FLOAT_PTR( imalp ) ;
get_pcor( ref_reg , voxcor_reg , pcorr_reg ) ;
get_coef( ref_reg , voxcor_reg , alpha_reg ) ;
for( ii=0 ; ii < numvox ; ii++ ){
pc = pcorr[ii] ; pcreg = pcorr_reg[ii] ;
if( fabs(pc) > fabs(pcreg) ){
pcorr[ii] = pcreg ; alpha[ii] = alpha_reg[ii] ;
}
}
mri_free( imalp_reg ) ; mri_free( imcor_reg ) ;
free_references( ref_reg ) ; free_voxel_corr( voxcor_reg ) ;
}
/*** June 1995: clip on intensities of first_flim ***/
if( CLIP_FRAC > 0.0 && opt.do_clip ){
float * far = MRI_FLOAT_PTR( opt.first_flim ) ;
for( vox=0 ; vox < numvox ; vox++ )
if( far[vox] == 0.0 ){ pcorr[vox] = alpha[vox] = 0.0 ; }
}
/*** store images for later processing ***/
ADDTO_IMARR( cor_ims , imcor ) ;
ADDTO_IMARR( alp_ims , imalp ) ;
/*** compute CNR and SIG, if desired ***/
if( do_cnr || do_sig ){
MRI_IMAGE * imcnr , * imsig ;
float * cnr , * sig ;
float rbot,rtop , scale , sbot ;
int ii , first , its ;
first = 1 ;
rbot = rtop = 0 ;
its = nref - 1 ; /* index of ideal */
for( ii=0 ; ii < opt.ntime ; ii++ ){
if( wtnew->ts[ii] > 0.0 ){
if( first ){
rtop = rbot = opt.refts->tsarr[its]->ts[ii] ;
first = 0 ;
} else {
rbot = MIN( opt.refts->tsarr[its]->ts[ii] , rbot ) ;
rtop = MAX( opt.refts->tsarr[its]->ts[ii] , rtop ) ;
}
}
}
scale = rtop-rbot ;
imsig = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
sig = MRI_FLOAT_PTR(imsig) ;
get_variance( voxcor , sig ) ;
sbot = 0.0 ;
for( vox=0 ; vox < numvox ; vox++ )
if( sig[vox] > 0.0 ){ sig[vox] = sqrt(sig[vox]) ; sbot += sig[vox] ; }
else sig[vox] = 0.0 ;
sbot = 0.001 * sbot / numvox ; /* for clipping cnr */
if( do_cnr ){
imcnr = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
cnr = MRI_FLOAT_PTR(imcnr) ;
for( vox=0 ; vox < numvox ; vox++ )
if( sig[vox] > sbot ) cnr[vox] = scale * alpha[vox] / sig[vox] ;
else cnr[vox] = 0.0 ;
ADDTO_IMARR( cnr_ims , imcnr ) ;
#ifdef DEBUG
{ char buf[64] ;
printf("ideal %d: sbot = %g\n",ii_idts,sbot) ;
sprintf(buf,"cnr.%02d",ii_idts) ; mri_write(buf,imcnr) ;
sprintf(buf,"sig.%02d",ii_idts) ; mri_write(buf,imsig) ;
sprintf(buf,"alp.%02d",ii_idts) ; mri_write(buf,imalp) ;
}
#endif
}
if( do_sig ) ADDTO_IMARR( sig_ims , imsig ) ;
else mri_free(imsig) ;
}
/** save fit coefficients for the -fitfile option **/
if( do_fit || do_subort ){
MRI_IMARR * fitim ;
MRI_IMAGE * tim ;
float ** fitar ;
int ir ;
INIT_IMARR(fitim) ; /* create array of fit coefficients */
/* (one for each ref vector) */
fitar = (float **) malloc( sizeof(float *) * nref ) ;
for( ir=0 ; ir < nref ; ir++ ){
tim = mri_new( opt.nxim , opt.nyim , MRI_float ) ; /* ir-th fit image */
fitar[ir] = MRI_FLOAT_PTR(tim) ; /* ir-th array */
ADDTO_IMARR(fitim,tim) ;
}
get_lsqfit( ref , voxcor , fitar ) ; /* compute all fit arrays at once */
fit_ims[ii_idts] = fitim ; /* add a whole new image array */
free( fitar ) ; /* don't need these pointers, are stored in fit_ims */
}
/*** cleanup for next ideal ***/
free_references( ref ) ; free_voxel_corr( voxcor ) ;
RWC_free_time_series( wtnew ) ;
} /*** end of loop over ideals (ii_idts) ***/
/** release images and other data (unless they are needed below) **/
if( !do_subort ) DESTROY_IMARR(opt.ims) ;
if( opt.rims != NULL ) DESTROY_IMARR(opt.rims) ;
opt.refts->tsarr[nref-1] = NULL ; /* this one is in opt.idts */
if( !do_subort ) DESTROY_TSARR(opt.refts) ;
DESTROY_TSARR(opt.idts) ;
mri_free( opt.first_flim ) ;
free( refvec ) ;
RWC_free_time_series( opt.weight ) ;
/*************** scan through all ideals and pick the "best" ***************/
if( nideal == 1 ){
imcor = cor_ims->imarr[0] ;
imalp = alp_ims->imarr[0] ;
} else {
imcor = mri_to_float( cor_ims->imarr[0] ) ; /* to be best correlation */
pcorr = MRI_FLOAT_PTR( imcor ) ;
best_im = mri_new( opt.nxim , opt.nxim , MRI_int ) ; /* to be index of best */
best = MRI_INT_PTR(best_im) ;
for( vox=0 ; vox < numvox ; vox++ ) best[vox] = 0 ; /* start at first ideal */
/** find best correlation image **/
for( ii_idts=1 ; ii_idts < nideal ; ii_idts++ ){
cornew = MRI_FLOAT_PTR( cor_ims->imarr[ii_idts] ) ; /* ii_idts'th correlation */
for( vox=0 ; vox < numvox ; vox++ ){
cnew = cornew[vox] ; cold = pcorr[vox] ;
if( fabs(cnew) > fabs(cold) ){
best[vox] = ii_idts ; pcorr[vox] = cnew ;
}
}
}
/** load best alpha image **/
imalp = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
alpha = MRI_FLOAT_PTR( imalp ) ;
for( vox=0 ; vox < numvox ; vox++ )
alpha[vox] = MRI_FLOAT_PTR( alp_ims->imarr[best[vox]] )[vox] ;
}
/*** write correlation and alpha out ***/
write_images( &opt , opt.fname_corr , imcor ,
opt.thresh_pcorr, opt.fname_fim , imalp ) ;
/*** write a report to the screen ***/
npass_pos = npass_neg = 0 ;
for( vox=0 ; vox < numvox ; vox++ )
if( fabs(pcorr[vox]) >= opt.thresh_pcorr ){
if( pcorr[vox] > 0 ) npass_pos++ ;
if( pcorr[vox] < 0 ) npass_neg++ ;
}
cormin = mri_min( imcor ) ; cormax = mri_max( imcor ) ;
alpmin = mri_min( imalp ) ; alpmax = mri_max( imalp ) ;
if( !opt.quiet ){
printf( "minimum activation = %11.4g maximum = %11.4g\n", alpmin,alpmax);
printf( "minimum correlation = %11.4g maximum = %11.4g\n", cormin,cormax);
printf( "number of voxels with correlation >= %5.3f is %d\n",
opt.thresh_pcorr , npass_pos ) ;
printf( "number of voxels with correlation <= -%5.3f is %d\n",
opt.thresh_pcorr , npass_neg ) ;
}
DESTROY_IMARR( cor_ims ) ;
DESTROY_IMARR( alp_ims ) ;
if( nideal > 1 ){ mri_free(imcor) ; mri_free(imalp) ; }
/*** write CNR out ***/
if( do_cnr ){
MRI_IMAGE * imcnr ;
float * cnr ;
if( nideal == 1 ){
imcnr = cnr_ims->imarr[0] ;
} else {
imcnr = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
cnr = MRI_FLOAT_PTR( imcnr ) ;
for( vox=0 ; vox < numvox ; vox++ )
cnr[vox] = MRI_FLOAT_PTR( cnr_ims->imarr[best[vox]] )[vox] ;
}
if( opt.flim ){
mri_write( opt.fname_cnr , imcnr ) ;
} else {
MRI_IMAGE * shim ;
shim = mri_to_short( 100.0 , imcnr ) ;
mri_write( opt.fname_cnr , shim ) ;
mri_free( shim ) ;
}
DESTROY_IMARR( cnr_ims ) ;
if( nideal > 1 ) mri_free(imcnr) ;
}
/*** write SIG out ***/
if( do_sig ){
MRI_IMAGE * imsig ;
float * sig ;
if( nideal == 1 ){
imsig = sig_ims->imarr[0] ;
} else {
imsig = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
sig = MRI_FLOAT_PTR( imsig ) ;
for( vox=0 ; vox < numvox ; vox++ )
sig[vox] = MRI_FLOAT_PTR( sig_ims->imarr[best[vox]] )[vox] ;
}
mri_write( opt.fname_sig , imsig ) ; /* always flim */
DESTROY_IMARR( sig_ims ) ;
if( nideal > 1 ) mri_free(imsig) ;
}
/*** write FIT out ***/
if( do_fit || do_subort ){
char root[128] , fname[128] ;
int ir , ib ;
MRI_IMAGE * tim ;
float * tar , * qar ;
float ortval ;
if( do_fit ){
strcpy(root,opt.fname_fit) ; ir = strlen(root) ;
if( root[ir-1] != '.' ){ root[ir] = '.' ; root[ir+1] = '\0' ; }
}
/* if have more than 1 ideal, must pick best fit and put into new image */
if( nideal > 1 ){
tim = mri_new( opt.nxim , opt.nyim , MRI_float ) ;
tar = MRI_FLOAT_PTR( tim ) ;
}
for( ir=0 ; ir < nref ; ir++ ){
if( nideal == 1 ){
tim = fit_ims[0]->imarr[ir] ; /* 1 ideal --> output the 1 fit */
tar = MRI_FLOAT_PTR( tim ) ; /* ptr to coefficients */
} else {
for( vox=0 ; vox < numvox ; vox++ )
tar[vox] = MRI_FLOAT_PTR( fit_ims[best[vox]]->imarr[ir] )[vox] ;
}
if( do_fit ){
sprintf(fname,"%s%02d",root,ir+1) ;
mri_write( fname , tim ) ; /* always flim */
}
if( do_subort && ir < nref-1 ){ /* subtract ort # ir from images */
for( itime=0 ; itime < opt.ntime ; itime++ ){ /* loop over time */
ortval = opt.refts->tsarr[ir]->ts[itime] ;
if( fabs(ortval) >= BLAST ) continue ; /* skip this one */
qar = MRI_FLOAT_PTR(opt.ims->imarr[itime]) ; /* current image */
for( vox=0 ; vox < numvox ; vox++ ) /* loop over space */
qar[vox] -= ortval * tar[vox] ; /* subtract ort */
} /* end of loop over time */
}
} /* end of loop over ref vectors */
for( ii_idts=0 ; ii_idts < nideal ; ii_idts++ ) /* don't need fits no more */
DESTROY_IMARR(fit_ims[ii_idts]) ;
free(fit_ims) ;
if( nideal > 1 ) mri_free(tim) ; /* don't need best fit image no more */
if( do_subort ){
strcpy(root,opt.fname_subort) ; ir = strlen(root) ;
if( root[ir-1] != '.' ){ root[ir] = '.' ; root[ir+1] = '\0' ; }
for( itime=0 ; itime < opt.ntime ; itime++ ){
sprintf(fname,"%s%04d",root,itime+1) ;
mri_write( fname , opt.ims->imarr[itime] ) ; /* always flim */
}
DESTROY_IMARR(opt.ims) ; /* no longer used */
DESTROY_TSARR(opt.refts) ;
}
}
/******** all done ********/
if( nideal > 1 ) mri_free( best_im ) ;
exit(0) ;
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*** write images out
(inputs must be floating mrilib format) ***/
void write_images( line_opt *opt , char *fname_imcor , MRI_IMAGE *imcor ,
float thresh_cor , char *fname_imalp , MRI_IMAGE *imalp )
{
int vox ,
numvox = imcor->nx * imcor->ny ;
float *pcorr = MRI_FLOAT_PTR(imcor) ,
*alpha = MRI_FLOAT_PTR(imalp) ;
MRI_IMAGE *shim ; /* image of shorts to actually write */
short *shar ; /* array inside shim */
float sfac , alpmin , alpmax ;
/*** threshold alpha on pcorr ***/
if( thresh_cor > 0.0 ){
for( vox=0 ; vox < numvox ; vox++ )
if( fabs(pcorr[vox]) < thresh_cor ) alpha[vox] = 0.0 ;
}
/*** write output images ***/
if( opt->flim ){ /* write floating point images [an advanced user :-)] */
if( fname_imcor != NULL ){ /* write correlation image */
mri_write( fname_imcor , imcor ) ;
}
if( fname_imalp != NULL ){ /* write activation image */
mri_write( fname_imalp , imalp ) ;
}
} else { /* write short images [a primitive user :-(] */
/*** scale and write correlation image ***/
if( fname_imcor != NULL ){
sfac = SCALE ;
shim = mri_to_short( sfac , imcor ) ; /* scale to shorts */
if( opt->scale_output ){
shar = MRI_SHORT_PTR( shim ) ; /* access array in shim */
shar[0] = 0 ; /* for display purposes */
shar[1] = -SCALE ;
shar[2] = SCALE ;
}
mri_write( fname_imcor , shim ) ; /* write short image */
mri_free( shim ) ; /* toss it away */
}
/*** scale and write activation image ***/
if( fname_imalp != NULL ){
alpmin = mri_min( imalp ) ; alpmin = fabs(alpmin) ; /* find */
alpmax = mri_max( imalp ) ; alpmax = fabs(alpmax) ; /* biggest */
alpmax = (alpmin < alpmax) ? (alpmax) : (alpmin) ; /* alpha */
/*** default scaling (normalization) ***/
if( opt->norm_fim ){ /* normalize alpha to +/-SCALE range */
if( alpmax==0.0 ){ /* no data range! */
sfac = SCALE ;
} else {
sfac = SCALE / alpmax ;
}
} else {
/*** user input scale factor ***/
sfac = opt->scale_fim ;
}
shim = mri_to_short( sfac , imalp ) ; /* do the scaling */
shar = MRI_SHORT_PTR( shim ) ; /* get array of shorts */
for( vox=0 ; vox < numvox ; vox++ )
if( shar[vox] > SCALE ) shar[vox] = SCALE ;
else if( shar[vox] < -SCALE ) shar[vox] = -SCALE ;
if( opt->scale_output ){
shar = MRI_SHORT_PTR( shim ) ; /* access array in shim */
shar[0] = 0 ; /* for display purposes */
shar[1] = -SCALE ;
shar[2] = SCALE ;
}
mri_write( fname_imalp , shim ) ; /* write short image */
mri_free( shim ) ; /* and toss it */
}
} /* endif (opt->flim) */
return ;
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*** read and interpret command line arguments ***/
#ifdef DEBUG
# define DBOPT(str) fprintf(stderr,"at option %s: %s\n",argv[nopt],str)
#else
# define DBOPT(str) /* nothing */
#endif
void get_line_opt( int argc , char *argv[] , line_opt *opt )
{
int nopt , ii , pp , rr , bad ;
int first_im , do_corr , num_im , num_time , polref ;
float ftemp , fscal ;
MRI_IMAGE *im ;
time_series *ideal , *vec ;
char err_note[128] , * regbase ;
/*** give help if needed ***/
if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) Syntax(NULL) ;
/*** setup opt with defaults ***/
opt->fname_corr = NULL ; /* no output names yet */
opt->fname_fim = NULL ;
opt->fname_cnr = NULL ;
opt->fname_fit = NULL ;
opt->fname_subort = NULL ;
opt->fname_sig = NULL ;
opt->thresh_pcorr = 0.70 ; /* default 70% threshold */
opt->thresh_report = 10.0 ; /* not currently used */
opt->nxim = 0 ;
opt->nyim = 0 ;
opt->ntime = 0 ;
opt->weight = NULL ; /* no weight vector yet */
opt->ims = NULL ;
opt->rims = NULL ;
opt->flim = FALSE ; /* default: don't write float images */
opt->scale_fim = 1.0 ; /* scale factor for writing short images */
/* (used only if norm_fim is FALSE) */
opt->norm_fim = TRUE ; /* default: normalize alpha image */
opt->scale_output = TRUE ; /* default: force scales into output corners */
opt->dfilt_code = DFILT_NONE ; /* default: do neither of these things */
opt->reg_bilinear = 0 ; /* default: bicubic */
opt->do_clip = 0 ;
opt->debug = 0 ;
opt->quiet = 0 ;
/*** initialize array of time series data ***/
INIT_TSARR(opt->refts) ;
INIT_TSARR(opt->idts) ;
/*** set defaults in local variables ***/
polref = 0 ; /* maximum order of polynomial reference vectors */
first_im = 1 ; /* first image to actually use */
num_im = 999999 ; /* maximum number of images to use */
do_corr = FALSE ; /* write correlation image out? */
ideal = NULL ; /* time_series of ideal response vector */
regbase = NULL ; /* pointer to name of registration image */
/*** read command line arguments and process them:
coding technique is to examine each argv, and if it matches
something we expect (e.g., -flim), process it, then skip to
the next loop through ('continue' statements in each strncmp if).
Note that the 'while' will increment the argument index (nopt) ***/
nopt = 1 ;
do{ /* a while loop, continuing until we reach the MRI filenames */
/** clip option **/
if( strncmp(argv[nopt],"-clip",5) == 0 ){
DBOPT("-clip") ;
opt->do_clip = 1 ;
continue ;
}
/** quiet option **/
if( strncmp(argv[nopt],"-q",2) == 0 ){
DBOPT("-q") ;
opt->quiet = 1 ;
continue ;
}
/** debug option **/
if( strncmp(argv[nopt],"-debug",5) == 0 ){
DBOPT("-debug") ;
opt->debug = 1 ;
continue ;
}
/** anti-motion filtering options **/
if( strncmp(argv[nopt],"-bilinear",6) == 0 ){
DBOPT("-bilinear") ;
opt->reg_bilinear = 1 ;
continue ;
}
if( strncmp(argv[nopt],"-regbase",6) == 0 ){
DBOPT("-regbase") ;
if( ++nopt >= argc ) Syntax("-regbase needs a filename!") ;
regbase = argv[nopt] ;
continue ;
}
#ifdef ALLOW_DFTIME
if( strstr(argv[nopt],"-dftime") != NULL ){
DBOPT("-dftime") ;
opt->dfilt_code = DFILT_TIME ;
if( strstr(argv[nopt],":0") != NULL ) opt->dfilt_code = DFILT_TIME0 ;
continue ;
}
#endif
if( strstr(argv[nopt],"-dfspace") != NULL && strstr(argv[nopt],"-dfspacetime") == NULL ){
DBOPT("-dfspace") ;
opt->dfilt_code = DFILT_SPACE ;
if( strstr(argv[nopt],":0") != NULL ) opt->dfilt_code = DFILT_SPACE0 ;
continue ;
}
#ifdef ALLOW_DFTIME
if( strstr(argv[nopt],"-dfboth") != NULL ){
DBOPT("-dfboth") ;
opt->dfilt_code = DFILT_BOTH ;
if( strstr(argv[nopt],":0") != NULL ) opt->dfilt_code = DFILT_BOTH0 ;
continue ;
}
#endif
#ifdef ALLOW_DFTIME
if( strstr(argv[nopt],"-dfspacetime") != NULL ){
DBOPT("-dfspacetime") ;
opt->dfilt_code = DFILT_SPACETIME ;
if( strstr(argv[nopt],":0") != NULL ) opt->dfilt_code = DFILT_SPACETIME0 ;
continue ;
}
#endif
/** -clean **/
if( strncmp(argv[nopt],"-clean",5) == 0 ){
DBOPT("-clean") ;
opt->scale_output = FALSE;
continue ;
}
/*** -flim ==> output floating point images (default is shorts) ***/
if( strncmp(argv[nopt],"-flim",5) == 0 ){
DBOPT("-flim") ;
opt->flim = TRUE ;
continue ;
}
/*** -non ==> don't normalize alpha image ***/
if( strncmp(argv[nopt],"-non",4) == 0 ){
DBOPT("-non") ;
opt->norm_fim = FALSE ;
continue ;
}
/*** -coef # ==> coefficient for normalizing alpha image ***/
if( strncmp(argv[nopt],"-coe",4) == 0 ){
DBOPT("-coef") ;
if( ++nopt >= argc ) Syntax("-coef needs an argument") ;
ftemp = strtod( argv[nopt] , NULL ) ;
if( ftemp <= 0 ){
sprintf( err_note , "illegal -coef value: %f" , ftemp ) ;
Syntax(err_note) ;
}
opt->scale_fim = ftemp ;
continue ;
}
/*** -list # ==> threshold for report listing (not used now) ***/
if( strncmp(argv[nopt],"-list",4) == 0 ){
DBOPT("-list") ;
if( ++nopt >= argc ) Syntax("-list needs an argument") ;
ftemp = strtod( argv[nopt] , NULL ) ;
#if 0
if( ftemp <= 0 ){
sprintf( err_note , "illegal -list value: %f" , ftemp ) ;
Syntax(err_note) ;
}
#else
fprintf(stderr,"WARNING: -list option is ignored in fim2!\n") ;
#endif
opt->thresh_report = ftemp ;
continue ;
}
/*** -polref # ==> order of polynomial reference functions ***/
if( strncmp(argv[nopt],"-polref",5) == 0 || strncmp(argv[nopt],"-polort",5) == 0 ){
DBOPT("-polref") ;
if( ++nopt >= argc ) Syntax("-polref needs an argument") ;
ftemp = strtod( argv[nopt] , NULL ) ;
if( ftemp > 3 ){
fprintf( stderr , "WARNING: large -polref value %f\n" , ftemp ) ;
}
polref = (int) ftemp ;
continue ;
}
/*** -pcnt # ==> percent of deviation from perfect correlation ***/
if( strncmp(argv[nopt],"-pcnt",4) == 0 ){
DBOPT("-pcnt") ;
if( ++nopt >= argc ) Syntax("-pcnt needs an argument") ;
ftemp = strtod( argv[nopt] , NULL ) ;
if( ftemp < 0.0 || ftemp > 100.0 ){
sprintf( err_note , "illegal -pcnt value: %f" , ftemp ) ;
Syntax(err_note) ;
}
opt->thresh_pcorr = 1.0 - 0.01 * ftemp ;
continue ;
}
/*** -pcthresh # ==> actual threshold on \rho ***/
if( strncmp(argv[nopt],"-pcthresh",5) == 0 ){
DBOPT("-pcthresh") ;
if( ++nopt >= argc ) Syntax("-pcthresh needs an argument") ;
ftemp = strtod( argv[nopt] , NULL ) ;
if( ftemp < 0.0 || ftemp > 1.0 ){
sprintf( err_note , "illegal -pcthresh value: %f" , ftemp ) ;
Syntax(err_note) ;
}
opt->thresh_pcorr = ftemp ;
continue ;
}
/*** -im1 # ==> index (1...) of first image to actually use ***/
if( strncmp(argv[nopt],"-im1",4) == 0 ){
DBOPT("-im1") ;
if( ++nopt >= argc ) Syntax("-im1 needs an argument") ;
ftemp = strtod( argv[nopt] , NULL ) ;
if( ftemp < 1.0 ){
sprintf( err_note , "illegal -im1 value: %f" , ftemp ) ;
Syntax(err_note) ;
}
first_im = (int)(ftemp+0.499) ;
continue ;
}
/*** -num # ==> maximum number of images to use from command line ***/
if( strncmp(argv[nopt],"-num",4) == 0 ){
DBOPT("-num") ;
if( ++nopt >= argc ) Syntax("-num needs an argument") ;
ftemp = strtod( argv[nopt] , NULL ) ;
if( ftemp < 2 ){
sprintf( err_note , "illegal -num value: %f" , ftemp ) ;
Syntax(err_note) ;
}
num_im = (int)(ftemp+0.499) ;
continue ;
}
#define OPENERS "[{/%"
#define CLOSERS "]}/%"
#define IS_OPENER(sss) (strlen((sss))==1 && strstr(OPENERS,(sss))!=NULL)
#define IS_CLOSER(sss) (strlen((sss))==1 && strstr(CLOSERS,(sss))!=NULL)
/*** -ort file ==> reference time series ***/
if( strncmp(argv[nopt],"-ort",4) == 0 ){
DBOPT("-ort") ;
if( ++nopt >= argc ) Syntax("-ort needs a filename") ;
/** July 1995: read a bunch of orts using [ a b c ... ] format **/
if( ! IS_OPENER(argv[nopt]) ){ /* one file */
vec = RWC_read_time_series( argv[nopt] ) ;
if( vec == NULL ){ fprintf( stderr , "cannot continue\a\n" ) ; exit(1) ; }
ADDTO_TSARR( opt->refts , vec ) ;
} else { /* many */
for( nopt++ ; !IS_CLOSER(argv[nopt]) && nopt<argc ; nopt++ ){
vec = RWC_read_time_series( argv[nopt] ) ;
if( vec == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
ADDTO_TSARR( opt->refts , vec ) ;
}
if( nopt == argc ) Syntax("-ort never finishes!") ;
}
continue ;
}
/*** -ideal file ==> ideal response vector time series ***/
if( strncmp(argv[nopt],"-ideal",5) == 0 ){
DBOPT("-ideal") ;
if( ++nopt >= argc ) Syntax("-ideal needs a filename") ;
/** July 1995: read a bunch of ideals using [ a b c ... ] format **/
if( ! IS_OPENER(argv[nopt]) ){ /* one file */
ideal = RWC_read_time_series( argv[nopt] ) ;
if( ideal == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
ADDTO_TSARR( opt->idts , ideal ) ;
} else { /* many */
for( nopt++ ; !IS_CLOSER(argv[nopt]) && nopt<argc ; nopt++ ){
ideal = RWC_read_time_series( argv[nopt] ) ;
if( ideal == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
ADDTO_TSARR( opt->idts , ideal ) ;
}
if( nopt == argc ) Syntax("-ideal never finishes!") ;
}
continue ;
}
/*** -weight file ==> weight vector time series ***/
if( strncmp(argv[nopt],"-weight",5) == 0 ){
DBOPT("-weight") ;
if( ++nopt >= argc ) Syntax("-weight needs a filename") ;
if( opt->weight != NULL ){
fprintf( stderr , "cannot have two weight vectors!\a\n" ) ;
exit(1) ;
}
opt->weight = RWC_read_time_series( argv[nopt] ) ;
if( opt->weight == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
DBOPT("read in OK") ;
continue ;
}
/*** -fimfile file ==> name of output file ***/
if( strncmp(argv[nopt],"-fimfile",5) == 0 ){
DBOPT("-fimfile") ;
if( ++nopt >= argc ) Syntax("-fimfile needs a filename") ;
if( opt->fname_fim != NULL ){
fprintf( stderr , "cannot have two fim output files!\a\n" ) ;
exit(1) ;
}
opt->fname_fim = argv[nopt] ;
DBOPT("stored as fimfile") ;
continue ;
}
/*** -corr ==> write correlation file as fimfile.CORR ***/
if( strncmp(argv[nopt],"-corr",5) == 0 ){
DBOPT("-corr") ;
do_corr = TRUE ;
continue ;
}
/*** -corfile file ==> write correlation file as this name ***/
if( strncmp(argv[nopt],"-corfile",5) == 0 ||
strncmp(argv[nopt],"-corrfile",6) == 0 ){
DBOPT("-corfile") ;
if( ++nopt >= argc ) Syntax("-corfile needs a filename") ;
if( opt->fname_corr != NULL ){
fprintf( stderr , "cannot have two corr output files!\a\n" ) ;
exit(1) ;
}
opt->fname_corr = argv[nopt] ;
do_corr = TRUE ;
DBOPT("stored as corfile") ;
continue ;
}
/*** -cnrfile file ==> write cnr file as this name ***/
if( strncmp(argv[nopt],"-cnrfile",5) == 0 ){
DBOPT("-cnrfile") ;
if( ++nopt >= argc ) Syntax("-cnrfile needs a filename") ;
if( opt->fname_cnr != NULL ){
fprintf( stderr , "cannot have two cnr output files!\a\n" ) ;
exit(1) ;
}
opt->fname_cnr = argv[nopt] ;
DBOPT("stored as cnrfile") ;
continue ;
}
/*** -sigfile file ==> write sig file as this name ***/
if( strncmp(argv[nopt],"-sigfile",5) == 0 ){
DBOPT("-sigfile") ;
if( ++nopt >= argc ) Syntax("-sigfile needs a filename") ;
if( opt->fname_sig != NULL ){
fprintf( stderr , "cannot have two sig output files!\a\n" ) ;
exit(1) ;
}
opt->fname_sig = argv[nopt] ;
DBOPT("stored as sigfile") ;
continue ;
}
/*** -fitfile file ==> write fit files as this name ***/
if( strncmp(argv[nopt],"-fitfile",5) == 0 ){
DBOPT("-fitfile") ;
if( ++nopt >= argc ) Syntax("-fitfile needs a filename") ;
if( opt->fname_fit != NULL ){
fprintf( stderr , "cannot have two fit output filenames!\a\n" ) ;
exit(1) ;
}
opt->fname_fit = argv[nopt] ;
DBOPT("stored as fitfile") ;
continue ;
}
/*** -subort file ==> compute time series with orts removed ***/
if( strncmp(argv[nopt],"-subort",5) == 0 ){
DBOPT("-subort") ;
if( ++nopt >= argc ) Syntax("-subort needs a filename") ;
if( opt->fname_subort != NULL ){
fprintf( stderr , "cannot have two subort output filenames!\a\n" ) ;
exit(1) ;
}
opt->fname_subort = argv[nopt] ;
DBOPT("stored as subortfile") ;
continue ;
}
/*** get to here, and start with a '-' ==> bad news city, Arizona ***/
if( strncmp(argv[nopt],"-",1) == 0 ){
sprintf( err_note , "unknown option %s" , argv[nopt] ) ;
Syntax(err_note) ;
}
/*** finished with switches. Anything that makes it here is a filename ***/
if( opt->idts->num == 0 ){ /* if ideal timeseries not given yet, is first */
DBOPT("will be ideal") ;
ideal = RWC_read_time_series( argv[nopt] ) ;
if( ideal == NULL ){ fprintf( stderr , "cannot continue\a\n" ); exit(1) ; }
ADDTO_TSARR( opt->idts , ideal ) ;
continue ;
}
/*** this point is the start of image filenames;
break out of the loop and process them ***/
DBOPT("1st image file") ;
break ; /* time to leave */
} while( ++nopt < argc ) ; /* loop until all args are read, or we break */
/***** when this loop ends, nopt = index of first image filename in argv[] *****/
/*** check the situation for reasonabilityositiness ***/
if( opt->idts->num == 0 ) Syntax("no ideal response vector is defined!") ;
/*** compute the number of image files ***/
num_time = argc - nopt ; /* # of arguments left */
if( opt->fname_fim == NULL ){
num_time -- ; /* one less, if need be */
opt->fname_fim = argv[argc-1] ; /* last name = output file */
}
/*** July 1995: read all images in now! ***/
if( num_time < 1 ) Syntax("No input files given on command line?!") ;
opt->ims = mri_read_many_files( num_time , argv+nopt ) ;
if( opt->ims == NULL ) Syntax("Cannot read all images!") ;
num_time = MIN( opt->ims->num , num_im ) ;
if( num_time < 2 ) Syntax("must have at least 2 images to correlate!") ;
opt->ntime = num_time ;
opt->nxim = opt->ims->imarr[0]->nx ;
opt->nyim = opt->ims->imarr[0]->ny ;
#ifdef DEBUG
fprintf(stderr,"num_time = %d nx = %d ny = %d\n",num_time,opt->nxim,opt->nyim) ;
#endif
/** convert images to float format **/
for( ii=0 ; ii < num_time ; ii++ ){
if( opt->ims->imarr[ii]->kind != MRI_float ){
im = mri_to_float( opt->ims->imarr[ii] ) ;
mri_free( opt->ims->imarr[ii] ) ;
opt->ims->imarr[ii] = im ;
}
if( opt->ims->imarr[ii]->nx != opt->nxim ||
opt->ims->imarr[ii]->ny != opt->nyim ){
fprintf(stderr,"** Image size mismatch at image # %d -- end of run!\a\n",ii+1) ;
exit(1) ;
}
}
if( first_im < 1 || first_im >= num_time ){
sprintf( err_note ,
"-im1 was %d, but only have %d images" , first_im , num_time ) ;
Syntax(err_note) ;
}
/*** replace earliest images with copies of "first_im", if needed ***/
for( ii=0 ; ii < first_im-1 ; ii++ ){
im = mri_to_float( opt->ims->imarr[first_im-1] ) ; /* copy data, not just pointer */
mri_free( opt->ims->imarr[ii] ) ;
opt->ims->imarr[ii] = im ;
}
if( do_corr && opt->fname_corr == NULL ){
#ifdef DEBUG
fprintf(stderr,"creating .CORR filename\n");
#endif
ii = strlen( opt->fname_fim ) + 7 ;
opt->fname_corr = (char *) malloc( sizeof(char) * ii ) ;
if( opt->fname_corr == NULL ) MALLOC_ERR(".CORR filename") ;
strcpy( opt->fname_corr , opt->fname_fim ) ;
strcat( opt->fname_corr , ".CORR" ) ;
}
/*** put the polynomial responses into opt->refts ***/
if( polref >= 0 ){
#ifdef DEBUG
fprintf(stderr,"creating polynomial references; polref=%d\n",polref);
#endif
for( pp=0 ; pp <= polref ; pp++ ){
vec = RWC_blank_time_series( num_time ) ;
#if 1
if( pp == 0 ){
for( ii=0 ; ii < num_time ; ii++ ) vec->ts[ii] = 1.0 ;
} else {
fscal = 1.0 / num_time ;
for( ii=0 ; ii < num_time ; ii++ ) vec->ts[ii] = pow(fscal*ii,pp) ;
}
#else
ftemp = 0.5 * num_time - 0.4321 ; /* 0.4321 ensures (ii-ftemp)!=0 */
fscal = 2.0 / num_time ; /* range of arg to pow is -1..1 */
for( ii=0 ; ii < num_time ; ii++ )
vec->ts[ii] = pow( fscal*(ii-ftemp) , pp ) ;
#endif
ADDTO_TSARR( opt->refts , vec ) ;
}
}
/*** if no weight vector input (via -weight filename), make one up ***/
if( opt->weight == NULL ){
#ifdef DEBUG
fprintf(stderr,"creating default weight\n");
#endif
vec = RWC_blank_time_series( num_time ) ;
for( ii=0 ; ii < num_time ; ii++ ) vec->ts[ii] = 1.0 ; /* weight = all ones */
opt->weight = vec ;
}
for( ii=0 ; ii < first_im-1 ; ii++ ) opt->weight->ts[ii] = 0.0 ; /* except skipped images */
/*** make space for the ideal vector in refts (but don't but it there yet) ***/
ADDTO_TSARR( opt->refts , NULL ) ;
/*** check all time series for being long enough ***/
bad = FALSE ;
for( ii=0 ; ii < opt->refts->num-1 ; ii++ ){ /* note not checking last one */
#ifdef DEBUG
fprintf(stderr,"checking ref %d for size\n",ii) ;
#endif
if( opt->refts->tsarr[ii]->len < num_time ){
fprintf( stderr ,
"reference vector %s has %d elements: too short!\a\n" ,
opt->refts->tsarr[ii]->fname , opt->refts->tsarr[ii]->len ) ;
bad = TRUE ;
}
}
for( ii=0 ; ii < opt->idts->num-1 ; ii++ ){
#ifdef DEBUG
fprintf(stderr,"checking ideal %d for size\n",ii) ;
#endif
if( opt->idts->tsarr[ii]->len < num_time ){
fprintf( stderr ,
"ideal vector %s has %d elements: too short!\a\n" ,
opt->idts->tsarr[ii]->fname , opt->idts->tsarr[ii]->len ) ;
bad = TRUE ;
}
}
if( opt->weight->len < num_time ){
fprintf( stderr ,
"weight vector %s has %d elements: too short!\a\n" ,
opt->weight->fname , opt->weight->len ) ;
bad = TRUE ;
}
if( bad ) exit(1) ;
/*** zero out weight at any time where a refts time_series is blasted ***/
#ifdef DEBUG
fprintf(stderr,"blasting away ... ") ;
#endif
for( ii=0 ; ii < num_time ; ii++ ){
bad = (opt->weight->ts[ii] >= BLAST) || (opt->weight->ts[ii] < 0.0) ;
for( rr=0 ; !bad && rr < opt->refts->num-1 ; rr++ )
bad = (opt->refts->tsarr[rr]->ts[ii] >= BLAST) ;
if( bad ){
opt->weight->ts[ii] = 0 ;
for( rr=0 ; rr < opt->refts->num-1 ; rr++ ) /* note not checking last one */
opt->refts->tsarr[rr]->ts[ii] = 0 ;
}
}
#ifdef DEBUG
fprintf(stderr,"\n") ;
#endif
/*** check to see if the edited vectors are still nonzero enough ***/
bad = FALSE ;
#ifdef DEBUG
fprintf(stderr,"checking weight for nonnegativity\n") ;
#endif
ftemp = RWC_norm_ts( num_time , opt->weight ) ;
if( ftemp < 1.e-9 ){
fprintf( stderr , "there is no time with nonzero weighting!\n" ) ;
bad = TRUE ;
}
for( rr=0 ; rr < opt->refts->num-1 ; rr++ ){ /* note not checking last one */
#ifdef DEBUG
fprintf(stderr,"checking ref %d for nonzeroness\n",rr) ;
#endif
ftemp = RWC_norm_ts( num_time , opt->refts->tsarr[rr] ) ;
if( ftemp < 1.e-9 ){
fprintf( stderr , "reference vector %d is all zeroes\n" , rr+1 ) ;
bad = TRUE ;
}
}
if( bad ) exit(1) ;
/*** June 1995: get first image with nonzero weighting ***/
if( regbase == NULL ){
for( ii=0 ; ii < num_time ; ii++ )
if( opt->weight->ts[ii] != 0.0 && ideal->ts[ii] < BLAST ) break ;
if( ii == num_time ){ fprintf(stderr,"FIRST_IM: scan error!\a\n");exit(1); }
opt->first_flim = mri_to_float( opt->ims->imarr[ii] ) ; /* copy it */
} else {
MRI_IMAGE * qim ;
qim = mri_read_just_one( regbase ) ;
if( qim == NULL ) Syntax("Couldn't read -regbase image!") ;
if( qim->kind == MRI_float ) opt->first_flim = qim ;
else {
opt->first_flim = mri_to_float( qim ) ;
mri_free( qim ) ;
}
if( opt->first_flim->nx != opt->nxim || opt->first_flim->ny != opt->nyim ){
fprintf(stderr,"-regbase: image size mismatch!\a\n") ; exit(1) ;
}
}
/*** Setup for differential filtering (registration) ***/
if( opt->dfilt_code != DFILT_NONE ){
int alcode ;
MRI_IMARR * tims ;
time_series * dxts , * dyts , * phits ;
switch( opt->dfilt_code ){
case DFILT_TIME: alcode = 0 ; break ;
case DFILT_TIME0: alcode = ALIGN_NOITER_CODE ; break ;
case DFILT_SPACETIME:
case DFILT_BOTH:
case DFILT_SPACE: alcode = ALIGN_REGISTER_CODE ; break ;
case DFILT_SPACETIME0:
case DFILT_BOTH0:
case DFILT_SPACE0: alcode = ALIGN_REGISTER_CODE | ALIGN_NOITER_CODE ; break ;
default:
Syntax("Internal error: opt->dfilt_code illegal!") ;
}
dxts = RWC_blank_time_series( num_time ) ; /* space for motion params */
dyts = RWC_blank_time_series( num_time ) ;
phits = RWC_blank_time_series( num_time ) ;
if( opt->reg_bilinear ) alcode |= ALIGN_BILINEAR_CODE ;
if( opt->debug ) alcode |= ALIGN_DEBUG_CODE ;
tims = mri_align_dfspace( opt->first_flim , NULL , opt->ims ,
alcode , dxts->ts , dyts->ts , phits->ts ) ;
switch( opt->dfilt_code ){
case DFILT_TIME:
case DFILT_TIME0:
ADDTO_TSARR( opt->refts , NULL ) ; /* make 3 blank spots at end */
ADDTO_TSARR( opt->refts , NULL ) ;
ADDTO_TSARR( opt->refts , NULL ) ;
/* move previously existing time series up by 3 */
for( ii=opt->refts->num-4 ; ii >=0 ; ii-- )
opt->refts->tsarr[ii+3] = opt->refts->tsarr[ii] ;
/* put dx,dy,phi at beginning of list */
opt->refts->tsarr[0] = dxts ;
opt->refts->tsarr[1] = dyts ;
opt->refts->tsarr[2] = phits ;
break ;
case DFILT_SPACE:
case DFILT_SPACE0:
DESTROY_IMARR( opt->ims ) ; /* put registered images in */
opt->ims = tims ; /* place of the input images */
#if 0
if( opt->debug ){
char fff[64] ;
for( ii=0 ; ii < IMARR_COUNT(opt->ims) ; ii++ ){
sprintf(fff,"rrr.%04d",ii+1) ;
mri_write( fff , IMARR_SUBIMAGE(opt->ims,ii) ) ;
}
}
#endif
RWC_free_time_series( dxts ) ;
RWC_free_time_series( dyts ) ;
RWC_free_time_series( phits ) ;
break ;
case DFILT_BOTH:
case DFILT_BOTH0:
ADDTO_TSARR( opt->refts , NULL ) ;
ADDTO_TSARR( opt->refts , NULL ) ;
ADDTO_TSARR( opt->refts , NULL ) ;
for( ii=opt->refts->num-4 ; ii >=0 ; ii-- )
opt->refts->tsarr[ii+3] = opt->refts->tsarr[ii] ;
opt->refts->tsarr[0] = dxts ;
opt->refts->tsarr[1] = dyts ;
opt->refts->tsarr[2] = phits ;
opt->rims = tims ; /* keep registered images as a separate data set */
break ;
case DFILT_SPACETIME:
case DFILT_SPACETIME0:
ADDTO_TSARR( opt->refts , NULL ) ;
ADDTO_TSARR( opt->refts , NULL ) ;
ADDTO_TSARR( opt->refts , NULL ) ;
for( ii=opt->refts->num-4 ; ii >=0 ; ii-- )
opt->refts->tsarr[ii+3] = opt->refts->tsarr[ii] ;
opt->refts->tsarr[0] = dxts ;
opt->refts->tsarr[1] = dyts ;
opt->refts->tsarr[2] = phits ;
DESTROY_IMARR( opt->ims ) ; /* put registered images in */
opt->ims = tims ; /* place of the input images */
break ;
}
}
/*** Done! ***/
return ;
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*** give some help and exit ***/
void Syntax( char *note )
{
static char *mesg[] = {
"Usage: fim2 [options] image_files ..." ,
"where 'image_files ...' is a sequence of MRI filenames," ,
" " ,
"options are:",
"-pcnt # correlation coeff. threshold will be 1 - 0.01 * #",
"-pcthresh # correlation coeff. threshold will be #",
"-im1 # index of image file to use as first in time series;",
" default is 1; previous images are filled with this",
" image to synchronize with the reference time series",
"-num # number of images to actually use, if more than this",
" many are specified on the command line; default is",
" to use all images",
"-non this option turns off the default normalization of",
" the output activation image; the user should provide",
" a scaling factor via '-coef #', or '1' will be used",
"-coef # the scaling factor used to convert the activation output",
" from floats to short ints (if -non is also present)",
" ",
"-ort fname fname = filename of a time series to which the image data",
" will be orthogonalized before correlations are computed;",
" any number of -ort options (from 0 on up) may be used",
"-ideal fname fname = filename of a time series to which the image data",
" is to be correlated; exactly one such time series is",
" required; if the -ideal option is not used, then the",
" first filename after all the options will be used",
" N.B.: This version of fim2 allows the specification of more than",
" one ideal time series file. Each one is separately correlated",
" with the image time series and the one most highly correlated",
" is selected for each pixel. Multiple ideals are specified",
" using more than one '-ideal fname' option, or by using the",
" form '-ideal [ fname1 fname2 ... ]' -- this latter method",
" allows the use of wildcarded ideal filenames.",
" The '[' character that indicates the start of a group of",
" ideals can actually be any ONE of these: " OPENERS ,
" and the ']' that ends the group can be: " CLOSERS ,
" ",
" [Format of ort and ideal time series files:",
" ASCII; one number per line;",
" Same number of lines as images in the time series;",
" Value over 33333 --> don't use this image in the analysis]",
" ",
"-polref # use polynomials of order 0..# as extra 'orts';",
"[or -polort #] default is 0 (yielding a constant vector).",
" Use # = -1 to suppress this feature.",
#if 0
"-weight fname fname = filename of a times series to be used as weights",
" in the correlation calculation; all time series",
" (orts, ideal, and image) will be weighted at time i",
" by the weight at that time; if not given, defaults to",
" 1 at all times (but any ort or ideal >= 33333 gives 0)",
#endif
" ",
"-fimfile fname fname = filename to save activation magnitudes in;",
" if not given, the last name on the command line will",
" be used",
"-corr if present, indicates to write correlation output to",
" image file 'fimfile.CORR' (next option is better)",
"-corfile fname fname = filename to save correlation image in;",
" if not present, and -corr is not present, correlation",
" image is not saved.",
"-cnrfile fname fname = filename to save contrast-to-noise image in;" ,
" if not present, will not be computed or saved;" ,
" CNR is scaled by 100 if images are output as shorts" ,
" and is written 'as-is' if output as floats (see -flim)." ,
" [CNR is defined here to be alpha/sigma, where" ,
" alpha = amplitude of normalized ideal in a pixel" ,
" sigma = standard deviation of pixel after removal" ,
" of orts and ideal" ,
" normalized ideal = ideal scaled so that trough-to-peak" ,
" height is one.]" ,
"-sigfile fname fname = filename to save standard deviation image in;" ,
" the standard deviation is of what is left after the" ,
" least squares removal of the -orts, -polrefs, and -ideal." ,
" N.B.: This is always output in the -flim format!" ,
"-fitfile fname Image files of the least squares fit coefficients of" ,
" all the -ort and -polref time series that" ,
" are projected out of the data time series before" ,
" the -ideal is fit. The actual filenames will" ,
" be fname.01 fname.02 ...." ,
" Their order is -orts, then -polrefs, and last -ideal." ,
" N.B.: These are always output in the -flim format!" ,
"-subort fname A new timeseries of images is written to disk, with",
" names of the form 'fname.0001', etc. These images",
" have the orts and polrefs (but not ideals) subtracted out.",
" N.B.: These are always output in the -flim format!" ,
"-flim If present, write outputs in mrilib 'float' format,",
" rather than scale+convert to integers",
" [The 'ftosh' program can later convert to short integers]",
"-clean if present, then output images won't have the +/- 10000",
" values forced into their corners for scaling purposes.",
"-clip if present, output correlations, etc., will be set to",
" zero in regions of low intensity.",
"-q if present, indicates 'quiet' operation.",
"-dfspace[:0] Indicates to use the 'dfspace' filter (a la imreg) to",
" register the images spatially before filtering." ,
#ifdef ALLOW_DFTIME
"-dftime[:0] Indicates to use the 'dftime' filter (a la imreg) to",
" produce 3 additional orts, in an attempt to reduce",
" motion artifacts.",
"-dfboth[:0] Indicates to use both -dftime and -dfspace (separately)," ,
" then take as the resulting correlation the smaller of the",
" two results in each pixel (the conservative approach: " ,
" to be above -pcthresh, both calculations must 'hit')." ,
" The resulting fim is the one corresponding to the chosen" ,
" correlation in each pixel." ,
"-dfspacetime:[0] Indicates to use -dfspace and then -dftime together",
" (not separately) on the time series of images." ,
#endif
"-regbase fname Indicates to read image in file 'fname' as the base",
" image for registration. If not given, the first image",
" in the time series that is used in the correlation",
" computations will be used. This is also the image",
" that is used to define 'low intensity' for the -clip option.",
" "
} ;
int nsize , ii ;
if( note != NULL && (nsize=strlen(note)) > 0 ){
fprintf(stderr,"\n") ;
for(ii=0;ii<nsize;ii++) fprintf(stderr,"*") ;
fprintf(stderr,"\a\n%s\n", note ) ;
for(ii=0;ii<nsize;ii++) fprintf(stderr,"*") ;
fprintf(stderr,"\ntry fim2 -help for instructions\a\n") ;
exit(-1) ;
}
for( ii=0 ; ii < sizeof(mesg)/sizeof(char *) ; ii++ ){
printf( " %s\n" , mesg[ii] );
}
exit(0) ;
}
/*----------------------------------------------------------------------------*/
time_series * edit_weight( time_series * wt , time_series * ideal )
{
time_series * wtnew ;
int ii , ntime ;
float ftemp ;
ntime = MIN(wt->len,ideal->len) ;
wtnew = RWC_blank_time_series( ntime ) ;
for( ii=0 ; ii < ntime ; ii++ )
wtnew->ts[ii] = (ideal->ts[ii] >= BLAST) ? (0.0) : (wt->ts[ii]) ;
ftemp = RWC_norm_ts( ntime , wtnew ) ;
if( ftemp < 1.e-9 ){
fprintf(stderr,"** bad ideal: no times with nonzero weight!\n") ;
RWC_free_time_series( wtnew ) ;
return NULL ;
}
return wtnew ;
}
syntax highlighted by Code2HTML, v. 0.9.1