Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
model_convgamma.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "NLfit_model.h"
00008
00009 static int refnum = 0 ;
00010 static int refnz = 0 ;
00011 static float * refts = NULL ;
00012 static int * refin = NULL ;
00013
00014 void gamma_model( float * , int , float ** , float * ) ;
00015
00016 #define ERREX(str) ( fprintf(stderr,"\n*** %s\a\n",str) , exit(1) )
00017
00018
00019
00020
00021
00022
00023 void conv_set_ref( int num , float * ref )
00024 {
00025 if( num > 0 && ref != NULL ){
00026 int ii ;
00027
00028
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++ )
00037 if( refts[ii] != 0 ) refin[refnz++] = ii ;
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 {
00048
00049 char * cp ;
00050 MRI_IMAGE * flim ;
00051 float one = 1.0 ;
00052
00053 cp = my_getenv("AFNI_CONVMODEL_REF") ;
00054 if( cp == NULL )
00055 ERREX("model_convgamma: Can't read AFNI_CONVMODEL_REF from environment") ;
00056
00057 flim = mri_read_1D(cp) ;
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) ) ;
00069 mri_free(flim) ;
00070 }
00071 return ;
00072 }
00073
00074
00075
00076
00077
00078 void conv_model( float * gs , int ts_length ,
00079 float ** x_array , float * ts_array )
00080 {
00081 int ii, jj,jbot,jtop , kk , nid_top,nid_bot ;
00082 float top , val ;
00083
00084 static int nid = 0 ;
00085 static float * fid = NULL ;
00086
00087
00088
00089 if( refnum <= 0 ) conv_set_ref( 0 , NULL ) ;
00090
00091
00092
00093 for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ;
00094
00095
00096
00097 if( nid < ts_length ){
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 ) ;
00104
00105 top = 0.0 ;
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 ;
00110 top *= 0.001 ;
00111 for( jj=0 ; jj < ts_length ; jj++ ){
00112 if( fabs(fid[jj]) < top ) fid[jj] = 0.0 ;
00113 }
00114 for( nid_bot=0 ; nid_bot < ts_length ; nid_bot++ )
00115 if( fid[nid_bot] != 0.0 ) break ;
00116 for( nid_top=ts_length-1 ; nid_top > nid_bot ; nid_top-- )
00117 if( fid[nid_top] != 0.0 ) break ;
00118
00119
00120
00121 for( ii=0 ; ii < refnz ; ii++ ){
00122 kk = refin[ii] ; if( kk >= ts_length ) break ;
00123 val = refts[kk] ;
00124
00125
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 }
00134
00135
00136
00137 DEFINE_MODEL_PROTOTYPE
00138
00139 MODEL_interface * initialize_model ()
00140 {
00141 MODEL_interface * mi = NULL;
00142
00143
00144
00145 mi = (MODEL_interface *) XtMalloc (sizeof(MODEL_interface));
00146
00147
00148
00149 strcpy (mi->label, "ConvGamma");
00150
00151
00152
00153 mi->model_type = MODEL_SIGNAL_TYPE;
00154
00155
00156
00157 mi->params = 4;
00158
00159
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
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
00174 mi->call_func = &conv_model;
00175
00176 return (mi);
00177 }
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194 void gamma_model
00195 (
00196 float * gs,
00197 int ts_length,
00198 float ** x_array,
00199 float * ts_array
00200 )
00201
00202 {
00203 int it;
00204 float t;
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
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 }