00118 {
00119    int ii, jj,jbot,jtop , kk , nid_top,nid_bot , jv ;
00120    float top , val , amp ;
00121 
00122    static int     nid = 0 ;     /* number of pts in fid array */
00123    static float * fid = NULL ;  /* impulse response function */
00124 
00125    /*** make sure there is a reference function to convolve with ***/
00126 
00127    if( refnum[0] <= 0 ) conv_set_ref( 0 , NULL ) ;
00128 
00129    /*** initialize the output ***/
00130 
00131    for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ;
00132 
00133    /*** initialize the impulse response ***/
00134 
00135    if( nid < ts_length ){              /* make some space for it */
00136       if( fid != NULL ) free(fid) ;
00137       nid = ts_length ;
00138       fid = (float *) malloc( sizeof(float) * nid ) ;
00139    }
00140 
00141    gamma_model( gs , ts_length , x_array , fid ) ;  /* compute impulse */
00142 
00143    /*** discard small values from the fid ***/
00144 
00145 #define EPS 0.001  /* definition of small */
00146 
00147 #undef FIND_TOP  /* don't need this, since our fid amplitude is set to 1.0 */
00148 #ifdef FIND_TOP
00149    top = 0.0 ;                                      /* find max value */
00150    for( jj=0 ; jj < ts_length ; jj++ ){
00151       val = fabs(fid[jj]) ; if( val > top ) top = val ;
00152    }
00153    if( top == 0.0 ) fid[0] = 1.0 ;                  /* very unlikely case */
00154    top *= EPS ;
00155 #else
00156    top = EPS ;
00157 #endif
00158    for( jj=0 ; jj < ts_length ; jj++ ){             /* discard small values */
00159       if( fabs(fid[jj]) < top ) fid[jj] = 0.0 ;
00160    }
00161    for( nid_bot=0 ; nid_bot < ts_length ; nid_bot++ )         /* find first nonzero */
00162       if( fid[nid_bot] != 0.0 ) break ;
00163    for( nid_top=ts_length-1 ; nid_top > nid_bot ; nid_top-- ) /* and last nonzero */
00164       if( fid[nid_top] != 0.0 ) break ;
00165 
00166    /*** loop over each reference ***/
00167 
00168    for( jv=0 ; jv < NREF ; jv++ ){
00169 
00170       amp = gs[jv+3] ; if( amp == 0.0 ) continue ;
00171 
00172       /*** loop over each nonzero point in the reference ***/
00173 
00174       for( ii=0 ; ii < refnz[jv] ; ii++ ){
00175          kk  = refin[jv][ii] ; if( kk >= ts_length ) break ;
00176          val = amp * refts[jv][kk] ;
00177 
00178          /*** for each point in the impulse ***/
00179 
00180          jtop = ts_length - kk ; if( jtop > nid_top ) jtop = nid_top+1 ;
00181          for( jj=nid_bot ; jj < jtop ; jj++ )
00182             ts_array[kk+jj] += val * fid[jj] ;
00183       }
00184    }
00185 
00186    return ;
00187 }

