ODIN
gensample.cpp
1 #include <odinseq/seqsim.h>
2 #include <odindata/data.h>
3 
4 
11 void inplane_reduction(Data<float,4>& map, int factor) {
12  if(factor>1) {
13  TinyVector<int,4> newshape=map.shape();
14  newshape(3)/=factor;
15  newshape(2)/=factor;
16  map.congrid(newshape);
17  }
18 }
19 
20 
21 
22 farray create_map(const Data<float,4>& inmap) {
23  Range all=Range::all();
24 
25  TinyVector<int,4> inshape=inmap.shape();
26 
27  Data<float,5> outmap(inshape(0),1,inshape(1),inshape(2),inshape(3));
28  outmap(all,0,all,all,all)=inmap;
29 
30  return outmap;
31 }
32 
33 
34 
35 
36 
38 
39 
40 void usage(const STD_string& missing_arg="") {
41  STD_cout << STD_endl;
42  STD_cout << "gensample: Generates a virtual sample suitable for sequence simulation in ODIN" << STD_endl;
43  STD_cout << configInfo() << STD_endl;
44  STD_cout << STD_endl;
45  STD_cout << "Usage and options:" << STD_endl;
46  STD_cout << " Create a sample from external maps:" << STD_endl;
47  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>] [-ft <comma-separated-list of frame-rates for dymaic phantoms>] -o <Sample-file>" << STD_endl;
48  STD_cout << " Create a point-spread function:" << STD_endl;
49  STD_cout << " gensample -psf -t1 <point T1> -t2 <point T2> [-dc <DiffusionCoeff>] -o <Sample-file>" << STD_endl;
50  STD_cout << " Create rectangle of uniform spin distribution about origin in all 3 spatial dimensions:" << STD_endl;
51  STD_cout << " gensample -uni -n <size> -f <FOV> [-t1 <uniform T1>] [-t2 <uniform T2>] -o <Sample-file>" << STD_endl;
52  STD_cout << " Create an isochromat distribution:" << STD_endl;
53  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;
54  STD_cout << " Create a homogenous disk with given radius:" << STD_endl;
55  STD_cout << " gensample -dsk -n <size> -f <FOV> -R <radius> [-t1 <uniform T1>] [-t2 <uniform T2>] [-dc <uniform DiffusionCoeff>] -o <Sample-file>" << STD_endl;
56  STD_cout << " Create coaxial cylinder:" << STD_endl;
57  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;
58  STD_cout << "Other options:" << STD_endl;
59  STD_cout << " " << LogBase::get_usage() << STD_endl;
60  STD_cout << " " << helpUsage() << STD_endl;
61  if(missing_arg!="") STD_cout << "Missing argument: -" << missing_arg << STD_endl;
62 }
63 
64 int main(int argc, char* argv[]) {
65  LogBase::set_log_levels(argc,argv);
66  Log<Para> odinlog("gensample","main");
67 
68 
69  if(hasHelpOption(argc,argv)) {usage(); return 0;}
70 
71  Range all=Range::all();
72  char optval[ODIN_MAXCHAR];
73 
74 
75  STD_string phantom_fname;
76  if(getCommandlineOption(argc,argv,"-o",optval,ODIN_MAXCHAR)) phantom_fname=optval; else {usage("o");exit(0);}
77 
78 
79  Sample sample("Virtual Sample",false);
80 
81  bool valid_mode=false;
82 
84 
85  if(isCommandlineOption(argc,argv,"-psf")) {
86  ODINLOG(odinlog,infoLog) << "Creating point spread function" << STD_endl;
87 
88  float t1=0.0;
89  float t2=0.0;
90  float D=0.0;
91 
92  if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval); else {usage("t1");exit(0);}
93  if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval); else {usage("t2");exit(0);}
94  if(getCommandlineOption(argc,argv,"-dc",optval,ODIN_MAXCHAR)) D=atof(optval);
95 
96  sample.resize(1,1,1,1,1);
97 
98  farray point(sample.get_extent());
99 
100  point=1.0; sample.set_spinDensity(point);
101  point=t1; sample.set_T1map(point);
102  point=t2; sample.set_T2map(point);
103  point=D; sample.set_DcoeffMap(point);
104 
105  sample.set_FOV(10.0);
106 
107 
108  valid_mode=true;
109  }
110 
112 
113  if(isCommandlineOption(argc,argv,"-iso")) {
114  ODINLOG(odinlog,infoLog) << "Creating isochromat distribution" << STD_endl;
115 
116  float t1,t2,df;
117  int n;
118  float D=0;
119 
120  if(getCommandlineOption(argc,argv,"-df",optval,ODIN_MAXCHAR)) df=atof(optval); else {usage("df");exit(0);}
121  if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval); else {usage("t1");exit(0);}
122  if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval); else {usage("t2");exit(0);}
123  if(getCommandlineOption(argc,argv,"-n" ,optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage("n");exit(0);}
124  if(getCommandlineOption(argc,argv,"-dc",optval,ODIN_MAXCHAR)) D=atof(optval);
125 
126  float total_width=4.0*df;
127 
128  sample.resize(1,n,1,1,1);
129  sample.set_freqrange(total_width);
130 
131  farray map(sample.get_extent());
132 
133  map=t1; sample.set_T1map(map);
134  map=t2; sample.set_T2map(map);
135  map=D; sample.set_DcoeffMap(map);
136 
137  float sum=0.0;
138 
139  for(int i=0; i<n; i++) {
140  float f=(secureDivision(i,n-1)-0.5)*total_width;
141  if(n==1) f=0.0;
142  map[i]=exp(-2.0*f*f/df/df);
143  sum+=map[i];
144  }
145 
146  sample.set_FOV(1.0e-6);
147 
148  map/=sum;
149  sample.set_spinDensity(map);
150 
151  valid_mode=true;
152  }
153 
155 
156  if(isCommandlineOption(argc,argv,"-cyl")) {
157  ODINLOG(odinlog,infoLog) << "Creating cylinder" << STD_endl;
158 
159  float t1=0.0;
160  float t2=0.0;
161  float fov,Xi,Xo,Ri,Ro;
162  int n;
163  int sl=-1;
164 
165  if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval);
166  if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval);
167  if(getCommandlineOption(argc,argv,"-sl",optval,ODIN_MAXCHAR)) sl=atoi(optval);
168  if(getCommandlineOption(argc,argv,"-n" ,optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage("n");exit(0);}
169  if(getCommandlineOption(argc,argv,"-f" ,optval,ODIN_MAXCHAR)) fov=atof(optval);else {usage("f");exit(0);}
170  if(getCommandlineOption(argc,argv,"-Xi" ,optval,ODIN_MAXCHAR)) Xi=atof(optval); else {usage("Xi");exit(0);}
171  if(getCommandlineOption(argc,argv,"-Xo" ,optval,ODIN_MAXCHAR)) Xo=atof(optval); else {usage("Xo");exit(0);}
172  if(getCommandlineOption(argc,argv,"-Ri" ,optval,ODIN_MAXCHAR)) Ri=atof(optval); else {usage("Ri");exit(0);}
173  if(getCommandlineOption(argc,argv,"-Ro" ,optval,ODIN_MAXCHAR)) Ro=atof(optval); else {usage("Ro");exit(0);}
174 
175  int nx=n;
176  int nz=n;
177  if(sl>=0) nz=1;
178  ODINLOG(odinlog,infoLog) << "nx/nz=" << nx << "/" << nz << STD_endl;
179 
180  ODINLOG(odinlog,infoLog) << "Initializing sample ..." << STD_endl;
181  sample.resize(1,1,nz,1,nx); // x-z-plane
182  sample.set_FOV(xAxis, fov);
183  sample.set_FOV(yAxis, 1.0);
184  sample.set_FOV(zAxis, fov);
185  if(sl>=0) sample.set_FOV(zAxis, fov/n);
186 
187  ODINLOG(odinlog,infoLog) << "Initializing maps ..." << STD_endl;
188 
189  farray map(sample.get_extent());
190  farray sdmap (sample.get_extent());
191  sdmap=1.0;
192 
193  ODINLOG(odinlog,infoLog) << "Setting T1/T2 maps " << STD_endl;
194  if(t1>0.0) {map=t1; sample.set_T1map(map);}
195  if(t2>0.0) {map=t2; sample.set_T2map(map);}
196 
197 
198  ODINLOG(odinlog,infoLog) << "Calculating sample ..." << STD_endl;
199  map=0.0;
200  ndim ii(n_sampleDim);
201  ii[frameDim]=0;
202  ii[freqDim]=0;
203  ii[yDim]=0;
204 
205  for(int iz=0; iz<nz; iz++) {
206  for(int ix=0; ix<nx; ix++) {
207  float x=(double(ix)/double(n)-0.5)*fov;
208  float z=(double(iz)/double(n)-0.5)*fov;
209  if(sl>=0) z=(double(sl)/double(n)-0.5)*fov;
210  float r=norm(x,z);
211  ii[zDim]=iz;
212  ii[xDim]=ix;
213  if(r>=Ri && r<=Ro) {
214  if(Xo==0.0 && Xi==0.0) map(ii)=0.0;
215  else {
216  float factor=0.5*(Xo-Xi)*Ri*Ri;
217  map(ii)=1.0-Xo/6.0 + secureDivision(factor*(x*x-z*z),r*r*r*r);
218  }
219  } else {
220  sdmap(ii)=0.0;
221  }
222  }
223  }
224 
225  ODINLOG(odinlog,infoLog) << "Setting spinDensity/ppmMap maps ..." << STD_endl;
226  sample.set_spinDensity(sdmap);
227  sample.set_ppmMap(map);
228 
229  valid_mode=true;
230  }
231 
233 
234  if(isCommandlineOption(argc,argv,"-dsk")) {
235  ODINLOG(odinlog,infoLog) << "Creating disk" << STD_endl;
236 
237  float t1=0.0;
238  float t2=0.0;
239  float fov,R;
240  int n;
241  int nz=1; //20;
242  float D=0;
243 
244  if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval);
245  if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval);
246  if(getCommandlineOption(argc,argv,"-dc",optval,ODIN_MAXCHAR)) D=atof(optval);
247  if(getCommandlineOption(argc,argv,"-n" ,optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage("n");exit(0);}
248  if(getCommandlineOption(argc,argv,"-f" ,optval,ODIN_MAXCHAR)) fov=atof(optval);else {usage("f");exit(0);}
249  if(getCommandlineOption(argc,argv,"-R" ,optval,ODIN_MAXCHAR)) R=atof(optval); else {usage("R");exit(0);}
250 
251 
252  ODINLOG(odinlog,infoLog) << "Initializing sample ..." << STD_endl;
253  sample.resize(1,1,nz,n,n);
254  sample.set_FOV(xAxis, fov);
255  sample.set_FOV(yAxis, fov);
256  sample.set_FOV(zAxis, nz);
257 
258  ODINLOG(odinlog,infoLog) << "Initializing maps ..." << STD_endl;
259  farray sdmap (sample.get_extent()); sdmap=0.0;
260  farray Dmap (sample.get_extent()); Dmap=0.0;
261 
262 
263  ODINLOG(odinlog,infoLog) << "Setting T1/T2 " << STD_endl;
264  if(t1>0.0) {sample.set_T1(t1);}
265  if(t2>0.0) {sample.set_T2(t2);}
266 
267 
268  ODINLOG(odinlog,infoLog) << "Calculating sample ..." << STD_endl;
269 
270  for(int iy=0; iy<n; iy++) {
271  for(int ix=0; ix<n; ix++) {
272  float x=(double(ix)/double(n)-0.5)*fov;
273  float y=(double(iy)/double(n)-0.5)*fov;
274  float r=norm(x,y);
275  if(r<R) {
276  for(int iz=0; iz<nz; iz++) {
277  sdmap(0,0,iz,iy,ix)=1.0;
278  Dmap(0,0,iz,iy,ix)=D;
279  }
280  }
281  }
282  }
283 
284  ODINLOG(odinlog,infoLog) << "Setting spinDensity/DcoeffMap maps ..." << STD_endl;
285  sample.set_spinDensity(sdmap);
286  if(D>0.0) sample.set_DcoeffMap(Dmap);
287 
288  valid_mode=true;
289  }
290 
292 
293  if(isCommandlineOption(argc,argv,"-uni")) {
294  ODINLOG(odinlog,infoLog) << "Creating uniform sample" << STD_endl;
295 
296  float t1=0.0;
297  float t2=0.0;
298  float fov;
299  int n;
300 
301  if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) t1=atof(optval);
302  if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) t2=atof(optval);
303  if(getCommandlineOption(argc,argv,"-n" ,optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage("n");exit(0);}
304  if(getCommandlineOption(argc,argv,"-f" ,optval,ODIN_MAXCHAR)) fov=atof(optval);else {usage("f");exit(0);}
305 
306  ODINLOG(odinlog,infoLog) << "Initializing sample ..." << STD_endl;
307  sample.resize(1,1,n,n,n);
308  sample.set_FOV(xAxis, fov);
309  sample.set_FOV(yAxis, fov);
310  sample.set_FOV(zAxis, fov);
311 
312  ODINLOG(odinlog,infoLog) << "Initializing maps ..." << STD_endl;
313  farray map(sample.get_extent());
314 
315  ODINLOG(odinlog,infoLog) << "Setting T1/T2/PD maps " << STD_endl;
316  if(t1>0.0) {map=t1; sample.set_T1map(map);}
317  if(t2>0.0) {map=t2; sample.set_T2map(map);}
318  if(t2>0.0) {map=1.0; sample.set_spinDensity(map);}
319 
320  valid_mode=true;
321  }
322 
324 
325  if(!valid_mode) {
326  ODINLOG(odinlog,infoLog) << "Combining arrays" << STD_endl;
327 
328  int inplane_reduction_factor=1;
329  int single_slice=-1;
330  dvector frame_durations;
331 
332  float FOVread;
333  float FOVphase;
334  float FOVslice;
335  STD_string S0map_fname;
336  STD_string T1map_fname;
337  STD_string T2map_fname;
338  STD_string ppmMap_fname;
339 
340  if(getCommandlineOption(argc,argv,"-fr",optval,ODIN_MAXCHAR)) FOVread=atof(optval); else {usage("fr");exit(0);}
341  if(getCommandlineOption(argc,argv,"-fp",optval,ODIN_MAXCHAR)) FOVphase=atof(optval); else {usage("fp");exit(0);}
342  if(getCommandlineOption(argc,argv,"-fs",optval,ODIN_MAXCHAR)) FOVslice=atof(optval); else {usage("fs");exit(0);}
343 
344  if(getCommandlineOption(argc,argv,"-s0",optval,ODIN_MAXCHAR)) S0map_fname=optval; else {usage("s0");exit(0);}
345 
346 
347  if(getCommandlineOption(argc,argv,"-t1",optval,ODIN_MAXCHAR)) T1map_fname=optval;
348  if(getCommandlineOption(argc,argv,"-t2",optval,ODIN_MAXCHAR)) T2map_fname=optval;
349  if(getCommandlineOption(argc,argv,"-pm",optval,ODIN_MAXCHAR)) ppmMap_fname=optval;
350 
351  if(getCommandlineOption(argc,argv,"-ir",optval,ODIN_MAXCHAR)) inplane_reduction_factor=atoi(optval);
352  if(getCommandlineOption(argc,argv,"-ss",optval,ODIN_MAXCHAR)) single_slice=atoi(optval);
353 
354 
355  int nframes_arg=0;
356  if(getCommandlineOption(argc,argv,"-ft",optval,ODIN_MAXCHAR)) {
357  svector toks=tokens(optval,',');
358  nframes_arg=toks.size();
359  ODINLOG(odinlog,infoLog) << "nframes_arg=" << nframes_arg << STD_endl;
360  frame_durations.resize(nframes_arg);
361  for(int i=0; i<nframes_arg; i++) {
362  frame_durations[i]=atof(toks[i].c_str());
363  }
364  }
365 
366  int nframes_max=0;
367 
368  Data<float,4> filedata_s0;
369  filedata_s0.autoread(S0map_fname);
370  nframes_max=STD_max(nframes_max, filedata_s0.extent(0));
371 
372  Data<float,4> filedata_t1;
373  if(T1map_fname!="") filedata_t1.autoread(T1map_fname);
374  nframes_max=STD_max(nframes_max, filedata_t1.extent(0));
375 
376  Data<float,4> filedata_t2;
377  if(T2map_fname!="") filedata_t2.autoread(T2map_fname);
378  nframes_max=STD_max(nframes_max, filedata_t2.extent(0));
379 
380  Data<float,4> filedata_ppmmap;
381  if(ppmMap_fname!="") filedata_ppmmap.autoread(ppmMap_fname);
382  nframes_max=STD_max(nframes_max, filedata_ppmmap.extent(0));
383 
384  ODINLOG(odinlog,infoLog) << "nframes_max=" << nframes_max << STD_endl;
385 
386  if(nframes_max==nframes_arg) {
387  sample.set_frame_durations(frame_durations);
388  } else {
389  ODINLOG(odinlog,infoLog) << "nframes_max(" << nframes_max << ") != nframes_arg(" << nframes_arg << ")" << STD_endl;
390  return -1;
391  }
392 
393 
394  if(single_slice>=0) {
395  FOVslice/=filedata_s0.extent(1);
396 
397  filedata_s0.reference(filedata_s0(all,Range(single_slice,single_slice),all,all));
398 
399  if(filedata_t1.size()) {
400  filedata_t1.reference(filedata_t1(all,Range(single_slice,single_slice),all,all));
401  }
402 
403  if(filedata_t2.size()) {
404  filedata_t2.reference(filedata_t2(all,Range(single_slice,single_slice),all,all));
405  }
406 
407  if(filedata_ppmmap.size()) {
408  filedata_ppmmap.reference(filedata_t2(all,Range(single_slice,single_slice),all,all));
409  }
410  }
411 
412 
413 /*
414  for(int i=0; i<product(fileshape); i++) {
415  TinyVector<int,4> index=filedata_s0.create_index(i);
416 
417  bool valid_voxel=true;
418 
419  if(filedata_s0(index)<=0.0) valid_voxel=false;
420  if(filedata_t1(index)<=0.0) valid_voxel=false;
421  if(filedata_t2(index)<=0.0) valid_voxel=false;
422 
423  if(!valid_voxel) {
424  filedata_s0(index)=0.0;
425  filedata_t1(index)=0.0;
426  filedata_t2(index)=0.0;
427  filedata_ppmmap(index)=0.0;
428  }
429 
430  }
431 */
432 
433  if(inplane_reduction_factor>1) {
434 
435  inplane_reduction(filedata_s0,inplane_reduction_factor);
436 
437  if(filedata_t1.size()) inplane_reduction(filedata_t1,inplane_reduction_factor);
438  if(filedata_t2.size()) inplane_reduction(filedata_t2,inplane_reduction_factor);
439  if(filedata_ppmmap.size()) inplane_reduction(filedata_ppmmap,inplane_reduction_factor);
440  }
441 
442 
443  TinyVector<int,4> s0shape=filedata_s0.shape();
444 
445  sample.resize(s0shape(0),1,s0shape(1),s0shape(2),s0shape(3));
446 
447  sample.set_FOV(xAxis, FOVread);
448  sample.set_FOV(yAxis, FOVphase);
449  sample.set_FOV(zAxis, FOVslice);
450 
451  sample.set_spinDensity(create_map(filedata_s0));
452  if(filedata_t1.size()) sample.set_T1map(create_map(filedata_t1));
453  else sample.set_T1(1000.0);
454  if(filedata_t2.size()) sample.set_T2map(create_map(filedata_t2));
455  else sample.set_T2(100.0);
456  if(filedata_ppmmap.size()) sample.set_ppmMap(create_map(filedata_ppmmap));
457  }
458 
459  ODINLOG(odinlog,infoLog) << "sample.get_extent()=" << STD_string(sample.get_extent()) << STD_endl;
460 
461  sample.write(phantom_fname);
462 
463  return 0;
464 }
Data::autoread
int autoread(const STD_string &filename, const FileReadOpts &opts=FileReadOpts(), Protocol *prot=0, ProgressMeter *progmeter=0)
Definition: odindata/data.h:757
Data::congrid
void congrid(const TinyVector< int, N_rank > &newshape, const TinyVector< float, N_rank > *subpixel_shift=0, bool left_to_right=false)
Definition: odindata/data.h:884
LogBase::set_log_levels
static bool set_log_levels(int argc, char *argv[], bool trigger_error=true)
LogBase::get_usage
static STD_string get_usage()
Log
Definition: tjlog.h:218
Sample
Virtual Sample for Simulation.
Definition: sample.h:55
ndim
Definition: tjarray.h:41
tjvector::resize
virtual tjvector< T > & resize(unsigned int newsize)
configInfo
const char * configInfo()
norm
double norm(double x, double y)
isCommandlineOption
int isCommandlineOption(int argc, char *argv[], const char *option, bool modify=true)
secureDivision
double secureDivision(double numerator, double denominator)
hasHelpOption
int hasHelpOption(int argc, char *argv[])
getCommandlineOption
int getCommandlineOption(int argc, char *argv[], const char *option, char *returnvalue, int maxchar, bool modify=true)
helpUsage
const char * helpUsage()
tokens
svector tokens(const STD_string &tokenstring, char custom_separator=0, char escape_begin='"', char escape_end='"')
Generated on Fri Oct 15 2021 11:36:56 for ODIN by doxygen 1.9.1