libodin-dev 1.8.8-1 / usr / include / odinseq / seqsim.h

/***************************************************************************
                          seqsim.h  -  description
                             -------------------
    begin                : Tue Jun 11 2002
    copyright            : (C) 2000-2014 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>


/**
  * @ingroup odinseq
  * \brief Time interval for simulation
  *
  * Data structure to hold the values for a single interval with constant fields
  * - dt:    Duration of the interval
  * - B1:    Complex RF field
  * - freq:  Transmit/receive frequency
  * - phase: Transmit/receive phase in deg
  * - rec:   Receiver (>0 means on)
  * - Gx:    Gradient in read direction
  * - Gy:    Gradient in phase direction
  * - Gz:    Gradient in slice direction
  */
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


/**
  * @ingroup odinseq_internals
  * Interface for Simulators
  */
class SeqSimAbstract : public virtual SeqClass {

 public:

  virtual ~SeqSimAbstract() {}

/**
  * Prepare a simulation (i.e. before successive calls to simulate() ) with the parameters:
  * - sample:  The virtual sample
  * - transmit_coil: Transmitter coil, 0 for none
  * - receive_coil:  Receiver coil, 0 for none
  * - progmeter:     Status indicator to trace progress, 0 for none
  */
  virtual void prepare_simulation(const Sample& sample, CoilSensitivity* transmit_coil=0, CoilSensitivity* receive_coil=0, ProgressMeter* progmeter=0) = 0;

/**
  * Simulation with
  * - simvals: The magnetic fields during a time interval
  * - gamma:   Gyromagnetic ration of the nucleus observed
  * Return value: Signal in each receiver channel
  */
  virtual cvector simulate(const SeqSimInterval& simvals, double gamma) = 0;

/**
  * Call this function after a simulation (i.e. after successive calls to simulate() )
  */
  virtual void finalize_simulation() = 0;

};


/////////////////////////////////////////////////////////////////////

/**
  * @addtogroup odinseq
  * @{
  */

/**
  * \brief MAGSI-based Magnetization Simulator
  *
  * This is a simulator to calculate the time evolution of a magnetization grid
  * in 4 dimension (frequency and three spatial dimensions). This simulation of the
  * Bloch-Torrey equations is performed by means of the MAGSI algorith
  * (c.f. Journal of Magnetic Resonance 180:29-38, 2006).
  * Simulation usually involves the following steps:
  * - Initialize the simulator by prepare_simulation() using a virtual sample.
  * - Iterative simulation by simulate() using a structure with piece-wise constant fields.
  * - Finish simulation by calling finalize_simulation().
  */
class SeqSimMagsi : public JcampDxBlock, public ThreadedLoop<SeqSimInterval,cvector,int>, public virtual SeqSimAbstract {


 public:
/**
  * Constructs a simulator labeled 'object_label'.
  */
  SeqSimMagsi(const STD_string& label="unnamedSeqSimMagsi");

/**
  * Copy constructor
  */
  SeqSimMagsi(const SeqSimMagsi& ssm);

/**
  * Destructor
  */
  ~SeqSimMagsi();

/**
  * Assignment operator
  */
  SeqSimMagsi& operator = (const SeqSimMagsi& ssm);

/**
  * Returns the overall size of the array
  */
  unsigned int  get_total_size() const {return Mx.total();}

/**
  * Set each magnetization to initial state 'initial_vector', which is (0,0,1) by default
  */
  SeqSimMagsi& reset_magnetization();

/**
  * Set the vector for the initial magnetization
  */
  SeqSimMagsi& set_initial_vector(float Mx, float My, float Mz);

/**
  * Returns the real part of the transverse magnetisation
  */
  const farray& get_Mx() const {return Mx;}

/**
  * Returns the imaginary part of the transverse magnetisation
  */
  const farray& get_My() const {return My;}

/**
  * Returns the longitudinal magnetisation
  */
  const farray& get_Mz() const {return Mz;}

/**
  * Returns the amplitude of the transverse magnetisation
  */
  const farray& get_Mamp() const {return Mamp;}

/**
  * Returns the phase of the transverse magnetisation
  */
  const farray& get_Mpha() const {return Mpha;}

/**
  * Updates all parameter relations
  */
  SeqSimMagsi& update();


/**
  * Specifies whether simulation should be performed everytime update() is called
  */
  SeqSimMagsi& set_online_simulation(bool onlineflag) { online=onlineflag; return *this;}


/**
  * Specifies whether intra-voxel magnetzation gradients are considered during simulation
  */
  SeqSimMagsi& set_intravoxel_simulation(bool ivflag) { magsi=ivflag; return *this;}

/**
  * Specifies the number of threads used during simulation
  */
  SeqSimMagsi& set_numof_threads(unsigned int n) { nthreads=n; return *this;}

/**
  * Specifies a rotation matrix for the spatial domain, i.e. the magnetization
  * array will be rotated in space using the specified rotation matrix.
  */
  SeqSimMagsi& set_spat_rotmatrix(const RotMatrix& rotmatrix);


/**
  * Returns whether simulation should be performed, i.e. whether the 'online' flag
  * is true or 'update' was activated.
  */
  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;



/**
  * Resize the array in the four dimensions according to the given sizes
  */
  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

/**
  * \brief Monte-Carlo-based Magnetization Simulator
  *
  * Monte-Carlo Simulator for diffusional averaging
  */
class SeqSimMonteCarlo : public ThreadedLoop<SeqSimInterval,cvector,RandomDist>, public virtual SeqSimAbstract {

 public:

  /**
  * Constructs a simulator labeled 'object_label' to simulate 'nparticles' diffusion trajectories using 'nthreads' threads.
  */
  SeqSimMonteCarlo(const STD_string& label="unnamedSeqSimMonteCarlo", unsigned int nparticles=10000, unsigned int nthreads=1);

/**
  * Copy constructor
  */
  SeqSimMonteCarlo(const SeqSimMonteCarlo& ssmc) {common_init(); SeqSimMonteCarlo::operator = (ssmc);}

/**
  * Assignment operator
  */
  SeqSimMonteCarlo& operator = (const SeqSimMonteCarlo& ssmc);

  
/**
  * Get spatial distribution of particles after simulation
  */
  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