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 Based on model_conv_diffgamma.c.
7
8 Written for C Connolly and Felix, requested on the AFNI message board,
9 regarding a "hemodynamic model incorporating an initial dip".
10
11 The model is from:
12
13 Fully Bayesian Spatio-Temporal Modeling of FMRI Data
14 IEEE Transactions on Medical Imaging,
15 Volume 23, Issue 2, February 2004, Pages 213-231
16 Woolrich, M.W., Jenkinson, M., Brady, J.M., Smith, S.M.
17
18 The model consists of 4 half-cosines, of duration m1, m2, m3, m4:
19
20 - the first half falling from 0 down to -c1
21 - the second half rising from -c1 up to 1
22 - the next first half falling from 1 down to -c2
23 - the next second half rising from -c2 back to 0
24
25 - plus an overall amplitude A, of course
26
27 This function will process parameters with amplitudes first:
28
29 A, c1, c2, m1, m2, m3, m4
30
31 Note that all parameters returned should be non-negative, except
32 possibly A. For a somewhat standard response to a brief event, one
33 might expect parameters near:
34
35 c1 = .25
36 c2 = .5
37 m1 = 1.0
38 m2 = 3.0
39 m3 = 4.0
40 m4 = 4.0
41
42 But of course, estimating the parameters is the point of this function.
43
44 R. Reynolds 6 Sep 2013
45 ******************************************************************************/
46
47 #include "NLfit_model.h"
48
49 static int refnum = 0 ; /* # pts in refts */
50 static int refnz = 0 ; /* # of nonzero pts */
51 static int g_debug = 0 ; /* debug level */
52 static int g_diter = -1 ; /* debug iteration number */
53 static float * refts = NULL ; /* reference time series */
54 static int * refin = NULL ; /* indexes of nonzero pts */
55
56 int signal_model( float * , int , float ** , float *, int );
57
58 static int disp_floats(char * mesg, float * p, int len);
59 static int model_help(void);
60
61 void conv_model( float * gs , int ts_length ,
62 float ** x_array , float * ts_array );
63
64 #define ERREX(str) ( fprintf(stderr,"\n*** %s\a\n",str) , exit(1) )
65
66 /*----------------------------------------------------------------------
67 Function to set the reference time series, with which the
68 model function is convolved to produce the simulated data.
69 ------------------------------------------------------------------------*/
70
conv_set_ref(int num,float * ref)71 void conv_set_ref( int num , float * ref )
72 {
73 if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/
74 int ii ;
75
76 /* get rid of old data */
77
78 if(refts != NULL){ free(refts); refts = NULL; free(refin); refin = NULL; }
79
80 refnum = num ;
81 refts = (float *) malloc( sizeof(float) * num ) ;
82 refin = (int *) malloc( sizeof(int) * num ) ;
83 memcpy( refts , ref , sizeof(float) * num ) ;
84 for( ii=0,refnz=0 ; ii < num ; ii++ ) /* build list of nonzero */
85 if( refts[ii] != 0.0 ) refin[refnz++] = ii ; /* points in refts */
86 if( refnz == 0 )
87 ERREX("model_conv_cosine4: All zero reference timeseries!") ;
88
89 if( g_debug ) {
90 fprintf(stderr,"+d conv_set_ref: num=%d nonzero=%d\n",num,refnz) ;
91 if( g_debug > 1 ) {
92 fprintf(stderr," TR locked stimuli :");
93 for( ii = 0; ii < refnz; ii++ ) fprintf(stderr," %d", refin[ii]);
94 fputc('\n',stderr);
95 }
96 }
97
98 return ;
99
100 } else { /*** if no inputs, read it from AFNI_CONVMODEL_REF 1D file ***/
101
102 char * cp ;
103 MRI_IMAGE * flim ;
104
105 cp = my_getenv("AFNI_CONVMODEL_REF") ; /* get name of reference file */
106 if( cp == NULL )
107 ERREX("model_conv_cosine4: need ref file as AFNI_CONVMODEL_REF") ;
108
109 flim = mri_read_1D(cp) ; /* 16 Nov 1999: replaces mri_read_ascii */
110 if( flim == NULL ){
111 char buf[256] ;
112 sprintf(buf,"model_conv_cosine4: Can't read timeseries file %s",cp) ;
113 ERREX(buf) ;
114 }
115
116 if( g_debug )
117 fprintf(stderr,"+d conv_set_ref: refts=%s nx=%d\n",cp,flim->ny) ;
118
119 conv_set_ref( flim->nx , MRI_FLOAT_PTR(flim) ) ; /* recursion! */
120 mri_free(flim) ;
121 }
122 return ;
123 }
124
125 /*-----------------------------------------------------------------------
126 Function to compute the simulated time series.
127 -------------------------------------------------------------------------*/
128
conv_model(float * gs,int ts_length,float ** x_array,float * ts_array)129 void conv_model( float * gs , int ts_length ,
130 float ** x_array , float * ts_array )
131 {
132 int ii, jj,jtop , kk;
133 int cur_debug = 0, irfdur=0;
134 float val;
135
136 static int iter = -1; /* iteration number */
137
138 static int nid = 0 ; /* number of pts in impulse */
139 static float * fid = NULL; /* impulse response function */
140
141 /*----- check for env vars -----*/
142
143 iter++ ;
144 if( iter == 0 ) {
145 double dval = AFNI_numenv("AFNI_MODEL_DITER");
146 if( dval >= 1.0 ) g_diter = (int)dval; /* zero is failure */
147 dval = AFNI_numenv("AFNI_MODEL_DEBUG");
148 if( dval >= 1.0 ) g_debug = (int)dval;
149 if(g_debug) fprintf(stderr,"\n+d TR = %f\n", x_array[1][1]-x_array[0][1]);
150 }
151
152 /*** make sure there is a reference function to convolve with ***/
153
154 if( refnum <= 0 ) conv_set_ref( 0 , NULL ) ;
155
156 /* to clean up, particularly as a parameter */
157 cur_debug = (iter == g_diter || (iter == 0 && g_debug > 1));
158 if( cur_debug ) disp_floats("+d params: ", gs, 8);
159
160 /*** initialize the output ***/
161
162 for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ;
163
164 /*** initialize the impulse response ***/
165
166 if( nid < ts_length ){ /* make some space for it */
167 if( fid ) free(fid) ;
168 nid = ts_length ;
169 fid = (float *) malloc( sizeof(float) * nid ) ;
170 }
171
172 /* compute first and second impulse functions */
173 irfdur = signal_model(gs, ts_length, x_array, fid, cur_debug);
174
175 /*** loop over each nonzero point in the reference ***/
176
177 /* TR-locked convolution */
178 for( ii=0 ; ii < refnz ; ii++ ){
179 kk = refin[ii] ; if( kk >= ts_length ) break ;
180 val = refts[kk] ;
181
182 /* for each point in the impulse, add its val times irf */
183 /* (top index offset is min(irfdur, ts_length-kk-1)) */
184 jtop = ts_length - kk ; if( jtop > irfdur ) jtop = irfdur ;
185 for( jj=0 ; jj < jtop ; jj++ )
186 ts_array[kk+jj] += val * fid[jj];
187 }
188
189 if( cur_debug ) disp_floats("+d conv : ", ts_array, ts_length);
190
191 return ;
192 }
193
disp_floats(char * mesg,float * p,int len)194 static int disp_floats(char * mesg, float * p, int len)
195 {
196 int c;
197 if( mesg ) fputs(mesg, stderr);
198 for( c = 0; c < len; c++ ) fprintf(stderr," %f ", p[c]);
199 fprintf(stderr,"\n\n");
200 return 0;
201 }
202
203
204 /*-----------------------------------------------------------------------*/
205
206 DEFINE_MODEL_PROTOTYPE
207
initialize_model()208 MODEL_interface * initialize_model ()
209 {
210 MODEL_interface * mi = NULL;
211
212 /*----- first, see if the user wants help -----*/
213 if( AFNI_yesenv("AFNI_MODEL_HELP_CONV_COSINE4") ||
214 AFNI_yesenv("AFNI_MODEL_HELP_ALL") ) model_help();
215
216 /*----- allocate memory space for model interface -----*/
217
218 mi = (MODEL_interface *) RwcMalloc (sizeof(MODEL_interface));
219
220 /*----- name of this model -----*/
221
222 strcpy (mi->label, "ConvCosine4");
223
224 /*----- this is a signal model -----*/
225
226 mi->model_type = MODEL_SIGNAL_TYPE;
227
228 /*----- number of parameters in the model -----*/
229
230 mi->params = 7;
231
232 /*----- parameter labels -----*/
233
234 strcpy (mi->plabel[0], "A"); /* might get fixed at 1.0 here */
235 strcpy (mi->plabel[1], "C1");
236 strcpy (mi->plabel[2], "C2");
237 strcpy (mi->plabel[3], "M1");
238 strcpy (mi->plabel[4], "M2");
239 strcpy (mi->plabel[5], "M3");
240 strcpy (mi->plabel[6], "M4");
241
242 /*----- minimum and maximum parameter constraints -----*/
243
244 /* first 3 are amplitudes, last 4 are durations */
245 /* (first amplitude may be fixed at 1.0 and computed externally) */
246 mi->min_constr[0] = -500.0; mi->max_constr[0] = 500.0;
247
248 mi->min_constr[1] = 0.0; mi->max_constr[1] = 1.0;
249 mi->min_constr[2] = 0.0; mi->max_constr[2] = 5.0;
250
251 mi->min_constr[3] = 0.0; mi->max_constr[3] = 5.0;
252 mi->min_constr[4] = 0.0; mi->max_constr[4] = 20.0;
253 mi->min_constr[5] = 0.0; mi->max_constr[5] = 20.0;
254 mi->min_constr[6] = 0.0; mi->max_constr[6] = 20.0;
255
256 /*----- function which implements the model -----*/
257 mi->call_func = conv_model;
258
259 return (mi);
260 }
261
262
263 /*----------------------------------------------------------------------*
264 * for use in signal model
265 *----------------------------------------------------------------------*/
266 #ifdef PI
267 #undef PI
268 #endif
269 #define PI 3.141592653589793238462643
270
271 /* F(t) is cosine on one interval:
272 *
273 * HEIGHT : height of single cosine (i.e. cos(t) gets scaled by HEIGHT)
274 * POFF : PI offset, either 0 or PI, depending on where curve starts
275 * TOFF : time offset = 0, M1, M1+M2 or M1+M2+M3
276 * DUR : duration, Mi
277 * NBASE : negative base value (minimum), i.e. C1 or C2
278 * t : current time value
279 */
280 #ifdef F
281 #undef F
282 #endif
283 #define F(HEIGHT,POFF,TOFF,DUR,NBASE,t) \
284 ((HEIGHT) * cos((POFF) + PI*((t)-(TOFF)) / (DUR)) + (HEIGHT) - (NBASE))
285
286 /*----------------------------------------------------------------------*/
287 /*
288 Routine to calculate the time series which results from using the
289 gamma variate drug response signal model with the specified
290 model parameters.
291
292 Note: gs[0,1] were reversed from those in model_convgamma.
293
294 Definition of model parameters (all >=0, except for gs[0])
295
296 gs[0] = A = multiplicitive constant
297 gs[1] = C1 = magnitude of initial undershoot
298 gs[2] = C2 = magnitude of post-undershoot
299 gs[3] = M1 = duration of fall from 0 to min pre-undershoot
300 gs[4] = M2 = duration of rise from pre-undershoot to peak
301 gs[5] = M3 = duration of fall from peak to min post-undershoot
302 gs[6] = M4 = duration of rise from post-undershoot back to 0
303
304 Return the non-zero duration of the resulting impulse response function,
305 as array length.
306 */
signal_model(float * gs,int ts_length,float ** x_array,float * ts_array,int debug)307 int signal_model
308 (
309 float * gs, /* parameters for signal model */
310 int ts_length, /* length of time series data */
311 float ** x_array, /* independent variable matrix */
312 float * ts_array, /* estimated signal model time series */
313 int debug /* make some noise */
314 )
315 {
316 int it, nt; /* time index, num non-zero points */
317 double t; /* time */
318 double A, C1, C2, M1, M2, M3, M4; /* input parameters */
319 double H1, H2, H3, H4; /* curve heights */
320 double T02, T03, T04, T05; /* time interval dividers */
321
322 /* assign parameters */
323 A = gs[0];
324 C1 = gs[1]; C2 = gs[2];
325 M1 = gs[3]; M2 = gs[4]; M3 = gs[5]; M4 = gs[6];
326
327 if( debug ) fprintf(stderr,"-d A=%.3f, C1=%.3f, C2=%.3f\n,"
328 "M1=%.3f, M2=%.3f, M3=%.3f, M4=%.3f\n",
329 A, C1, C2, M1, M2, M3, M4);
330
331 /* make a list of t0 values, skipping T01=0.0 */
332 T02 = M1;
333 T03 = M1 + M2;
334 T04 = M1 + M2 + M3;
335 T05 = M1 + M2 + M3 + M4;
336
337 /* and heights */
338 H1 = C1/2.0;
339 H2 = (1.0 + C1)/2.0;
340 H3 = (1.0 + C2)/2.0;
341 H4 = C2/2.0;
342
343 /* do not assume time array is ordered? */
344 nt = 0; /* largest time point applied */
345 for( it=0; it < ts_length; it++){
346 t = x_array[it][1];
347
348 /* if not in bounds, skip this time */
349 if( t < 0 || t > T05 ) { ts_array[it] = 0.0; continue; }
350
351 /* track number of usable indices */
352 nt = it + 1;
353
354 /* set values per interval (do not assume time is ordered or regular?) */
355 if ( t < T02 ) ts_array[it] = F(H1, 0.0, 0.0, M1, C1, t);
356 else if( t < T03 ) ts_array[it] = F(H2, PI, T02, M2, C1, t);
357 else if( t < T04 ) ts_array[it] = F(H3, 0.0, T03, M3, C2, t);
358 else /* t <= T05 */ ts_array[it] = F(H4, PI, T04, M4, C2, t);
359 }
360
361 for( it=0; it < nt; it++) ts_array[it] *= A; /* scale by amplitude */
362
363 if( debug )
364 disp_floats("+d signal model : ", ts_array, ts_length);
365
366 return nt;
367 }
368
369 /*----------------------------------------------------------------------*/
model_help(void)370 static int model_help(void)
371 {
372 printf(
373 "----------------------------------------------------------------------\n"
374 "ConvCosine4 - sum of 4 half-cosine curves\n"
375 "\n"
376 " This model is from the paper:\n"
377 "\n"
378 " Fully Bayesian Spatio-Temporal Modeling of FMRI Data\n"
379 " IEEE Transactions on Medical Imaging,\n"
380 " Volume 23, Issue 2, February 2004, Pages 213-231\n"
381 " Woolrich, M.W., Jenkinson, M., Brady, J.M., Smith, S.M.\n"
382 "\n"
383 " The model consists of 4 half-cosines, of duration m1, m2, m3, m4,\n"
384 " respectively, all scaled by an amplitude A:\n"
385 "\n"
386 " - the first half falling from 0 down to -c1\n"
387 " - the second half rising from -c1 up to 1\n"
388 " - the next first half falling from 1 down to -c2\n"
389 " - the next second half rising from -c2 back to 0\n"
390 "\n"
391 " 1 | /\\\n"
392 " 0___ | / \\\n"
393 " -c1 | \\/ \\ / - each of 4 'segments' is a half cosine, either\n"
394 " | \\/ -c2 falling (first half) or rising (second half)\n"
395 " +------------\n"
396 " || | | |\n"
397 " HI J K L\n"
398 "\n"
399 " time H = 0 (starting time of response)\n"
400 " I = m1 (duration of falling segment)\n"
401 " J = m1+m2 (plus duration m2 of rising segment)\n"
402 " K = m1+m2+m3 (plus duration m3 of falling segment)\n"
403 " L = m1+m2+m3+m4 (plus duration m4 of rising segment)\n"
404 "\n"
405 " This function will process parameters with amplitudes first:\n"
406 "\n"
407 " A, c1, c2, m1, m2, m3, m4\n"
408 "\n"
409 " Note that all parameters returned should be positive, except\n"
410 " possibly A. For a somewhat standard response to a brief event,\n"
411 " one might expect parameters near or maybe 50%% longer:\n"
412 "\n"
413 " c1 = .25\n"
414 " c2 = .5\n"
415 " m1 = 1.0\n"
416 " m2 = 3.0\n"
417 " m3 = 4.0\n"
418 " m4 = 4.0\n"
419 "\n"
420 "--------------------------------------------------\n"
421 "To use this model function:\n"
422 "\n"
423 " 1. Generate a stimulus event file. This is like a stimulus timing\n"
424 " file but is just a 1D file of TR-locked events, where a 1 denotes\n"
425 " a stimulus event. Note that a 0/1 file is not necessary. If there\n"
426 " is a reason to use a range of stim event magnitudes, they will be\n"
427 " convolved into the model, if that is what the user wishes\n"
428 "\n"
429 " This should be a vertical text file of stimulus events or magnitudes.\n"
430 " A simple example would be just having an event at the first TR. The\n"
431 " file could be a single line containing '1' in such a case, e.g.\n"
432 "\n"
433 " echo 1 > cos4_stim_file.1D\n"
434 "\n"
435 " If there is interest, using a non-TR-locked stimulus timing file\n"
436 " could be added.\n"
437 "\n"
438 " 2. Set the environment variable AFNI_CONVMODEL_REF to the name of the\n"
439 " stimulus event file from step 1.\n"
440 "\n"
441 " tcsh: setenv AFNI_CONVMODEL_REF cos4_stim_file.1D\n"
442 " bash: export AFNI_CONVMODEL_REF=cos4_stim_file.1D\n"
443 "\n"
444 " 3. Run 3dNLfim as usual, e.g.,\n"
445 "\n"
446 " 3dNLfim -input epi_data+tlrc \\\n"
447 " -noise Quadratic \\\n"
448 " -signal ConvCosine4 \\\n"
449 " -sconstr 0 0.8 1.2 \\\n"
450 " -sconstr 1 0.2 0.3 \\\n"
451 " -sconstr 2 0.3 0.7 \\\n"
452 " -sconstr 3 0.5 1.5 \\\n"
453 " -sconstr 4 2 4 \\\n"
454 " -sconstr 5 3 5 \\\n"
455 " -sconstr 6 3 5 \\\n"
456 " -BOTH \\\n"
457 " -ignore 0 \\\n"
458 " -nrand 10000 \\\n"
459 " -nbest 10 \\\n"
460 " -jobs 4 \\\n"
461 " -bucket 0 buck.cos4 \\\n"
462 " -snfit snfit.cos4\n"
463 "\n"
464 "--------------------------------------------------\n"
465 " For C Connolly and Felix, requested on the AFNI message board.\n"
466 "\n"
467 " R. Reynolds 10 Sep 2013\n"
468 "----------------------------------------------------------------------\n"
469 );
470
471 return 0 ;
472 }
473
474