simple generator

go.sh

@@ -84,4 +84,4 @@ MPI_options=
 #MPI_options="--bind-to-core"
 MPI_size=1
 cd $path_bin
-mpirun -np $MPI_size $MPI_options ./$mri_bin #mri.config
+time mpirun -np $MPI_size $MPI_options ./$mri_bin #mri.config

src/generator-rf.cc

@@ -1,21 +1,52 @@
 #include "./generator-rf.h"
 #include <random>
+#include <cmath>
 #include <stdexcept>
 namespace mri {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  void GeneratorRF::EvaluateSignal(std::vector<double> &signal_rf) {
+    signal_rf.resize(total_samples_);
+    for (long i = 0; i < total_samples_; ++i) {
+      signal_rf[i] = 0.0;
+      double time = timesteps_[i];
+      for (unsigned int voxel = 0; voxel < voxel_num_; ++voxel) {
+        signal_rf[i] +=
+          amplitudes_[voxel]
+          * std::sin(
+                     omega_*time + phases_[voxel]
+                     )
+          * std::exp(-time/
+                     (T2s_decays_norm_[voxel]*T2s_scale_
+                      )
+                     )
+          // ) + ( 0.0 
+          //                                           + np.random.rand()/noise_ratio
+
+          ;
+      }
+    }
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
   void GeneratorRF::RandomizeVoxels(){
     std::random_device rd;
-    std::mt19937_64 gn;
-    gn.seed(rd());
+    std::mt19937_64 om;
+    om.seed(rd());
     double lbound = 0.0, ubound = 1.0;
     std::uniform_real_distribution<double> rand(lbound, ubound);
-    for (int i = 0; i < 10; ++i;){
-      std::cout<< rand(gn) <<" ";
+    amplitudes_.clear();
+    T2s_decays_norm_.clear();
+    amplitudes_.resize(voxel_num_);
+    T2s_decays_norm_.resize(voxel_num_);
+    for (double &amp:amplitudes_){
+      amp = rand(om);
+    }
+    for (double &decay:T2s_decays_norm_){
+      decay = rand(om);
     }
-    std::cout<<std::endl;
-  
   };
 
   // ********************************************************************** //
@@ -29,8 +60,8 @@ namespace mri {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  void GeneratorRF::SetT2sDecays(std::vector<double> T2s_decays){
-    T2s_decays_ = T2s_decays;
+  void GeneratorRF::SetT2sDecaysNorm(std::vector<double> T2s_decays_norm){
+    T2s_decays_norm_ = T2s_decays_norm;
   };
   // ********************************************************************** //
   // ********************************************************************** //
@@ -42,10 +73,10 @@ namespace mri {
     period_ms_ = period_ms;
     samples_per_period_ = samples_per_period;
     omega_ = 2.0*PI_/period_ms_;
-    long total_samples = total_periods * samples_per_period;
-    timesteps_.resize(total_samples);
+    total_samples_ = total_periods * samples_per_period;
+    timesteps_.resize(total_samples_);
     double a_timestep = period_ms/samples_per_period;
-    for (long i = 0; i < total_samples; ++i) {
+    for (long i = 0; i < total_samples_; ++i) {
       timesteps_[i] = i*a_timestep;      
     }
     // for (double step: timesteps_) std::cout << step<<' ';
@@ -56,10 +87,18 @@ namespace mri {
   // ********************************************************************** //
   // ********************************************************************** //
   void GeneratorRF::SetVoxels(double voxel_num, double phase_init,
-                              double phase_range) {
+                              double phase_range, double T2s_scale) {
+    T2s_scale_ = T2s_scale;
+    if (voxel_num < 2)
+      throw std::invalid_argument("We need at least two voxels for phase gradient");    
     voxel_num_ = voxel_num;
     phase_init_ = phase_init;
     phase_range_ = phase_range;
+    double phase_step = phase_range_/(voxel_num_ - 1);
+    phases_.resize(voxel_num);
+    for (unsigned int i =0; i< voxel_num_; ++i) {
+      phases_[i] = phase_init_ + i*phase_step;
+    }
   };
   // ********************************************************************** //
   // ********************************************************************** //

src/generator-rf.h

@@ -6,14 +6,18 @@ namespace mri {
   class GeneratorRF {
   public:
     const double PI_=3.14159265358979323846;
+    void EvaluateSignal(std::vector<double> &signal_rf);
     void RandomizeVoxels();
     void SetAmplitudes(std::vector<double> amplitudes);
-    void SetT2sDecays(std::vector<double> T2s_decays);
+    void SetT2sDecaysNorm(std::vector<double> T2s_decays_norm);
     void SetTimer(long total_periods, double period_ms, int samples_per_period);
-    void SetVoxels(double voxel_num, double phase_init_, double phase_range_);
+    void SetVoxels(double voxel_num, double phase_init_, double phase_range_, double T2s_scale);
   private:
     std::vector<double> amplitudes_;
-    std::vector<double> T2s_decays_;
+    std::vector<double> phases_;
+    std::vector<double> T2s_decays_norm_;
+    double T2s_scale_ = 0.01; // #ms # need to be 10ms
+    /* double T2s_min_ = T2s_scale_/1000.0; */
     
     unsigned int voxel_num_ = 3;
     double phase_range_ = PI_/2;
@@ -26,6 +30,8 @@ namespace mri {
     double omega_ = 2.0*PI_/period_ms_; 
     long total_periods_ = 1000;
     int samples_per_period_ = 10;
+    long total_samples_ = total_periods_ * samples_per_period_;
+
     std::vector<double> timesteps_;
 
   };

src/main.cc

@@ -23,14 +23,15 @@ double noise_ratio = 1e8;
 
 // B0=1.5T freq=64Mhz, period = 15.6 ns
 double period_ms = 15.6/1000/1000; //ms
-double total_periods = 1000;
-int rf_samples_per_period = 10;
+double total_periods = 10;
+int rf_samples_per_period = 20;
 // double period_ms = 10; //ms
 // double total_periods = 2;
 // int rf_samples_per_period = 3;
 
-double T2s_scale = 0.01; //ms # need to be 10ms
-double T2s_min = T2s_scale/1000.0; //
+//double T2s_scale = 0.01; //ms # need to be 10ms
+//double T2s_min = T2s_scale/1000.0; //
+double T2s_scale = period_ms*total_periods; //ms # need to be 10ms
 
   
 int main(int argc, char *argv[]) {
@@ -38,12 +39,18 @@ int main(int argc, char *argv[]) {
   int rank;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
   try {
+    std::cout<<"Runing .... "<<std::endl;
     jade::SubPopulation sub_population;  // Optimizer.
     mri::GeneratorRF generator_rf;
-    generator_rf.SetVoxels(voxel_num, phase_init, phase_range);
+    generator_rf.SetVoxels(voxel_num, phase_init, phase_range, T2s_scale);
     generator_rf.SetTimer(total_periods, period_ms, rf_samples_per_period);
-    std::cout<<"Run"<<std::endl;
-
+    generator_rf.RandomizeVoxels();
+    std::vector<double> signal_rf;
+    generator_rf.EvaluateSignal(signal_rf);
+    // for (auto value:signal_rf){
+    //   std::cout << value <<std::endl;
+    // }
+    
   } catch( const std::invalid_argument& ia ) {
     // Will catch if  multi_layer_mie_ fails or other errors.
     std::cerr << "Invalid argument: " << ia.what() << std::endl;