1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*----------------------------------------------------------------------
8 * Fit parameters t0, k, r and b in the equation:
9 *
10 * f(t) = k * (t-t0)^r * e^(-(t-t0)/b)
11 *
12 * t0 : time offset before response begins
13 * k : amplitude
14 * r : rise parameter (exponent on polynomial component)
15 * b : fall parameter (exponent on decay component)
16 *
17 * The above reference function is then convolved with the time series
18 * specified by AFNI_CONVMODEL_REF, which might just be a TR-locked binary
19 * time series specifying onset events. In such a case, the returned time
20 * series would be the sum of f(t), starting at each of the "events" in
21 * AFNI_CONVMODEL_REF.
22 *
23 * Convolving with AFNI_CONVMODEL_REF could also handle events of different
24 * durations (spanning multiple consecutive onset time points), as well as
25 * events with varying (but known) relative magnitudes, akin to ampiltude
26 * modulation.
27 *----------------------------------------------------------------------*/
28
29 #include "NLfit_model.h"
30
31 static int refnum = 0 ; /* # pts in refts */
32 static int refnz = 0 ; /* # of nonzero pts */
33 static float * refts = NULL ; /* reference time series */
34 static int * refin = NULL ; /* indexes of nonzero pts */
35
36 void gamma_model( float * , int , float ** , float * ) ;
37
38 #define ERREX(str) ( fprintf(stderr,"\n*** %s\a\n",str) , exit(1) )
39
40 /*----------------------------------------------------------------------
41 Function to set the reference time series, with which the
42 model function is convolved to produce the simulated data.
43 ------------------------------------------------------------------------*/
44
conv_set_ref(int num,float * ref)45 void conv_set_ref( int num , float * ref )
46 {
47 if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/
48 int ii ;
49
50 /* get rid of old data */
51
52 if( refts != NULL ){ free(refts); refts = NULL; free(refin); refin = NULL; }
53
54 refnum = num ;
55 refts = (float *) malloc( sizeof(float) * num ) ;
56 refin = (int *) malloc( sizeof(int) * num ) ;
57 memcpy( refts , ref , sizeof(float) * num ) ;
58 for( ii=0,refnz=0 ; ii < num ; ii++ ) /* build list of nonzero */
59 if( refts[ii] != 0 ) refin[refnz++] = ii ; /* points in refts */
60 if( refnz == 0 )
61 ERREX("model_convgamma: All zero reference timeseries!") ;
62
63 #if 0
64 fprintf(stderr,"conv_set_ref: num=%d nonzero=%d\n",num,refnz) ;
65 #endif
66
67 return ;
68
69 } else { /*** if no inputs, do something special ***/
70
71 char * cp ;
72 MRI_IMAGE * flim ;
73 float one = 1.0 ;
74
75 cp = my_getenv("AFNI_CONVMODEL_REF") ; /* get name of reference file */
76 if( cp == NULL )
77 ERREX("model_convgamma: Can't read AFNI_CONVMODEL_REF from environment") ;
78
79 flim = mri_read_1D(cp) ; /* 16 Nov 1999: replaces mri_read_ascii */
80 if( flim == NULL ){
81 char buf[256] ;
82 sprintf(buf,"model_convgamma: Can't read timeseries file %s",cp) ;
83 ERREX(buf) ;
84 }
85
86 #if 0
87 fprintf(stderr,"conv_set_ref: refts=%s nx=%d\n",cp,flim->ny) ;
88 #endif
89
90 conv_set_ref( flim->nx , MRI_FLOAT_PTR(flim) ) ; /* recursion! */
91 mri_free(flim) ;
92 }
93 return ;
94 }
95
96 /*-----------------------------------------------------------------------
97 Function to compute the simulated time series.
98 -------------------------------------------------------------------------*/
99
conv_model(float * gs,int ts_length,float ** x_array,float * ts_array)100 void conv_model( float * gs , int ts_length ,
101 float ** x_array , float * ts_array )
102 {
103 int ii, jj,jbot,jtop , kk , nid_top,nid_bot ;
104 float top , val ;
105
106 static int nid = 0 ; /* number of pts in impulse */
107 static float * fid = NULL ; /* impulse response function */
108
109 /*** make sure there is a reference function to convolve with ***/
110
111 if( refnum <= 0 ) conv_set_ref( 0 , NULL ) ;
112
113 /*** initialize the output ***/
114
115 for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ;
116
117 /*** initialize the impulse response ***/
118
119 if( nid < ts_length ){ /* make some space for it */
120 if( fid != NULL ) free(fid) ;
121 nid = ts_length ;
122 fid = (float *) malloc( sizeof(float) * nid ) ;
123 }
124
125 gamma_model( gs , ts_length , x_array , fid ) ; /* compute impulse */
126
127 top = 0.0 ; /* find max value */
128 for( jj=0 ; jj < ts_length ; jj++ ){
129 val = fabs(fid[jj]) ; if( val > top ) top = val ;
130 }
131 if( top == 0.0 ) fid[0] = 1.0 ; /* very unlikely case */
132 top *= 0.001 ;
133 for( jj=0 ; jj < ts_length ; jj++ ){ /* discard small values */
134 if( fabs(fid[jj]) < top ) fid[jj] = 0.0 ;
135 }
136 for( nid_bot=0 ; nid_bot < ts_length ; nid_bot++ ) /* find first nonzero */
137 if( fid[nid_bot] != 0.0 ) break ;
138 for( nid_top=ts_length-1 ; nid_top > nid_bot ; nid_top-- ) /* and last nonzero */
139 if( fid[nid_top] != 0.0 ) break ;
140
141 /*** loop over each nonzero point in the reference ***/
142
143 for( ii=0 ; ii < refnz ; ii++ ){
144 kk = refin[ii] ; if( kk >= ts_length ) break ;
145 val = refts[kk] ;
146
147 /*** for each point in the impulse ***/
148
149 jtop = ts_length - kk ; if( jtop > nid_top ) jtop = nid_top+1 ;
150 for( jj=nid_bot ; jj < jtop ; jj++ )
151 ts_array[kk+jj] += val * fid[jj] ;
152 }
153
154 return ;
155 }
156
157 /*-----------------------------------------------------------------------*/
158
159 DEFINE_MODEL_PROTOTYPE
160
initialize_model()161 MODEL_interface * initialize_model ()
162 {
163 MODEL_interface * mi = NULL;
164
165 /*----- allocate memory space for model interface -----*/
166
167 mi = (MODEL_interface *) RwcMalloc (sizeof(MODEL_interface));
168
169 /*----- name of this model -----*/
170
171 strcpy (mi->label, "ConvGamma");
172
173 /*----- this is a signal model -----*/
174
175 mi->model_type = MODEL_SIGNAL_TYPE;
176
177 /*----- number of parameters in the model -----*/
178
179 mi->params = 4;
180
181 /*----- parameter labels -----*/
182
183 strcpy (mi->plabel[0], "t0");
184 strcpy (mi->plabel[1], "amp");
185 strcpy (mi->plabel[2], "r");
186 strcpy (mi->plabel[3], "b");
187
188 /*----- minimum and maximum parameter constraints -----*/
189
190 mi->min_constr[0] = 0.0; mi->max_constr[0] = 10.0;
191 mi->min_constr[1] = 0.0; mi->max_constr[1] = 200.0;
192 mi->min_constr[2] = 1.0; mi->max_constr[2] = 19.0;
193 mi->min_constr[3] = 0.1; mi->max_constr[3] = 5.0;
194
195 /*----- function which implements the model -----*/
196 mi->call_func = &conv_model;
197
198 return (mi);
199 }
200
201 /*----------------------------------------------------------------------*/
202 /*
203 Routine to calculate the time series which results from using the
204 gamma variate drug response signal model with the specified
205 model parameters.
206
207 Definition of model parameters:
208
209 gs[0] = time delay of response (t0)
210 gs[1] = multiplicative constant (k)
211 gs[2] = rise rate exponent (r)
212 gs[3] = decay rate constant (b)
213
214 f(t) = k * t^r * e^(-t/b)
215
216 */
217
gamma_model(float * gs,int ts_length,float ** x_array,float * ts_array)218 void gamma_model
219 (
220 float * gs, /* parameters for signal model */
221 int ts_length, /* length of time series data */
222 float ** x_array, /* independent variable matrix */
223 float * ts_array /* estimated signal model time series */
224 )
225
226 {
227 int it; /* time index */
228 float t; /* time */
229 float gsi , fac ;
230
231 if( gs[3] <= 0.0 || gs[2] <= 0.0 ){
232 for( it=0 ; it < ts_length ; it++ ) ts_array[it] = 0.0 ;
233 return ;
234 }
235
236 /* fac is chosen to make the peak value equal to gs[1] */
237
238 gsi = 1.0 / gs[3] ;
239 fac = gs[1] * exp( gs[2] * ( 1.0 - log(gs[2]*gs[3]) ) ) ;
240 for( it=0; it < ts_length; it++){
241 t = x_array[it][1] - gs[0] ;
242 ts_array[it] = (t <= 0.0) ? 0.0
243 : fac * exp( log(t) * gs[2] - t * gsi ) ;
244 }
245 return ;
246 }
247