00001 #include <odinseq/seqsim.h>
00002 #include <odindata/data.h>
00003
00004
00011 void inplane_reduction(Data<float,4>& map, int factor) {
00012 if(factor>1) {
00013 TinyVector<int,4> newshape=map.shape();
00014 newshape(3)/=factor;
00015 newshape(2)/=factor;
00016 map.congrid(newshape);
00017 }
00018 }
00019
00020
00021
00022 farray create_map(const Data<float,4>& inmap) {
00023 Range all=Range::all();
00024
00025 TinyVector<int,4> inshape=inmap.shape();
00026
00027 Data<float,5> outmap(inshape(0),1,inshape(1),inshape(2),inshape(3));
00028 outmap(all,0,all,all,all)=inmap;
00029
00030 return outmap;
00031 }
00032
00033
00034
00035
00036
00038
00039
00040 void usage(const STD_string& missing_arg="") {
00041 STD_cout << STD_endl;
00042 STD_cout << "gensample: Generates a virtual sample suitable for sequence simulation in ODIN" << STD_endl;
00043 STD_cout << STD_endl;
00044 STD_cout << "Usage and options:" << STD_endl;
00045 STD_cout << " Create a sample from external maps:" << STD_endl;
00046 STD_cout << " gensample -fr <FOVread> -fp <FOVphase> -fs <FOVslice> -s0 <spin-density> [-t1 <T1map>] [-t2 <T2map>] [-pm <ppmMap>] [-ir <inplane-reduction>] [-ss <single-slice number>] [-fr <comma-separated-list of frame-rates>] -o <Sample-file>" << STD_endl;
00047 STD_cout << " Create a point-spread function:" << STD_endl;
00048 STD_cout << " gensample -psf -t1 <point T1> -t2 <point T2> [-dc <DiffusionCoeff>] -o <Sample-file>" << STD_endl;
00049 STD_cout << " Create rectangle of uniform spin distribution about origin in all 3 spatial dimensions:" << STD_endl;
00050 STD_cout << " gensample -uni -n <size> -f <FOV> [-t1 <uniform T1>] [-t2 <uniform T2>] -o <Sample-file>" << STD_endl;
00051 STD_cout << " Create an isochromat distribution:" << STD_endl;
00052 STD_cout << " gensample -iso -n <size> -df <freq-distribution-width[kHz]> -t1 <point T1> -t2 <point T2> [-dc <point Dcoeff>] -o <Sample-file>" << STD_endl;
00053 STD_cout << " Create a homogenous disk with given radius:" << STD_endl;
00054 STD_cout << " gensample -dsk -n <size> -f <FOV> -R <radius> [-t1 <uniform T1>] [-t2 <uniform T2>] [-dc <uniform DiffusionCoeff>] -o <Sample-file>" << STD_endl;
00055 STD_cout << " Create coaxial cylinder:" << STD_endl;
00056 STD_cout << " gensample -cyl -n <size> -f <FOV> -Xi <inner-suscept[ppm]> -Xo <outer-suscept[ppm]> -Ri <inner-radius> -Ro <outer-radius> [-t1 <uniform T1> -t2 <uniform T2> -sl <single-line>] -o <Sample-file>" << STD_endl;
00057 STD_cout << "Other options:" << STD_endl;
00058 STD_cout << " " << LogBase::get_usage() << STD_endl;
00059 STD_cout << " " << helpUsage() << STD_endl;
00060 if(missing_arg!="") STD_cout << "Missing argument: -" << missing_arg << STD_endl;
00061 }
00062
00063 int main(int argc, char* argv[]) {
00064 LogBase::set_log_levels(argc,argv);
00065 Log<Para> odinlog("gensample","main");
00066
00067
00068 if(hasHelpOption(argc,argv)) {usage(); return 0;}
00069
00070 Range all=Range::all();
00071 char optval[ODIN_MAXCHAR];
00072
00073
00074 STD_string phantom_fname;
00075 if(getCommandlineOption(argc,argv,"-o",optval,ODIN_MAXCHAR)) phantom_fname=optval; else {usage("o");exit(0);}
00076
00077
00078 Sample sample("Virtual Sample",false);
00079
00080 bool valid_mode=false;
00081
00083
00084 if(isCommandlineOption(argc,argv,"-psf")) {
00085 ODINLOG(odinlog,infoLog) << "Creating point spread function" << STD_endl;
00086
00087 float t1=0.0;
00088 float t2=0.0;
00089 float D=0.0;
00090
00091 if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval); else {usage("t1");exit(0);}
00092 if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval); else {usage("t2");exit(0);}
00093 if(getCommandlineOption(argc,argv,"-dc",optval,ODIN_MAXCHAR)) D=atof(optval);
00094
00095 sample.resize(1,1,1,1,1);
00096
00097 farray point(sample.get_extent());
00098
00099 point=1.0; sample.set_spinDensity(point);
00100 point=t1; sample.set_T1map(point);
00101 point=t2; sample.set_T2map(point);
00102 point=D; sample.set_DcoeffMap(point);
00103
00104 sample.set_FOV(10.0);
00105
00106
00107 valid_mode=true;
00108 }
00109
00111
00112 if(isCommandlineOption(argc,argv,"-iso")) {
00113 ODINLOG(odinlog,infoLog) << "Creating isochromat distribution" << STD_endl;
00114
00115 float t1,t2,df;
00116 int n;
00117 float D=0;
00118
00119 if(getCommandlineOption(argc,argv,"-df",optval,ODIN_MAXCHAR)) df=atof(optval); else {usage("df");exit(0);}
00120 if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval); else {usage("t1");exit(0);}
00121 if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval); else {usage("t2");exit(0);}
00122 if(getCommandlineOption(argc,argv,"-n" ,optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage("n");exit(0);}
00123 if(getCommandlineOption(argc,argv,"-dc",optval,ODIN_MAXCHAR)) D=atof(optval);
00124
00125 float total_width=4.0*df;
00126
00127 sample.resize(1,n,1,1,1);
00128 sample.set_freqrange(total_width);
00129
00130 farray map(sample.get_extent());
00131
00132 map=t1; sample.set_T1map(map);
00133 map=t2; sample.set_T2map(map);
00134 map=D; sample.set_DcoeffMap(map);
00135
00136 float sum=0.0;
00137
00138 for(int i=0; i<n; i++) {
00139 float f=(secureDivision(i,n-1)-0.5)*total_width;
00140 if(n==1) f=0.0;
00141 map[i]=exp(-2.0*f*f/df/df);
00142 sum+=map[i];
00143 }
00144
00145 sample.set_FOV(1.0e-6);
00146
00147 map/=sum;
00148 sample.set_spinDensity(map);
00149
00150 valid_mode=true;
00151 }
00152
00154
00155 if(isCommandlineOption(argc,argv,"-cyl")) {
00156 ODINLOG(odinlog,infoLog) << "Creating cylinder" << STD_endl;
00157
00158 float t1=0.0;
00159 float t2=0.0;
00160 float fov,Xi,Xo,Ri,Ro;
00161 int n;
00162 int sl=-1;
00163
00164 if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval);
00165 if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval);
00166 if(getCommandlineOption(argc,argv,"-sl",optval,ODIN_MAXCHAR)) sl=atoi(optval);
00167 if(getCommandlineOption(argc,argv,"-n" ,optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage("n");exit(0);}
00168 if(getCommandlineOption(argc,argv,"-f" ,optval,ODIN_MAXCHAR)) fov=atof(optval);else {usage("f");exit(0);}
00169 if(getCommandlineOption(argc,argv,"-Xi" ,optval,ODIN_MAXCHAR)) Xi=atof(optval); else {usage("Xi");exit(0);}
00170 if(getCommandlineOption(argc,argv,"-Xo" ,optval,ODIN_MAXCHAR)) Xo=atof(optval); else {usage("Xo");exit(0);}
00171 if(getCommandlineOption(argc,argv,"-Ri" ,optval,ODIN_MAXCHAR)) Ri=atof(optval); else {usage("Ri");exit(0);}
00172 if(getCommandlineOption(argc,argv,"-Ro" ,optval,ODIN_MAXCHAR)) Ro=atof(optval); else {usage("Ro");exit(0);}
00173
00174 int nx=n;
00175 int nz=n;
00176 if(sl>=0) nz=1;
00177 ODINLOG(odinlog,infoLog) << "nx/nz=" << nx << "/" << nz << STD_endl;
00178
00179 ODINLOG(odinlog,infoLog) << "Initializing sample ..." << STD_endl;
00180 sample.resize(1,1,nz,1,nx);
00181 sample.set_FOV(xAxis, fov);
00182 sample.set_FOV(yAxis, 1.0);
00183 sample.set_FOV(zAxis, fov);
00184 if(sl>=0) sample.set_FOV(zAxis, fov/n);
00185
00186 ODINLOG(odinlog,infoLog) << "Initializing maps ..." << STD_endl;
00187
00188 farray map(sample.get_extent());
00189 farray sdmap (sample.get_extent());
00190 sdmap=1.0;
00191
00192 ODINLOG(odinlog,infoLog) << "Setting T1/T2 maps " << STD_endl;
00193 if(t1>0.0) {map=t1; sample.set_T1map(map);}
00194 if(t2>0.0) {map=t2; sample.set_T2map(map);}
00195
00196
00197 ODINLOG(odinlog,infoLog) << "Calculating sample ..." << STD_endl;
00198 map=0.0;
00199 ndim ii(n_sampleDim);
00200 ii[frameDim]=0;
00201 ii[freqDim]=0;
00202 ii[yDim]=0;
00203
00204 for(int iz=0; iz<nz; iz++) {
00205 for(int ix=0; ix<nx; ix++) {
00206 float x=(double(ix)/double(n)-0.5)*fov;
00207 float z=(double(iz)/double(n)-0.5)*fov;
00208 if(sl>=0) z=(double(sl)/double(n)-0.5)*fov;
00209 float r=norm(x,z);
00210 ii[zDim]=iz;
00211 ii[xDim]=ix;
00212 if(r>=Ri && r<=Ro) {
00213 if(Xo==0.0 && Xi==0.0) map(ii)=0.0;
00214 else {
00215 float factor=0.5*(Xo-Xi)*Ri*Ri;
00216 map(ii)=1.0-Xo/6.0 + secureDivision(factor*(x*x-z*z),r*r*r*r);
00217 }
00218 } else {
00219 sdmap(ii)=0.0;
00220 }
00221 }
00222 }
00223
00224 ODINLOG(odinlog,infoLog) << "Setting spinDensity/ppmMap maps ..." << STD_endl;
00225 sample.set_spinDensity(sdmap);
00226 sample.set_ppmMap(map);
00227
00228 valid_mode=true;
00229 }
00230
00232
00233 if(isCommandlineOption(argc,argv,"-dsk")) {
00234 ODINLOG(odinlog,infoLog) << "Creating disk" << STD_endl;
00235
00236 float t1=0.0;
00237 float t2=0.0;
00238 float fov,R;
00239 int n;
00240 int nz=1;
00241 float D=0;
00242
00243 if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval);
00244 if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval);
00245 if(getCommandlineOption(argc,argv,"-dc",optval,ODIN_MAXCHAR)) D=atof(optval);
00246 if(getCommandlineOption(argc,argv,"-n" ,optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage("n");exit(0);}
00247 if(getCommandlineOption(argc,argv,"-f" ,optval,ODIN_MAXCHAR)) fov=atof(optval);else {usage("f");exit(0);}
00248 if(getCommandlineOption(argc,argv,"-R" ,optval,ODIN_MAXCHAR)) R=atof(optval); else {usage("R");exit(0);}
00249
00250
00251 ODINLOG(odinlog,infoLog) << "Initializing sample ..." << STD_endl;
00252 sample.resize(1,1,nz,n,n);
00253 sample.set_FOV(xAxis, fov);
00254 sample.set_FOV(yAxis, fov);
00255 sample.set_FOV(zAxis, nz);
00256
00257 ODINLOG(odinlog,infoLog) << "Initializing maps ..." << STD_endl;
00258 farray sdmap (sample.get_extent()); sdmap=0.0;
00259 farray Dmap (sample.get_extent()); Dmap=0.0;
00260
00261
00262 ODINLOG(odinlog,infoLog) << "Setting T1/T2 " << STD_endl;
00263 if(t1>0.0) {sample.set_T1(t1);}
00264 if(t2>0.0) {sample.set_T2(t2);}
00265
00266
00267 ODINLOG(odinlog,infoLog) << "Calculating sample ..." << STD_endl;
00268
00269 for(int iy=0; iy<n; iy++) {
00270 for(int ix=0; ix<n; ix++) {
00271 float x=(double(ix)/double(n)-0.5)*fov;
00272 float y=(double(iy)/double(n)-0.5)*fov;
00273 float r=norm(x,y);
00274 if(r<R) {
00275 for(int iz=0; iz<nz; iz++) {
00276 sdmap(0,0,iz,iy,ix)=1.0;
00277 Dmap(0,0,iz,iy,ix)=D;
00278 }
00279 }
00280 }
00281 }
00282
00283 ODINLOG(odinlog,infoLog) << "Setting spinDensity/DcoeffMap maps ..." << STD_endl;
00284 sample.set_spinDensity(sdmap);
00285 sample.set_DcoeffMap(Dmap);
00286
00287 valid_mode=true;
00288 }
00289
00291
00292 if(isCommandlineOption(argc,argv,"-uni")) {
00293 ODINLOG(odinlog,infoLog) << "Creating uniform sample" << STD_endl;
00294
00295 float t1=0.0;
00296 float t2=0.0;
00297 float fov;
00298 int n;
00299
00300 if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval);
00301 if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval);
00302 if(getCommandlineOption(argc,argv,"-n" ,optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage("n");exit(0);}
00303 if(getCommandlineOption(argc,argv,"-f" ,optval,ODIN_MAXCHAR)) fov=atof(optval);else {usage("f");exit(0);}
00304
00305 ODINLOG(odinlog,infoLog) << "Initializing sample ..." << STD_endl;
00306 sample.resize(1,1,n,n,n);
00307 sample.set_FOV(xAxis, fov);
00308 sample.set_FOV(yAxis, fov);
00309 sample.set_FOV(zAxis, fov);
00310
00311 ODINLOG(odinlog,infoLog) << "Initializing maps ..." << STD_endl;
00312 farray map(sample.get_extent());
00313
00314 ODINLOG(odinlog,infoLog) << "Setting T1/T2/PD maps " << STD_endl;
00315 if(t1>0.0) {map=t1; sample.set_T1map(map);}
00316 if(t2>0.0) {map=t2; sample.set_T2map(map);}
00317 if(t2>0.0) {map=1.0; sample.set_spinDensity(map);}
00318
00319 valid_mode=true;
00320 }
00321
00323
00324 if(!valid_mode) {
00325 ODINLOG(odinlog,infoLog) << "Combining arrays" << STD_endl;
00326
00327 int inplane_reduction_factor=1;
00328 int single_slice=-1;
00329 dvector frame_durations;
00330
00331 float FOVread;
00332 float FOVphase;
00333 float FOVslice;
00334 STD_string S0map_fname;
00335 STD_string T1map_fname;
00336 STD_string T2map_fname;
00337 STD_string ppmMap_fname;
00338
00339 if(getCommandlineOption(argc,argv,"-fr",optval,ODIN_MAXCHAR)) FOVread=atof(optval); else {usage("fr");exit(0);}
00340 if(getCommandlineOption(argc,argv,"-fp",optval,ODIN_MAXCHAR)) FOVphase=atof(optval); else {usage("fp");exit(0);}
00341 if(getCommandlineOption(argc,argv,"-fs",optval,ODIN_MAXCHAR)) FOVslice=atof(optval); else {usage("fs");exit(0);}
00342
00343 if(getCommandlineOption(argc,argv,"-s0",optval,ODIN_MAXCHAR)) S0map_fname=optval; else {usage("s0");exit(0);}
00344
00345
00346 if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) T1map_fname=optval;
00347 if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) T2map_fname=optval;
00348 if(getCommandlineOption(argc,argv,"-pm",optval,ODIN_MAXCHAR)) ppmMap_fname=optval;
00349
00350 if(getCommandlineOption(argc,argv,"-ir",optval,ODIN_MAXCHAR)) inplane_reduction_factor=atoi(optval);
00351 if(getCommandlineOption(argc,argv,"-ss",optval,ODIN_MAXCHAR)) single_slice=atoi(optval);
00352
00353
00354 int nframes_arg=0;
00355 if(getCommandlineOption(argc,argv,"-fr",optval,ODIN_MAXCHAR)) {
00356 svector toks=tokens(optval,',');
00357 nframes_arg=toks.size();
00358 ODINLOG(odinlog,infoLog) << "nframes_arg=" << nframes_arg << STD_endl;
00359 frame_durations.resize(nframes_arg);
00360 for(unsigned int i=0; i<nframes_arg; i++) {
00361 frame_durations[i]=atof(toks[i].c_str());
00362 }
00363 }
00364
00365 int nframes_max=0;
00366
00367 Data<float,4> filedata_s0;
00368 filedata_s0.autoread(S0map_fname);
00369 nframes_max=STD_max(nframes_max, filedata_s0.extent(0));
00370
00371 Data<float,4> filedata_t1;
00372 if(T1map_fname!="") filedata_t1.autoread(T1map_fname);
00373 nframes_max=STD_max(nframes_max, filedata_t1.extent(0));
00374
00375 Data<float,4> filedata_t2;
00376 if(T2map_fname!="") filedata_t2.autoread(T2map_fname);
00377 nframes_max=STD_max(nframes_max, filedata_t2.extent(0));
00378
00379 Data<float,4> filedata_ppmmap;
00380 if(ppmMap_fname!="") filedata_ppmmap.autoread(ppmMap_fname);
00381 nframes_max=STD_max(nframes_max, filedata_ppmmap.extent(0));
00382
00383 ODINLOG(odinlog,infoLog) << "nframes_max=" << nframes_max << STD_endl;
00384
00385 if(nframes_max==nframes_arg) {
00386 sample.set_frame_durations(frame_durations);
00387 } else {
00388 ODINLOG(odinlog,infoLog) << "nframes_max(" << nframes_max << ") != nframes_arg(" << nframes_arg << ")" << STD_endl;
00389 return -1;
00390 }
00391
00392
00393 if(single_slice>=0) {
00394 FOVslice/=filedata_s0.extent(1);
00395
00396 filedata_s0.reference(filedata_s0(all,Range(single_slice,single_slice),all,all));
00397
00398 if(filedata_t1.size()) {
00399 filedata_t1.reference(filedata_t1(all,Range(single_slice,single_slice),all,all));
00400 }
00401
00402 if(filedata_t2.size()) {
00403 filedata_t2.reference(filedata_t2(all,Range(single_slice,single_slice),all,all));
00404 }
00405
00406 if(filedata_ppmmap.size()) {
00407 filedata_ppmmap.reference(filedata_t2(all,Range(single_slice,single_slice),all,all));
00408 }
00409 }
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432 if(inplane_reduction_factor>1) {
00433
00434 inplane_reduction(filedata_s0,inplane_reduction_factor);
00435
00436 if(filedata_t1.size()) inplane_reduction(filedata_t1,inplane_reduction_factor);
00437 if(filedata_t2.size()) inplane_reduction(filedata_t2,inplane_reduction_factor);
00438 if(filedata_ppmmap.size()) inplane_reduction(filedata_ppmmap,inplane_reduction_factor);
00439 }
00440
00441
00442 TinyVector<int,4> s0shape=filedata_s0.shape();
00443
00444 sample.resize(s0shape(0),1,s0shape(1),s0shape(2),s0shape(3));
00445
00446 sample.set_FOV(xAxis, FOVread);
00447 sample.set_FOV(yAxis, FOVphase);
00448 sample.set_FOV(zAxis, FOVslice);
00449
00450 sample.set_spinDensity(create_map(filedata_s0));
00451 if(filedata_t1.size()) sample.set_T1map(create_map(filedata_t1));
00452 else sample.set_T1(1000.0);
00453 if(filedata_t2.size()) sample.set_T2map(create_map(filedata_t2));
00454 else sample.set_T2(100.0);
00455 if(filedata_ppmmap.size()) sample.set_ppmMap(create_map(filedata_ppmmap));
00456 }
00457
00458 ODINLOG(odinlog,infoLog) << "sample.get_extent()=" << STD_string(sample.get_extent()) << STD_endl;
00459
00460 sample.write(phantom_fname);
00461
00462 return 0;
00463 }
