small fov bug fixed

Author: aghaeifar

src/file_utils.cpp

@@ -327,7 +327,7 @@ bool file_utils::read_phantom(std::string phantom_filename, std::vector<float> &
     std::copy(fov.begin() , fov.end() , param->fov);
     std::copy(dims.begin(), dims.end(), param->fieldmap_size);
     BOOST_LOG_TRIVIAL(info) << "Size = " << dims[0] << " x " << dims[1] << " x " << dims[2] << std::endl;
-    BOOST_LOG_TRIVIAL(info) << "FoV = " << fov[0]*1e6 << " x " << fov[1]*1e6 << " x " << fov[2]*1e6 << " um^3" << std::endl;
+    BOOST_LOG_TRIVIAL(info) << "FoV = " << param->fov[0]*1e6 << " x " << param->fov[1]*1e6 << " x " << param->fov[2]*1e6 << " um^3" << std::endl;
 
     int n_tissue = find_max(mask) + 1;
     if (n_tissue > param->n_tissue_type)

src/kernels.cu

@@ -127,16 +127,16 @@ void sim(const simulation_parameters *param, const float *pFieldMap, const uint8
                 rnd_wlk = dist_random_walk_xyz(gen_r) * diffusivity_scale;
                 xyz_new[i] = xyz_old[i] + rnd_wlk; // new spin position after random-walk
                 if (xyz_new[i] < 0)
-                    xyz_new[i] += (param->enCrossFOV ? param->fov[i] : 2*ABS(rnd_wlk)); // rnd_wlk is negative here
+                    xyz_new[i] += (param->enCrossFOV ? param->fov[i] : 2*std::abs(rnd_wlk)); // rnd_wlk is negative here
                 else if (xyz_new[i] >= param->fov[i])
-                    xyz_new[i] -= (param->enCrossFOV ? param->fov[i] : 2*ABS(rnd_wlk)); // rnd_wlk is positive here
+                    xyz_new[i] -= (param->enCrossFOV ? param->fov[i] : 2*std::abs(rnd_wlk)); // rnd_wlk is positive here
             }
 
             // ------ subscripts to linear indices ------
             ind = sub2ind(xyz_new[0]*param->scale2grid[0], xyz_new[1]*param->scale2grid[1], xyz_new[2]*param->scale2grid[2], param->fieldmap_size[0], param->fieldmap_size[1], param->fieldmap_size[2]);
             if(ind >= param->matrix_length || ind < 0)
             {
-                printf("Error:spin=%d, ind=%" PRId64 ", %d,\n\tscale=(%.10f, %.10f, %.10f)\n\txyz=(%.10f, %.10f, %.10f)\n\tscale*xyz(%.10f, %.10f, %.10f)\n\tFoV=(%.10f, %.10f, %.10f)\n\t(%d, %d, %d)\n",spin_no, ind, current_timepoint, param->scale2grid[0], param->scale2grid[1], param->scale2grid[2], xyz_new[0], xyz_new[1], xyz_new[2], param->scale2grid[0]*xyz_new[0], param->scale2grid[1]*xyz_new[1], param->scale2grid[2]*xyz_new[2], param->fov[0], param->fov[1], param->fov[2], xyz_new[0] >= param->fov[0], xyz_new[1] >= param->fov[1], xyz_new[2] >= param->fov[2]);
+                printf("Error:spin=%d, ind=%" PRId64 ", %d,\n\tscale2grid=(%.10f, %.10f, %.10f)\n\txyz=(%.10f, %.10f, %.10f)\n\txyz_old=(%.10f, %.10f, %.10f)\n\tscale2grid*xyz(%.10f, %.10f, %.10f)\n\tPhantomSize=(%d, %d, %d)\n\tFoV=(%.10f, %.10f, %.10f)\n\tError?(%d, %d, %d)\n",spin_no, ind, current_timepoint, param->scale2grid[0], param->scale2grid[1], param->scale2grid[2], xyz_new[0], xyz_new[1], xyz_new[2], xyz_old[0], xyz_old[1], xyz_old[2], param->scale2grid[0]*xyz_new[0], param->scale2grid[1]*xyz_new[1], param->scale2grid[2]*xyz_new[2], param->fieldmap_size[0], param->fieldmap_size[1], param->fieldmap_size[2], param->fov[0], param->fov[1], param->fov[2], xyz_new[0] >= param->fov[0], xyz_new[1] >= param->fov[1], xyz_new[2] >= param->fov[2]);
                 return;
             }
             // ------ accumulate phase ------

src/simulation_parameters.h

@@ -102,7 +102,7 @@ typedef struct simulation_parameters
         ss<<"TE = "; for(int i=0; i<n_TE; i++) ss<<TE_us[i]*timestep_us/1000.<<' '; ss<<"ms.\n";
         ss<<"T2 = "; for(int i=0; i<n_tissue_type; i++) ss<<T2_ms[i]<<' '; ss<<"ms.\n";
         ss<<"T1 = "; for(int i=0; i<n_tissue_type; i++) ss<<T1_ms[i]<<' '; ss<<"ms.\n";
-        ss<<"Diffusivity = "; for(int i=0; i<n_tissue_type; i++) ss<<diffusivity[i]<<' '; ss<<"\n";
+        ss<<"sqrt(2*diffusivity*timestep) = "; for(int i=0; i<n_tissue_type; i++) ss<<diffusivity[i]<<' '; ss<<"\n";
         ss<<"Cross Tissue Probability =\n"; for(int i=0; i<n_tissue_type; i++) {for(int j=0; j<n_tissue_type; j++) ss<<pXY[j+i*n_tissue_type]<<' '; ss<<'\n';};
         ss<<"RF flip-angle   = "; for(int i=0; i<n_RF; i++) ss<<RF_FA_deg[i]<<' '; ss<<"deg.\n";
         ss<<"RF phase        = "; for(int i=0; i<n_RF; i++) ss<<RF_PH_deg[i]<<' '; ss<<"deg.\n";

src/spinwalk.cu

@@ -168,7 +168,6 @@ bool run(simulation_parameters param, std::map<std::string, std::vector<std::str
         // tqdm bar;
         using namespace indicators;
         ProgressBar bar{option::ShowPercentage{true}, option::Start{"["}, option::Fill{"="}, option::Lead{">"}, option::End{"]"}};
-        bar.set_progress(0);
         simulation_parameters param_local;
         memcpy(&param_local, &param, sizeof(simulation_parameters));
         std::vector<uint32_t> v(param_local.n_spins);
@@ -179,6 +178,12 @@ bool run(simulation_parameters param, std::map<std::string, std::vector<std::str
             {
                 param_local.fov[i] = fov_scale[sl] * param.fov[i];
                 param_local.scale2grid[i]    = param_local.fieldmap_size[i] / param_local.fov[i];
+                for (int n = 0; n < param.n_tissue_type; i++)
+                    if(param_local.fov[i] < 100 * param.diffusivity[n])
+                    {
+                        BOOST_LOG_TRIVIAL(error) << "FoV (= " << param_local.fov[i] << ") must be 100x larger than sqrt(2*diffusivity*timestep) (= " << param.diffusivity[n] << ")!";
+                        return false;
+                    }
             }           
 
             if(param.no_gpu)

src/version.h

@@ -6,7 +6,7 @@
 
 #define SPINWALK_VERSION_MAJOR 1
 #define SPINWALK_VERSION_MINOR 13
-#define SPINWALK_VERSION_PATCH 5
+#define SPINWALK_VERSION_PATCH 6
 
 //---------------------------------------------------------------------------------------------
 //