DICOM Files Made for Halcyon QA - Snooker Cue

Copyright 2019 Nick Harding of Queen's Centre for Oncology, Castle Hill Hospital.

Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License here.

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

The following cell imports the relevant libraries to run the code below.

In [1]:
import math
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pydicom
from pydicom.sequence import Sequence
from pydicom.dataset import Dataset
import os.path
import tkinter as tk
from tkinter import filedialog

%matplotlib inline

The original Snooker Cue test is described in Van Esch et al. This test is in regular use on our TrueBeams and iX machines for the quick routine assessment of various factors in RapidArc delivery. However, this test needs adapting for the Halcyon due to the faster gantry and MLC speeds. Offsetting the snooker cue by 5cm left and 10cm down as in the original test results in the MLCs not having to be driven at maximum speed in the fourth quadrant (270° - 180°). To created a modified plan for Halycon we need to understand the maths relating the MLC position to the gantry angle and snooker cue position.

The formula describing the central MLC position is given in the following code:

In [1]:
def get_centre(a=5.0,d=9.0,ga=0.0,iso=100.0):
    """Returns the centre of MLCs position for the snooker cue test
    
    Keyword arguments
    a is the offset left in cm (default 5.0)
    d is the offset down in cm (default 9.0)
    ga is the gantry angle in radians as per Halcyon IEC scale (default = 0.0)
    Iso is the isocentre distance in cm (Default 100.0cm)
    
    """
    
    h = math.sqrt(a**2 + d**2)
    delta = ((3/2 * np.pi)- np.arctan(d/a))
    
    
    return (iso * np.sin(delta - ga) * h)/(iso - (np.cos(delta-ga)*h))

The image below visualises the parameters described in the formula above:sc.png See Appendix A to visualise the MLC positions as a function of gantry angle

Using this formula it is possible to determine the MLC positions for given gantry angles in a segment. Eclipse allows RapidArc plans with two avoidance segments which means each arc segment can take three images of the snooker cue. As per the original paper the snooker cue is imaged in four separate quadrants (dose is delivered over 0.4°). The MLCs are held in the previous image position until the last point to make them drive at maximum speed. Each plan will hence need eight control points:

  1. Initial Position
  2. 0.4° further with increment in meterset (gives image)
  3. MLCs held until required GA to force max MLC speed
  4. Move MLCs to next segment position
  5. 0.4° further with increment in meterset (gives image)
  6. MLCs held until required GA to force max MLC speed
  7. Move MLCs to next segment position
  8. 0.4° further with increment in meterset (gives image)

With this information an existing Halcyon RapidArc plan can be modified. First, various parameters to be input into the DICOM plan need to be defined:

In [3]:
#this cell defines a function used below
def get_prev_ga(mlc_move, next_ga):
    """returns the gantry angle for max MLC speed
    
    mlc_move is the distance the mlcs are required to move in cm
    next_ga is the next gantry angle to start imaging at
    """
    
    abs_move = abs(mlc_move)
    t_move = abs_move / LS
    g_move = t_move * GS
    if DIRECTION == 'CC':
        return next_ga + g_move
    else:
        return next_ga - g_move
    
In [4]:
#this cell sets up various parameters to allow the dicom file to be created
#CONSTANTS
DIRECTION = 'CC' #'CC' = counter clockwise, 'CW' = clockwise 
LS = 5.0 # MLC max speed in cm/s
GS = 24.0 # max gantry speed in deg/s
DR = 740.0 # Dose rate in MU/min. Will need changing depending on cal conditions
NO_OF_CP = 8 # these plans need eight cps as described above
TOTAL_MU = 22 # number of MU for plan to deliver. 22 MU works for 740 DR. Suggest scaling this value with DR e.g.
#24 MU for 800 MU/min

#this list is three gantry angles in degress that will give images at positions where the images will be discrete
#the angles should also be chosen to cover as much of the quadrant as possible to ensure testing of the MLCs at the full
#range of gantry angles.
#The angles used here are for the example are for the quadrant 180° - 90°
#Update the angles as required.
#ensure they are entered correctly depending on whether CW or CC is chosen for direction above.
three_GA = [176.9,166.7,101.3]
three_MLC = []
#get the three MLC positions for required GA
for ga in three_GA:
    three_MLC.append(round(get_centre(ga = np.deg2rad(ga)),1))

#g_inc absolute value is half the value in degrees over which the SC will be imaged
#Typically this is 0.4° /2
#The sign depends on the arc direction 
g_inc = 0.0
if DIRECTION == 'CC':
    g_inc = -0.2
elif DIRECTION == 'CW':
    g_inc = 0.2
else:
    raise Exception(f"DIRECTION should be CW or CC. Value of DIRECTION was {DIRECTION}")
