seqsim.h
1 /***************************************************************************
2  seqsim.h - description
3  -------------------
4  begin : Tue Jun 11 2002
5  copyright : (C) 2000-2015 by Thies H. Jochimsen
6  email : thies@jochimsen.de
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #ifndef SEQSIM_H
19 #define SEQSIM_H
20 
21 #include <tjutils/tjthread.h>
22 #include <tjutils/tjnumeric.h> // for RandomDist
23 
24 #include <odinpara/sample.h>
25 #include <odinseq/seqclass.h>
26 
27 
43  SeqSimInterval() : dt(0.0), B1(0.0), freq(0.0), phase(0.0), rec(0.0), Gx(0.0), Gy(0.0), Gz(0.0) {}
44  float dt;
45  STD_complex B1;
46  float freq;
47  float phase;
48  float rec;
49  float Gx;
50  float Gy;
51  float Gz;
52 };
53 
55 
56 class ProgressMeter; // forward declaration
57 
58 
63 class SeqSimAbstract : public virtual SeqClass {
64 
65  public:
66 
67  virtual ~SeqSimAbstract() {}
68 
76  virtual void prepare_simulation(const Sample& sample, CoilSensitivity* transmit_coil=0, CoilSensitivity* receive_coil=0, ProgressMeter* progmeter=0) = 0;
77 
84  virtual cvector simulate(const SeqSimInterval& simvals, double gamma) = 0;
85 
89  virtual void finalize_simulation() = 0;
90 
91 };
92 
93 
95 
113 class SeqSimMagsi : public LDRblock, public ThreadedLoop<SeqSimInterval,cvector,int>, public virtual SeqSimAbstract {
114 
115 
116  public:
120  SeqSimMagsi(const STD_string& label="unnamedSeqSimMagsi");
121 
125  SeqSimMagsi(const SeqSimMagsi& ssm);
126 
130  ~SeqSimMagsi();
131 
135  SeqSimMagsi& operator = (const SeqSimMagsi& ssm);
136 
140  unsigned int get_total_size() const {return Mx.total();}
141 
145  SeqSimMagsi& reset_magnetization();
146 
150  SeqSimMagsi& set_initial_vector(float Mx, float My, float Mz);
151 
155  const farray& get_Mx() const {return Mx;}
156 
160  const farray& get_My() const {return My;}
161 
165  const farray& get_Mz() const {return Mz;}
166 
170  const farray& get_Mamp() const {return Mamp;}
171 
175  const farray& get_Mpha() const {return Mpha;}
176 
180  SeqSimMagsi& update();
181 
182 
186  SeqSimMagsi& set_online_simulation(bool onlineflag) { online=onlineflag; return *this;}
187 
188 
192  SeqSimMagsi& set_intravoxel_simulation(bool ivflag) { magsi=ivflag; return *this;}
193 
197  SeqSimMagsi& set_numof_threads(unsigned int n) { nthreads=n; return *this;}
198 
203  SeqSimMagsi& set_spat_rotmatrix(const RotMatrix& rotmatrix);
204 
205 
210  bool do_simulation();
211 
212 
213  // implementing virtual functions of SeqSimAbstract
214  void prepare_simulation(const Sample& sample, CoilSensitivity* transmit_coil=0, CoilSensitivity* receive_coil=0, ProgressMeter* progmeter=0);
215  cvector simulate(const SeqSimInterval& simvals, double gamma);
216  void finalize_simulation();
217 
218  // implementing virtual functions of ThreadedLoop
219  bool kernel(const SeqSimInterval& simvals, cvector& signal, int&, unsigned int begin, unsigned int end);
220 
221  private:
222  friend class SeqTimecourse;
223 
224 
225 
229  SeqSimMagsi& resize(unsigned int xsize, unsigned int ysize, unsigned int zsize, unsigned int freqsize=1);
230 
231 
232 
233  void common_init();
234 
235  int append_all_members();
236 
237  SeqSimMagsi& MampMpha2MxMy();
238  SeqSimMagsi& MxMy2MampMpha();
239 
240  void update_axes();
241 
242  void set_axes_cache(const Sample& sample);
243 
244  LDRfloatArr Mx;
245  LDRfloatArr My;
246  LDRfloatArr Mz;
247  LDRfloatArr Mamp;
248  LDRfloatArr Mpha;
249 
250  LDRbool online;
251  LDRaction update_now;
252  LDRtriple initial_vector;
253 
254 
255  bool iactive;
256  bool magsi;
257  unsigned int nthreads;
258  RotMatrix* spat_rotmatrix;
259 
260 
261  double gamma_cache;
262 
263  double elapsed_time; // within current time frame
264  unsigned int time_index_cache;
265  unsigned int numof_time_intervals_cache;
266  double* time_intervals_cache;
267 
268 
269  // cache for update_axes()
270  float x_low;
271  float x_upp;
272  float y_low;
273  float y_upp;
274  float z_low;
275  float z_upp;
276  float freq_low; // in rad/s
277  float freq_upp; // in rad/s
278 
279 
280  // intra-voxel magn gradients
281  float *dMx[4];
282  float *dMy[4];
283  float *dMz[4];
284  float *dppm[3]; // readonly
285 
286 
287 
288  // use raw pointers to avoid slower []-operator of STD_vector
289  unsigned int oneframe_size_cache; // size of one frame
290  float* xpos_cache;
291  float* ypos_cache;
292  float* zpos_cache;
293  float* freqoffset_cache; // in rad*kHz
294  unsigned int nframes_ppm_cache;
295  float* ppm_cache;
296  unsigned int nframes_spin_density_cache;
297  float* spin_density_cache;
298  STD_complex* B1map_transm_cache;
299  unsigned int num_rec_channel_cache;
300  STD_complex** B1map_receiv_cache;
301  unsigned int nframes_Dcoeff_cache;
302  float* Dcoeff_cache;
303  bool sim_diffusion;
304  unsigned int nframes_r1_cache;
305  float* r1_cache;
306  unsigned int nframes_r2_cache;
307  float* r2_cache;
308  bool* has_relax_cache;
309  float L[4];
310  float B0_ppm;
311  bool simcache_up2date;
312  void outdate_simcache();
313 };
314 
315 
316 
318 
319 
320 #ifdef STANDALONE_PLUGIN // exclude from Siemens DLLs
321 
327 class SeqSimMonteCarlo : public ThreadedLoop<SeqSimInterval,cvector,RandomDist>, public virtual SeqSimAbstract {
328 
329  public:
330 
334  SeqSimMonteCarlo(const STD_string& label="unnamedSeqSimMonteCarlo", unsigned int nparticles=10000, unsigned int nthreads=1);
335 
339  SeqSimMonteCarlo(const SeqSimMonteCarlo& ssmc) {common_init(); SeqSimMonteCarlo::operator = (ssmc);}
340 
344  SeqSimMonteCarlo& operator = (const SeqSimMonteCarlo& ssmc);
345 
346 
350  farray get_spatial_dist() const;
351 
352 
353  // implementing virtual functions of SeqSimAbstract
354  void prepare_simulation(const Sample& sample, CoilSensitivity* transmit_coil=0, CoilSensitivity* receive_coil=0, ProgressMeter* progmeter=0);
355  cvector simulate(const SeqSimInterval& simvals, double gamma);
356  void finalize_simulation();
357 
358  // implementing virtual functions of ThreadedLoop
359  bool kernel(const SeqSimInterval& simvals, cvector& signal, RandomDist& local_rng, unsigned int begin, unsigned int end);
360 
361  private:
362  struct Particle {
363  float pos[3];
364  float Mx, My, Mz;
365  };
366 
367  void common_init();
368  void clear_cache();
369 
370  unsigned int linear_index(const float pos[3]) const;
371 
372  STD_vector<Particle> particle;
373  unsigned int numof_threads;
374 
375  RandomDist rng; // seed only once per simulator
376 
377  double gamma_cache;
378 
379  unsigned int size_cache[3];
380 
381  // use raw pointers to avoid slower []-operator of STD_vector
382  float* Dcoeff_cache;
383  float* ppmMap_cache;
384  float* R1map_cache;
385  float* R2map_cache;
386  float* spinDensity_cache;
387 
388  float pixelspacing_cache[3];
389  float B0_ppm_cache;
390 
391 
392 };
393 
394 #endif
395 
396 
397 
398 
399 
401 
402 
406 #endif
MAGSI-based Magnetization Simulator.
Definition: seqsim.h:113
unsigned int get_total_size() const
Definition: seqsim.h:140
const farray & get_My() const
Definition: seqsim.h:160
Rotation Matrix.
Definition: geometry.h:67
SeqSimMagsi & set_intravoxel_simulation(bool ivflag)
Definition: seqsim.h:192
const farray & get_Mamp() const
Definition: seqsim.h:170
SeqSimMagsi & set_numof_threads(unsigned int n)
Definition: seqsim.h:197
SeqSimMagsi & set_online_simulation(bool onlineflag)
Definition: seqsim.h:186
Time interval for simulation.
Definition: seqsim.h:42
const farray & get_Mx() const
Definition: seqsim.h:155
const farray & get_Mz() const
Definition: seqsim.h:165
const farray & get_Mpha() const
Definition: seqsim.h:175
Virtual Sample for Simulation.
Definition: sample.h:55