00080 {
00081    int ii, jj,jbot,jtop , kk , nid_top,nid_bot ;
00082    float top , val ;
00083 
00084    static int     nid = 0 ;     /* number of pts in impulse */
00085    static float * fid = NULL ;  /* impulse response function */
00086 
00087    /*** make sure there is a reference function to convolve with ***/
00088 
00089    if( refnum <= 0 ) conv_set_ref( 0 , NULL ) ;
00090 
00091    /*** initialize the output ***/
00092 
00093    for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ;
00094 
00095    /*** initialize the impulse response ***/
00096 
00097    if( nid < ts_length ){              /* make some space for it */
00098       if( fid != NULL ) free(fid) ;
00099       nid = ts_length ;
00100       fid = (float *) malloc( sizeof(float) * nid ) ;
00101    }
00102 
00103    gamma_model( gs , ts_length , x_array , fid ) ;  /* compute impulse */
00104 
00105    top = 0.0 ;                                      /* find max value */
00106    for( jj=0 ; jj < ts_length ; jj++ ){
00107       val = fabs(fid[jj]) ; if( val > top ) top = val ;
00108    }
00109    if( top == 0.0 ) fid[0] = 1.0 ;                  /* very unlikely case */
00110    top *= 0.001 ;
00111    for( jj=0 ; jj < ts_length ; jj++ ){             /* discard small values */
00112       if( fabs(fid[jj]) < top ) fid[jj] = 0.0 ;
00113    }
00114    for( nid_bot=0 ; nid_bot < ts_length ; nid_bot++ )         /* find first nonzero */
00115       if( fid[nid_bot] != 0.0 ) break ;
00116    for( nid_top=ts_length-1 ; nid_top > nid_bot ; nid_top-- ) /* and last nonzero */
00117       if( fid[nid_top] != 0.0 ) break ;
00118 
00119    /*** loop over each nonzero point in the reference ***/
00120 
00121    for( ii=0 ; ii < refnz ; ii++ ){
00122       kk  = refin[ii] ; if( kk >= ts_length ) break ;
00123       val = refts[kk] ;
00124 
00125       /*** for each point in the impulse ***/
00126 
00127       jtop = ts_length - kk ; if( jtop > nid_top ) jtop = nid_top+1 ;
00128       for( jj=nid_bot ; jj < jtop ; jj++ )
00129          ts_array[kk+jj] += val * fid[jj] ;
00130    }
00131 
00132    return ;
00133 }

00024 {
00025    if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/
00026       int ii ;
00027 
00028       /* get rid of old data */
00029 
00030       if( refts != NULL ){ free(refts); refts = NULL; free(refin); refin = NULL; }
00031 
00032       refnum = num ;
00033       refts  = (float *) malloc( sizeof(float) * num ) ;
00034       refin  = (int *)   malloc( sizeof(int)   * num ) ;
00035       memcpy( refts , ref , sizeof(float) * num ) ;
00036       for( ii=0,refnz=0 ; ii < num ; ii++ )              /* build list of nonzero */
00037          if( refts[ii] != 0 ) refin[refnz++] = ii ;      /* points in refts */
00038       if( refnz == 0 )
00039          ERREX("model_convgamma: All zero reference timeseries!") ;
00040 
00041 #if 0
00042 fprintf(stderr,"conv_set_ref: num=%d nonzero=%d\n",num,refnz) ;
00043 #endif
00044 
00045       return ;
00046 
00047    } else { /*** if no inputs, do something special ***/
00048 
00049      char * cp ;
00050      MRI_IMAGE * flim ;
00051      float one = 1.0 ;
00052 
00053      cp = my_getenv("AFNI_CONVMODEL_REF") ;  /* get name of reference file */
00054      if( cp == NULL )
00055         ERREX("model_convgamma: Can't read AFNI_CONVMODEL_REF from environment") ;
00056 
00057      flim = mri_read_1D(cp) ;            /* 16 Nov 1999: replaces mri_read_ascii */
00058      if( flim == NULL ){
00059         char buf[256] ;
00060         sprintf(buf,"model_convgamma: Can't read timeseries file %s",cp) ;
00061         ERREX(buf) ;
00062      }
00063 
00064 #if 0
00065 fprintf(stderr,"conv_set_ref: refts=%s  nx=%d\n",cp,flim->ny) ;
00066 #endif
00067 
00068      conv_set_ref( flim->nx , MRI_FLOAT_PTR(flim) ) ;  /* recursion! */
00069      mri_free(flim) ;
00070    }
00071    return ;
00072 }

00202 {
00203   int it;                           /* time index */
00204   float t;                          /* time */
00205   float gsi , fac ;
00206 
00207   if( gs[3] <= 0.0 || gs[2] <= 0.0 ){
00208      for( it=0 ; it < ts_length ; it++ ) ts_array[it] = 0.0 ;
00209      return ;
00210   }
00211 
00212   /* fac is chosen to make the peak value equal to gs[1] */
00213 
00214   gsi = 1.0 / gs[3] ;
00215   fac = gs[1] * exp( gs[2] * ( 1.0 - log(gs[2]*gs[3]) ) ) ;
00216   for( it=0;  it < ts_length;  it++){
00217      t = x_array[it][1] - gs[0] ;
00218      ts_array[it] = (t <= 0.0) ? 0.0
00219                                : fac * exp( log(t) * gs[2] - t * gsi ) ;
00220   }
00221   return ;
00222 }

00140 {
00141   MODEL_interface * mi = NULL;
00142 
00143   /*----- allocate memory space for model interface -----*/
00144 
00145   mi = (MODEL_interface *) XtMalloc (sizeof(MODEL_interface));
00146 
00147   /*----- name of this model -----*/
00148 
00149   strcpy (mi->label, "ConvGamma");
00150 
00151   /*----- this is a signal model -----*/
00152 
00153   mi->model_type = MODEL_SIGNAL_TYPE;
00154 
00155   /*----- number of parameters in the model -----*/
00156 
00157   mi->params = 4;
00158 
00159   /*----- parameter labels -----*/
00160 
00161   strcpy (mi->plabel[0], "t0");
00162   strcpy (mi->plabel[1], "amp");
00163   strcpy (mi->plabel[2], "r");
00164   strcpy (mi->plabel[3], "b");
00165 
00166   /*----- minimum and maximum parameter constraints -----*/
00167 
00168   mi->min_constr[0] =     0.0;    mi->max_constr[0] =    10.0;
00169   mi->min_constr[1] =     0.0;    mi->max_constr[1] =   200.0;
00170   mi->min_constr[2] =     1.0;    mi->max_constr[2] =    19.0;
00171   mi->min_constr[3] =     0.1;    mi->max_constr[3] =     5.0;
00172 
00173   /*----- function which implements the model -----*/
00174   mi->call_func = &conv_model;
00175 
00176   return (mi);
00177 }