1 #include "mrilib.h"
2 #include "zgaussian.c"
3
4 #include <time.h>
5 #include <unistd.h>
6
7 #undef LAMBDA
8 #define LAMBDA(a,b) ((b+a)*(1.0+a*b)/(1.0+2.0*a*b+b*b))
9
main(int argc,char * argv[])10 int main( int argc , char *argv[] )
11 {
12 int nxy=256 , nxyq , nz=1 , ntr=200 , ndiv=16 , nbb=3 , iarg ;
13 float TR=2.5f , blur=5.0f , base=30.0f , bsig=3.0f ;
14 float amx=0.8f , bmx=0.5f , amn=0.0f , bmn=-0.5f ;
15 float taumx=3.0f , sigmx=5.0f , sigmn=1.0f ;
16 MRI_IMAGE *rim=NULL ;
17 char *prefix="Simm" ;
18
19 THD_3dim_dataset *dset ;
20 THD_ivec3 ixyz , oxyz ; THD_fvec3 dxyz ;
21 int dvx,dvy,xx,yy,zz,tt , kk , kold ;
22 float *tsar , val , wal , a0,a1 , b0,b1 , dsig , aa,bb,lam ;
23 MRI_IMAGE *armim ; float *armar ;
24 long seed ;
25
26 /*-- the pitiful help --*/
27
28 if( argc == 1 || strcmp(argv[1],"-help") == 0 ){
29 printf(
30 "This program generates a simulated FMRI time series dataset\n"
31 "for the nefarious and inscrutable purposes of RW Cox and G Chen.\n"
32 "At least one option must be given (or you get this help message)!\n"
33 "\n"
34 "Parameters to determine structure of output dataset\n"
35 "---------------------------------------------------\n"
36 "-nxy xx = size of 2D image to produce [default=256 ]\n"
37 "-nz = number of slices to produce [default=1 ]\n"
38 "-div xx = size of image subdivisions to use [default=16 ]\n"
39 "-nTR xx = number of time points to produce [default=200 ]\n"
40 "-TR xx = TR to use for the output dataset [default=2.5 ]\n"
41 "-pre pp = output dataset prefix [default=Simm]\n"
42 "\n"
43 "Parameters to determine baseline in each voxel time series\n"
44 " [set bas==0 and nbb==0 to skip baseline production]\n"
45 "----------------------------------------------------------\n"
46 "-bas xx = mean value to use for constant baseline [default=30 ]\n"
47 "-nbb xx = number of baseline polynomials to produce [default=3 ]\n"
48 "-vbb xx = stdev of baseline polynomial coefficients [default=3 ]\n"
49 "-blr xx = FWHM of blurring to be applied to baseline [default=5 ]\n"
50 "\n"
51 "Parameters to determine random noise added to each time series\n"
52 " [set mxa==mna mxb==mnb smx==smn for spatial homogeneity]\n"
53 "--------------------------------------------------------------\n"
54 "-mxa xx = max value of ARMA11 parameter 'a' to use [default=+0.8]\n"
55 "-mna xx = min value of ARMA11 parameter 'a' to use [default= 0.0]\n"
56 "-mxb xx = max value of ARMA11 parameter 'b' to use [default=+0.5]\n"
57 "-mnb xx = max value of ARMA11 parameter 'b' to use [default=-0.5]\n"
58 "-smx xx = max value of time series stdev to use [default=5 ]\n"
59 "-smn xx = min value of time series stdev to use [default=1 ]\n"
60 "\n"
61 "Parameters to add a signal component to each time series\n"
62 " [if you want to test 3dREMLfit with a signal model]\n"
63 "--------------------------------------------------------\n"
64 "-1D xx = 1D file to use for response model [no default ]\n"
65 "-tmx xx = max value of tau (signal amplitude) to use [default=3 ]\n"
66 "\n"
67 "There is very little error checking. If you do something stooopid\n"
68 "(e.g., '-mxa 2.0'), then that's just tough for you.\n"
69 ) ;
70 exit(0) ;
71 }
72
73 /*-- process command line options --*/
74
75 #undef OOP
76 #define OOP(ostr,ovar) \
77 if( strcmp(argv[iarg],(ostr)) == 0 ){ \
78 iarg++; if( iarg >= argc )ERROR_exit("need arg after %s",argv[iarg-1]); \
79 (ovar) = strtod(argv[iarg],NULL) ; iarg++ ; continue ; \
80 } else
81
82 iarg = 1 ;
83 while( iarg < argc ){
84
85
86 if( strcmp(argv[iarg],"-tdof") == 0 ){ /* hidden option */
87 float tdof=0.0f , atd=0.0f ;
88 iarg++ ; if( iarg >= argc ) ERROR_exit("need argument after -tdof!") ;
89 tdof = (float)strtod(argv[iarg],NULL) ; atd = fabsf(tdof) ;
90 if( atd > 0.0f && atd < 4.0f ) ERROR_exit("illegal value after -tdof") ;
91 if( tdof != 0.0f ){
92 INFO_message("Transforming Gaussian deviates to t(DOF=%g)",tdof) ;
93 mri_genARMA11_set_tdof(tdof) ;
94 }
95 iarg++ ; continue ;
96 }
97
98 if( strncmp(argv[iarg],"-prefix",4) == 0 ){
99 iarg++; if( iarg >= argc ) ERROR_exit("need arg after %s",argv[iarg-1]);
100 prefix = strdup(argv[iarg]) ;
101 iarg++ ; continue ;
102 }
103
104 if( strcmp(argv[iarg],"-1D") == 0 ){
105 iarg++; if( iarg >= argc ) ERROR_exit("need arg after %s",argv[iarg-1]);
106 rim = mri_read_1D(argv[iarg]) ;
107 if( rim == NULL ) ERROR_exit("Can't read -1D file '%s'",argv[iarg]) ;
108 iarg++ ; continue ;
109 }
110
111 OOP("-nxy",nxy ) ;
112 OOP("-div",ndiv ) ;
113 OOP("-nTR",ntr ) ;
114 OOP("-TR" ,TR ) ;
115 OOP("-bas",base ) ;
116 OOP("-nbb",nbb ) ;
117 OOP("-vbb",bsig ) ;
118 OOP("-blr",blur ) ;
119 OOP("-mxa",amx ) ;
120 OOP("-mxb",bmx ) ;
121 OOP("-mna",amn ) ;
122 OOP("-mnb",bmn ) ;
123 OOP("-tmx",taumx) ;
124 OOP("-smx",sigmx) ;
125 OOP("-smn",sigmn) ;
126 OOP("-nz" ,nz ) ;
127
128 ERROR_exit("Unknown option '%s'",argv[iarg]) ;
129 }
130
131 if( rim != NULL && rim->nx < ntr )
132 ERROR_exit("-1D file length %d than nTR=%d",rim->nx,ntr) ;
133
134 /*-- create output dataset ---*/
135
136 dset = EDIT_empty_copy(NULL) ;
137
138 LOAD_IVEC3(ixyz,nxy,nxy,nz) ;
139 LOAD_IVEC3(oxyz,ORI_L2R_TYPE,ORI_P2A_TYPE,ORI_I2S_TYPE) ;
140 LOAD_FVEC3(dxyz,-1.0,-1.0,1.0) ;
141 EDIT_dset_items( dset ,
142 ADN_nxyz , ixyz ,
143 ADN_xyzorient , oxyz ,
144 ADN_xyzdel , dxyz ,
145 ADN_prefix , prefix ,
146 ADN_nvals , ntr ,
147 ADN_ntt , ntr ,
148 ADN_ttdel , TR ,
149 ADN_ttorg , 0.0 ,
150 ADN_nsl , 0 ,
151 ADN_tunits , UNITS_SEC_TYPE ,
152 ADN_datum_all , MRI_float ,
153 ADN_none ) ;
154 for( tt=0 ; tt < ntr ; tt++ ) /* attach zero filled volumes */
155 EDIT_substitute_brick(dset,tt,MRI_float,NULL) ;
156
157 tsar = (float *)malloc(sizeof(float)*ntr) ;
158 seed = (long)GSEED ; srand48(seed) ; INFO_message("seed = %ld",seed) ;
159
160 /* build baseline model into dataset */
161
162 nxyq = nxy * nxy ;
163
164 if( base != 0.0f || (nbb > 0 && bsig > 0.0f) ){
165 INFO_message("3dSimARMA11: build baseline model") ;
166
167 for( zz=0 ; zz < nz ; zz++ ){
168 for( yy=0 ; yy < nxy ; yy++ ){
169 for( xx=0 ; xx < nxy ; xx++ ){
170 val = zgaussian() ; val *= val * base ;
171 for( tt=0 ; tt < ntr ; tt++ ) tsar[tt] = val ;
172 wal = 1.99998f / (ntr-1.0f) ;
173 for( kk=1 ; kk <= nbb ; kk++ ){
174 val = bsig * zgaussian() ;
175 for( tt=0 ; tt < ntr ; tt++ )
176 tsar[tt] += val * Plegendre( wal*tt-0.99999 , kk ) ;
177 }
178 THD_insert_series( xx+yy*nxy+zz*nxyq , dset , ntr , MRI_float , tsar , 1 ) ;
179 }}}
180 } else {
181 blur = 0.0f ;
182 }
183
184 /* smooth baseline model */
185
186 if( blur > 0.0f ){
187 INFO_message("blur baseline model") ;
188 val = FWHM_TO_SIGMA(blur) ;
189 for( tt=0 ; tt < ntr ; tt++ )
190 EDIT_blur_volume( nxy,nxy,nz , 1.0f,1.0f,1.0f ,
191 MRI_float , DSET_ARRAY(dset,tt) , val ) ;
192 }
193
194 /* add signal model */
195
196 if( rim == NULL )
197 INFO_message("no -1D file ==> 'signal' strength will be zero") ;
198
199 if( taumx > 0.0f && rim != NULL ){
200 float *rar=MRI_FLOAT_PTR(rim) , dtau , tau , beta ;
201 float *bar ; int nbar ;
202
203 INFO_message("adding signal model") ;
204
205 nbar = (nxy-1)/ndiv + 1 ; dtau = taumx / (nbar-1) ;
206 bar = (float *)malloc(sizeof(float)*nbar*nbar) ;
207
208 for( zz=0 ; zz < nz ; zz++ ){
209 for( yy=0 ; yy < nbar ; yy++ )
210 for( xx=0 ; xx < nbar ; xx++ )
211 bar[xx+yy*nbar] = zgaussian() * xx*dtau ;
212 for( xx=0 ; xx < nxy ; xx++ ){
213 dvx = xx/ndiv ;
214 for( yy=0 ; yy < nxy ; yy++ ){
215 dvy = yy/ndiv ;
216 beta = bar[dvx+dvy*nbar] ;
217 if( beta != 0.0f ){
218 THD_extract_array( xx+yy*nxy+zz*nxyq , dset , 1 , tsar ) ;
219 for( tt=0 ; tt < ntr ; tt++ ) tsar[tt] += beta*rar[tt] ;
220 THD_insert_series( xx+yy*nxy+zz*nxyq , dset , ntr , MRI_float , tsar , 1 ) ;
221 }
222 }
223 }
224 }
225
226 free(bar) ;
227 }
228
229 /* add noise model */
230
231 fprintf(stderr,"++ adding ARMA(1,1) model") ;
232 kk = (nxy+nxy-2)/(2*ndiv) ; kk = MAX(kk,1) ; a0 = amn ; a1 = (amx-amn)/kk ;
233 kk = (nxy-1) /(2*ndiv) ; kk = MAX(kk,1) ; b0 = 0.5f*(bmn+bmx) ; b1 = (bmx-bmn)/(2.0f*kk) ;
234 kk = 1 + (nxy-1)/ndiv ; dsig = (sigmx-sigmn)/kk ; kold = -1 ;
235 for( zz=0 ; zz < nz ; zz++,kold=-1 ){
236 fprintf(stderr,".") ;
237 for( yy=0 ; yy < nxy ; yy++ ){
238 kk = yy/ndiv ; val = (1+kk)*dsig + sigmn ; if( val <= 0.0 ) val = 0.01 ;
239 for( xx=0 ; xx < nxy ; xx++ ){
240 dvx = (xx+yy)/(2*ndiv) ; dvy = (xx-yy)/(2*ndiv) ;
241 aa = a0+a1*dvx ; bb = b0+b1*dvy ;
242 if( aa < -0.888 ) aa = -0.888 ; else if( aa > 0.888 ) aa = 0.888 ;
243 if( bb < -0.888 ) bb = -0.888 ; else if( bb > 0.888 ) bb = 0.888 ;
244 lam = LAMBDA(aa,bb) ;
245 armim = mri_genARMA11( ntr , 1 , aa,lam , val ) ;
246 if( armim == NULL ){
247 WARNING_message("ARMA11 fails: xx=%d yy=%d a=%g b=%g lam=%g sig=%g",
248 xx,yy,aa,bb,lam,val);
249 continue ;
250 }
251 armar = MRI_FLOAT_PTR(armim) ;
252 THD_extract_array( xx+yy*nxy+zz*nxyq , dset , 1 , tsar ) ;
253 for( tt=0 ; tt < ntr ; tt++ ) tsar[tt] += armar[tt] ;
254 THD_insert_series( xx+yy*nxy+zz*nxyq , dset , ntr , MRI_float , tsar , 1 ) ;
255 mri_free(armim) ;
256 }
257 }
258 }
259 fprintf(stderr,"\n") ;
260
261 /* write output */
262
263 putenv("AFNI_DECONFLICT=OVERWRITE") ;
264 DSET_write(dset) ; WROTE_DSET(dset) ; exit(0) ;
265 }
266