import matplotlib.pyplot as plt

from pylinac.core.io import retrieve_demo_file
from pylinac.core.gamma import gamma_geometric
from pylinac.core.image import DicomImage

# get the files from the cloud
ref_file = retrieve_demo_file(name='gamma/HN_VMAT_Reference_1mmPx.dcm')
eval_file = retrieve_demo_file(name='gamma/HN_VMAT_Evaluated_1mmPx.dcm')
# load them as DICOM images; we want to use the raw pixels
ref_img = DicomImage(ref_file, raw_pixels=True)
eval_img = DicomImage(eval_file, raw_pixels=True)
# convert the images to float arrays; uints can result in overflow errors
ref_array = ref_img.array.astype(float)
eval_array = eval_img.array.astype(float)
# take a random sample profile through the middle of the images
ref_prof = ref_img[:, 90]
eval_prof = eval_img[:, 90]

# compute the gamma
gamma_map = gamma_geometric(ref_prof, eval_prof, dose_threshold=0)

# plot the results
fig, ax = plt.subplots(figsize=(10, 7))
ax.plot(ref_prof, 'k+', label='Original Reference')
ax.plot(eval_prof, 'bo', label='Evaluation')
ax.set_ylabel('Pixel Value')
ax.set_xlabel('Pixel Position')
ax.legend(loc='upper left')
g_ax = ax.twinx()
g_ax.plot(gamma_map, 'rx', label='Gamma')
g_ax.legend(loc='upper right')
g_ax.set_ylabel('Gamma')
fig.suptitle(f'1D Gamma Analysis; max \N{Greek Small Letter Gamma}:{gamma_map.max():.2f}; avg \N{Greek Small Letter Gamma}: {gamma_map.mean():.2f}; pass rate: {100*gamma_map[gamma_map <= 1].size/gamma_map.size:.2f}%')
plt.show()