#make the eight required gantry angles for the eight control points
req_GA = [
    three_GA[0] - g_inc,
    three_GA[0] + g_inc,
    np.NaN,
    three_GA[1] - g_inc,
    three_GA[1] + g_inc,
    np.NaN,
    three_GA[2] - g_inc,
    three_GA[2] + g_inc
]
#the NaN values are placed as the gantry angle to hold needs calculating.
#This depends on knowing what the next list values are.
#We will calculate those values now:
req_GA[2] = get_prev_ga(three_MLC[1]-three_MLC[0],req_GA[3])
req_GA[5] = get_prev_ga(three_MLC[2]-three_MLC[1],req_GA[6])

#req_pos is an array of the required positions for the MLCs to be centered around (in cm)
req_pos = [three_MLC[0]] * 3 + [three_MLC[1]] * 3 + [three_MLC[2]] * 2

#make a bankA list to hold the position of the MLCs for the A bank. NB the DICOM file wants MLC positions
#in mm hence the * 10 conversion.
#the A bank will be offset by 5mm from the central position to give a 1cm wide slit
bankA = []
for pos in req_pos:
    bankA.append((10.0 * pos) - 5.0)

#list for the cumulative meterset of each CP
cp_cumulative_meterset =(
        [0.0000]
        + [round(1/3,4)] * 3
        + [round(2/3,4)] * 3
        + [1.0000]
        )

#This code has now made all the required parameters to create the SC DICOM plan.

We will now open an existing Halcyon DICOM plan to edit it and return a snooker cue plan. The Halcyon RP DICOM files have various hidden tags in them which the machine reads to determine the calibration conditions, hence the need to edit an existing plan rather than create a new DICOM RP plan from scratch. This means it is important to open a plan that is created with the right calibration for your machine. These plans should be shipped with your Halcyon. Whilst any Halcyon RA plan could be used for this process, the plan we used is the RA T3 Leaf Speed plan labelled RP.Phantom.T3_LS.dcm in our system. This plan has two fields in it. The code below will edit the RA beam and largely ignore the open beam. The code will need editing if the plan has different beams present.

In [5]:
#This function is required in the code below.
def writeMLC(distal,posA,posB):
    """returns a list of MLC positions
    
    if distal is True then a return a list of 28 leaves in posA
    and 28 leaves in posB
    if distal is False i.e. proximal return 29 leaves in each
    Pos A should be in mm
    Pos B should be in mm
    """
    if distal:
        return [posA]*28 + [posB]*28
    else:
        return [posA]*29 + [posB]*29
In [6]:
#the following two lines enable a GUI to open files with
root = tk.Tk()
root.withdraw()

#open a file dialog to open the required RP file to edit
ds =pydicom.dcmread(filedialog.askopenfilename())

#set MU for all beams in plan
for rbs in ds.FractionGroupSequence[0].ReferencedBeamSequence:
    rbs.BeamMeterset = str(TOTAL_MU)

#now work on RA plan:
ds.BeamSequence[1].NumberOfControlPoints = str(NO_OF_CP)
cp_sequence = Sequence()
ds.BeamSequence[1].ControlPointSequence = cp_sequence

#set up cp 0 to ensure banks are set in correct order and initial position correct
# Control Point Sequence: Control Point 0
cp0 = Dataset()
cp0.ControlPointIndex = "0"
cp0.NominalBeamEnergy = "6"
cp0.DoseRateSet = "740"

# Beam Limiting Device Position Sequence
beam_limiting_device_position_sequence = Sequence()
cp0.BeamLimitingDevicePositionSequence = beam_limiting_device_position_sequence

# Beam Limiting Device Position Sequence: Beam Limiting Device Position 1
beam_limiting_device_position1 = Dataset()
beam_limiting_device_position1.RTBeamLimitingDeviceType = 'X'
beam_limiting_device_position1.LeafJawPositions = ['-140', '140']
beam_limiting_device_position_sequence.append(beam_limiting_device_position1)

# Beam Limiting Device Position Sequence: Beam Limiting Device Position 2
beam_limiting_device_position2 = Dataset()
beam_limiting_device_position2.RTBeamLimitingDeviceType = 'Y'
beam_limiting_device_position2.LeafJawPositions = ['-140', '140']
beam_limiting_device_position_sequence.append(beam_limiting_device_position2)

# Beam Limiting Device Position Sequence: Beam Limiting Device Position 3
beam_limiting_device_position3 = Dataset()
beam_limiting_device_position3.RTBeamLimitingDeviceType = 'MLCX1'
beam_limiting_device_position3.LeafJawPositions = writeMLC(True,bankA[0],bankA[0]+10)
beam_limiting_device_position_sequence.append(beam_limiting_device_position3)

# Beam Limiting Device Position Sequence: Beam Limiting Device Position 4
beam_limiting_device_position4 = Dataset()
beam_limiting_device_position4.RTBeamLimitingDeviceType = 'MLCX2'
beam_limiting_device_position4.LeafJawPositions = writeMLC(False,bankA[0],bankA[0]+10)
beam_limiting_device_position_sequence.append(beam_limiting_device_position4)

