diff --git a/config/config_default.ini b/config/config_default.ini index da8d0e1..9002d56 100644 --- a/config/config_default.ini +++ b/config/config_default.ini @@ -63,12 +63,12 @@ DEPHASING_T[1] = DEPHASING_T[2] = ; Apply gradient in T/m for each axis. Gradients are active for one DWELL_TIME GRADIENT_XYZ[0] = 0.0 0.0 0.0 -GRADIENT_XYZ[1] = 0.0 0.0 0.0 -GRADIENT_XYZ[2] = 0.0 0.0 0.0 +GRADIENT_XYZ[1] = +GRADIENT_XYZ[2] = ; Time (sec) to apply gradient. GRADIENT_T[0] = 0.0 -GRADIENT_T[1] = 0.0 -GRADIENT_T[2] = 0.0 +GRADIENT_T[1] = +GRADIENT_T[2] = DWELL_TIME = 50e-6 diff --git a/src/file_utils.h b/src/file_utils.h index e80ba33..26305cd 100644 --- a/src/file_utils.h +++ b/src/file_utils.h @@ -110,6 +110,9 @@ bool read_config(std::string config_filename, simulation_parameters& param, std: for (uint16_t i = 0; ifirst + "[" + std::to_string(i) + "]"); i++) { it->second[i] = ini.get("files").get(it->first + "[" + std::to_string(i) + "]"); + if (std::all_of(it->second[i].begin(), it->second[i].end(), isspace)) + break; // if it is empty, break the loop + std::filesystem::path f(it->second[i]); if (f.is_relative()) { // if it is relative, make it absolute @@ -131,7 +134,6 @@ bool read_config(std::string config_filename, simulation_parameters& param, std: } } - // ============== reading section SCAN_PARAMETERS ============== if(ini.has("SCAN_PARAMETERS")) { @@ -192,16 +194,22 @@ bool read_config(std::string config_filename, simulation_parameters& param, std: std::cout << ERR_MSG << "RF Times must be in ascending order, starts with 0 and must not have duplicates values" << std::endl; return false; } - + // ---------------- dephasing (start times, Flip angles ) ---------------- // Dephase start times - for(i=0; i> *(param.gradient_xyz+ind++) && ind < (i+1)*3); + } + catch(const std::exception& e){break;} + } if(i != param.n_gradient) { @@ -234,7 +252,7 @@ bool read_config(std::string config_filename, simulation_parameters& param, std: return false; } - // check Dephase start time conditions + // check Gradient start time conditions if (std::is_sorted(param.gradient_T, param.gradient_T + i) == false || std::adjacent_find(param.gradient_T, param.gradient_T + i) != param.gradient_T + i ) { diff --git a/src/kernels.cu b/src/kernels.cu index a77acbf..2008a04 100644 --- a/src/kernels.cu +++ b/src/kernels.cu @@ -14,16 +14,14 @@ #include "rotation.cuh" #include "helper_cuda.h" -#define GAMMA 267515315. // rad/s.T - //--------------------------------------------------------------------------------------------- // //--------------------------------------------------------------------------------------------- -__global__ void cu_sim(const simulation_parameters *param, const float *pFieldMap, const bool *pMask, const float *M0, const float *XYZ0, float *M1, float *XYZ1, uint32_t spin_no) +__global__ void cu_sim(const simulation_parameters *param, const float *pFieldMap, const bool *pMask, const float *M0, const float *XYZ0, float *M1, float *XYZ1) { - spin_no = blockIdx.x * blockDim.x + threadIdx.x ; + uint32_t spin_no = blockIdx.x * blockDim.x + threadIdx.x ; if (spin_no >= param->n_spins) return; @@ -69,8 +67,8 @@ __global__ void cu_sim(const simulation_parameters *param, const float *pFieldMa // ------ loop over timepoints ------ uint64_t ind=0, ind_old=param->matrix_length+1; uint32_t current_timepoint = 0, old_timepoint = 0; - uint16_t current_rf = 1, current_te = 0, counter_dephasing = 0; - float accumulated_phase = 0.f, dephase_deg = 0.f; + uint16_t current_rf = 1, current_te = 0, counter_dephasing = 0, counter_gradient = 0; + float accumulated_phase = 0.f; while (current_timepoint < param->n_timepoints) // param->n_timepoints is the total number of timepoints (= TR/dwelltime) { // ------ generate random walks and wrap around the boundries ------ @@ -100,18 +98,24 @@ __global__ void cu_sim(const simulation_parameters *param, const float *pFieldMa // ------ apply dephasing if there is any ------ if(counter_dephasing < param->n_dephasing && param->dephasing_T[counter_dephasing] == current_timepoint) { - float dephase_deg = (float)spin_no * param->dephasing[counter_dephasing] / (float)param->n_spins; // assign dephasing linearly to spins - accumulated_phase += dephase_deg / (param->B0 * GAMMA * param->dt * RAD2DEG); // scale dephasing to Tesla per dt + accumulated_phase += (float)spin_no * param->dephasing[counter_dephasing] / (float)param->n_spins; // assign dephasing linearly to spins counter_dephasing++; } + + // ------ apply gradient if there is any ------ + if(counter_gradient < param->n_gradient && param->gradient_T[counter_gradient] == current_timepoint) + { + const float *Gxyz = param->gradient_xyz + 3*counter_gradient; + accumulated_phase += (Gxyz[0]*xyz_new[0] + Gxyz[1]*xyz_new[1] + Gxyz[2]*xyz_new[2]) * param->dt*GAMMA*RAD2DEG; // Gx * x + Gy * y + Gz * z + counter_gradient++; + } // ------ apply other RF pulse if there is any ------ if(current_rf < param->n_RF && param->RF_ST[current_rf] == current_timepoint) { // dephase - dephase_deg = accumulated_phase * param->B0 * GAMMA * param->dt * RAD2DEG; // convert accumulated phase to degree - zrot(dephase_deg, m0, m1); + zrot(accumulated_phase, m0, m1); // relax time_elapsed = (current_timepoint - old_timepoint) * param->dt; relax(exp(-time_elapsed/param->T1), exp(-time_elapsed/param->T2), m1); @@ -126,8 +130,7 @@ __global__ void cu_sim(const simulation_parameters *param, const float *pFieldMa if (is_lastscan && current_te < param->n_TE && param->TE[current_te] == current_timepoint) { // dephase - dephase_deg = accumulated_phase * param->B0 * GAMMA * param->dt * RAD2DEG; // convert accumulated phase to degree - zrot(dephase_deg, m0, m1); + zrot(accumulated_phase, m0, m1); // relax time_elapsed = (current_timepoint - old_timepoint) * param->dt; relax(exp(-time_elapsed/param->T1), exp(-time_elapsed/param->T2), m1); @@ -146,8 +149,7 @@ __global__ void cu_sim(const simulation_parameters *param, const float *pFieldMa current_timepoint++; } // dephase - dephase_deg = accumulated_phase * param->B0 * GAMMA * param->dt * RAD2DEG; // convert accumulated phase to degree - zrot(dephase_deg, m0, m1); + zrot(accumulated_phase, m0, m1); // relax time_elapsed = (current_timepoint - old_timepoint) * param->dt; relax(exp(-time_elapsed/param->T1), exp(-time_elapsed/param->T2), m1); @@ -166,18 +168,27 @@ __global__ void cu_sim(const simulation_parameters *param, const float *pFieldMa // //--------------------------------------------------------------------------------------------- -__global__ void cu_scalePos(float *scaled_xyz, float *initial_xyz, float scale, uint32_t size) +__global__ void cu_scalePos(float *scaled_xyz, float *initial_xyz, float scale, uint64_t size) { - int n = blockIdx.x * blockDim.x + threadIdx.x ; + uint64_t n = blockIdx.x * blockDim.x + threadIdx.x ; if(n < size) { - uint32_t ind = 3*n; + uint64_t ind = 3*n; scaled_xyz[ind+0] = initial_xyz[ind+0] * scale; scaled_xyz[ind+1] = initial_xyz[ind+1] * scale; scaled_xyz[ind+2] = initial_xyz[ind+2] * scale; } } +//--------------------------------------------------------------------------------------------- +// CUDA kernel to perform array multiplication with a constant +//--------------------------------------------------------------------------------------------- +__global__ void cu_scaleArray(float *array, float scale, uint64_t size) +{ + uint64_t n = blockIdx.x * blockDim.x + threadIdx.x ; + if(n < size) + array[n] *= scale; +} //--------------------------------------------------------------------------------------------- // generate random initial position diff --git a/src/kernels.cuh b/src/kernels.cuh index 044242e..9b96d4d 100644 --- a/src/kernels.cuh +++ b/src/kernels.cuh @@ -16,10 +16,13 @@ #define GAMMA 267515315. // rad/s.T -__global__ void cu_sim(const simulation_parameters *param, const float *pFieldMap, const bool *pMask, const float *M0, const float *XYZ0, float *M1, float *XYZ1, uint32_t spin_no = 0); +__global__ void cu_sim(const simulation_parameters *param, const float *pFieldMap, const bool *pMask, const float *M0, const float *XYZ0, float *M1, float *XYZ1); // scale position to mimic the different volume size -__global__ void cu_scalePos(float *scaled_xyz, float *initial_xyz, float scale, uint32_t size); +__global__ void cu_scalePos(float *scaled_xyz, float *initial_xyz, float scale, uint64_t size); + +// CUDA kernel to perform array multiplication with a constant +__global__ void cu_scaleArray(float *array, float scale, uint64_t size); // generate random initial position __global__ void cu_randPosGen(float *spin_position_xyz, simulation_parameters *param, bool *pMask, uint32_t spin_no = 0); @@ -46,7 +49,7 @@ uint32_t is_masked(std::vector &XYZ0, std::vector &mask, simulation_par { uint64_t index = sub2ind(ROUND(XYZ0[3*i] * scale2grid[0]), ROUND(XYZ0[3*i+1] * scale2grid[1]), ROUND(XYZ0[3*i+2] * scale2grid[2]), param->fieldmap_size[0], param->fieldmap_size[1]); if (index >= mask.size()) - std::cout << ERR_MSG << "Index out of range: " << index << " >= " << mask.size() << std::endl; + std::cout << ERR_MSG << "For "<= " << mask.size() << std::endl; mask_counter[i] = mask[index]; } return std::accumulate(mask_counter.begin(), mask_counter.end(), 0); diff --git a/src/miscellaneous.h b/src/miscellaneous.h index c824a36..de2765b 100644 --- a/src/miscellaneous.h +++ b/src/miscellaneous.h @@ -77,8 +77,8 @@ typedef struct simulation_parameters std::cout<<"RF time = "; for(int i=0; i>>(d_pFieldMap[d], param.B0*GAMMA*param.dt*RAD2DEG, fieldmap.size()); if(hasXYZ0 == false) { // generate initial spatial position for spins, based on sample_length_ref printf("GPU %d) Generating random initial position for spins... ", d); @@ -132,7 +134,8 @@ bool simulate(simulation_parameters param, std::map