There are a few key stages required to make this notebook work. Firstly the iCOM stream from the Linac needs to be recorded to disk, this iCOM stream then needs to be archived and grouped by patient.
Once these iCOM streams are running the plan needs to be delivered on the machine being listened to. From this point the iCOM stream is then able to be compared directly to the Monaco tel file.
This comparison is done by calculating an MU Density upon which the difference is reported on using the Gamma comparison tool.
A Windows service is created to listen to the iCOM stream. This service is made from a .bat file, an example of which can be seen at:
https://github.com/CCA-Physics/physics-server/blob/8a0954e5/RCCC/icom/harry_listening.bat
To create this service the nssm tool is used, with an example of its usage available at:
https://github.com/CCA-Physics/physics-server/blob/8a0954e5/RCCC/icom/services-install.bat
Warning: Take Note
A force closing of the service which is listening to the iCOM stream done in such a way that it is not able to properly close down the listening socket will cause the Linac being listened to raise an interlock which won't interupt the beam, but it will not let a new beam be delivered until the machine is logged out of and re-logged back in.
This can happen when the service is force killed, the host machine is force shutdown, or the network connection is abruptly disconnected. Normal shutdown of the service or machine should not have this effect.
For the iCOM stream to be easily indexed in the future deliveries are stored by patient id and name. This is done by setting up a .bat file to run on machine boot. An example .bat file that achieves this can be seen at:
https://github.com/CCA-Physics/physics-server/blob/8a0954e5/RCCC/icom/patient_archiving.bat
The resulting files are then loaded into a PyMedPhys Delivery object from which an MU Density can be calculated and used as a comparison and reporting tool.
These steps will be further expanded on below, prior to the lines of code that implement them.
In [1]:
import pathlib # for filepath path tooling
import lzma # to decompress the iCOM file
import numpy as np # for array tooling
import matplotlib.pyplot as plt # for plotting
In [2]:
# Makes it so that any changes in pymedphys is automatically
# propagated into the notebook without needing a kernel reset.
from IPython.lib.deepreload import reload
%load_ext autoreload
%autoreload 2
In [3]:
import pymedphys
In [4]:
patient_id = '015112'
In [5]:
icom_directory = pathlib.Path(r'\\physics-server\iComLogFiles\patients')
monaco_directory = pathlib.Path(r'\\monacoda\FocalData\RCCC\1~Clinical')
In [6]:
output_directory = pathlib.Path(r'S:\Physics\Patient Specific Logfile Fluence')
pdf_directory = pathlib.Path(r'P:\Scanned Documents\RT\PhysChecks\Logfile PDFs')
In [7]:
GRID = pymedphys.mudensity.grid()
COORDS = (GRID["jaw"], GRID["mlc"])
GAMMA_OPTIONS = {
'dose_percent_threshold': 2, # Not actually comparing dose though
'distance_mm_threshold': 0.5,
'local_gamma': True,
'quiet': True,
'max_gamma': 2,
}
In [8]:
all_tel_paths = list(monaco_directory.glob(f'*~{patient_id}/plan/*/tel.1'))
all_tel_paths
Out[8]:
In [9]:
plan_names_to_choose_from = [
path.parent.name for path in all_tel_paths
]
plan_names_to_choose_from
Out[9]:
In [10]:
icom_deliveries = list(icom_directory.glob(f'{patient_id}_*/*.xz'))
icom_deliveries
Out[10]:
In [11]:
icom_files_to_choose_from = [
path.stem for path in icom_deliveries
]
icom_files_to_choose_from
Out[11]:
In [12]:
monaco_plan_name = 'LeftIlium1' # plan directory name
icom_delivery = '20200213_133208' # iCOM timestamp
In [13]:
tel_path = list(monaco_directory.glob(f'*~{patient_id}/plan/{monaco_plan_name}/tel.1'))[-1]
tel_path
Out[13]:
In [14]:
icom_path = list(icom_directory.glob(f'{patient_id}_*/{icom_delivery}.xz'))[-1]
icom_path
Out[14]:
In [15]:
with lzma.open(icom_path, 'r') as f:
icom_stream = f.read()
Within PyMedPhys there is a Delivery object. This object can be created from a range of sources such as RT Plan DICOM files, Mosaiq SQL queries, iCOM data, trf log files, as well as Monaco tel files.
From this Delivery object many tasks can be undergone. The available methods and attributes on the Delivery object are given below:
In [16]:
# Print out available methods and attributes on the Delivery object
[command for command in dir(pymedphys.Delivery) if not command.startswith('_')]
Out[16]:
In [17]:
delivery_icom = pymedphys.Delivery.from_icom(icom_stream)
delivery_tel = pymedphys.Delivery.from_monaco(tel_path)
In [18]:
mudensity_icom = delivery_icom.mudensity()
mudensity_tel = delivery_tel.mudensity()
In [19]:
def to_tuple(array):
return tuple(map(tuple, array))
gamma = pymedphys.gamma(
COORDS,
to_tuple(mudensity_tel),
COORDS,
to_tuple(mudensity_icom),
**GAMMA_OPTIONS
)
In [20]:
def plot_gamma_hist(gamma, percent, dist):
valid_gamma = gamma[~np.isnan(gamma)]
plt.hist(valid_gamma, 50, density=True)
pass_ratio = np.sum(valid_gamma <= 1) / len(valid_gamma)
plt.title(
"Local Gamma ({0}%/{1}mm) | Percent Pass: {2:.2f} % | Max Gamma: {3:.2f}".format(
percent, dist, pass_ratio * 100, np.max(valid_gamma)
)
)
In [21]:
def plot_and_save_results(
mudensity_tel,
mudensity_icom,
gamma,
png_filepath,
pdf_filepath,
header_text="",
footer_text="",
):
diff = mudensity_icom - mudensity_tel
largest_item = np.max(np.abs(diff))
widths = [1, 1]
heights = [0.3, 1, 1, 1, 0.1]
gs_kw = dict(width_ratios=widths, height_ratios=heights)
fig, axs = plt.subplots(5, 2, figsize=(10, 16), gridspec_kw=gs_kw)
gs = axs[0, 0].get_gridspec()
for ax in axs[0, 0:]:
ax.remove()
for ax in axs[1, 0:]:
ax.remove()
for ax in axs[4, 0:]:
ax.remove()
axheader = fig.add_subplot(gs[0, :])
axhist = fig.add_subplot(gs[1, :])
axfooter = fig.add_subplot(gs[4, :])
axheader.axis("off")
axfooter.axis("off")
axheader.text(0, 0, header_text, ha="left", wrap=True, fontsize=30)
axfooter.text(0, 1, footer_text, ha="left", va="top", wrap=True, fontsize=6)
plt.sca(axs[2, 0])
pymedphys.mudensity.display(GRID, mudensity_tel)
axs[2, 0].set_title("Monaco Plan MU Density")
plt.sca(axs[2, 1])
pymedphys.mudensity.display(GRID, mudensity_icom)
axs[2, 1].set_title("Recorded iCOM MU Density")
plt.sca(axs[3, 0])
pymedphys.mudensity.display(
GRID, diff, cmap="seismic", vmin=-largest_item, vmax=largest_item
)
plt.title("iCOM - Monaco")
plt.sca(axs[3, 1])
pymedphys.mudensity.display(GRID, gamma, cmap="coolwarm", vmin=0, vmax=2)
plt.title(
"Local Gamma | "
f"{GAMMA_OPTIONS['dose_percent_threshold']}%/"
f"{GAMMA_OPTIONS['distance_mm_threshold']}mm")
plt.sca(axhist)
plot_gamma_hist(
gamma,
GAMMA_OPTIONS['dose_percent_threshold'],
GAMMA_OPTIONS['distance_mm_threshold'])
return fig
In [22]:
results_dir = output_directory.joinpath(patient_id, tel_path.parent.name, icom_path.stem)
results_dir.mkdir(exist_ok=True, parents=True)
header_text = (
f"Patient ID: {patient_id}\n"
f"Plan Name: {tel_path.parent.name}\n"
)
footer_text = (
f"tel.1 file path: {str(tel_path)}\n"
f"icom file path: {str(icom_path)}\n"
f"results path: {str(results_dir)}"
)
png_filepath = str(results_dir.joinpath("result.png").resolve())
pdf_filepath = str(pdf_directory.joinpath(f"{patient_id}.pdf").resolve())
fig = plot_and_save_results(
mudensity_tel, mudensity_icom,
gamma, png_filepath, pdf_filepath,
header_text=header_text, footer_text=footer_text
)
fig.tight_layout()
plt.savefig(png_filepath, dpi=300)
plt.show()
To create a pdf, the just created png file can be converted to pdf. To do this the tool imagemagick needs to be installed on your system. If you install this now you will need to reset your Jupyter server in a new command prompt so that the magick command is available within your path.
In [23]:
!magick convert "{png_filepath}" "{pdf_filepath}"