cp0.GantryAngle = str(req_GA[0])
cp0.GantryRotationDirection = DIRECTION
cp0.BeamLimitingDeviceAngle = "0"
cp0.BeamLimitingDeviceRotationDirection = 'NONE'
cp0.PatientSupportAngle = "0"
cp0.PatientSupportRotationDirection = 'NONE'
cp0.TableTopEccentricAngle = "0"
cp0.TableTopEccentricRotationDirection = 'NONE'
cp0.TableTopVerticalPosition = ''
cp0.TableTopLongitudinalPosition = ''
cp0.TableTopLateralPosition = ''
cp0.IsocenterPosition = ['-5.3710938e-1', '-5.3710938e-1', '0']
cp0.CumulativeMetersetWeight = str(cp_cumulative_meterset[0])
cp0.TableTopPitchAngle = 0.0
cp0.TableTopPitchRotationDirection = 'NONE'
cp0.TableTopRollAngle = 0.0
cp0.TableTopRollRotationDirection = 'NONE'

# Referenced Dose Reference Sequence
refd_dose_ref_sequence = Sequence()
cp0.ReferencedDoseReferenceSequence = refd_dose_ref_sequence

# Referenced Dose Reference Sequence: Referenced Dose Reference 1
refd_dose_ref1 = Dataset()
refd_dose_ref1.CumulativeDoseReferenceCoefficient = "0"
refd_dose_ref1.ReferencedDoseReferenceNumber = "1"
refd_dose_ref_sequence.append(refd_dose_ref1)
cp_sequence.append(cp0)

#a lot of code to set up cp0. Some may be superfluous but it should help make sure the plan is deliverable.

#next code to set up cp index 1 - 7 inclusive
for i in range(1,len(bankA)):
    cp = Dataset()
    cp.ControlPointIndex = str(i)
    # Beam Limiting Device Position Sequence
    beam_limiting_device_position_sequence = Sequence()
    cp.BeamLimitingDevicePositionSequence = beam_limiting_device_position_sequence
    
    # Beam Limiting Device Position Sequence: Beam Limiting Device Position 1
    beam_limiting_device_position1 = Dataset()
    beam_limiting_device_position1.RTBeamLimitingDeviceType = 'MLCX1'
    beam_limiting_device_position1.LeafJawPositions = writeMLC(True,bankA[i],bankA[i]+10)
    beam_limiting_device_position_sequence.append(beam_limiting_device_position1)
    
    # Beam Limiting Device Position Sequence: Beam Limiting Device Position 2
    beam_limiting_device_position2 = Dataset()
    beam_limiting_device_position2.RTBeamLimitingDeviceType = 'MLCX2'
    beam_limiting_device_position2.LeafJawPositions = writeMLC(False,bankA[i],bankA[i]+10)

    beam_limiting_device_position_sequence.append(beam_limiting_device_position2)
    
    #gantry angle
    cp.GantryAngle = str(req_GA[i])
    #meterset
    cp.CumulativeMetersetWeight = str(cp_cumulative_meterset[i])
    cp_sequence.append(cp)

#can now save file.
#This will open a file directory window to choose where to save

output_dir = filedialog.askdirectory()
save_name = f"sc_{int(req_GA[0])}to{int(req_GA[7])}"
f_suffix = ".dcm"
output_name = os.path.join(output_dir,save_name,f_suffix)
output_name = f"{output_dir}/{save_name}.dcm"
ds.save_as(output_name)
print (f"File saved as {output_name}")
File saved as H:/Python/sc_177to101.dcm

Some of the code in the above cell was generated using the codify tool included with pydicom.

The created DICOM file can now either be annoymised and imported into Eclipse or alternatively delivered directly on the Halcyon in MachineQA mode. Note that the above code will only generate code for one quadrant. Each of the four quadrants will need creating separately.

An example screenshot can be seen below which is taken from Portal Dosimetry. This shows the test passing on the Halcyon demonstarting good interplay between MU, gantry angle and MLC movement. sc_pd.PNG

Conclusion

The Snooker Cue test has been adapted for the Halcyon and has been incorporated in the routine QA of this device.

Appendix A Visualisation of MLC Sinusoid

In [7]:
#Run this cell to visualise the sinusoid described by the MLCs
#x and y are empty lists for plotting
#x holds the gantry angle in deg
x=[]
#y will will hold the MLC central postion
y=[]

#make a numpy range running from pi (180) to -pi (-180) in 5deg intervals
ga_rad = np.arange(np.pi,-np.pi,-np.deg2rad(5.0))
with plt.style.context('seaborn-white'):
    for ga in ga_rad:
        centre = get_centre(ga=ga)
        ga_deg = np.rad2deg(ga)
        if ga_deg < 0:
            ga_deg += 360.0
        x.append(ga_deg)
        y.append(centre)
    sns.scatterplot(x,y)
    plt.xlim(0,360)
    plt.xlabel('Gantry Angle (°)')
    plt.ylabel('Centre position (cm)')
    plt.show()