seqsim.h

/***************************************************************************
                          seqsim.h  -  description
                             -------------------
    begin                : Tue Jun 11 2002
    copyright            : (C) 2000-2021 by Thies H. Jochimsen
    email                : thies@jochimsen.de
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef SEQSIM_H
#define SEQSIM_H
 
#include <tjutils/tjthread.h>
#include <tjutils/tjnumeric.h> // for RandomDist
 
#include <odinpara/sample.h>
#include <odinseq/seqclass.h>
 
 
struct SeqSimInterval {
  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) {}
  float dt;
  STD_complex B1;
  float freq;
  float phase;
  float rec;
  float Gx;
  float Gy;
  float Gz;
};
 
 
class ProgressMeter; // forward declaration
 
 
class SeqSimAbstract : public virtual SeqClass {
 
 public:
 
  virtual ~SeqSimAbstract() {}
 
  virtual void prepare_simulation(const Sample& sample, CoilSensitivity* transmit_coil=0, CoilSensitivity* receive_coil=0, ProgressMeter* progmeter=0) = 0;
 
  virtual cvector simulate(const SeqSimInterval& simvals, double gamma) = 0;
 
  virtual void finalize_simulation() = 0;
 
};
 
 
 
class SeqSimMagsi : public LDRblock, public ThreadedLoop<SeqSimInterval,cvector,int>, public virtual SeqSimAbstract {
 
 
 public:
  SeqSimMagsi(const STD_string& label="unnamedSeqSimMagsi");
 
  SeqSimMagsi(const SeqSimMagsi& ssm);
 
  ~SeqSimMagsi();
 
  SeqSimMagsi& operator = (const SeqSimMagsi& ssm);
 
  unsigned int  get_total_size() const {return Mx.total();}
 
  SeqSimMagsi& reset_magnetization();
 
  SeqSimMagsi& set_initial_vector(float Mx, float My, float Mz);
 
  const farray& get_Mx() const {return Mx;}
 
  const farray& get_My() const {return My;}
 
  const farray& get_Mz() const {return Mz;}
 
  const farray& get_Mamp() const {return Mamp;}
 
  const farray& get_Mpha() const {return Mpha;}
 
  SeqSimMagsi& update();
 
 
  SeqSimMagsi& set_online_simulation(bool onlineflag) { online=onlineflag; return *this;}
 
 
  SeqSimMagsi& set_intravoxel_simulation(bool ivflag) { magsi=ivflag; return *this;}
 
  SeqSimMagsi& set_numof_threads(unsigned int n) { nthreads=n; return *this;}
 
  SeqSimMagsi& set_spat_rotmatrix(const RotMatrix& rotmatrix);
 
 
  bool do_simulation();
 
 
  // implementing virtual functions of SeqSimAbstract
  void prepare_simulation(const Sample& sample, CoilSensitivity* transmit_coil=0, CoilSensitivity* receive_coil=0, ProgressMeter* progmeter=0);
  cvector simulate(const SeqSimInterval& simvals, double gamma);
  void finalize_simulation();
 
  // implementing virtual functions of ThreadedLoop
  bool kernel(const SeqSimInterval& simvals, cvector& signal, int&, unsigned int begin, unsigned int end);
 
 private:
  friend class SeqTimecourse;
 
 
 
  SeqSimMagsi& resize(unsigned int xsize, unsigned int ysize, unsigned int zsize, unsigned int freqsize=1);
 
 
 
  void common_init();
 
  int append_all_members();
 
  SeqSimMagsi& MampMpha2MxMy();
  SeqSimMagsi& MxMy2MampMpha();
 
  void update_axes();
 
  void set_axes_cache(const Sample& sample);
 
  LDRfloatArr Mx;
  LDRfloatArr My;
  LDRfloatArr Mz;
  LDRfloatArr Mamp;
  LDRfloatArr Mpha;
 
  LDRbool   online;
  LDRaction update_now;
  LDRtriple initial_vector;
 
 
  bool iactive;
  bool magsi;
  unsigned int nthreads;
  RotMatrix* spat_rotmatrix;
 
 
  double gamma_cache;
 
  double elapsed_time; // within current time frame
  unsigned int time_index_cache;
  unsigned int numof_time_intervals_cache;
  double* time_intervals_cache;
 
 
  // cache for update_axes()
  float x_low;
  float x_upp;
  float y_low;
  float y_upp;
  float z_low;
  float z_upp;
  float freq_low; // in rad/s
  float freq_upp; // in rad/s
 
 
  // intra-voxel magn gradients
  float *dMx[4];
  float *dMy[4];
  float *dMz[4];
  float *dppm[3];  // readonly
 
 
 
  // use raw pointers to avoid slower []-operator of STD_vector
  unsigned int oneframe_size_cache; // size of one frame
  float* xpos_cache;
  float* ypos_cache;
  float* zpos_cache;
  float* freqoffset_cache; // in rad*kHz
  unsigned int nframes_ppm_cache;
  float* ppm_cache;
  unsigned int nframes_spin_density_cache;
  float* spin_density_cache;
  STD_complex* B1map_transm_cache;
  unsigned int num_rec_channel_cache;
  STD_complex** B1map_receiv_cache;
  unsigned int nframes_Dcoeff_cache;
  float* Dcoeff_cache;
  bool sim_diffusion;
  unsigned int nframes_r1_cache;
  float* r1_cache;
  unsigned int nframes_r2_cache;
  float* r2_cache;
  bool* has_relax_cache;
  float L[4];
  float B0_ppm;
  bool simcache_up2date;
  void outdate_simcache();
};
 
 
 
 
 
#ifdef STANDALONE_PLUGIN
 
class SeqSimMonteCarlo : public ThreadedLoop<SeqSimInterval,cvector,RandomDist>, public virtual SeqSimAbstract {
 
 public:
 
  SeqSimMonteCarlo(const STD_string& label="unnamedSeqSimMonteCarlo", unsigned int nparticles=10000, unsigned int nthreads=1);
 
  SeqSimMonteCarlo(const SeqSimMonteCarlo& ssmc) {common_init(); SeqSimMonteCarlo::operator = (ssmc);}
 
  SeqSimMonteCarlo& operator = (const SeqSimMonteCarlo& ssmc);
 
  
  farray get_spatial_dist() const;
 
 
  // implementing virtual functions of SeqSimAbstract
  void prepare_simulation(const Sample& sample, CoilSensitivity* transmit_coil=0, CoilSensitivity* receive_coil=0, ProgressMeter* progmeter=0);
  cvector simulate(const SeqSimInterval& simvals, double gamma);
  void finalize_simulation();
 
  // implementing virtual functions of ThreadedLoop
  bool kernel(const SeqSimInterval& simvals, cvector& signal, RandomDist& local_rng, unsigned int begin, unsigned int end);
 
 private:
  struct Particle {
    float pos[3];
    float Mx, My, Mz;
  };
 
  void common_init();
  void clear_cache();
 
  unsigned int linear_index(const float pos[3]) const;
 
  STD_vector<Particle> particle;
  unsigned int numof_threads;
 
  RandomDist rng; // seed only once per simulator
 
  double gamma_cache;
 
  unsigned int size_cache[3];
 
  // use raw pointers to avoid slower []-operator of STD_vector
  float* Dcoeff_cache;
  float* ppmMap_cache;
  float* R1map_cache;
  float* R2map_cache;
  float* spinDensity_cache;
 
  float pixelspacing_cache[3];
  float B0_ppm_cache;
 
 
};
 
#endif
 
 
 
 
 
 
 
#endif