@@ -94,7 +94,7 @@ void sim(const simulation_parameters *param, const float *pFieldMap, const uint8
94
94
// tissue type
95
95
uint8_t ts, ts_old;
96
96
auto indx = sub2ind (xyz1[0 ]*param->scale2grid [0 ], xyz1[1 ]*param->scale2grid [1 ], xyz1[2 ]*param->scale2grid [2 ], param->fieldmap_size [0 ], param->fieldmap_size [1 ], param->fieldmap_size [2 ]);
97
- ts_old = pMask[indx];
97
+ ts = ts_old = pMask[indx];
98
98
double diffusivity_scale = param->diffusivity [ts_old];
99
99
100
100
bool is_lastscan = false ;
@@ -122,10 +122,9 @@ void sim(const simulation_parameters *param, const float *pFieldMap, const uint8
122
122
while (current_timepoint < param->n_timepoints ) // param->n_timepoints is the total number of timepoints (= TR/dwelltime)
123
123
{
124
124
// ------ generate random walks and wrap around the boundries ------
125
- double rnd_wlk;
126
125
for (uint8_t i=0 ; i<3 ; i++)
127
126
{
128
- rnd_wlk = dist_random_walk_xyz (gen_r) * diffusivity_scale;
127
+ double rnd_wlk = dist_random_walk_xyz (gen_r) * diffusivity_scale;
129
128
xyz_new[i] = xyz_old[i] + rnd_wlk; // new spin position after random-walk
130
129
if (xyz_new[i] < 0 )
131
130
xyz_new[i] += (param->enCrossFOV ? param->fov [i] : 2 *std::abs (rnd_wlk)); // rnd_wlk is negative here
@@ -144,8 +143,9 @@ void sim(const simulation_parameters *param, const float *pFieldMap, const uint8
144
143
if (ind != ind_old) // fewer access to the global memory which is slow. Helpful for large samples!
145
144
{
146
145
// cross-tissue diffusion
147
- ts = pMask[ind];
146
+ ts = pMask[ind];
148
147
if (ts != ts_old)
148
+ {
149
149
if (dist_cross_tissue (gen_u) >= param->pXY [ts_old*param->n_tissue_type + ts])
150
150
{
151
151
if (itr++ > param->max_iterations )
@@ -155,16 +155,18 @@ void sim(const simulation_parameters *param, const float *pFieldMap, const uint8
155
155
}
156
156
continue ;
157
157
}
158
- else
159
- ts_old = ts;
160
- itr = 0 ;
161
- field = pFieldMap != nullptr ? pFieldMap[ind_old = ind]:0 .f ;
162
- ind = pMask[ind]; // the index of the tissue type
163
- T1 = param->T1_ms [ind] * 1e-3 ; // ms -> s
164
- T2 = param->T2_ms [ind] * 1e-3 ; // ms -> s
165
- diffusivity_scale = param->diffusivity [ind];
166
- }
158
+ ts_old = ts;
159
+ }
160
+
161
+ ind_old = ind;
162
+ field = pFieldMap != nullptr ? pFieldMap[ind]:0 .f ;
163
+ T1 = param->T1_ms [ts_old] * 1e-3 ; // ms -> s
164
+ T2 = param->T2_ms [ts_old] * 1e-3 ; // ms -> s
165
+ diffusivity_scale = param->diffusivity [ts_old];
166
+ }
167
+
167
168
accumulated_phase += field;
169
+ itr = 0 ;
168
170
169
171
// ------ apply ideal dephasing if there is any ------
170
172
if (counter_dephasing < param->n_dephasing && param->dephasing_us [counter_dephasing] == current_timepoint)
@@ -203,7 +205,7 @@ void sim(const simulation_parameters *param, const float *pFieldMap, const uint8
203
205
// save echo and copy m1 to m0 for the next iteration
204
206
for (uint32_t i=0 , shift=3 *param->n_TE *spin_no + 3 *current_te; i<3 ; i++)
205
207
M1[shift + i] = m0[i] = m1[i];
206
- T[spin_no*param->n_TE + current_te] = ts ;
208
+ T[spin_no*param->n_TE + current_te] = ts_old ;
207
209
208
210
accumulated_phase = 0 ; // reset phase since we have applied it in the previous step
209
211
old_timepoint = current_timepoint;
0 commit comments