5 points with noise

phase-encoding-simple.py

@@ -27,6 +27,9 @@ voxel_T2s_decay = np.random.rand(voxel_num)*(T2s_scale-T2s_min) + T2s_min
 voxel_all = np.append(voxel_amplitudes,voxel_T2s_decay/T2s_scale)
 if voxel_num == 5:
     voxel_all = [0.822628,0.691376,0.282906,0.226013,0.90703,0.144985,0.328563,0.440353,0.662462,0.720518]
+    voxel_all = [0.592606,0.135168,0.365712,0.667536,0.437378,0.918822,0.943879,0.590338,0.685997,0.658839]
+#first estimate    0.551777 0.190833 0.271438 0.814036 0.347389 0.926153 0.908453 0.581414 0.666012 0.673226
+
 #voxel_amplitudes = [0.4, 0.8, 0]
 #voxel_amplitudes = [0.9, 0.092893218813452, 0.5]
 #voxel_amplitudes = [0.6, 0.517157287525381, 0.4]
@@ -111,8 +114,11 @@ for i in range(voxel_num):
     bounds.append(T2s_minmax)
 
 x0 = np.full(2*voxel_num,0.5)
+x0 = [0.551777,0.190833,0.271438,0.814036,0.347389,0.926153,0.908453,0.581414,0.666012,0.673226]
+
+# result = minimize(hyper_fit, x0, bounds = bounds, method='L-BFGS-B', options={'ftol': 1e-16,'gtol': 1e-16, 'disp': True})
+result = minimize(fitness, x0, bounds = bounds, method='L-BFGS-B', options={'ftol': 1e-16,'gtol': 1e-16, 'disp': True})
 
-result = minimize(hyper_fit, x0, bounds = bounds, method='L-BFGS-B', options={'ftol': 1e-16,'gtol': 1e-16, 'disp': True})
 # result = differential_evolution(hyper_fit, bounds, polish=True
 #                                     , maxiter = voxel_num*2*500
 #                                 , disp = True

src/generator-rf.cc

@@ -8,6 +8,12 @@ namespace mri {
   // ********************************************************************** //
   void GeneratorRF::EvaluateDemodulated(std::vector<double> &mag_sin,
                                         std::vector<double> &mag_cos) {
+    std::random_device rd;
+    std::mt19937_64 om;
+    om.seed(rd());
+    double lbound = 0.0, ubound = 1.0;
+    std::uniform_real_distribution<double> rand(lbound, ubound);
+
     mag_sin.resize(total_samples_);
     mag_cos.resize(total_samples_);
     for (unsigned long i = 0; i < total_samples_; ++i) {
@@ -23,6 +29,8 @@ namespace mri {
                      (T2s_decays_norm_[voxel]*T2s_scale_
                       )
                      )
+          + (rand(om)-0.5)*noise_amplitude_
+
           ;
         mag_cos[i] +=
           amplitudes_[voxel]
@@ -32,6 +40,8 @@ namespace mri {
                      (T2s_decays_norm_[voxel]*T2s_scale_
                       )
                      )
+          + (rand(om)-0.5)*noise_amplitude_
+          
           ;
        }
     }
@@ -55,8 +65,8 @@ namespace mri {
                       )
                      )
           // ) + ( 0.0 
-          //                                           + np.random.rand()/noise_ratio
-
+          //+ np.random.rand()/noise_ratio
+          
           ;
       }
     }
@@ -192,7 +202,9 @@ namespace mri {
     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;
+      //phases_[i] = phase_init_ + i*phase_step;
+      phases_[i] = phase_init_ + std::sin(1.0*i/voxel_num_);
+      std::cout<<phases_[i]<<std::endl;
     }
   };
   // ********************************************************************** //

src/generator-rf.h

@@ -18,6 +18,7 @@ namespace mri {
     void SetAmplitudes(std::vector<double> amplitudes);
     void SetBase(const std::vector<double> &mag_sin,
                  const std::vector<double> &mag_cos);
+    void SetNoiseAmplitude(double noise){noise_amplitude_ = noise;};
     void SetT2sDecaysNorm(std::vector<double> T2s_decays_norm);
 
     // amplitudes -> values[:size/2]
@@ -40,7 +41,7 @@ namespace mri {
     double phase_range_ = PI_/2;
     double phase_init_ = 0.0;
     
-    double noise_ratio_ = 1e8;
+    double noise_amplitude_ = 1e-8;
 
     // B0=1.5T freq=64Mhz, period = 15.6 ns
     double period_ms_ = 15.6/1000/1000; //ms

src/main.cc

@@ -14,7 +14,8 @@
 
 const double PI_=3.14159265358979323846;
 const double eps_=1e-11;
-long total_generations_ = 5000;
+long total_generations_ = 15000;
+double noise_amp = 1e-8;
 
 // int voxel_num = 10;
 // double vox[] = {0.652282,0.286431,0.727236,0.271936,0.526648,0.47129,0.631995,0.808049,0.322672,0.0329494,0.283877,0.930383,0.0681125,0.269938,0.205617,0.810485,0.615148,0.717859,0.846448,0.264384};
@@ -22,16 +23,22 @@ long total_generations_ = 5000;
 int voxel_num = 5;
 double vox[] = {0.822628,0.691376,0.282906,0.226013,0.90703,0.144985,0.328563,0.440353,0.662462,0.720518};
 
+// int voxel_num = 4;
+// double vox[] = {0.822628,0.691376,0.282906,0.226013,0.90703,0.144985,0.328563,0.440353};
+
+// int voxel_num = 3;
+// //double vox[] = {0.822628,0.691376,0.282906,0.226013,0.90703,0.144985};
+// double vox[] = {0.822628,0.691376,0.282906,0.226013,0.30703,0.244985};
+
 double total_periods = 20;
 int rf_samples_per_period = 1;
 
 std::vector<double> init_vox (vox, vox + sizeof(vox) / sizeof(double) );
 
 //double phase_range = 0;
-double phase_range = PI_/2;
-double phase_init = 0.0;
+double phase_range = PI_/3.0;
+double phase_init = PI_/8.0;
 
-double noise_ratio = 1e8;
 
 
 // B0=1.5T freq=64Mhz, period = 15.6 ns
@@ -111,7 +118,9 @@ int main(int argc, char *argv[]) {
 // ********************************************************************** //
 void InitGenerators() {
     generator_rf_.InitVoxels(voxel_num, phase_init, phase_range, T2s_scale);
+    generator_rf_.SetNoiseAmplitude(noise_amp);
     fitness_rf_.InitVoxels(voxel_num, phase_init, phase_range, T2s_scale);
+    fitness_rf_.SetNoiseAmplitude(0.0);
     
     generator_rf_.InitTimer(total_periods, period_ms, rf_samples_per_period);
     fitness_rf_.InitTimer(total_periods, period_ms, rf_samples_per_period);