00043 {
00044    int jv ;
00045 
00046    if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/
00047       int ii ;
00048 
00049       for( jv=0 ; jv < NREF ; jv++ ){
00050 
00051          /* get rid of old data? */
00052 
00053          if( refts[jv] != NULL ){
00054             free(refts[jv]); refts[jv] = NULL;
00055             free(refin[jv]); refin[jv] = NULL;
00056          }
00057 
00058          /* copy new data */
00059 
00060          refnum[jv] = num ;
00061          refts[jv]  = (float *) malloc( sizeof(float) * num ) ;
00062          refin[jv]  = (int *)   malloc( sizeof(int)   * num ) ;
00063          memcpy( refts[jv] , ref[jv] , sizeof(float) * num ) ;
00064 
00065          /* build a list of nonzero entries in this column */
00066 
00067          for( ii=0,refnz[jv]=0 ; ii < num ; ii++ )
00068             if( refts[jv][ii] != 0.0 ) refin[jv][refnz[jv]++] = ii ;
00069 
00070          if( refnz[jv] == 0 )
00071             ERREX(__FILE__ ": All zero reference timeseries column!") ;
00072       }
00073       return ;
00074 
00075    } else { /*** if no inputs, do something special ***/
00076 
00077      char * cp ;
00078      MRI_IMAGE * flim ;
00079      int jv , nx ;
00080      float * ref[NREF] ;
00081 
00082      cp = my_getenv("AFNI_CONVMODEL_REF") ;  /* get name of reference file */
00083      if( cp == NULL )
00084         ERREX(__FILE__ ": Can't read AFNI_CONVMODEL_REF from environment") ;
00085 
00086      flim = mri_read_1D(cp) ;            /* 16 Nov 1999: replaces mri_read_ascii */
00087      if( flim == NULL ){
00088         char buf[256] ;
00089         sprintf(buf,__FILE__ ": Can't read timeseries file %s",cp) ;
00090         ERREX(buf) ;
00091      } else {
00092         fprintf(stderr,__FILE__ ": Read reference file %s\n",cp) ;
00093      }
00094 
00095      if( flim->ny < NREF )
00096         ERREX(__FILE__ ": reference file has too few columns!") ;
00097      else if( flim->nv > NREF )
00098         fprintf(stderr,__FILE__ " WARNING: reference file has too many columns!\n") ;
00099 
00100      nx = flim->nx ;
00101      for( jv=0 ; jv < NREF ; jv++ )
00102         ref[jv] = MRI_FLOAT_PTR(flim) + jv*nx ;
00103 
00104      conv_set_ref( nx , ref ) ;  /* recursion! */
00105      mri_free(flim) ;
00106    }
00107    return ;
00108 }

00265 {
00266   int it;                           /* time index */
00267   float t;                          /* time */
00268   float gsi , fac ;
00269 
00270   if( gs[2] <= 0.0 || gs[1] <= 0.0 ){
00271      ts_array[0] = 1.0 ;
00272      for( it=1 ; it < ts_length ; it++ ) ts_array[it] = 0.0 ;
00273      return ;
00274   }
00275 
00276   /* fac is chosen to make the peak value equal to 1 (at t = gs[1]*gs[2]) */
00277 
00278   gsi = 1.0 / gs[2] ;
00279   fac = exp( gs[1] * ( 1.0 - log(gs[1]*gs[2]) ) ) ;
00280   for( it=0;  it < ts_length;  it++){
00281      t = x_array[it][1] - gs[0] ;
00282      ts_array[it] = (t <= 0.0) ? 0.0
00283                                : fac * exp( log(t) * gs[1] - t * gsi ) ;
00284   }
00285   return ;
00286 }

00194 {
00195   MODEL_interface * mi = NULL;
00196   int jv ;
00197   char buf[32] ;
00198 
00199   /*----- allocate memory space for model interface -----*/
00200 
00201   mi = (MODEL_interface *) XtMalloc (sizeof(MODEL_interface));
00202 
00203   /*----- name of this model -----*/
00204 
00205   sprintf( mi->label , "ConvGamma%da" , NREF ) ;
00206 
00207   /*----- this is a signal model -----*/
00208 
00209   mi->model_type = MODEL_SIGNAL_TYPE;
00210 
00211   /*----- number of parameters in the model -----*/
00212 
00213   mi->params = 3 + NREF ;
00214 
00215   /*----- parameter labels -----*/
00216 
00217   strcpy (mi->plabel[0], "t0");  /* impulse response is proportional to */
00218   strcpy (mi->plabel[1], "r");   /*       r  -(t-t0)/b  */
00219   strcpy (mi->plabel[2], "b");   /* (t-t0)  e           */
00220 
00221   for( jv=0 ; jv < NREF ; jv++ )
00222      sprintf( mi->plabel[jv+3] , "amp_%d" , jv+1 ) ;
00223 
00224   /*----- minimum and maximum parameter constraints -----*/
00225 
00226   mi->min_constr[0] = 0.0;  mi->max_constr[0] = 10.0;  /* delay */
00227   mi->min_constr[1] = 1.0;  mi->max_constr[1] = 19.0;  /* power */
00228   mi->min_constr[2] = 0.1;  mi->max_constr[2] =  5.0;  /* time scale */
00229 
00230   for( jv=0 ; jv < NREF ; jv++ ){
00231      mi->min_constr[jv+3] = 0.0;  mi->max_constr[jv+3] = 200.0;
00232   }
00233 
00234   /*----- function which implements the model -----*/
00235   mi->call_func = &conv_model;
00236 
00237   return (mi);
00238 }