seqsim.h

/***************************************************************************
                          seqsim.h  -  description
                             -------------------
    begin                : Tue Jun 11 2002
    copyright            : (C) 2001 by Thies H. Jochimsen
    email                : jochimse@cns.mpg.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 JcampDxBlock, public ThreadedLoop<SeqSimInterval,cvector>, 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, 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);

  JDXfloatArr Mx;
  JDXfloatArr My;
  JDXfloatArr Mz;
  JDXfloatArr Mamp;
  JDXfloatArr Mpha;

  JDXbool   online;
  JDXaction update_now;
  JDXtriple 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 // exclude from Siemens DLLs

class SeqSimMonteCarlo : public virtual SeqSimAbstract {

 public:

  SeqSimMonteCarlo(const STD_string& label="unnamedSeqSimMonteCarlo", unsigned int nparticles=10000);

  SeqSimMonteCarlo(const SeqSimMonteCarlo& ssmc) {common_init(); SeqSimMonteCarlo::operator = (ssmc);}

  SeqSimMonteCarlo& operator = (const SeqSimMonteCarlo& ssmc);


  // 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();


 private:
  struct Particle {
    float xpos, ypos, zpos;
    float Mx, My, Mz;
  };

  void common_init();
  void clear_cache();

  unsigned int linear_index(float xpos, float ypos, float zpos) const;

  STD_vector<Particle> particle;

  RandomDist rng; // seed only once per simulator

  unsigned int xsize;
  unsigned int ysize;
  unsigned int zsize;

  // 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 xpixelspacing_cache;
  float ypixelspacing_cache;
  float zpixelspacing_cache;
  float B0_ppm_cache;


};

#endif







#endif

---