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_expMEMRI.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, "expMEMRI");
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 = 3;
63
64 /*----- parameter labels -----*/
65 strcpy (mi->plabel[0], "a");
66 strcpy (mi->plabel[1], "b");
67 strcpy (mi->plabel[2], "c scalefactor");
68
69 /*----- minimum and maximum parameter constraints -----*/
70 mi->min_constr[0] = 0.0; mi->max_constr[0] = 100.0;
71 mi->min_constr[1] = -100.00; mi->max_constr[1] = 0.00;
72 mi->min_constr[2] = 0.0; mi->max_constr[2] = 100.0;
73
74 /*----- function which implements the model -----*/
75 mi->call_func = &signal_model;
76
77 #if 0
78 envp = my_getenv("AFNI_MODEL_EXPMEMRI_T_FILE");
79 if( !envp )
80 {
81 fprintf(stderr,"\n** NLfim: need env var AFNI_EXPMEMRI_T_FILE\n");
82 fprintf(stderr," (a 1D file of time values for each sub-brick)\n");
83 return(NULL);
84 }
85
86 im = mri_read_1D(envp);
87 if( !im )
88 {
89 fprintf(stderr,"** failed to open time file %s for model\n", envp);
90 return(NULL);
91 }
92
93 /* nx == 1 and ny > 1, take the transpose */
94 if( im->nx == 1 && im->ny > 1 )
95 {
96 MRI_IMAGE * flim = mri_transpose(im);
97 mri_free(im);
98 im = flim;
99 if( !im ) { fprintf(stderr,"** time transposing failure\n"); return(NULL); }
100 fprintf(stderr,"taking transpose of time file, new len = %d\n",im->nx);
101 }
102
103 rt = MRI_FLOAT_PTR(im); /* do not free this */
104 rt_len = im->nx;
105 fprintf(stderr,"Times at each sub-brick :\n");
106 for(ii=0;ii<im->nx;ii++)
107 fprintf(stderr, "%f ", rt[ii]);
108 fprintf(stderr, "\n");
109
110 #endif
111
112 /*----- return pointer to the model interface -----*/
113 return (mi);
114 }
115
116
117 /*---------------------------------------------------------------------------*/
118 /*
119 Routine to calculate the time series which results from using the
120 an exponential signal model with the specified
121 model parameters.
122
123 Definition of model parameters:
124
125 y = c / (1 + exp(a+bt))
126
127 gs[0] = a
128 gs[1] = b
129 gs[2] = c
130
131 */
132
signal_model(float * gs,int ts_length,float ** x_array,float * ts_array)133 void signal_model
134 (
135 float * gs, /* parameters for signal model */
136 int ts_length, /* length of time series data */
137 float ** x_array, /* independent variable matrix */
138 float * ts_array /* estimated signal model time series */
139 )
140
141 {
142 int it; /* time index */
143 float t; /* time */
144 float fval; /* time series value at time t */
145
146
147 for (it = 0; it < ts_length; it++)
148 {
149 t = x_array[it][1];
150 fval = gs[2] / (1.0+exp(gs[0] + gs[1]*t));
151 ts_array[it] = fval;
152 }
153
154 }
155
156
157
158
159