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 This file contains routines to initialize and implement the
9 differential exponential signal model for MEMRI applications
10
11 File: model_expMEMRI2.c
12 Author: D. Glen based on model_exp.c by Z.Saad and
13 model_diffexp.c by B. Douglas Ward
14 */
15
16
17 /*---------------------------------------------------------------------------*/
18
19 #include <math.h>
20 #include "NLfit_model.h"
21
22 void signal_model
23 (
24 float * gs, /* parameters for signal model */
25 int ts_length, /* length of time series data */
26 float ** x_array, /* independent variable matrix */
27 float * ts_array /* estimated signal model time series */
28 );
29
30 /*static float *rt;
31 static int rt_len;*/
32
33 /*---------------------------------------------------------------------------*/
34 /*
35 Routine to initialize the signal model by defining the number of parameters
36 in the signal model, the name of the signal model, and the default values
37 for the minimum and maximum parameter constraints.
38 */
39
40 DEFINE_MODEL_PROTOTYPE
41
initialize_model()42 MODEL_interface * initialize_model ()
43 {
44 MODEL_interface * mi = NULL;
45 MRI_IMAGE * im;
46 char * envp;
47 int ii;
48
49 /*----- allocate memory space for model interface -----*/
50 mi = (MODEL_interface *) RwcMalloc (sizeof(MODEL_interface));
51
52
53 /*----- define interface for the differential - exponential model -----*/
54
55 /*----- name of this model -----*/
56 strcpy (mi->label, "expMEMRI3");
57
58 /*----- this is a signal model -----*/
59 mi->model_type = MODEL_SIGNAL_TYPE;
60
61 /*----- number of parameters in the model -----*/
62 mi->params = 5;
63
64 /*----- parameter labels -----*/
65 strcpy (mi->plabel[0], "a - wash-in delay");
66 strcpy (mi->plabel[1], "b - wash-in rate");
67 strcpy (mi->plabel[2], "c - wash-in scalefactor");
68 strcpy (mi->plabel[3], "d - wash-out delay");
69 strcpy (mi->plabel[4], "e - wash-out rate");
70 /* strcpy (mi->plabel[5], "f - wash-out scalefactor");*/
71
72 /*----- minimum and maximum parameter constraints -----*/
73 mi->min_constr[0] = 0.0; mi->max_constr[0] = 100.0;
74 mi->min_constr[1] = -100.00; mi->max_constr[1] = 0.00;
75 mi->min_constr[2] = 0.0; mi->max_constr[2] = 100.0;
76 mi->min_constr[3] = 0.0; mi->max_constr[3] = 500.0;
77 mi->min_constr[4] = -100.00; mi->max_constr[4] = 0.00;
78 /* mi->min_constr[5] = 0.0; mi->max_constr[5] = 100.0;*/
79
80 /*----- function which implements the model -----*/
81 mi->call_func = &signal_model;
82
83 #if 0
84 envp = my_getenv("AFNI_MODEL_EXPMEMRI_T_FILE");
85 if( !envp )
86 {
87 fprintf(stderr,"\n** NLfim: need env var AFNI_EXPMEMRI_T_FILE\n");
88 fprintf(stderr," (a 1D file of time values for each sub-brick)\n");
89 return(NULL);
90 }
91
92 im = mri_read_1D(envp);
93 if( !im )
94 {
95 fprintf(stderr,"** failed to open time file %s for model\n", envp);
96 return(NULL);
97 }
98
99 /* nx == 1 and ny > 1, take the transpose */
100 if( im->nx == 1 && im->ny > 1 )
101 {
102 MRI_IMAGE * flim = mri_transpose(im);
103 mri_free(im);
104 im = flim;
105 if( !im ) { fprintf(stderr,"** time transposing failure\n"); return(NULL); }
106 fprintf(stderr,"taking transpose of time file, new len = %d\n",im->nx);
107 }
108
109 rt = MRI_FLOAT_PTR(im); /* do not free this */
110 rt_len = im->nx;
111 fprintf(stderr,"Times at each sub-brick :\n");
112 for(ii=0;ii<im->nx;ii++)
113 fprintf(stderr, "%f ", rt[ii]);
114 fprintf(stderr, "\n");
115
116 #endif
117
118 /*----- return pointer to the model interface -----*/
119 return (mi);
120 }
121
122
123 /*---------------------------------------------------------------------------*/
124 /*
125 Routine to calculate the time series which results from using the
126 an exponential signal model with the specified
127 model parameters.
128
129 Definition of model parameters:
130
131 y = c / (1 + exp(a+bt)) - c/(1+exp(d+et)) (baseline in noise model +g)
132
133 gs[0] = a
134 gs[1] = b
135
136 */
137
signal_model(float * gs,int ts_length,float ** x_array,float * ts_array)138 void signal_model
139 (
140 float * gs, /* parameters for signal model */
141 int ts_length, /* length of time series data */
142 float ** x_array, /* independent variable matrix */
143 float * ts_array /* estimated signal model time series */
144 )
145
146 {
147 int it; /* time index */
148 float t; /* time */
149 float fval; /* time series value at time t */
150
151
152 for (it = 0; it < ts_length; it++)
153 {
154 t = x_array[it][1];
155 fval = gs[2] / (1.0+exp(gs[0] + gs[1]*t)) - \
156 gs[2] / (1.0+exp(gs[3] + gs[4]*t));
157 ts_array[it] = fval;
158 }
159
160 }
161
162
163
164
165