Hi Simon,
I am trying to use the pymedphys.gamma function on a numpy array. I have done this by following the steps used to get gamma from DICOM images. However, I have run into an issue when I do this. I am trying to generate gamma analysis maps from 3D numpy arrays, however, the array that is being pulled from the reference and evaluation is outputted in the wrong axis and in the wrong orientation. I was wondering how I would go about re-orienting the array so that the data can be outputted correctly. I have tried to fix the array myself but I seem to be getting nowhere. I have attached an image of how the data is currently being outputted, it looks like a dose profile on its side at the moment.
Please note, I also noticed that the axes_reference and axes_evaluation are ordered so that the axes are pulled in the order z,y,x whereas I ordered my arrays x,y,z but when I tried to change this in the code, the gamma analysis stopped printing as well?
Thanks,
Nevin
Attached below is my code
x = (np.linspace(-1.9875,1.9875,x_voxel))
y = (np.linspace(-1.9875,1.9875,y_voxel))
z = (np.linspace(0.0125,29.9875,z_voxel))
coords = (x,y,z)
axes_reference = coords
axes_evaluation = coords
dose_reference = reduced_dose
dose_evaluation = sens_combined
(z_ref, y_ref, x_ref) = axes_reference
(z_eval, y_eval, x_eval) = axes_evaluation
gamma_options = {
'dose_percent_threshold': 3,
'distance_mm_threshold': 2,
'lower_percent_dose_cutoff': 2.8,
'interp_fraction': 10, # Should be 10 or more for more accurate results
'max_gamma': 3,
'random_subset': None,
'local_gamma': True,
'ram_available': 2**29 # 1/2 GB
}
gamma = pymedphys.gamma(
axes_reference, dose_reference,
axes_evaluation, dose_evaluation,
**gamma_options)
valid_gamma = gamma[~np.isnan(gamma)]
num_bins = (
gamma_options['interp_fraction'] * gamma_options['max_gamma'])
bins = np.linspace(0, gamma_options['max_gamma'], num_bins + 1)
plt.figure(2)
plt.hist(valid_gamma, bins, density=True)
plt.xlim([0, gamma_options['max_gamma']])
pass_ratio = np.sum(valid_gamma <= 1) / len(valid_gamma)
plt.title(f"Local Gamma ({gamma_options['dose_percent_threshold']}%/{gamma_options['distance_mm_threshold']}mm) | Percent Pass: {pass_ratio*100:.2f} %")
plt.xlabel("Gamma")
plt.ylabel("Probability Density")
plt.show()
# plt.savefig('gamma_hist.png', dpi=300)
########### GAMMA AT VARIOUS DEPTHS ##############
max_ref_dose = np.max(dose_reference)
lower_dose_cutoff = gamma_options['lower_percent_dose_cutoff'] / 100 * max_ref_dose
relevant_slice = (
np.max(dose_reference, axis=(1, 0)) >
lower_dose_cutoff)
slice_start = np.max([
np.where(relevant_slice)[0][0],
0])
slice_end = np.min([
np.where(relevant_slice)[0][-1],
len(z_ref)])
z_vals = z_ref[slice(slice_start, slice_end, 5)]
eval_slices = [
dose_evaluation[np.where(z_i == z_eval)[0][0], :, :]
for z_i in z_vals
]
ref_slices = [
dose_reference[np.where(z_i == z_ref)[0][0], :, :]
for z_i in z_vals
]
gamma_slices = [
gamma[np.where(z_i == z_ref)[0][0], :, :]
for z_i in z_vals
]
diffs = [
eval_slice - ref_slice
for eval_slice, ref_slice
in zip(eval_slices, ref_slices)
]
max_diff = np.max(np.abs(diffs))
for i, (eval_slice, ref_slice, diff, gamma_slice) in enumerate(zip(eval_slices, ref_slices, diffs, gamma_slices)):
fig, ax = plt.subplots(figsize=(13,10), nrows=2, ncols=2)
c00 = ax[0,0].contourf(
x_eval, y_eval, eval_slice, 100,
vmin=0, vmax=max_ref_dose)
ax[0,0].set_title("Evaluation")
fig.colorbar(c00, ax=ax[0,0], label='Dose (Gy)')
ax[0,0].invert_yaxis()
ax[0,0].set_xlabel('x (mm)')
ax[0,0].set_ylabel('z (mm)')
c01 = ax[0,1].contourf(
x_ref, y_ref, ref_slice, 100,
vmin=0, vmax=max_ref_dose)
ax[0,1].set_title("Reference")
fig.colorbar(c01, ax=ax[0,1], label='Dose (Gy)')
ax[0,1].invert_yaxis()
ax[0,1].set_xlabel('x (mm)')
ax[0,1].set_ylabel('z (mm)')
c10 = ax[1,0].contourf(
x_ref, y_ref, diff, 100,
vmin=-max_diff, vmax=max_diff, cmap=plt.get_cmap('seismic'))
ax[1,0].set_title("Dose difference")
fig.colorbar(c10, ax=ax[1,0], label='[Dose Eval] - [Dose Ref] (Gy)')
ax[1,0].invert_yaxis()
ax[1,0].set_xlabel('x (mm)')
ax[1,0].set_ylabel('z (mm)')
c11 = ax[1,1].contourf(
x_ref, y_ref, gamma_slice, 100,
vmin=0, vmax=2, cmap=plt.get_cmap('coolwarm'))
ax[1,1].set_title(
f"Local Gamma ({gamma_options['dose_percent_threshold']} % / {gamma_options['distance_mm_threshold']} mm)")
fig.colorbar(c11, ax=ax[1,1], label='Gamma Value')
ax[1,1].invert_yaxis()
ax[1,1].set_xlabel('x (mm)')
ax[1,1].set_ylabel('z (mm)')
plt.show()
print("\n")
,
,
.