UNIVERSITY OF CAMBRIDGE
DOCTORAL THESIS
Incorporating Range Uncertainty into
Proton Therapy Treatment Planning
Author:
Stacey Elizabeth McGowan
Supervisors:
Prof. Neil G. Burnet
Dr Simon J. Thomas
Selwyn College
This dissertation is submitted for the degree of
Doctor of Philosophy
in the Department of Oncology
Tuesday 14th April, 2015

In memory of my Grandmother
Janet Jenkins
6th July 1923 - 17th August 2014

Declaration
This dissertation is the result of my own work and includes nothing which is the outcome
of work done in collaboration except as specified in the text. It is not substantially the same
as any that I have submitted, or, is being concurrently submitted for a degree or diploma
or other qualification at the University of Cambridge or any other University or similar
institution except as specified in the text. I further state that no substantial part of my
dissertation has already been submitted, or, is being concurrently submitted for any such
degree, diploma or other qualification at the University of Cambridge or any other University
or similar institution except as declared in specified in the text. It does not exceed the
prescribed word limit for the relevant Degree Committee.
Stacey Elizabeth McGowan
April 2015

Acknowledgements
I would like to acknowledge and thank both my supervisors, Professor Neil Burnet and Dr
Simon Thomas, for all their guidance and support. I would also like to acknowledge the
support provided by Dr Hayley Woffendin.
I would like to thank everyone in Medical Physics, all the staff in the Radiotherapy and
Oncology Departments at Addenbrooke’s Hospital, the members of VoxTox and everyone
who worked with me at the Paul Scherrer Institute in Switzerland, especially Professor Tony
Lomax and Dr Francesca Albertini.
Thank you to the Medical Research Council for funding my research and financially sup-
porting my project in Switzerland. Also thank you to the NHS, ACT, IPEM, PARTNER
and CNAO who have all either contributed financially to this work or to my professional
and academic development during my time at Cambridge. Thank you to ETH, Orsay Proton
Centre and the Surrey Ion Beam centre for their contributions to this work in the form of
data, lectures, presentations and interesting conversation.
I want to mention the support from my friends and family and thank them. I also want to
thank Baden Rowing Club for their friendships and for making me feel so welcome during
my time working in Switzerland.
Finally, a special thank you to the Trebilcock family and their amazing contribution to the
Addenbrooke’s Charitable Trust, supporting research in this area.

Abstract
This dissertation addresses the issue of robustness in proton therapy treatment planning
for cancer treatment. Proton therapy is considered to be advantageous in treating most
childhood cancers and certain adult cancers, including those of the skull base, spine and
head and neck. Protons, unlike X-rays, have a finite range highly dependent on the electron
density of the material they are traversing, resulting in a steep dose gradient at the distal edge
of the Bragg peak. These characteristics, together with advancements in computation and
technology have led to the ability to plan and deliver treatments with greater conformality,
sparing normal tissue and organs at risk.
Radiotherapy treatment plans aim to meet set dosimetric constraints, and meet them at every
fraction. Plan robustness is a measure of deviation between the delivered dose distribution
and the planned dose distribution. Due to the same characteristics that make protons advan-
tageous, conventional means of using margins to create a Planning Target Volume (PTV) to
ensure plan robustness are inadequate. Additional to this, without a PTV, a new method of
analysing plan quality is required in proton therapy.
My original contribution to the knowledge in this area is the demonstration of how site-
and centre- specific robustness constraints can be established. Robustness constraints can
be used both for proton plan analysis and to identify patients that require plans of greater
individualisation. I have also used the daily volumetric imaging from patients previously
treated with conventional radiotherapy to quantify range uncertainty from inter- and intra-
fraction motion. These new methods of both quantifying and analysing the change in proton
range in the patient can aid in the choice of beam directions, provide input into a multi- cri-
teria optimisation algorithm or can be used as criteria to determine when adaptive planning
may be required. This greater understanding in range uncertainty better informs the planner
on how best to balance the trade-off between plan conformality and robustness in proton
therapy.
10
This research is directly relevant to furthering the knowledge base in light of HM Govern-
ment pledging £250 million to build two proton centres in England, to treat NHS patients
from 2018. Use of methods described in this dissertation will aid in the establishment of
clear and pre-defined protocols for treating patients in the future.
Table of contents
Table of contents 11
List of figures 13
List of tables 17
Abbreviations 21
1 Introduction 1
1.1 The Cancer Burden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 NHS England . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Proton Therapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Accelerating Protons for Therapy . . . . . . . . . . . . . . . . . . 6
1.2.2 Delivering Proton Therapy . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Treatment Planning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.4 Overview of PhD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2 Uncertainties in Proton Therapy 17
2.1 Uncertainties in Proton Planning . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.1 Sources of Range Uncertainty . . . . . . . . . . . . . . . . . . . . 17
2.1.2 Range Calculation Uncertainties in the Treatment Planning System 19
2.1.3 Discrepancies Between Planned Dose and Delivered Dose . . . . . 21
2.2 Managing Uncertainties Including Positional Discrepancies and Motion . . 23
2.2.1 Motion Mitigation . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2.2 Margins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.3 Probablistic Planning and Robust Optimisation . . . . . . . . . . . 26
2.2.4 Treatment Analysis and Validation . . . . . . . . . . . . . . . . . . 30
2.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
12 Table of contents
3 A Method for Quantifying Random Range Uncertainty 35
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 TomoTherapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3 Calculating WEPL Changes in the CTV from TomoTherapy Shift & MVCT
Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3.1 Coding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3.2 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4 Quantifying Site-specific Range Uncertainty 53
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2 Methods and Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.3.1 Prostate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.3.2 Head and Neck . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.3.3 Population Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5 Defining Robustness Protocols 81
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.1.1 Methods and Materials . . . . . . . . . . . . . . . . . . . . . . . . 82
5.1.2 Retrospective Robustness Analysis . . . . . . . . . . . . . . . . . . 85
5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.2.1 Retrospective Robustness Analysis . . . . . . . . . . . . . . . . . . 87
5.2.2 How to use the Robustness Database in Practice . . . . . . . . . . . 89
5.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6 Summary 99
6.1 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
6.2 MVCT and Proton Therapy . . . . . . . . . . . . . . . . . . . . . . . . . . 104
6.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
References 109
Appendix A Published Work and Presentations 121
A.0.1 Other Presentations . . . . . . . . . . . . . . . . . . . . . . . . . . 125
Table of contents 13
Appendix B Code 129

List of figures
1.1 An illustration of the proton dose fall-off with depth known as the Bragg
Peak. The target drawn above the Bragg peak represents how dose can be
deposited in the tumour at depth in the patient. . . . . . . . . . . . . . . . . 4
1.2 The world showing current and planned proton therapy centres (Centre lo-
cations provided at (Proton Therapy Today - The online magazine for proton
therapy, 2014), image created using Google Maps) . . . . . . . . . . . . . 5
1.3 Photograph of the cyclotron at Orsay, which feeds 3 treatment rooms with a
high energy of 230MeV . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Photograph of the synchrotron at CNAO, Italy used to accelerate high en-
ergy protons and carbon ions for patient treatments . . . . . . . . . . . . . 8
1.5 Schematic of Bragg peak delivery along a single profile through a target.
(a) A flat spread-out Bragg peak (SOBP) is achieved by placing spots with
increasing weights throughout the target to produce a uniform field, as used
in passive scattering and single-field uniform dose. (b) Only the most distal
single pristine Bragg peak (BP) is used for distal-edge tracking. (c) Opti-
mally weighted spots are positioned throughout the volume to achieve fields
with non-uniform doses for three dimensional intensity-modulated particle
therapy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.6 Photo of a patient specific collimator . . . . . . . . . . . . . . . . . . . . . 11
1.7 Photo of a scattering beam compensator . . . . . . . . . . . . . . . . . . . 11
2.1 Schematic of sources of range uncertainties in proton therapy . . . . . . . . 18
2.2 Illustrating two parameters for optimisation, conformity and robustness, and
how a Pareto front can be used to investigate clinically achievable plans with
the same Pareto-optimised solution, but with different values of importance
on robustness and conformity. . . . . . . . . . . . . . . . . . . . . . . . . 30
16 List of figures
3.1 TomoTherapy at Addenbrooke’s. Provides daily imaging with positional
correction on all patients (Burnet et al., 2010). . . . . . . . . . . . . . . . . 38
3.2 Example of how CT data (left image) can be converted into WEPL (right
image) at a beam angle of 90◦ using calibration curves and ray-tracing . . . 39
3.3 Image taken from the Addenbrooke’s Radiotherapy Department’s work in-
struction for carrying out MVCT HU QA using the Cheese phantom on
TomoTherapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.4 Average MVCT IVDT from several months on both TomoTherapy units at
Addenbrooke’s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.5 Calibration curve of MVCT HU to electron density relative to water . . . . 45
3.6 Calibration curve of MVCT HU to proton stopping power relative to water . 47
3.7 Image of truncated patient on TomoTherapy showing the the scanning circle
(white circle). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.1 Heat plots of %volume range error histograms for prostate patients 4, 7, 21,
34, and 36 (in order top to bottom) with fraction using 0 and 90 degree beam
angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.2 Heat plots of %volume range error histograms for prostate patients 58, 72,
78 and 95 (in order top to bottom) with fraction using 0 and 90 degree beam
angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.3 Plots of the standard deviation for patients 4, 7, 21, 34, and 36 (in order top
to bottom) of range error with fraction and beam angle. . . . . . . . . . . . 60
4.4 Plots of the standard deviation for patients 58, 72, 78 and 95 (in order top to
bottom) of range error with fraction and beam angle. . . . . . . . . . . . . 61
4.5 Screen shot of a single slice of Patient 78 showing the high dose CTV and
approximate angles used to calulate WEPL . . . . . . . . . . . . . . . . . 63
4.6 Heat plots of %volume range error histograms for all head and neck patients
with treatment fraction at 0 and 25 degree beam angles. . . . . . . . . . . . 65
4.7 Heat plots of %volume range error histograms for all head and neck patients
with treatment fraction at 70 and 110 degree beam angles. . . . . . . . . . 66
4.8 Heat plots of %volume range error histograms for all head and neck patients
with treatment fraction at 180 degree beam angle. . . . . . . . . . . . . . . 67
4.9 Plots of the standard deviation of range error for each patient with fraction
and beam angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
List of figures 17
4.10 The fractional means of WEPL difference for each prostate patient was plot-
ted as a histogram for both angles (top is 0◦ and bottom 90◦). The total mean
can be taken from these plots for each angle as well as the standard deviation
which is the systematic component of the range error for these sites. . . . . 70
4.11 The fractional means of WEPL difference for each head and neck patient
was plotted as a histogram for three angles (top is 0◦, middle is 25◦ and
bottom 70◦). The total mean can be taken from these plots for each angle
as well as the standard deviation which is the systematic component of the
range error for these sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.12 The fractional means of WEPL difference for each head and neck patient
was plotted as a histogram for three angles (top is 110◦ and bottom 180◦).
The total mean can be taken from these plots for each angle as well as the
standard deviation which is the systematic component of the range error for
these sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.13 The fractional means of WEPL difference for all prostate patients (top) and
all the head and neck (bottom) were plotted as a histogram for all angles.
The total mean can be taken from these plots as well as the standard deviation. 73
4.14 The binned prostate data for each beam angle from each fraction for all
patients was summed and plotted as a histogram normalised to one. Each
colour represents a beam angle, and the black histogram is of all binned data. 74
4.15 The binned Head and neck data for each beam angle from each fraction for
all patients was summed and plotted as a histogram normalised to one. Each
colour represents a beam angle, and the black histogram is of all binned data. 75
4.16 Example of tooth filling artefact in kV CT image (left) and in MVCT image
(right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.17 Possible weight loss leading to negative Σ for Head and Neck. Images left
to right of patients 52, 91, 92 . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.1 The ebDD represents in each voxel a dose-error-bar that brackets the cal-
culated deviations from the nominal dose distribution. Here are the screen-
shots of the nominal plan (using typical PSI beam arrangements referenced
as A(4f) in this chapter) and the ebDD calculated for ±3% HU from one
patient. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
18 List of figures
5.2 The A(4f) plan consists of four non-coplanar fields, two posterior oblique
and two anterior lateral oblique fields. Clockwise from top left the beam
angles and couch twists used are (290◦, 340◦), (70◦, 20◦), (110◦, 20◦) and
(250◦, 340◦), where the first angle is the gantry orientation and the second
is the table orientation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.3 Robustness results, represented as percentage dose uncertainty for 16 skull-
base IMPT plans from 8 patients to both the range and set-up uncertainties
modelled in (a) in the brainstem, (b) in the chiasm and (c) in the CTV. On
each box plot the central mark is the median, the edges of the box are the
25th and 75th percentiles, the whiskers extend to the most extreme data
points not considered outliers, and outliers are plotted individually. Outliers
are considered points ±2.7σ . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.4 Volume Histogram of (a) systematic range uncertainties and (b) random set-
up uncertainties for an underdose scenario in the CTV for five plans of dif-
ferent beam orientations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.5 Dose distribution for the A(4f) (top left) and B(4f) (top right) plans and
ebDD of systematic range uncertainties for the worst-case underdose sce-
nario for the A(4f) (bottom left) and B(4f)(bottom right) treatment plans
showing the CTV in green. . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.6 Dose distributions from the A(4f) (top left) and B(4f) (top right) plans and
the corresponding ebDD from B(4f) (bottom left) and A(4f) (bottom right)
treatment plans showing the brainstem in green and the target in yellow. . . 94
5.7 Volume histogram of (a) the systematic range uncertainties, (b) the random
set-up uncertainties for an overdose scenario in the brainstem and (c) the
random set-up uncertainties for overdose scenario in the chiasm, for five
plans of different beam orientations. . . . . . . . . . . . . . . . . . . . . . 95
A.1 REVIEW ARTICLE: Treatment planning optimisation in proton therapy . . 122
A.2 Accepted ESTRO33 ePoster 2014 . . . . . . . . . . . . . . . . . . . . . . 123
A.3 Winning popular press poster describing my research to the general public
at the School of Clinical Medicine Research day 2012 . . . . . . . . . . . . 126
A.4 Copy of my certificate for presenting a popular press poster at the clinical
school Research day . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
List of tables
1.1 Proton beam delivery techniques, production methods and planning tech-
niques. The further down the table, the more conformal the technique . . . 12
3.1 Calculating a mean β to use in the Bethe-Bloch equation to establish an
energy independent proton stopping power relative to water . . . . . . . . . 43
3.2 calculation of proton stopping power relative to water for each tissue using
the energy specific β in the Bethe-Bloch equation. The I-values have been
interpolated from the Gammex data tables . . . . . . . . . . . . . . . . . . 44
3.3 Table of density inserts used to calibrate the MVCT and their corresponding
physical densities, estimated electron densities and proton stopping powers
relative to water as calculated from the Bethe-Bloch equation using the β
value to a 150MeV proton beam . . . . . . . . . . . . . . . . . . . . . . . 46
4.1 Prostate data showing the target and number of fractions used. The final col-
umn shows the percentage of fractions that did not suffer target truncation
at daily imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2 Table of Intra-fraction changes in WEPL for a single prostate patient, 58, at
beam angles 0◦ and 90◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.3 Table of Head and Neck data showing the targets and number of fractions
used. The final column shows the percentage of fractions that did not suffer
target truncation at imaging . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.4 Table of prostate range error in mm for all patients at each beam angle. . . . 74
4.5 Table of head and neck range error in mm for all patients at each beam angle 75
5.1 Descriptions of each plan . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.2 The use of retrospective plan analysis to establish planning parameters of
robustness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
20 List of tables
5.3 Robustness values for each parameter set in the site specific robustness
database in Table 5.2 for plans A(4f) and B(4f). Bold italic values are those
that are out of tolerance . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.4 This table shows robustness results from the original A(4f) plan that patient
X was treated with and those from the re-plan created for this work B(4f).
Tolerances have not been set for these OARs as they had for the brainstem
and chiasm in the robustness database, but they have been supplied here to
use to further compare plans. . . . . . . . . . . . . . . . . . . . . . . . . . 97
Abbreviations
4DCT Four Dimensional Comptued Tomography
BP Bragg Peak
CBCT Cone Beam Computed Tomography
CNAO Centro Nazionale di Adroterapia Oncologica
CNS Central Nervous System
CPAC Compact Particle Acceleration Corporation
CRUK Cancer Research UK
CT Computed Tomography
CTV Clinical Target Volume
DDA Dose Difference Analysis
DECT Dual Energy Computed Tomography
DET Distal Edge Tracking
DH Department of Health
DICOM Digital Imaging and Communications in Medicine
DV H Dose Volume Histogram
DWA Dielectric wall accelerator
ebDD error-bar Dose Distribution
ebDD error-bar Dose Distribution
22 Abbreviations
ebV H error-bar Volume Histogram
FOV Field of View
GTV Gross Tumour Volume
HU Hounsfield Unit
ICRU International Committee on Radiation Units and Measurements
IGRT Image Guided Radiotherapy Treatment
IMPT Intensity-Modulated Particle Therapy
IMRT Intensity-Modulated Radiotherapy
IV DT Image Value Density Table
kVCT Kilovoltage Computed Tomography
Linac Linear Accelerator
MCO Multi-Criteria Optimisation
MVCT Megavoltage Computed Tomography
NICE National Institute of Health and Care Excellence
NRAG National Radiotherapy Advisory Group
OAR Organ at Risk
PET Positron Electron Emission
PRV Planning organ at Risk Volume
PSI Paul Scherrer Institute
PTV Planning Target Volume
QALY Quality-Adjusted Life Years
RBE Relative Biological Effectiveness
RF Radio Frequency
Abbreviations 23
ROI Region of Interest
SFUD Single Field Uniform Dose
SID Source to Isocentre Distance
SOBP spread-out Bragg peak
T PS Treatment Planning System
VOI Volume of Interest
WEPL Water Equivalent Path-Length

Chapter 1
Introduction
Throughout the first two chapters I have reproduced, with permission, the literature review
which I wrote, entitled Treatment planning optimisation in proton therapy published in the
British Journal of Radiology (McGowan et al., 2013).
1.1 The Cancer Burden
Cancer is a major worldwide health problem with approximately 11 million people being
diagnosed with cancer and 7 million people dying from cancer each year (Dosanjh et al.,
2008). Cancer Research UK (CRUK) reports that in the UK there were 331,487 cancer
cases in 2011, 274,233 of which presented in England alone. Of these, 7490 cases were
brain, other central nervous system (CNS) and intracranial tumours (Cancer Research UK,
2014c). Based on the same set of statistics, childhood cancers (<14 years) account for
approximately 1% of cancer cases in the UK with 1,574 new cases of cancer diagnosed in
children between 2009 and 2011. Brain, other CNS and intracranial tumours are the most
common cause of childhood cancer death. (Cancer Research UK, 2014a).
Despite increased post-diagnosis life expectancy due to developments in early diagnosis
and treatment, incidences of cancer will increase in the coming years due to the increas-
ing lifespan of the general population. Of the adult patients diagnosed with a malignant
tumour, CRUK estimates that 50% will survive their cancer for at least 10 years (Cancer
Research UK, 2014b). For 40% of these patients, radiotherapy will play a role in their treat-
ment, whether on its own or as part of a combined therapy. This corresponds to around
100,000 patients receiving radical radiotherapy each year in the UK; for 60% of patients it
is given with curative intent. Radiotherapy includes external beam radiotherapy and also
2 Introduction
brachytherapy. External beam radiotherapy can include the use of high energy photons,
electrons, protons or carbon ions. In this thesis photon based, or X-ray based radiotherapy,
will be referred to as conventional radiotherapy whereas radiotherapy using protons will be
referred to as proton therapy. Other forms of particle therapy are beyond the scope of this
thesis, the focus of which is on developments in proton therapy in line with the prospect of
bringing proton therapy to England.
1.1.1 NHS England
In 2007, the National Radiotherapy Advisory Group (NRAG) recommended for the De-
partment of Health to give NHS patients access to proton therapy (Department of Health,
2014). A business case was proposed to establish English proton therapy centres and in the
meantime overseas access was established. In the UK, patients have been able to access
proton therapy abroad under the auspices of the National Health Service Proton Overseas
Programme since 2008 (NHS Specialised Services, 2011). Between 2008 and 2012 160
patients, 107 of whom were children, were treated with proton therapy overseas. In 2009
NRAG determined that there was a need to treat up to 1,487 patients from England and the
devolved administrations with proton therapy each year. Four hundred of these patients are
classed as high priority patients, consisting of a casemix of complex cases including adults,
children and very young children (Department of Health, 2014).
The decision to establish two proton therapy centres in the UK was made based on the costs,
benefits and risks of each option investigated. Two centres, one in Manchester and one in
London, has been determined to be the most cost effective option for the provision of the
NHS service. This solution is based on the incremental cost effectiveness ratio of each
option compared to the baseline option, misleadingly named, ‘do nothing’.
The ‘do nothing’ option refers to the cost to treat these patients with conventional radio-
therapy within the UK. Conventional radiotherapy is priced in fractions, with each fraction
costing between £85 and £242. This works out on average, £7,500 to treat a child, and
£4,600 to treat an adult (Department of Health, 2014). The measure of benefit to a patient
from carrying out an intervening treatment is measured in Quality-Adjusted Life Years, or
QALYs . QALYs take into account both the quantity and quality of life created by the in-
tervening healthcare, in this case radiotherapy. It is the product of life expectancy and a
measure of the quality of the remaining life-years. A QALY places a weight on the time the
patient has at different levels of healthiness. For example a year of perfect health is worth 1
1.1 The Cancer Burden 3
and a year of less than perfect health is worth less than 1, with death considered to be equiv-
alent to zero (Phillips, 2009). The estimate of QALYs per patient gained over their lifetime
from radiotherapy treatment is between 9.4 and 14 years (Department of Health, 2014).
This is the ‘do nothing’ baseline: on average from radiotherapy treatment one quality- ad-
justed life year costs £517. However, it should be noted that the concept of the QALY has
limitations; it is based on survival probabilities and does not take into account the context of
each given patient’s life or chronic side effects from treatment such as deafness or a drop in
IQ (Phillips, 2009). However, the QALY is a useful metric for comparing treatment options
in terms of benefit received and the cost incurred.
Models have been used to determine the difference in QALYs between treating patients
deemed to benefit from proton therapy based on current evidence. It was estimated that
the high priority patients could benefit from an extra 2.5 QALYs from being treated with
proton therapy over conventional radiotherapy and other patients referred for protons could
gain 1.8 extra QALYs from this treatment (Department of Health, 2014). Based on this,
the number of patients treated per year and the 30-year lifetime of the machines the cost
is £19,187 per QALY. Due to that fact that 80% of the estimated costs are fixed, it is the
limited throughput of patients that drives the cost up. To put this into perspective, NICE’s
(National Institute for Health and Care Excellence) notional threshold on per QALY cost is
£30,000. Costs for drugs such as trastuzumab (Herceptin) and imatinib (Glivec), which are
both supplied on the NHS are above this threshold (Towse, 2002). It is estimated that each
of the two centres will treat 675 patients each a year. The low throughput takes into account
the complicated casemix of patients to be treated. There is no facility in the world that treats
an equally complex casemix of patients as those proposed for England.
On the 10th February 2012, the Her Majesty’s Treasury and the Department of Health for-
mally confirmed full approval of the Strategic Outline Case for the National Proton Beam
Therapy Development Programme, pledging £250 million (Department of Health, 2014). It
has decided that all patients treated in England with proton therapy, will be treated within
clearly pre-defined protocols (Department of Health, 2014). With the aim to provide fur-
ther evidence of the effectiveness of proton therapy every patient will also be enrolled into
a prospective programme of evaluation with outcome tracking. It is hoped that the centres
will start to treat patients in 2018 (Burnet et al., 2012; www.dh.gov.uk, 2012), such that
developments discussed within this thesis will be directly relevant.
4 Introduction
1.2 Proton Therapy
The aim of radiotherapy is to create and deliver the ideal treatment plan, where the tar-
get volume receives 100% of the prescribed dose and normal tissue receives 0% (Bortfeld,
1999). It is, however, impossible to achieve this perfect balance. Instead, multiple trade-
offs are required to achieve a clinically acceptable plan, so the problem becomes one of
optimisation.
Fig. 1.1 An illustration of the proton dose fall-off with depth known as the Bragg Peak. The
target drawn above the Bragg peak represents how dose can be deposited in the tumour at
depth in the patient.
The nature of proton therapy makes the aim of cure without complications potentially more
achievable due to the highly localised deposition of dose in the characteristic Bragg peak
(Bragg, 1904) (see Figure 1.1). This relates predominantly to the ability to deliver high
doses of radiation close to normal tissue structures which are dose-limiting in conventional
X-ray treatments, and to the finite range of protons, which results in a reduced integral dose
to surrounding normal tissues. However, for a plan to be optimal it must also be deliverable.
The term robustness is used within this thesis to define how well a plan can maintain its
dosimetric integrity at every fraction throughout the patient’s treatment. It is the finite range
and the steep dose gradients achievable in proton beam therapy that lead to its inherent lack
of robustness when compared to conventional radiotherapy.
From a clinical perspective, the exact role of proton therapy has yet to be defined. However,
1.2 Proton Therapy 5
for childhood cancers proton therapy delivers a lower dose to tissues around the tumour
and lower whole body dose. This results in less growth disturbance and lower risk of sec-
ondary malignancies, compared to X-rays. There is also the suggestion that the use of
proton therapy can reduce impairment of neuropsychological and IQ development (Mer-
chant et al., 2008). In adults, proton therapy seems particularly effective in the treatment
of radio-resistant tumours close to critical structures such as the brain stem and spinal cord.
For example, outstanding results have been published for the use of proton therapy in the
treatment of chordoma and chondrosarcoma (Ares et al., 2009). The current evidence for
the use of proton therapy at different sites has been extensively reviewed (Allen et al., 2012;
Ares et al., 2009; De Ruysscher et al., 2012; Jones, 2008). The reader is referred to the
literature (Durante and Loeffler, 2010; Lodge et al., 2007; Ramaekers et al., 2011; van de
Water et al., 2011) for more clinical data. Nevertheless, substantial opportunity for further
clinical research development and evaluation remains (McGowan et al., 2013).
Presently, proton therapy has not been investigated to either the desired or required extent
through clinical trials. This is due to several factors including limited samples sizes, ethical
implications (Glimelius and Montelius, 2007; Sheehan et al., 2014) and, most significantly,
the lack of experience as compared to that in conventional radiotherapy. A significant ad-
vantage of conventional X-ray therapy over proton therapy is the wealth of experience and
knowledge available. A key area for optimising the treatment planning process is in gaining
experience in planning proton treatments. The number of proton facilities available world-
wide is rapidly increasing (Figure1.2), yet there is a substantial shortage of oncologists,
dosimetrists and physicists with the required expertise (Dosanjh et al., 2008).
1.2.1 Accelerating Protons for Therapy
The aim of radiotherapy is to be able to deliver the desired dose anywhere in a patient and
for protons, the penetration depth is energy dependent. Therefore to be able to adequately
treat to depths of 26-38cm in a patient, beam energies of 200-250 MeV are required. The
energy also needs to be variable during treatment delivery to allow for uniform coverage
of the treatment volume at all depths. It is also required that the treatment delivery time
be as short as possible so that more patients can be treated, for patient comfort, and to
reduce uncertainties caused by motion. Depending on the delivery system used, to achieve
a 2Gy dose uniformly to a target volume of 1 litre in 1 minute intensities of about 3-6×1010
particles per second are required (ICRU 78, 2007). The current methods for accelerating
protons for radiotherapy are briefly described.
6 Introduction
Fig. 1.2 The world showing current and planned proton therapy centres (Centre locations
provided at (Proton Therapy Today - The online magazine for proton therapy, 2014), image
created using Google Maps)
Cyclotron
The cyclotron, (Figure 1.3), was invented by Ernest Lawrence in 1929 (Lawrence, 1934). It
uses a large magnet to confine the motion of the particles to that of a spiral path. As par-
ticles enter into the opposite half of the accelerator, known as a ‘dee’, they are accelerated
by a radio frequency (RF) electric field. A large potential difference is applied across the
gap between the dees, giving particles a ‘kick’ every half orbit. As they gain energy the
radius of their orbit increases spiralling outward until they are extracted at the maximum ra-
dius, and energy, capable for the dimensions of the cyclotron. For non-relativistic particles,
Bvq = mv
2
r , where B is the magnetic field, q is the charge of the particle, v is the velocity
of the particle, m is the mass of the particle and r is the radius of orbit. In considering
the frequency, f, of the RF field required, the relationship between frequency, velocity and
radius can be used. Since ω = 2π f = vr then, f =
Bq
2πm . This shows that for a particle with
a constant mass the frequency is independent of the radius and so a fixed RF can be used.
This only works for non-relativistic cases as the particles nears the speed of light its mass
will increase and this convenience no longer stands (Kirby, 2010)). There are two types of
accelerators to work around this problem. One is a Synchrocyclotron that uses a variable
RF. The other is an isochronous cyclotron, which increases the magnetic field with radius.
1.2 Proton Therapy 7
Fig. 1.3 Photograph of the cyclotron at Orsay, which feeds 3 treatment rooms with a high
energy of 230MeV
Synchrotron
The synchrotron has a fixed radius so that a bunch of particles can be confined in an evacu-
ated pipeline (see Figure 1.4). The particles travel in straight segments with bending mag-
nets in between. The particles are accelerated using an RF oscillator that supplies an energy
boost each time a particle passes through the accelerating cavity. The RF increases as the
particles increase in energy and velocity and is synchronised with an increase in the mag-
netic field in the bending magnets. The synchrotron offers a flexible energy variation system
where the beam energy can quickly be changed instead of using a physical range modifier
in the treatment delivery system. A limitation though is the number of particles per bunch
and therefore a limit on intensity.
New Technology
As all the methods of particle acceleration mentioned above require large accelerators to be
built they cannot be used on a gantry as with conventional radiotherapy Linacs (radiotherapy
linear accelerators). This means that one accelerator will feed particles to several treatment
rooms allowing only one treatment room to operate at a time. Several new technologies are
being developed by different groups to overcome these problems.
■ The Dielectric wall accelerator (DWA) - The Compact Particle Acceleration Cor-
poration (CPAC) is developing and commercializing a highly compact fixed beam
8 Introduction
Fig. 1.4 Photograph of the synchrotron at CNAO, Italy used to accelerate high energy pro-
tons and carbon ions for patient treatments
intensity-modulated proton therapy (IMPT) system powered by the dielectric-wall
accelerator. DWA uses high gradient linear accelerator technology to produce a com-
pact accelerator capable of producing proton energies up to 200MeV and suitable for
gantry mounting (Caporaso et al., 2009).
■ Laser plasma particle accelerators can create a very large and localised electric field
using a laser to knock protons from a thin target. It is theorized that by using a given
target thickness, therapeutic energies can be reached using an accelerator that can be
mounted on a gantry (DeLaney and Kooy, 2008).
■ Non-Scaling Fixed Field Alternating Gradient accelerator known as PAMELA in the
UK. This design aims to incorporate the best qualities of a cyclotron (fixed magnetic
field) for rapid acceleration and of a synchrotron (alternating gradient) for variable
energy extraction. This group envisage that this type of accelerator will be able to
be gantry mounted and will be suitable for delivering spot scanning and rescanning
techniques (Peach et al., 2009).
1.2 Proton Therapy 9
1.2.2 Delivering Proton Therapy
The choice in how the protons are delivered can have a large impact on the ability to both
produce conformal dose distributions, and to produce plans that are robust to uncertainties.
There are three main treatment delivery techniques used clinically (Table 1.1):
1. passive scattering (Haberer et al., 1993),
2. uniform, and
3. active scanning (DeLaney and Kooy, 2008; Pedroni et al., 1995).
These techniques are used to broaden the narrow proton beam created by the accelerator
into one that can achieve uniform dose coverage of the target at all depths. This is achieved
for passive scattering and uniform scanning through the delivery of so-called spread-out
Bragg peaks (SOBP) (Figure 1.5). Orthogonal to the beam direction, the beam is spread
using carefully designed scatterers (for passive scattering) or by continually deflecting the
proton beam in a regular pattern orthogonal to the beam direction with constant intensity
(uniform scanning). For both approaches, three-dimensional (3D) conformation of the final
dose to the target is achieved through the additional use of patient- and field-specific colli-
mators (which conform the dose in directions orthogonal to the beam - see Figure 1.6) and
compensators (which conform the dose in the beam direction - see Figure 1.7) inserted in
the beam nozzle (DeLaney and Kooy, 2008). Active scanning also uses magnets to scan
the proton beam across the target volume. In contrast to uniform scanning, however, active
scanning allows the fluence (dose) applied at each Bragg position to be continuously varied.
Active scanning can offer an advantage to the patient by allowing for greater flexibility in
the delivered dose and a reduction in integral dose to healthy tissues. This allows for the
delivery of intensity-modulated particle therapy (IMPT) (Lomax, 1999), which is analogous
with Intensity- Modulated Radiotherapy (IMRT) in conventional radiotherapy. Although
there is in principle, a continuum of solutions to the IMPT problem, at its extremes IMPT
can be divided into two flavours: distal-edge tracking (DET) (Deasy, 1998), where Bragg
peaks are placed only at the edges of the target volume, and 3D IMPT (Lomax, 1999), where
Bragg peaks are optimally distributed and weighted throughout the target volume (Figure
1.5). IMPT allows for delivery of single inhomogeneous but optimised fields to produce a
final homogeneous dose distribution in the target volume. This permits the planner to be
more flexible in the placement of residual dose to healthy tissues. However, although IMPT
offers greater optimisation of dose delivery at the planning stage, it has the potential to be es-
pecially sensitive to range uncertainties (Nill et al., 2004; Oelfke and Bortfeld, 2000). Table
10 Introduction
1.1 summarises all these modes of producing and manipulating proton dose distributions.
Fig. 1.5 Schematic of Bragg peak delivery along a single profile through a target. (a) A
flat spread-out Bragg peak (SOBP) is achieved by placing spots with increasing weights
throughout the target to produce a uniform field, as used in passive scattering and single-
field uniform dose. (b) Only the most distal single pristine Bragg peak (BP) is used for
distal-edge tracking. (c) Optimally weighted spots are positioned throughout the volume
to achieve fields with non-uniform doses for three dimensional intensity-modulated particle
therapy.
1.2 Proton Therapy 11
Fig. 1.6 Photo of a patient specific collimator
Fig. 1.7 Photo of a scattering beam compensator
12 Introduction
Table 1.1 Proton beam delivery techniques, production methods and planning techniques.
The further down the table, the more conformal the technique
Methods of producing a
clinical proton beam to treat
entire target volume Descriptions
Passive scattering Works on the principle that high atomic number materials, such
as lead, scatter the beam with minimum energy loss and low
atomic number materials, such as plastic, decrease proton energy
with minimum scatter. Combining these materials to produce
patient-specific collimators and compensators results in a
conformal treatment beam with a spread-out Bragg peak.
Uniform scanning This is similar to passive scattering with the difference that the
beam is spread in the lateral direction through magnetically deflecting
the beam with constant fluence instead of using a scattering foil.
Different spot weights are produced using a compensator, as in
passive scattering
Active Scanning This uses magnetic fields to deflect the path of each proton beam
towards the planned position in the target volume. Individual Bragg
peaks are distributed within the target volume and the cumulative effect
produces an effective SOBP without the need for compensators. This is
achieved by either continuous magnetic scanning or spot scanning.
The latter is analogous to the step-and-shoot mode in IMRT, i.e. a non-
continuous delivery of dose, where the exact position is determined before
the dose is delivered
Methods of achieving
adequate dose distributions Descriptions
Single Field Uniform Dose (SFUD) Single individually optimised proton fields that each deliver a homogeneous
dose to a volume. If necessary, these can be combined by simple addition
Field patching The sharp distal edge dose gradient can be matched up to the lateral edges
of another “patch” field to produce a continuous dose distribution. Where
possible, equivalent opposite fields are also used to reduce the potential for
dose variation at the abutting edges. Multiple fields in patch work can be
used to achieve multiple dose gradients inside a treatment volume. Field
patching is a 3D extension of matching lateral field edges. Therefore, if
multiple fields are used, each one can deliver a homogeneous dose to part of
the volume
IMPT IMPT is analogous to IMRT, and is a mode of treatment delivery achievable
only with active scanning beams. IMPT uses narrow proton beams which are
magnetically moved over the volume in the transverse plane while the
energy and intensity are altered to control dose to a point and sculpt the dose
at depth. Unlike SFUD treatments, IMPT can deliver a number of non-
uniform fields to produce the desired dose distribution.
"Flavours" of IMPT Descriptions
3D IMPT This is most similar to IMRT. Bragg peaks are placed throughout the entire
volume and their weights optimally adjusted.
DET DET is a method by which pristine Bragg peaks of optimal weights are
distributed only along the distal edge of the target and not throughout
the target volume.
1.3 Treatment Planning 13
1.3 Treatment Planning
The developments in technology have allowed for the ability to deliver treatments of greater
conformality and complexity in radiotherapy. This led to the need for treatment planning
using inverse planning and optimisation. In current inverse planning systems for IMRT the
dose distribution is determined using a computerised optimisation based on dose objectives
for targets and other volumes which have been assigned an importance level (Thieke et al.,
2007). To determine the plan quality, a number is assigned based on the deviation from ob-
jective dose for each volume and the optimisation result is the plan with the lowest number.
This is a trial and error process, given that the resulting plan may not be clinically acceptable
and the importance levels assigned to each volume may need adjusting and the optimisation
re-run. This can be time-consuming and the best-quality plan may not have been achieved
as the planner cannot try every combination of parameters. There also exists a problem that,
if upper and lower constraints are met, the optimisation process will not further improve
doses to these volumes. This means that IMRT and IMPT cannot necessarily be exploited
to their full potential owing to limitations with inverse planning (Thieke et al., 2007).
Besides meeting planning constraints the treatment plan is also required to meet the aims
of the planner at each and every fraction. The ability to plan and deliver treatments of
greater conformality and complexity led to the need for uniform recording, reporting and
prescribing to account for geometric uncertainties. An uncertainty refers to any variation
between the planned and delivered dose distribution and they arise from the fact that treat-
ment plans are based on a static view of a living, and therefore, a very non-static patient.
The variation between the planned dose distribution and the delivered dose distribution is
a measure of a plan’s robustness: the smaller the variation the greater the plan robustness.
Residual uncertainties are unavoidable despite efforts in variation reduction such as patient
immobilisation and image guided radiotherapy treatment (IGRT) . International Committee
on Radiation Units and Measurements (ICRU) Reports 50 and 62 (ICRU 50, 1993; ICRU
62, 1999) provided a solution by defining target volumes in conventional radiotherapy in-
cluding the Gross Tumour Volume (GTV) , the visible tumour, the Clinical Target Volume
(CTV) , which includes the GTV and microscopic growth, and the Planning Target Volume
(PTV) which encompasses both the GTV and CTV as well as accounting for any deviations
from our static model. Analogous to the PTV, a Planning organ at Risk Volume (PRV) is
used for Organs at Risk (OAR) (ICRU 83, 2010).
In regard to conventional radiotherapy, Stroom et al. (1999) and van Herk et al. (2000)
14 Introduction
suggested that uncertainties fall into two categories, namely preparation and treatment exe-
cution uncertainties. Treatment execution uncertainties are random and act to blur the ideal
dose distribution by a Gaussian distribution with a width dependent on the photon penum-
bra and the standard deviation of day-to-day variations in the CTV location. Preparation
uncertainties are present at every fraction of the treatment and are therefore systematic in
nature. Using the shift invariance assumption, systematic uncertainties result simply in a
shift of the dose distribution with the same direction and magnitude.
Van Herk also addressed the problem of how to combine both these types of uncertainty in
an analytical model to determine the CTV-PTV margin. It was shown that it was incorrect
to combine the standard deviations of the uncertainties linearly due to the lack of correlation
between systematic and random uncertainties. It was acknowledged that the systematic un-
certainty is stochastic across the patient population. Van Herk’s PTV margin for systematic
uncertainties is defined so that the CTV has a 90% chance of being in the PTV. This, now
well known, margin recipe (equation 1.1) allows that in 90% of patients treated to meet the
ICRU recommendations, the CTV receives at least 95% of the prescribed dose.
PTV margin = 2.5Σ+1.64(
√
(σ2p +σ
2))−1.64σp (1.1)
where, Σ is the standard deviation of the systematic uncertainty, σ if the standard deviation
of the random uncertinaity and σp is the unblurred beam penumbra width. In this paper Van
Herk (van Herk et al., 2000) suggests redefining the PTV as,
"The volume defined in treatment room coordinates to which the prescribed
dose must be delivered in order to obtain a clinically acceptable and specified
probability that the prescribed dose is actually received by the CTV, which has
an uncertain location.”
In this way uncertainties are condensed into the CTV-to-PTV margin. The static PTV repre-
sents the moving CTV and is therefore a useful method of evaluating a plan, by using Dose
Volume Histograms (DVHs) to report the minimum dose to the CTV with a pre-specified
confidence level. The PTV then becomes a surrogate target for both treatment plan optimi-
sation and evaluation.
Planning with a PTV has some limitations; for example it is problematic in cases where the
CTV is close to the skin. It is important that there is confidence in quantifying and identify-
ing uncertainties used to grow the PTV margin and that they are specific for treatment site,
1.3 Treatment Planning 15
modality used, and treatment regime to ensure optimal treatment. In using Image Guided
Radiotherapy (IGRT) greater information about the systematic uncertainty can be acquired,
allowing for corrections to be made and a reduction in the margins used. IGRT is the use
of imaging at pre-treatment and delivery to improve or verifiy the accuracy of radiother-
apy. IGRT includes methods of simple visual field alignment checks, to the more complex
volumetric imaging that allows direct visualisation of the target volume and surrounding
anatomy (NRIG, 2012). In the case of online IGRT (where the images are analysed and
used on set), information about the random component can be gained and corrected for by
applying rigid shifts to the treatment table, aligning either soft tissue or bony landmarks in
the daily image to the planning image. This again, allows for a reduction in the CTV to PTV
margin. However, if reduced too much the margin may become inadequate and the tumour
control probability or other endpoint may worsen, such as biochemicall failure as seen by
(Engels et al., 2009). If the margin is too large, greater incidence and severity in normal
tissue complication may be seen.
Furthermore, the PTV concept is still present in newer approaches where inhomogeneous
doses are prescribed. Such cases are in conflict with several of the underlying assumptions
taken in van Herk’s original margin design. As a consequence, the indirect estimation of the
confidence level becomes unreliable. Intensity modulated particle therapy (IMPT) is one
such example, where complex and inhomogeneous dose distributions are delivered, often
with steep dose gradients. Along with the extra degree of freedom in the proton range,
this renders the PTV an inadequate planning tool. This is partly due to the shift invariance
assumption (also known as static-cloud) as it cannot be applied to the proton case. This
assumption relies on a rigid isocentric shift of the dose distribution, in the same direction
and magnitude as the geometric uncertainty, and is representative of the delivered dose when
experiencing the same geometric uncertainty. For example if a shift, S is applied, this means
that the dose value at a voxel Di, is now D(i−S). The shift invariance relies on assuming the
patient contour is unchanged when the dose isocentre is moved and also on the radiation
passing through tissue of the same density before and after the dose isocentre shift (Nguyen
et al., 2009). The sensitivity of the proton range to density heterogeneities in the traversing
material means this assumption can no longer be made. Even if there has been no shift by
the patient and no change in the tumour size, shape or positioning there is still an uncertainty
due to changes in the radiological path-length, resulting from patient changes such as weight
gain or loss.
16 Introduction
1.4 Overview of PhD
With the primary focus on including range uncertainty into proton therapy treatment plan-
ning, the purpose of this PhD is to explore how range uncertainties may be quantified and
to develop a solution for evaluating plan quality in terms of robustness that can be included
into the treatment planning workflow. The structure of the thesis is as follows:
Chapter 2 The current knowledge of range uncertainties, including source, magnitude and
management, as well as their impact on proton dose is discussed. Methods of how
to incorporate range uncertainties into the plan optimisation and evaluation process is
examined.
Chapter 3 A method has been suggested to quantify the inter-fraction range uncertainty in
proton therapy by calculating range changes throughout treatment using daily volu-
metric imaging. The process of creating a calibration curve of HU to relative proton
stopping power for the TomoTherapy MVCT has been laid out and assumptions and
limitations considered.
Chapter 4 Clinical patient data is used to measure random range uncertainty from inter-
fraction motion for a smaple of patients and two treatment sites. This includes analysis
of individual patient data and of the samples as a whole. This information may be used
within plan optimisation or evaluation, or both.
Chapter 5 A method to introduce robustness analysis into plan evaluation is presented. In
particular, a method to define a site-specific, and centre-specific, robustness database
by retrospectively analysing clinical IMPT plans is shown. Such a database could
be used during the plan evaluation and during the plan calculation phase to aid the
planner to select the optimal plan.
Chapter 6 Discussions
Chapter 2
Uncertainties in Proton Therapy
2.1 Uncertainties in Proton Planning
Sources of uncertainty present in conventional radiotherapy also apply to proton therapy.
Most geometrical uncertainties can be managed in the same way, including variation in de-
lineation, set-up uncertainties, imaging inaccuracies and patient motion. However, there
exists an important and fundamental difference in proton therapy, namely the proton range.
Protons have a finite range highly dependent on the electron density of the material they
are traversing, resulting in a steep dose gradient at the distal edge of the Bragg peak. Posi-
tioning of these dose gradients is critical to successful planning and treatment. Therefore,
an uncertainty of even a few millimetres can lead to under-dosage in the target volume or
over-dosage of an organ at risk (OAR).
Several authors have addressed the problem of range uncertainties in proton therapy, and the
purpose of this section is to analyse their conclusions. Understanding the causes and mag-
nitude of range uncertainties and incorporating them into the planning process is essential
for optimised proton planning.
2.1.1 Sources of Range Uncertainty
The main factors leading to range uncertainty are shown in (Figure 2.1). Because the main
advantage of using protons in cancer treatment is their finite range, this advantage can be
fully exploited only if the proton range in the patient can be precisely predicted (Lomax,
2008a). It has been suggested that range uncertainties can be between 1 and 15mm for lung
tumours (España and Paganetti, 2010), but larger changes are possible owing to anatomical
18 Uncertainties in Proton Therapy
changes in the patient (e.g. weight loss or gain and differential filling of anatomical cavities).
Uncertainties are normally compensated for in conventional radiotherapy by introducing
safety margins around the treatment volume and around OARs to produce a planning target
volume (PTV) and planning OAR volume(PRV), as recommended by the ICRU Reports 50,
62 and 83 (ICRU 50, 1993; ICRU 62, 1999; ICRU 83, 2010). A similar method has been
recommended by the ICRU for protons (ICRU 78, 2007). The larger the safety margin, the
less conformal the resulting dose distribution. Therefore, to achieve an optimum proton
treatment plan, the range prediction needs to be as accurate as possible.
The variables that give rise to uncertainties in the range prediction (Figure 2.1) can be di-
vided into two main groups: those causing uncertainties in the range calculation in the treat-
ment planning system (TPS) , and those leading to discrepancies between planning dose and
delivered dose.
Fig. 2.1 Schematic of sources of range uncertainties in proton therapy
2.1 Uncertainties in Proton Planning 19
2.1.2 Range Calculation Uncertainties in the Treatment Planning Sys-
tem
Inaccuracies Arising from the Planning CT
With respect to range calculation uncertainties, these can arise from inaccurate data exported
to the TPS. CT is used to acquire patient image data and the Hounsfield units (HUs) are then
converted into proton-stopping powers so dose calculations can be made. Uncertainties arise
in proton range calculation from CT-based plans owing to inaccuracy in the HU to proton
stopping power conversion and inaccuracies in the HU values themselves (Chvetsov and
Paige, 2010). Inaccuracies in the HU values are caused by noise, CT artefacts and beam
hardening.
Schaffner and Pedroni (1998) and España and Paganetti (2010) investigated how the con-
version of HUs to stopping power affects the range calculation in comparison with the real
treatment range. España and Paganetti (2010) tested different conversion methods, in-
cluding the traditional stoichiometric calibration method. Positron emission therapy (PET)
imaging was used to determine the range of the proton beam in a phantom and compare
it with the calculated range. It was noted that the back-to-back photons imaged with PET
are generated from electron– positron annihilations caused by inelastic nuclear interactions
between protons and target nuclei and not from atomic interactions, which primarily leads
to dose deposition. España and Paganetti (2010) and Schaffner and Pedroni (1998) both
concluded the same result: that the uncertainty caused by conversion was, ±1%. There is
research being undertaken into the development of proton CT. This would remove the un-
certainty associated with the conversion from HUs to proton stopping powers (Schulte et al.,
2004).
Noise is a stochastic uncertainty that either adds or subtracts from the HU value; this type
of uncertainty is important only if the proton beam is sensitive enough to be affected by
these changes in HU value. It was concluded that uncertainties caused by noise have a
similar contribution to conversion uncertainties ±1% (Schaffner and Pedroni, 1998). Beam
hardening had a greater effect on the assigned HU value. This was dependent on the position
and density of the tissue and added uncertainties of the order of±1.8% and±1.1% for bone
and soft tissue, respectively (Schaffner and Pedroni, 1998).
It is essential for CT-scanner-specific calibrations to be carried out (ICRU 78, 2007; Schaffner
and Pedroni, 1998). This is because, even though proton stopping power is independent of
20 Uncertainties in Proton Therapy
the imaging X-ray energy and the position of the target, HUs are dependent on the X-ray
spectrum and target position. Each scanner will produce a different X-ray spectrum gener-
ated with a different tube potential and current, therefore requiring individual calibration.
Moyers et al. (2010) carried out the measurement of the relative linear stopping powers of
21 different tissue substitutes. These were then scanned using both kilo-voltage and mega-
voltage CT and the relationship between stopping power and HU value was determined.
Lomax (2008b) has described how the combined effects of the proton’s sharp distal fall-
off, finite range and multiple Coulomb scattering can have an impact on the sensitivity of
plans to density heterogeneities in the patient. At the Paul Scherrer Institute in Switzerland
the CT scanner has been stoichiometircally calibrated to have an accuracy of 1% for soft
tissue and 2% for bony tissue, but, owing to inherent uncertainties such as beam hardening,
reconstruction artefacts and reconstruction algorithms (Lomax, 2008b), it has been stated
that an uncertainty in HU value of 3% is more realistic. To investigate the effect of this
uncertainty, each plan was recalculated with the HU value increased and decreased by 3%
to simulate an undershoot and overshoot scenario. The results were obtained for a simple
prostate case and a skull base case. The distal-edge tracking (DET) approach was found
to incur a systematic over- and under-dosage of the clinical treatment volume (CTV) of 5%
when a±3% HU uncertainty was introduced, whereas the 3D IMPT dose–volume histogram
(DVH) for the nominal, under- and overshoot plans showed little difference suggesting DET
is less robust than 3D IMPT.
Inaccuracies arising from the Dose Calculation Algorithm
Lomax (2008a) investigated the limitations in analytical dose calculations and the effects of
uncertainties in density calculation from CT data in DET and 3D IMPT dose distributions.
To investigate uncertainties arising from using an analytical dose calculation algorithm the
dose distributions for Version 1 and Version 2 DET and 3D IMPT skull base plans were
compared with the same plans calculated using Monte Carlo models (Version 1 plans use
less stringent constraints on OARs than Version 2 plans, therefore, target coverage has a
greater priority in the optimisation in Version 1 plans). In contrast to using an analytical
mathematical algorithm to calculate dose distributions, Monte Carlo models use a probabil-
ity distribution to model interactions and the production of secondary particles in a medium
for a given energy. For the skull base case the two calculation methods were in close agree-
ment for an acceptance level of ±10%, but decreased for lower acceptance levels. At all
acceptance levels, 3D IMPT plans showed better agreement in the PTV than DET plans,
2.1 Uncertainties in Proton Planning 21
with 87% of points in the 3D IMPT plan agreeing to within ±3% and only 80% of points
agreeing to within ±3% in the DET plan for Version 1 plans. However, for the Version 2
plans, agreement falls to 77% and 70% for 3D IMPT and DET, respectively. In the OARs
it was found that the Monte Carlo calculation predicted smaller doses in the optic nerve
and brain stem than the analytical calculation, especially in the Version 2 plans. The trend
showed that for both cases 3D IMPT showed smaller differences than the DET plans. It was
concluded that DET was relativity sensitive to calculation uncertainties; in comparison, 3D
IMPT was more robust with respect to both types of uncertainty.
In planning with photons there is a wealth of experience in prescribing doses to PTVs and
dose constraints to OARs. This has meant that when planning with protons there was a need
to prescribe doses relative to those used previously with photons by applying the concept of
radiobiological equivalent dose, Gy(RBE). ICRU report 78 (ICRU 78, 2007) offers a method
for dealing with the Relative Biological Effectiveness (RBE) of protons when it comes to
determining the best estimate of RBE to be used in any calculation relating proton and x-ray
doses. A generic RBE value for protons for all therapeutic energies (60-250 MeV) of 1.10
is used. This value is derived from RBE values for a variety of in vivo systems using animal
tissues and organs. 1.10 is the mean value for all proton energies, dose levels and tissue
systems (ICRU 78, 2007). There is, however, a related range uncertainty owing to the RBE
of proton beams, which is beyond the scope of this thesis. For more information on RBE the
reader is referred to the literature (Carabe et al., 2012; Grassberger et al., 2011; Matsuura
et al., 2010; Robertson et al., 1975).
2.1.3 Discrepancies Between Planned Dose and Delivered Dose
Despite algorithms being able to calculate range in the presence of density heterogeneities,
range uncertainties can be introduced by geometric changes in the position of density het-
erogeneities relative to the proton beam by set-up uncertainties and patient motion (Deasy,
1998).
Motion
In many treatment sites, organ motion has to be considered and incorporated into the plan-
ning and delivery workflow. Organ motion includes inter-fraction motion (between each
fraction) and intra-fraction motion (during the treatment delivery). Patient motion influ-
ences the position of the inter-fractionally moving organ and intervention in the form of
immobilisation and image guidance for precise bony or soft tissue anatomy setup at the be-
22 Uncertainties in Proton Therapy
ginning of each fraction is required to optimise treatment delivery and minimise CTV to
PTV margins. These are well-established techniques in radiotherapy (Burnet et al., 2010).
Intra-fraction motion includes both inter-field (i.e. between each field) and intra-field (i.e.
during the field) motion. Intra-fraction motion changes during the delivery over the time
period of seconds and minutes and includes anatomical changes such as bowel movements,
respiration and heartbeats. Studies have been carried out to determine the magnitude of
inter- and intra-fraction motion at different tumour sites (Booth and Zavgorodni, 1999; Lan-
gen and Jones, 2001; Rimmer et al., 2008).
With proton therapy, geometric changes caused by motion can also result in density changes,
and therefore a change in radiological path length, along the beam path. In X-ray therapy,
the dose distribution changes by only a few per cent, owing to density changes. However,
their influence in proton therapy can result in severe under-dosage of the CTV and over-
dosage of OARs and normal tissues distal to the target (Lambert et al., 2005). This effect
is the same for both scattering and scanning delivery techniques. In addition, for active
scanning, the major effect of intra-field motion is interplay, which relates to motion, usually
respiratory motion, with a frequency similar to that of the scanned beam, and which can lead
to over- and under-dosage in the target volume (Knopf et al., 2011; Phillips et al., 1992; Seco
et al., 2009; Zenklusen et al., 2010).
Lomax (2008a) investigated the effects of inter-fraction, intra-fraction and inter-field motion
for both 3D and DET IMPT treatment plans. It was reported that, for a 5mm shift in the dose
distribution, an under-dosage of up to 20% can occur in the CTV when plan optimisation for
maximum OAR sparing is used. Treatment deliveries involving high dose gradients that rely
on matching contributing fields are very sensitive to any changes in position between deliv-
eries of each field. An important conclusion from this paper is that for certain IMPT plans a
simple PTV margin cannot be applied to compensate for inter-fraction motion. Further in-
vestigation into the management of uncertainties and more assessment for IMPT treatments
at the treatment planning level are needed. Lomax (2008a) also described the greater sensi-
tivity protons have to density heterogeneities owing to their physical characteristics and that
these could accentuate motion uncertainties. Without a PTV, there is no method for record-
ing dose to a moving CTV or for evaluating plans through the use of DVH’s and for IMPT
plans, no other compensation method to ensure that the CTV is covered. After comparing
the effects of geometrical and density heterogeneities for both DET and 3D IMPT it was
found that DET plans are very sensitive to motion uncertainties and considerable changes in
2.2 Managing Uncertainties Including Positional Discrepancies and Motion 23
the dose distribution were found. The reasoning behind this is that internal dose gradients
in the individual fields in DET can cause large variations in dose within the target volume
when mismatched. This was observed for inter-field motions, but would also be an issue
when intra-field motions were present.
Simulations have also been carried out by Lambert et al. (2005), albeit only in homogeneous
geometries, to investigate how interplay affects the dose distribution in the CTV. Lambert
et al took the ICRU 50 (ICRU 50, 1993) recommendations for PTV dose homogeneity of
95– 107% as threshold. Their results showed that in extreme cases up to 100% of the target
volume received doses below that recommended by ICRU 50 and with a minimum dose as
low as 34% of the prescribed dose. These results were backed up by simulations carried out
by Grözinger et al. (2006) and experimental work by Bert et al. (2008). Bert et al. (2008)
carried out the first patient simulation that confirmed under-dosage using four-dimensional
computed tomograph (4DCT) lung data. Despite using margins that consider the effect of
the changing radiobiological path length, adequate CTV coverage could not be achieved.
2.2 Managing Uncertainties Including Positional Discrep-
ancies and Motion
The management of uncertainties is critical to successful radiotherapy. Current methods
of reducing geometric uncertainties include immobilisation, in-room imaging, image guid-
ance, planning from 4D CT and gated radiotherapy for respiratory motion. These are all
methods that are routine in X-ray therapy and are now being applied to proton therapy (Bert
and Durante, 2011).
2.2.1 Motion Mitigation
Two main methods of motion mitigation have been developed for active scanning to miti-
gate the effects of interplay: rescanning (repainting) and beam tracking (Bert and Durante,
2011).
Rescanning
Typically in proton therapy, multiple dose painting is required to deliver the prescribed dose
distribution to each layer of the target volume. In IMPT the dose is delivered using multiple
24 Uncertainties in Proton Therapy
field directions and with a range of different proton energies to produce a uniform dose dis-
tribution throughout the target volume. To increase dose conformity, steep dose gradients
are used at the target border and field edges. This increases the complexity of the fluence
maps per field and therefore makes the plan less robust to uncertainties. Intra-fraction mo-
tion can lead to an under and overdose pattern that is dependent on the motion parameters
(initial phase, period and amplitude) and the speed of the scanning process or direction of
scanning. By rescanning the PTV several times per treatment fraction an averaging effect
of the over- and under-dose pattern can be achieved. As long as the intra-fraction motion
changes between each rescan, so that there is no synchronisation between delivery and or-
gan motion, a homogeneous dose distribution to the CTV can be achieved with a “blurred”
dose distribution in the region of the margins (Bert and Durante, 2011). There are two main
types of rescanning:
■ rescanning by energy slice, also known as slice-by-slice, level painting or non-volumetric
rescanning.
■ rescanning of volume, also known as volumetric rescanning or uniform repainting.
The problem of organ motion and rescanning synchronism has also been tackled by several
groups, including Furukawa et al. (2007) and Seco et al. (2009). Solutions include:
■ using random modulations e.g. a change in scan speed
■ repainting energy slices in different orders (random repainting)
■ random delays between repaints (time delay)
■ change in scan paths between two rescans
■ use of data from motion monitoring systems in combination with modulation dose
rate, either as phase-controlled rescanning or as breath-sampled rescanning.
Data from phase-controlled rescanning and breath sampled rescanning show that uniform
spreading of rescans over the motion of the breathing cycle leads to more robust treatment
delivery, requiring fewer rescans than the other methods for the same level of homogeneity
in the CTV. There are also two methods of delivering rescanned treatments (Zenklusen et al.,
2010):
■ Scaled rescanning – this is delivering each rescan with a proportionally reduced dose
per scan and is the most typical method.
2.2 Managing Uncertainties Including Positional Discrepancies and Motion 25
■ Isolayered rescanning – where a specified number of protons per spot are delivered
at each scan, which leads to a different spot position in each rescan. This is because
some spots will have received sufficient dose from previously completed scans.
Beam Tracking
In beam tracking the motion of the CTV is monitored and compensated for so that a PTV
margin is reduced (Bert and Durante, 2011). Owing to the need to compensate for both
CTV motion and the change in radiological path length in proton therapy, the German na-
tional heavy-ion physics laboratory, the Gesellschaft für Schwerionen, in Darmstadt, has
established a method of using a motor-driven compensation system for changes in radiolog-
ical path-length (Bert et al., 2007; Saito et al., 2009) and raster scanning (Furukawa et al.,
2007).
2.2.2 Margins
The widely used CTV–PTV margin recipe, derived by van Herk (2004), is used to ensure
that 90% of patients have CTV coverage of at least 95% of the prescribed dose. The static
PTV represents the moving CTV and is therefore a useful method of evaluating a plan, by
using DVHs to report the minimum dose to the CTV. In proton therapy, if conventional PTV
margins do not produce plans robust to uncertainties, another method of ensuring confidence
in a planned dose distribution representing the delivered dose distribution is required. This
can be achieved by increasing margins or by smearing the proton range with a compensator
in the case of passive scattering. However, this will lead to a reduction in dose conformity
in the plan. The concept of a safety margin also partly fails for set-up uncertainties because
they not only shift the dose distribution but also change the range where there are density
changes in the beam path, leading to a distorted dose distribution (Grözinger et al., 2008).
In proton therapy, any concept used to compensate for organ motion-generated uncertainties
must include both geometric motion and the influence it has on the beam range, as this can
have a severe dosimetric impact (Grözinger et al., 2008).
PTV margin sizes in proton therapy have been investigated by Thomas (2006). In many
comparative studies of achievable dose distributions between protons and X-rays, the CTV
and PTV margin sizes are the same for both modalities. The CTV and OAR volumes are the
same for any treatment modality. However, the size of the margin between CTV to PTV and
OAR to PRV is modality dependant. Thomas (2006) quantified margins required for PTVs
26 Uncertainties in Proton Therapy
and PRVs in proton therapy for anterior single, anterior–posterior parallel opposed and four-
field brick proton beam arrangements. Each type of systematic and random uncertainty was
considered and the effect they had in geometry and range was discussed, as was the modality
dependency. The key message is that isodoses defined by lateral edges, for instance, can be
treated in the same way as for X-rays, whereas isodoses defined by the distal edge will
have a different uncertainty arising from inaccurate electron densities derived from CT data.
Margins were calculated for a head and neck (H&N) and a prostate case. For the H&N case
the CTV–PTV margins in all three planes for the single field and parallel opposed fields
were smaller than those required by X-rays; however, they were greater for the four-field
brick plan. This same pattern was observed for the prostate case. This is because a smaller
margin size is needed in the anterior–posterior direction for single and parallel opposed
(3mm compared with 10.5mm for the H&N case), as set-up uncertainties and motion in this
direction will not affect the dose distribution. Beam specific margins were further explored
by Rietzel and Bert (2010) and Park et al. (2012). Park et al. (2012) showed a method
of creating beam specific PTV margins by ray-tracing and shifting ray lines to account for
tissue misalignment. They found that plans created using them provided better coverage in
the presence of setup and range uncertainties.
Albertini et al. (2011b) addressed the issue of whether safety margins in proton planning
are necessary . In this paper two types of treatment delivery were investigated: single field
uniform dose (SFUD) plans and IMPT plans. Plans included:
(i) 3 field SFUD plan to the PTV.
(ii) 2 field IMPT plan delivering uniform fields to the PTV.
(iii) 4 field Non-uniform field IMPT plan, with strict constraints on OARs, planned to the
PTV.
(iv) Non-uniform field IMPT plan, with strict constraints on OARs, planned to the CTV.
The robustness of each plan to random set-up and systematic range uncertainties was com-
pared. Robustness was determined using the concept of error-bar dose distributions, by
shifting dose distributions relative to the expected uncertainties and then displaying a final
error-bar dose distribution and error-bar volume histogram; therefore, it was representative
of the possible discrepancies in dose between planning and delivery. The results from this
work show that, for uniform field deliveries [(i) and (ii)], the use of margins improved the
plan robustness, where <5% of the CTV contained dose uncertainties >10%. However, in
highly complex non-uniform IMPT plans [(iii) and (iv)], margins improved robustness only
2.2 Managing Uncertainties Including Positional Discrepancies and Motion 27
slightly. For 5% of the CTV, uncertainties of up to 55% were observed. For plans (iii) and
(iv), steep dose gradients existed within the target, leading to uncertainties within the target.
Because margins can help CTV coverage only at the edges of the target volume and not in
the centre, they have little effect on plan robustness when steep dose gradients exist within
the target volume (Albertini et al., 2011b). Interestingly, it was found that complex IMPT
plans were robust to OARs. It was concluded that there is a need for more sophisticated
methods for taking into account uncertainties in highly modulated IMPT plans, such as in-
cluding them in the optimisation. This work assumed that the set-up uncertainty was indeed
random and considered only the systematic component of the range uncertainty. It did not
include the effect on the dose distribution from motion uncertainties.
2.2.3 Probablistic Planning and Robust Optimisation
An alternative to using margins is being developed by incorporating uncertainties directly
into the optimisation algorithm, as proposed by Unkelbach et al. (2007) and Pflugfelder
(2008) to achieve a high quality plan with a high probability of being deliverable (Chu,
2006). Due to the number of spots in 3D IMPT required to fill the target an optimisation
algorithm is required to determine the optimal solution for their individual weights. Intensity
modulation optimisation is highly degenerate (Albertini, 2011; Webb, 2003). There exist
many dosimetrically equivalent solutions to the given problem of delivering a homogenous
dose to the target whilst ensuring the dose to OARs are within predefined constraints. Some
of these solutions may offer greater plan robustness and therefore a greater chance of cure
whilst limiting the chance of normal tissue complications. This “degeneracy” of solutions
can be used to reduce the sensitivity of the plan to uncertainties if they are incorporated into
the optimisation process (Pflugfelder, 2008).
Probabilistic planning assumes prior knowledge of the probability distribution of the uncer-
tainty; in most cases, the distribution was assumed to be normal. The probabilistic approach
works on the principle that the total delivered dose to a voxel is random due to its depen-
dency on both the dose delivered at each fraction and the location of the voxel during that
fraction (Chu, 2006). It can be assumed that the total dose to a voxel is the sum of the
N random variables, where N is the number of fractions. The total dose can therefore be
approximated as a Gaussian distribution by following the central limit theorem. An approx-
imation can be made then that the N random variables, which give the dose per fraction, are
independently and identically distributed (Chu, 2006). The use of this approximation allows
for the solution to be less computationally expensive by computing the mean and variance of
28 Uncertainties in Proton Therapy
the dose from a single fraction to form a probability distribution (Chu, 2006). However, this
approximation relies on the assumption that patient changes over time, for example weight
gain or loss, or ‘relaxation creep’, are negligible to the uncertainty problem. To address this
Bangert et al. (2013) use a method of analytical probabilistic modelling (APM) to include
complex correlation models of the underlying uncertainties. Robust optimisation effectively
optimises the worst case that may occur, assuming that the proton range may vary within
some known interval (Unkelbach et al., 2007). An example of a robust optimisation algo-
rithm is the worst-case dose distribution introduced by Lomax et al. (2004) and is a method
of combining multiple dose distributions into a single one. Uncertainties can be applied to
the CT data, for example a percentage change in HU or several isocentric shifts, and the
dose calculated. For a voxel inside the target volume the minimum dose to this voxel is
stored and for a voxel outside of the target volume the maximum dose is stored. Minimis-
ing these values results in the optimal worst-case plan. Each voxel is treated independently
so the worst-case dose distribution is an “unphysical” one. Although the worst-case dose
distribution is unphysical it can be used as a lower bound for the worst quality treatment
plan. A best-case plan can be seen as the upper bound in the same respect, but to the best
achievable plan quality. In this optimisation it is assumed that the range uncertainties for
each Bragg peak are correlated, so that at the target the range uncertainty is accumulated and
effects within the target are ignored. Both methods, probabilistic and robust planning, pro-
duce plans less sensitive to uncertainties, and mathematically each approach yields similar
results as shown by (Chu, 2006).
Unkelbach et al. (2007) investigated how incorporating uncertainties into the IMPT optimi-
sation yields qualitatively different treatment plans compared to conventional plans which
do not account for uncertainty. Both probabilistic and robust methods were used to include
uncertainties into the plan optimisation to allow for a comparison to be made. The dose
from each beamlet was calculated at multiple ranges; three ranges between 2 and 5mm
were used, and these were the nominal, maximum and minimum range uncertainties (Lo-
max et al., 2004). The set-up uncertainty was modelled as a shift of the target inside a
sphere with a radius equal to the maximum setup uncertainty. Both approaches yielded sim-
ilar results and both fell into a multi-criteria problem - that plans had greater robustness to
uncertainties, but at a cost to conformality. They concluded that both methods led to treat-
ment plans that were less sensitive to range variations and that both plans were of a similar
quality. This was achieved by using the lateral edge instead of the distal edge to shape the
dose distribution at the transition between the OAR and the tumour for both methods.
2.2 Managing Uncertainties Including Positional Discrepancies and Motion 29
Pflugfelder (2008) applied the worst-case optimisation method to real patient data where
the target surrounded the spinal cord. Uncertainties of the beamlet ranges of ±5mm were
considered. The range uncertainty was sampled at three positions: the nominal range, the
maximum range and the minimum range. Using this optimisation process, the plan’s sensi-
tivity to range uncertainties was decreased. However, this was achieved by compromising
the dosimetric quality of the plan by:
(i) using the lateral instead of the distal edges of the Bragg peaks to shape the dose gradients
between the target and the OAR
(ii) adding a safety margin automatically at the distal field edge for each treatment beam
(iii) flattening the dose profile in depth for each treatment beam compared with the nominal
plan
For the tumour, the largest deviations between delivered and prescribed doses corresponded
to an under-dosage close to the OAR, whereas the deviations in most other parts of the
tumour were small. The maximum doses delivered to OAR voxels reached approximately
80% of the prescribed dose (Pflugfelder, 2008). The DVHs from each method showed
very similar results. Pflugfelder (2008) showed comparable results to those of Unkelbach
et al. (2007) using their worst-case optimisation function (Pflugfelder, 2008). When using
the worst-case optimisation incorporating both setup and range uncertainties, the resulting
change to the dose distribution was such that the distal dose gradient of the treatment beam
was smoothed.
Chen et al. (2012) introduced a method of including robustness into a multi-criteria opti-
misation (MCO) framework for IMPT. The concept of multi-objective Pareto optimisation
(often known MCO) has been introduced into radiotherapy treatment planning to overcome
these problems (Cotrutz et al., 2001). A Pareto optimal treatment plan is not a single plan
but a database of plans where each one represents a Pareto optimal solution which cannot
be improved without worsening at least one other parameter (Thieke et al., 2007). In this
case the planner can navigate through the pre-calculated database of Pareto optimal plans
and visualise in real time the trade-offs for each case (Craft et al., 2005). MCO allows for
the plan which strikes the best balance between different objectives to be selected from a
Pareto front.
Chen et al. (2012) suggested that the trade-off between robustness and dosimetric quality
in IMPT can be investigated using MCO. An example is shown in Figure 2.2. A planning
30 Uncertainties in Proton Therapy
objective for robustness represents the worst case realised for any uncertainty scenario and
so can provide a measure of plan robustness. The MCO method used by this group was
based on a linear projection solver using minimax optimisation, meaning that the objective
was minimised for the worst-case uncertainty scenario and was used to investigate plan ro-
bustness to uncertainties. The computing time required for the robustness optimisation for
each Pareto optimal plan was 5 minutes. Range uncertainties were modelled using an over-
shoot and undershoot scenario of ±3% (as used by Lomax (2008b)). The set-up uncertainty
was modelled using discrete rigid shifts of the patient with respect to the isocentre. A base
of skull tumour was planned twice: once using robust optimisation (but without MCO) and
once with margins. A chordoma case was planned with MCO as an ideal example of the
trade-off between brain stem sparing, CTV coverage and robustness. In the skull base case
the DVH of the CTV for both the robust plan and the margin plan showed similar coverage
for the nominal plan. However, once plans with uncertainties were introduced, the cover-
age was worse for the margin case. In consideration of brain stem sparing the robust plan
ensured that in all uncertainty cases the dose was limited to 60Gy; this was not seen in the
margin plan. In the chordoma case the Pareto surface of objectives allowed the planner to
have greater control when deciding between a robust plan, a conformal plan or somewhere
in between.
2.2.4 Treatment Analysis and Validation
For all areas in radiotherapy, validation and quality assurance are routinely carried and are a
necessity for safe, effective treatment. Phantom work and technical assessment are carried
out when a new technology or method is being implemented, such as motion mitigation
techniques. However, patient-specific quality assurance is also required, and this starts at
the decision to accept a given plan design. In cases where margins are in use, the PTV
is used as a tool for reporting the minimum dose to the CTV. If the PTV method is to be
discarded then a new method for recording dose to a moving CTV, and for evaluating plans
through the use of DVHs, is needed. As discussed above, Albertini et al. (2011a) have de-
vised a possible solution to this problem by determining plan robustness for set-up and range
uncertainties using an error-bar dose distribution method. This method was later validated
experimentally using a customised anthropomorphic phantom based on a diagnostic head
phantom, and GafChromic H EBT2 batch F100070903B film, to carry out patient-specific
quality assurance under realistic conditions and with deliberately introduced uncertainties.
A large part of ensuring optimum treatment delivery to the patient is verifying the accuracy
2.2 Managing Uncertainties Including Positional Discrepancies and Motion 31
Fig. 2.2 Illustrating two parameters for optimisation, conformity and robustness, and how
a Pareto front can be used to investigate clinically achievable plans with the same Pareto-
optimised solution, but with different values of importance on robustness and conformity.
of the dose delivered and comparing it with that predicted by the planning system. In cases
where uniform fields are being delivered it is enough to experimentally measure the individ-
ual field doses, and in the case of proton therapy, the dose in homogeneous water (Albertini
et al., 2011a).
The highly modulated nature of IMPT means that ever more complex methods of plan ver-
ification are being developed to cope with the advancing treatment technology to ensure
safe patient treatment. Safai et al. (2004) described how a scintillation dosimetry system
was used to verify such proton treatments. The problem with this method was that, even
though it could accurately measure the dose delivered, the dose from the TPS needed to
be recalculated in a homogeneous material to be compared with the measured results (Al-
bertini et al., 2011a). This approach of patient-specific plan verification is not adequate
when fields of inhomogeneous dose distributions are being applied. Each inhomogeneous
dose distribution will, in combination, achieve a uniform dose in the target volume, but any
32 Uncertainties in Proton Therapy
slight misalignment of a steep dose gradient could lead to a severe under- or over-dosage
in comparison with the plan (Albertini et al., 2011a). This effect could be worsened by the
presence of density heterogeneities in the patient (Albertini et al., 2011a). By using the cus-
tomised phantom, Albertini et al. (2011a) have been able to measure the accuracy of dose
and the effect of spatial and range uncertainties by deliberately introducing uncertainties.
Setup uncertainties were introduced by moving the phantom by known amounts, and range
uncertainties simulated by modifying the HU values by ±3% (Lomax, 2008a). From this
work it was found that 3D IMPT plans were more robust than DET plans. The theory is that
there are fewer Bragg peaks used in DET, whereas in 3D IMPT more spots are used (about
180 for DET compared with about 1500 for 3D IMPT in the example used), each with a
lower weighting than those in DET, so that that any misalignments would have a greater
impact on the DET plan.
2.3 Discussion
To produce optimal plans using protons, knowledge of the proton range in the patient is
essential. There are many factors that can contribute to range uncertainty. HU uncertain-
ties can contribute to approximately ±3% uncertainty in range even after site-specific CT
scanner calibrations have been carried out. Owing to the steep dose gradients that can be
achieved at the edges of, and within, the target volume, precise field matching is required
to prevent over- or under-dosage in the target. Simulation and patient data examples have
been used to show that the use of a PTV margin does not ensure uniform dose coverage
to the CTV for complex IMPT plans. Despite these results, there is still a need for facility
based and treatment protocol-specific simulations to be carried out to achieve quantitative
assessments for different tumour sites (Bert and Durante, 2011). Parameters such as dose
regime, scan path pattern and beam extraction rate need to be included because they will
affect the resulting dose distribution.
Methods of improving plan robustness to range uncertainties beyond the use of margins are
being developed for use in complex IMPT plans to ensure CTV coverage. However, all
of these methods will decrease the achievable conformity of a plan. Robust optimisation
using the worst-case scenario will produce a robust plan, but at the cost of sacrificing con-
formity by placing a lateral edge instead of the sharp distal edge to shape the dose between
the target and the OAR. The use of robust MCO gives the planner greater control over ob-
jective weights. This is important because there is little evidence or known experience in
2.3 Discussion 33
determining what type of plan, conformal or robust, is most beneficial for the population of
patients.
The difficulty with deliberately introducing uncertainties into either the optimisation or the
plan evaluation is that all probability distributions are assumed known. Albertini et al.
(2011b) and Chen et al. (2012) only considered range uncertainties created in the range
calculation arm in Figure 2.1. The assumption is then that any change in radiological
path length caused by set-up discrepancies or motion, whether during a fraction or between
fractions, has a negligible effect on the range. Motion uncertainties will be dependent on
anatomical location and setup uncertainty will be protocol specific, based on what IGRT
schedules are in place and the type of immobilisation used. More research is needed in
quantifying range uncertainties for different anatomical locations so that they can be imple-
mented into the robust optimisation of dose modelling. This would include investigating
range uncertainties associated with different IGRT schemes and types of motion.
In the absence of a PTV, novel methods to evaluate dose coverage have been developed such
as the error-bar DVH and the use of worst-case and best-case optimisation to give upper and
lower bounds on a DVH. A challenge for all these methods will be how to carry out adequate
and efficient patient-specific verification. The data available from non-patient heterogene-
ity studies have shown that rescanning techniques produce improved results compared with
beam tracking, but further investigations using real patient data are required (Bert and Du-
rante, 2011). There is no information in the literature regarding patient-specific validation
for rescanning and beam tracking methods because they are not yet in routine clinical use.
However, the group at PSI did publish the first experimental results of motion mitigation by
discrete spot (Schätti et al., 2013) and continuous line (Schätti et al., 2014) rescanning in
2013 and 2014 respectively.
The question remains as to whether adaptive therapy can become an integral part of a proton
therapy treatment, to allow re-optimisation during the course of a patient’s treatment and if
so how this can be achieved. Adaptive radiotherapy is the adjustment during the treatment
course of the parameters initially chosen at planning, in order to re-optimise the treatment
as a direct result of unavoidable changes in the patient. In proton therapy, the use of adap-
tive therapy may prove to be valuable for ensuring the delivered dose matches the planned
dose at each fraction. In this context, access to high quality daily volumetric imaging, fast
dose recalculation algorithms and computing power will be required. This may provide the
opportunity for greater individualisation of the patient’s treatment. Use of a multi-objective
34 Uncertainties in Proton Therapy
Pareto optimisation function could allow plans to be recalculated when necessary using dif-
ferent weights for robustness and conformity, incorporating image guidance data to optimise
the plan.
It has been shown that a PTV may not be the best solution for ensuring target coverage in
the case of complex IMPT plans. This is the result of the sensitivity of the proton range
to density heterogeneities and the achievable steep dose gradients inside the volume. Many
groups have been working towards possible solutions for achieving dosimetric quality and
plan evaluation without a PTV. The solutions described here show great promise for opti-
mising proton therapy planning. Based on the literature discussed here the purpose for the
following work in this thesis was determined. The following chapters explore the effect of
start conditions on plan robustness and look to develop methods of plan robustness analysis,
and to quantify range uncertainties for different anatomical locations due to patient motion.
Chapter 3
A Method for Quantifying Random
Range Uncertainty
Techniques such as probabilistic planning, robust optimisation, and methods of simulat-
ing uncertainties to analyse plan robustness are becoming more widely used. For example
these include, the worst case optimisation (Lomax et al., 2004; Pflugfelder, 2008), used in a
research capacity, and the RayStation TPS 4.5 which supports robust optimisation commer-
cially for both photon and proton planning (Raysearch Laboratory, 2014). It is becoming
evident that there is a need for a greater understanding of the magnitude and probability of
range uncertainties for the population.
In this chapter a method for acquiring centre- and site-specific prior knowledge of random
range uncertainties for the population is provided. In Chapter 4 patient examples are given.
3.1 Introduction
It has been previously established that there are many sources of uncertainty in proton ther-
apy. These include both the systematic and random components arising from the uncertainty
in HU, dose calculation algorithms, daily setup, target delineation and organ motion (see fig-
ure 2.1). Within conventional radiotherapy, IGRT provides information about the systematic
uncertainty, allowing for this uncertainty to be corrected and for the possibility of reducing
margins. Using IGRT online can allow for inter-fraction changes to be partly compensated
for by applying rigid couch shifts before each fraction. The same technology is available in
proton therapy. Though margins may not be an adequate way of providing CTV coverage
(to a pre-designed probability) for complex proton plans, the same methods can be imple-
36 A Method for Quantifying Random Range Uncertainty
mented for calculating the prior information required as an input into the robust optimisation
and for models designed to analyse plan robustness.
The range uncertainty problem may be relatively new to radiotherapy, but the uncertainty
problem is not. Previously the ICRU recommendations for uniform prescribing and report-
ing (ICRU 50, 1993; ICRU 62, 1999; ICRU 83, 2010), as well as the margin recipes (Stroom
et al., 1999; van Herk et al., 2000), came as a solution to solve the sensitivity of treatment
plans to geometric uncertainties. The margins are grown based on known systematic and
random uncertainties inserted into the recipe. In deriving his margin recipe, van Herk and
colleagues (van Herk et al., 2000) started from a well-defined probabilistic criterion; that
the CTV-to-PTV margin should be grown to ensure the CTV minimum dose is greater than
or equal to the planned PTV minimum dose, for 90% of patient treatments (van Herk et al.,
2000). The limitations of the PTV in complex proton cases have been discussed and sev-
eral authors (Chen et al., 2012; Fredriksson et al., 2011; Liu et al., 2012; Pflugfelder, 2008;
Unkelbach et al., 2007) have addressed this head on by developing probabilistic and robust
optimisation algorithms for use in radiotherapy. To clarify, the probabilistic approach con-
siders all possible range variation scenarios as a term in the objective function by weighting
each with its probability of occurring. This is represented by optimizing the expected value
of the objective function. The robust approach does not perform a weighting of each dif-
ferent scenario, but instead minimises the maximum deviation of delivered and prescribed
dose, which can occur for all allowable ranges (Unkelbach et al., 2007). Through using
these methods, the treatment planner is no longer responsible for adding margins. Instead,
the TPS builds its own “margins” into the dose distribution, governed by probabilistic plan-
ning criteria or a lower boundary with a predefined probability of occurring. (Gordon et al.,
2010). These approaches appear to be the correct way to deal with uncertainties when plan-
ning complex proton therapy treatments and it should be noted that treatment plans obtained
by the probabilistic or robust methods cannot be reproduced by a safety margin approach
(Unkelbach et al., 2007).
Though finding, potentially, a limited place for PTVs in complex proton therapy the quan-
tification of range uncertainties is required as input into robust optimisation (Gordon et al.,
2010; Unkelbach et al., 2007). It has been shown how systematic setup uncertainties can be
measured for each centre and site as described by Bolsi and colleagues (Bolsi et al., 2008),
but there is still a challenge in quantifying the random range uncertainty, and the associated
probability arising from patient motion. In most cases robust optimisation and probabilistic
3.1 Introduction 37
planning approaches had been modelled relying on rigid isocentric shifts in CT data sets or
from uniform percentage changes in HUs. In Unkelbach’s paper (Unkelbach et al., 2007)
values of±5mm and±2.5mm were used to represent the systematic range and setup uncer-
tainties respectively, random errors were not taken into account. These assumptions on the
amount of range uncertainty are due to the lack of measured data to support more sophis-
ticated models (Unkelbach et al., 2007). Unkelbach et al. (2007) suggested that for future
studies tighter bounds on the actual magnitude of range uncertainty could be added with the
aim to derive more precise uncertainty models for specific tumour sites and planning pro-
tocols. Gordon and colleagues (Gordon et al., 2010) suggest estimates could be obtained
from a population of patients or, in adaptive therapy, from patient-specific data. PTV-based
and probabilistic/robust plans are equally dependent on sufficiently accurate knowledge of
the distribution of the geometric uncertainty (Gordon et al., 2010). In proton therapy this
statement is equally true of the range uncertainty. In both conventional radiotherapy and
in proton therapy, if the uncertainty distribution model is inaccurate, both types of plan,
whether PTV or probabilistic, will produce underdosing of the target and/or overdosing of
OAR(Gordon et al., 2010).
In this work the focus is on determining the random range error from inter-fraction motion
after online IGRT based corrections have been applied. The clinical Mega-Voltage Com-
puted Tomography (MVCT) data has been used offline along with the shifts recorded by
TomoTherapy for each fraction to provide the relevant information. Using the shifted 3D
MVCT data and structure sets (outlined by the doctor at planning) the Water Equivalent
Path-Length (WEPL) within the CTV can be calculated. By taking the difference in WEPL
between the first and consecutive fractions an indication of inter-fraction range uncertainty
for a given set of beam directions can be obtained. When carried out for a large sample of
patients an indication of the population probability distribution of WEPL changes that oc-
cur between fractions may be gained. This information would be specific to the anatomical
location, immobilisation and method of IGRT. Using this method, probability distributions
can be generated and used in either the robust optimisation or at plan analysis, or in both. By
calculating this information online, action levels for carrying out dose difference analysis
can be established allowing for adaptive planning.
38 A Method for Quantifying Random Range Uncertainty
3.2 TomoTherapy
TomoTherapy, shown in Figure 3.1 is an integrated helical IMRT and IGRT system consist-
ing of a compact in-line 6MV accelerator mounted on a CT ring gantry with a detector array
(xenon gas ion chamber) opposite. It provides volumetric, fan-beam megavoltage (MV) CT
imaging on the treatment couch (Langen et al., 2005). Since 2007 at Addenbrooke’s Hospi-
tal patients have been imaged daily on TomoTherapy using MVCT with rigid couch shifts
and roll applied to correct for set-up and inter- fraction motion.
Fig. 3.1 TomoTherapy at Addenbrooke’s. Provides daily imaging with positional correction
on all patients (Burnet et al., 2010).
3.3 Calculating WEPL Changes in the CTV from TomoTher-
apy Shift & MVCT Data
In proton therapy it is essential for dose calculation that accurate conversion of HU to rel-
ative proton stopping power is carried out. In this work proton dose will not be calculated,
yet it will still be required to establish the relationship between the CT data and the water-
equivalent path length for protons to account for density variations between fractions. The
method used to determine this relationship between HU and relative proton stopping power
is explained so that WEPL calculations may be made from clinical TomoTherapy MVCT
data.
3.3 Calculating WEPL Changes in the CTV from TomoTherapy Shift & MVCT Data 39
The WEPL of each pixel in the beam direction is calculated using an effective depth algo-
rithm. This means that high density pixels will have a larger path-length than that for water
and low density pixels will have a shorter path-length. Therefore, confidence in the con-
version between HU and proton stopping power is required. In this way the trajectory of a
proton through the CT is transformed from the CT system into a water-equivalent system in
the beam’s-eye view. This then transforms the original target volume shape into a distorted
one in the water-equivalent system. Figure 3.2 is an example to illustrate the process in a
prostate MVCT scan. A 90◦ beam direction was used when converting from CT system
(right image) to the water equivalent one (left image). The orange area in this image is the
mask used to represent the prostate.
Fig. 3.2 Example of how CT data (left image) can be converted into WEPL (right image) at
a beam angle of 90◦ using calibration curves and ray-tracing
There is no direct functional dependence between HU and WEPL. The HUs reflect the at-
tenuation of X-rays, mainly by the Compton scattering, whereas ion energy loss processes
dominate the proton path length. The actual chemical composition of the material (or tissue)
affects both quantities, but with a different functional dependence on the Atomic number,
Z, of the contributing atoms. Hence the desired relation cannot be retrieved from the knowl-
edge of the HU alone. Until work investigating proton tomography (Schulte et al., 2004) is
complete the option left is to carry out CT-scanner and site-specific CT calibration.
40 A Method for Quantifying Random Range Uncertainty
A Quick Note on Stopping Power
Charged particles interact primarily by ionisation and excitation and deposit dose directly.
These interactions are mediated by the Coulomb force between the electric fields of the in-
coming proton and the orbital electrons and nuclei of atoms in the medium (Khan, 2010).
Collisions between the proton and orbital electrons result in excitation and ionisation of the
atom by raising an electron to a higher shell or, if provided with sufficient energy, the elec-
tron will leave the atom. If this happens the electron may have sufficient kinetic energy to
go on to cause further ionisations. These electrons are known as δ rays which are respon-
sible for the microscopic structure of energy deposition around the primary track. Atomic
interactions occur so frequently that energy loss appears continuous. The rate of loss of
energy per travelled distance for charged particles is known as stopping power (S), which is
the change in energy with distance, S = dEdx . However, this equation does not hold unless we
assume that all the dose is deposited at the site of interaction and ignore dose carried away
by δ rays. The stopping power is proportional to the square of the particle’s charge and
inversely proportional to the square of its velocity, q
2
v2 (Khan, 2010). As the proton slows
down, the rate of energy loss increases and therefore so does the amount of dose deposited.
The dose increases until near the proton range where it rises steeply at the Bragg peak be-
fore falling steeply to zero. This is described using the Bethe-Bloch formula (equation 3.1),
which can be used to determine the mean range of the particle. As defined in ICRU 78
(2007) the distal fall off is the distance in gcm−1 in which the dose decreases from 80%
to 20% of the maximum value when measured along the central axis in water. The range
is defined as the dpth in gcm−1 along the the beam central axis is water to the distal 90%
point of the maximum dose value. There is also a mass stopping power, Sρ , known as total
stopping power and is expressed in MeV cm2g−1 which takes into account the mass density
of the medium and includes energy losses from both collisional and nuclear interactions.
Calibrating CT Scanners
There exists a significant source of uncertainty in the HU conversion arising from the use of
tissue substitutes, that have different chemical compositions to that of real tissue, to calibrate
scanners. This leads to different values of electron density as the oxygen, carbon, hydrogen
and calcium contents of these substitutes are changed to allow them to be usable often for
radiobiological purposes (Schneider et al., 1996). To address this challenge, many centres,
including PSI in Switzerland and Orsay in France, have adopted a method of stoichiometric
3.3 Calculating WEPL Changes in the CTV from TomoTherapy Shift & MVCT Data 41
calibration as described by Schneider et al. (1996) to calibrate their CT scanners.
In a stoichiometric calibration both real tissue and tissue substitutes are used to predict
the HUs for human tissues (Schneider et al., 1996). Schneider et al. (1996) verified this
method of calibration for protons using a sheep’s head and radiographic film measurements.
The relative stopping power through the sheep’s head was calculated using both the con-
version established from the real tissue samples and the substitute tissues. The comparison
in measurements and calculations showed that the stoichiometric calibration was of greater
precision than the tissue substitutes for proton beam therapy. Schneider et al. (1996) also
advises that this method could improve on the precision for X-ray radiotherapy treatment
planning.
Calibrating the TomoTherapy Imaging Beam
The TomoTherapy imaging beam also supplies a second use beyond daily correction: it
provides a method of planned adaptive radiotherapy. At Addenbrooke’s, Dose Difference
Analysis (DDA) is routinely carried out for TomoTherapy patients to create a comparison
of the planned dose and the delivered dose distributions for a fraction of treatment using the
TomoTherapy planned adaptive module. The delivered dose is calculated by applying the
delivered sinogram (describing MLC opening with time and gantry angle) to the MVCT im-
age taken for the specified fraction. The dose distribution from the fraction and the planned
distribution are then compared by viewing a dose difference map or by analysing DVHs.
DDAs are performed routinely on specified fractions for several types of treatment, and
also if a change in dose distribution is suspected, for example if there are set up issues or a
change in weight. Out of tolerance DDAs are discussed with the clinician, and may indicate
that a re-plan is necessary. As the MVCT images from TomoTherapy are used for DDAs,
the TomoTherapy MVCT images require a conversion of HU to relative electron density
for dose calculation. TomoTherapy use a method of converting to physical density, rather
than relative electron density, using Image Value to Density Tables (IVDT) to allow for
dose calculation. The energy relationship in both Compton Scattering and the Photoelec-
tric Effect, as well as their atomic number dependency, accounts for a decrease in contrast
in MVCT images compared to kVCT images. In the megavoltage energy range, Compton
interactions are dominant even in high Z materials. Due to this difference in physical inter-
action probabilities a linear HU to electron density curve is expected for MV images. To
use the MVCT images from TomoTherapy for dose computations, an MVCT HU number
42 A Method for Quantifying Random Range Uncertainty
to physical density calibration curve must be periodically obtained (Langen et al., 2010).
Fig. 3.3 Image taken from the Addenbrooke’s Radiotherapy Department’s work instruction
for carrying out MVCT HU QA using the Cheese phantom on TomoTherapy
Without a stoichiometric phantom available for calibrating the MVCT to proton stopping
powers relative to water another method was used. To calculate the WEPL from the MVCT
data sets a calibration curve of HU to relative proton stopping power to water is required.
The Addenbrooke’s data were calculated by calibrating the scanner to electron density using
the TomoTherapy ‘Cheese’ phantom with density plugs inserted as shown in Figure 3.3.
This process is carried out routinely as part of the monthly MVCT imaging QA. The density
plugs provided have values specified in g/m2. To convert the TomoTherapy HU to relative
proton stopping power the approximate values for the electron density relative to water for
the density plugs directly from Gammex were interpolated for values required. These values
were then converted to relative proton stopping power using used the ratio shown in equation
3.2 (Schaffner and Pedroni, 1998)) and the mean beta value shown in the Table 3.1, where
II and Iw are the ionisation potentials of tissue and water respectively, and Iw has been taken
as 75eV (ICRU 50, 1993). The approximate mean excitation values, I-values, were supplied
by Gammex (Gammex, 2009).
−dE
dx
= (
1
4 ·π · ε0 )
2 · (4π2
2e2ρe
mev2
) ·
[
ln
(
2mec2β 2
I · (1−β 2)
)
−β 2
]
(3.1)
3.3 Calculating WEPL Changes in the CTV from TomoTherapy Shift & MVCT Data 43
hence,
Srel =
ρe
ρew
·
[
ln
(
2mec2β 2
I·(1−β 2)
)
−β 2
]
[
ln
(
2mec2β 2
Iw·(1−β 2)
)
−β 2
] (3.2)
The absolute stopping power is clearly energy dependent from the Bethe-Bloch formula
(β = vc ). Using Bortfeld and Schlegel (1996)’s range-energy relationship shown in equation
3.3, four proton beam energies were chosen, within (or a bit above) the therapetic range, to
test if the relative stopping power can be assumed independent of proton energy allowing
for one curve to be used for all therapeutic energies (Hünemohr et al., 2014; Yang et al.,
2012).
R = α ·EP0 (3.3)
where R is the range, E0 is the initial particle energy; P is taken as 1.77 (ICRU 78, 2007)
and α is proportional to the square root of the effective atomic mass divided by the density
of the medium, which has been taken as 0.0022 for water (Jette and Chen, 2011).
In Table 3.2, β is calculated for four proton energies. For each insert the percentage dif-
ference between the relative proton stopping power calculated for each energy is shown in
the far right column. As cortical bone has an I-value furthest from that of water, this is
the tissue that will show the largest β dependency. The percentage difference in relative
stopping power for the cortical bone insert was less than 0.6% for the two most extreme
energies used, 90 MeV and 310MeV. With an error on the β value less than 0.6%, energy
independency can be assumed for this work. The relative proton stopping powers have been
calculated using the β associated with a 150MeV proton beam, as this beam energy can de-
liver a suitable treatment depth. It is already known from previous studies (such as Yadav
et al. (2010) and Duchateau et al. (2010)) that the IVDT can change over time, especially
with target change, which if not taken into account, could result in an error in the dose cal-
culation. Clinically we take a new IVDT each month to be compared to the one used in the
adaptive module. However for this case each patient was treated at a different time period
so an average IVDT has been used. The data plotted in Figure 3.4 is the average IVDT
taken from both TomoTherapy units at Addenbrooke’s taken over several months. Using all
44 A Method for Quantifying Random Range Uncertainty
Table 3.1 Calculating a mean β to use in the Bethe-Bloch equation to establish an energy
independent proton stopping power relative to water
Energy (MeV) 90 150 260 310
Proton Range in water (mm) 63.3 156.35 413.92 565.11
β 0.4092 0.5067 0.6221 0.6596
Table 3.2 calculation of proton stopping power relative to water for each tissue using the
energy specific β in the Bethe-Bloch equation. The I-values have been interpolated from
the Gammex data tables
Tissue Sub I value (eV) 90 MeV 150 MeV 260 MeV 310 MeV %diff
Lung 73.8 0.481 0.481 0.481 0.481 -0.03
Water 75 1.00 1.00 1.00 1.00 0.00
inner bone 80.1 1.081 1.081 1.082 1.082 0.11
Bone mineral 80.2 1.091 1.091 1.091 1.092 0.11
30% CaCO3 80.8 1.259 1.258 1.259 1.259 0.12
50% CaCO3 93.1 1.432 1.431 1.433 1.434 0.36
Cortical bone 104.5 1.622 1.622 1.625 1.626 0.56
of the above, the following Table 3.3 and the plots in figures 3.5 and 3.6 were created. The
equation of the line in figure 3.6 is the equation required in the equivalent depth algorithm
to calculate the WEPL through the MVCT images for a proton beam. (If we wanted to
know the WEPL change for conventional radiotherapy the line equation from 3.5 would be
required.)
3.3 Calculating WEPL Changes in the CTV from TomoTherapy Shift & MVCT Data 45
Fig. 3.4 Average MVCT IVDT from several months on both TomoTherapy units at Adden-
brooke’s
Fig. 3.5 Calibration curve of MVCT HU to electron density relative to water
46 A Method for Quantifying Random Range Uncertainty
Table 3.3 Table of density inserts used to calibrate the MVCT and their corresponding phys-
ical densities, estimated electron densities and proton stopping powers relative to water as
calculated from the Bethe-Bloch equation using the β value to a 150MeV proton beam
Tissue Sub Density e- Density SPw,150MeV HU
Air 0 0 0 -970
Lung 0.49 0.48 0.481 -479
Water 1 1.00 1.00 9.5
inner bone 1.139 1.09 1.081 118
Bone mineral 1.152 1.10 1.091 127
30% CaCO3 1.334 1.27 1.258 292
50% CaCO3 1.562 1.47 1.431 474
Cortical bone 1.824 1.69 1.622 673
3.3 Calculating WEPL Changes in the CTV from TomoTherapy Shift & MVCT Data 47
Fig. 3.6 Calibration curve of MVCT HU to proton stopping power relative to water
A shift in the HU representing water will have a larger impact than a shift in the HU repre-
senting bone since a typical patient image contains more water equivalent density material
than bone-like materials (Langen et al., 2010). For this reason real water has been used in
the HU calibration. There is a large variability in the Addenbrooke’s MVCT calibrations
between the two scanners used to treat and image the patients in this study. There also exist
errors related in interpolated I-values and electron densities of the inserts used in the Cheese
phantom. The values supplied by Gammex required an interpolation, and as the relationship
was hard to predict, a linear model was used. The I-value themselves constitute the largest
error in stopping power in the Bethe-Bloch equation, with up to a 10% mean excitation en-
ergy variation, even for water (Yang et al., 2012), the clinical impact of which is discussed
in Besemer et al. (2013). Taking into account many of the uncertainties in calibrating HU
to relative proton stopping power, it is estimated that even using a stoichiometric calibration
errors of ±3.5% still exist (Yang et al., 2012). There is currently research investigating the
use of both kV-kV and kV-MV Dual Energy CT (DECT), where kV-MV DECT has the
advantage of being less sensitive to uncertainties in the patient CT imaging (Yang et al.,
2012). The method described to calibrate the Addenbrooke’s MVCTs to proton stopping
power is not ideal for accurate measurements of WEPL and highlights the importance of
using a stoichiometric calibration and purpose built phantoms. For the patient data used,
48 A Method for Quantifying Random Range Uncertainty
patients have been treated across both scanners and in different months throughout 2012 -
2014. For an exact study, appropriate curves are required for each scanner and each time
period. The HU for each density insert has been acquired as part of the department’s routine
TomoTherapy QA. Due to beam hardening effects it would be better practice to image each
insert at a time at the centre of the phantom. In this work it has been investigated whether
this method of WEPL calculation could be used and whether it can yield a better estimate
of random range error due to inter- and intra- fraction motion. It is envisioned that such
a method could be applied in a proton centre using MVCT, or potentially CBCT. Due to
the nature of this calibration and the assumptions made, the data presented in Chapter 4 is
relative and not absolute and would be specific to Addenbrooke’s Hospital.
3.3.1 Coding
I carried out all coding (described below) using MatLab which is a high level computing
language that has pre-built functions which allows for fast programming (Matlab UK, 2011).
It also has capabilities for reading and writing DICOM data and reading information from
DICOM headers, which makes it a good language to use for this project. DICOM stands
for Digital Imaging and Communications in Medicine and is a network protocol for medical
data exchange. DICOM data includes a variety of image sources such as CT, MRI and PET.
DICOM data can also be in the form of radiotherapy objects such as treatment plans, dose
cubes and delineated structures. DICOM defines how all this data is stored and shared.
The list of the key functions, the description of their roles and the assumptions used as
well as the code can be found documented in Appendix B. Here a short description of the
process used to generate the probability distributions of inter-fraction range changes from
patient MVCT data.
Daily Shifts
Beyond having an effective depth algorithm and a calibration to be used within it, the daily
rigid couch shifts used on set for each fraction are required. For each treatment fraction on
TomoTherapy the radiographer will set the patient up for treatment and perform an MVCT
scan. The radiographer will perform image registration of the MVCT image to the planning
kVCT making any of their own adjustments online and determine the final shift required in
each axis (x, y, z) and for the roll. These shifts are translated to changes in the couch position
3.3 Calculating WEPL Changes in the CTV from TomoTherapy Shift & MVCT Data 49
and are applied. The patient will then be treated in this position, further details of the
treatment approach are given in Burnet et al. (2010). By applying these shifts to the MVCT
image the patient position at the time of treatment can be simulated. The correction for roll
has not been included in the calculation. As part of the VoxTox (voxtox, 2014) study shift
data was extracted directly from TomoTherapy for all patients. This data was embedded
into the DICOM header of each MVCT scan eliminating the human error involved in using
numbers quoted on set.
To apply a shift to the 3D MVCT for a given fraction two coordinate systems of the same
size of the MVCT matrix are created. The first is the coordinate system as retrieved from
the DICOM header; the second is this coordinate system with the shifts added to it in each
direction. For example if say pixel (1,1,1) existed at the coordinates (0,0,-12) and rigid
couch shifts (0.5,1.5,-1) were applied, the new coordinate system for pixel (1,1,1) would
be (0.5,1.5,-13). This is done to all pixels, in pixel dimension rather than in mm, and a 3D
cubic interpolation is carried out to determine the new values of pixel data. In the x and
y directions missing patient information is masked to correct for truncation artefacts. The
pixels representing air are then all assigned to the same value.
Truncation Artefacts
If a patient is large compared to the field of view (FOV) or positioned with part of their
body outside the FOV of the CT then truncation artefacts in the transmission data can occur
(McGowan et al., 2012). In radiotherapy truncation of either the region of interest, or of the
patient outline is a problem for dose calculation. On TomoTherapy the scanning circle has
a diameter of 38.6cm and truncation due to the scanning circle size is common, especially
of the shoulders in head and neck plans and of the body in prostate plans. An example of
scanning circle truncation of shown in Figure 3.7. Primarily the MVCT is used for daily
registration, but in the case of dose re-computation for adaptive planning truncation needs to
be corrected for. Solutions to correct for truncation artefacts include solutions such as using
the kVCT data to fill in the missing data may be used. Alternatively, the body outline can
be used to mask this missing data with a uniform pixel value. In both cases the daily shift
will need to be applied to either the skin outline or the kVCT, so it will align correctly with
the shifted MVCT. In this work a mixture of both methods was used. The kVCT was read
in with the MVCT data and the shifts were applied to both. Knowing the scanning circle
size, a mask was created for both CTs. The kVCT contributed the data that existed outside
of the scanning circle and the MVCT provided the data within the scanning circle.
50 A Method for Quantifying Random Range Uncertainty
Fig. 3.7 Image of truncated patient on TomoTherapy showing the the scanning circle (white
circle).
WEPL
To calculate the change in WEPL between fractions, each fraction is compared to a refer-
ence image. Due to the daily images and planning images having been obtained via different
imaging modalities the reference image used was the MVCT of the first fraction in the treat-
ment. The structure sets created at planning were used to determine the region of interest
(ROI). The structure set of the ROI was read into the program and the coordinates of the
structure were aligned to the corresponding shifted MVCT image slices. The resolution of
the MVCT and the structure set differ in slice thickness so to determine which structure
set slices were required the in-built MatLab function ismember was used to pinpoint the
elements of the two matrix sets where the z coordinates were in correspondence.
The WEPL was calculated for a given beam direction, and beam parameters, using the
shifted MVCT data and the effective depth algorithm previously described. To ensure the
quality of the effective depth algorithm several verification tests have been carried out; these
included testing the WEPL is accurate by setting the source to the surface of the CT image
and the beam angle perpendicular to the pixels. This way the value of the pixel can be
checked to determine if it was as expected and that the next pixel had a value equal to its
own plus that of the previous pixel. An easy way to check the summation process is to look
at water and see that the pixel value increases by a distance of one pixel each time.
3.3 Calculating WEPL Changes in the CTV from TomoTherapy Shift & MVCT Data 51
Following the calculation of the WEPL, masks were used to retrieve the WEPL data from
only pixels that exist inside the structure set for each fraction. Each subsequent fraction is
then compared to the reference. All dimensions are converted into distances and volumes as
required. Statistical parameters, such as the mean and standard deviation (in mm) of each
data set were calculated; these represent the systematic and range errors associated with the
treatment, both as a whole and at each fraction for that given beam angle.
Assumptions
Several assumptions have been made; these include
■ The target has been delineated correctly.
■ The target is in the correct location.
■ The target has not rotated.
■ Image matching and the shifts applied at treatment were correct.
■ The target has not been deformed.
■ Air outside the patient can be represented as a vacuum.
■ The Source to Isocentre Distance (SID) has been set to seven metres in this study. This
corresponds with the distance in the proton centre at CNAO, Pavia as proveided to us
by the physicists there, this is used simply as an example. At PSI, at dose calculation
the SID is assumed infinite.
3.3.2 Analysis
For each set of patient data the WEPL for each pixel in the volume of interest for every
fraction. The WEPL value stored in each pixel of the reference is then subtracted from that
in each fraction. This gives the change in WEPL between the first and consecutive fractions.
The data is stored in this raw form and also as binned data. A positive difference in range
reflects undershoot scenario when compared to first fraction and for negative differences an
overshoot scenario. The results can be displayed as histograms for each individual fraction.
52 A Method for Quantifying Random Range Uncertainty
The aim of this work is not to analyse individual fractions or even patients (though this
information may be useful in analysing plan quality throughout the treatment delivery), but
to analyse a sample of patient’s treatments. The mean and standard deviation, σ of WEPL
changes in mm are calculated for each beam angle across all patients and all fractions for
both treatment site. The means of each individual fraction are also plotted as a histogram
and the standard deviation calculated, Σ. It is these values that are being investigated and
are required for a greater understanding of random range uncertainty in proton therapy.
3.4 Summary
There is a need to be able to quantify the random range uncertainty in proton therapy. In
this chapter a method has been suggested by calculating the change in WEPL throughout
treatment using daily volumetric imaging with rigid couch shifts applied. The process of
creating a calibration curve of HU to relative proton stopping power for the TomoTherapy
MVCT has been laid out. Though it is not a stoichiometric calibration, it is suitable for
measuring changes in WEPL between fractions using daily MVCT images. Finally, an
outline of the methods used to carry out this calculation have been described.
Chapter 4
Quantifying Site-specific Range
Uncertainty
4.1 Introduction
The method described in chapter 3 is tested using clinical patient data to calculate the change
in water equivalent path-length due to inter-fraction motion. Using a sample of patients and
two treatment sites, the limitations of this method are explored. Finally, calculations for the
overall mean, systematic and random range error for the patient sample are calculated.
4.2 Methods and Materials
Nine prostate and seven head and neck patients have been selected; these patients have been
anonymised and tokenised for their relevant trials to uphold patient confidentiality. The
standard protocol for prostate is to image daily using MVCT with a slice thickness of 6mm,
and for head and neck it is 3mm. For the prostate cases, six of the patients were on the
PIVOTAL trial. All patients, prostate and head and neck, are in the VoxTox study (voxtox,
2014). PIVOTAL is a phase II trial randomising between prostate and prostate + pelvic
IMRT in patients with locally advanced prostate cancer (Harris et al., 2011). VoxTox is a
Cancer Research UK funded observational study at the University of Cambridge collecting
toxicity data for patients undergoing image guided IMRT to the head and neck, prostate
and CNS. For both treatment sites daily imaging has been carried out using MVCT on
TomoTherapy and rigid shifts applied after soft tissue matching has been carried out by
the treatment radiographers. This means that the range errors calculated in this chapter are
54 Quantifying Site-specific Range Uncertainty
residual errors remaining after daily IGRT.
These specific patients have been selected because they were all treated on TomoTherapy at
Addenbrooke’s and have all given verbal and written consent for their images to be used for
research. The VoxTox study was approved ethically by the NRES East of England Essex
Ethics Committee in 2013. In the VoxTox study the head and neck cases required scans
to have been long enough to ensure both parotid glands are imaged entirely on the MVCT.
An increase in scan length by approximately 2cm was needed for around 80% of the cases.
For the head and neck patients the estimated maximum dose using the standard protocol is
11.2 mSv. The estimated maximum dose using the study protocol is 12.2 mSv (Bates et al.,
2013). This dose is 0.2% of the maximum dose from their radiotherapy. The combined
additional total dose corresponds to a risk of radiation-induced cancer during the rest of the
lifetime of the standard population of approximately 1 in 20,750 (ICRP, 2007). The prostate
patients were all treated for both prostate and nodes and no extra scan length was applied.
In this work the change in WEPL within a Region of Interest (ROI) is carried out. The
ROIs include the prostate without margins and the high dose CTV for the prostate and head
and neck patients respectively. The advantage of the increased scan length in the head and
neck data set, and the need for nodal irradiation in the prostate data set is that, despite the
investigation only being focused on either the high dose CTV or prostate itself, there is a
greater chance that the daily imaging will cover the entire ROI irrespective of the shifts that
will need to be applied to align the patient. This will hopefully yield minimal truncation
of the ROI in the scanning direction (superior-inferior direction). This is necessary for this
study, as the results cannot be obtained if there is truncation of the ROI in the MVCT image
at each fraction.
In the case of the prostate data two beam angles have been used to investigate inter-fractional
changes in WEPL; a coplanar anterior and a lateral beam, 0◦ and 90◦. In the case of the head
and neck coplanar beam angles of 0◦, 25◦, 70◦, 110◦, and 180◦ have been used. The use of
several beam angles was to gain an insight into how direction and the amount of traversing
material affects the magnitude of range error. In the prostate case these two beam angles
have been chosen, as they are clinically relevant beam angles. Though prostate is not a site
that will be treated in the new proton centres in England, these data sets have been used
in this work. For one of these patients there was opportunity to investigate intra-fraction
changes. Permission was granted by the patient and ethically approved by NRES to image
twice on some fractions, both before and after treatment as part of a study to determine
4.3 Results 55
intra-fraction motion (Thomas et al., 2013). However, using prostate cases within this work
has been due partly to the data that was available. Irrespective of cancer site and whether
this site will be eligible for treatment on the NHS, this data can be used to investigate how
this method of quantifying range uncertainty could be used for any location. It also offers as
an opportunity to investigate how range uncertainty can differ between site locations. The
benefit of treating prostate cancer with protons over conventional radiotherapy has been a
controversial topic in the proton therapy world and papers have been written on the topic
(Slater et al., 2004; Vees et al., 2014; Zietman et al., 2010).
Scanning circle truncation artefacts are present for both prostate and head and neck treat-
ment sites and are an issue for several beam angles. In the case of prostate often a few
centimetres (or less) of the body may suffer from truncation on several slices. This part of
the body is fairly homogeneous and it has been assumed to have a uniform density. For the
prostate patients a HU of -100 was used in the mask created from the kVCT. This was an
estimate of the average HU value in this area based on several of the patient MVCT images.
Truncation of the shoulders in the head and neck is more problematic due to the density
changes in this region. A blanket HU number substitution could not be used and there was
not a calibration created for the kVCT scanner. Due to the large number of slices the head
and neck high dose CTV extends over, it was deemed acceptable to remove the last few
MVCT slices so that the WEPL changes could be analysed without being compromised by
truncation artefact.
In this chapter it is attempted to characterise the magnitude and probability of inter-fraction
range error for prostate and head and neck treatments for the beam directions investigated.
Though these sites and beam angles are not necessarily those that would be used clinically,
the aim of this work is to demonstrate how using daily MVCT images could be used to
quantify this error.
4.3 Results
4.3.1 Prostate
All prostate data sets were run for two beam angles, 0◦ and 90◦. Despite using patient data
sets that required prostate and nodal irradiation some fractions were found to be unusable
due to missing prostate data in the scanning direction once the shifts were applied. This
is a form of truncation artefact of the prostate itself. This is a real issue for IGRT studies,
56 Quantifying Site-specific Range Uncertainty
including for VoxTox (Bates et al., 2013). In Table 4.1 the full information about each
prostate data set is supplied. This includes the patient token, the name of the structure
investigated, the total number of fractions treated, the total number of fractions used in
the analysis, what fraction the comparison was made to and finally, the percentage of the
fractions that were used in the analysis. Though the names of the structures differ, all the
volumes used were the prostate without margins applied. The different name was because
for patients in the PIVOTAL trial it was required that the prostate be named ‘Tp’. Despite
attempting to use patient data that would not suffer from ROI truncation in the scanning
direction many of the prostate fractions could not be analysed due to part of the ROI being
missed off the scanning image.
The WEPL for each fraction is compared to the WEPL of a reference fraction and the
difference is stored. This data has been presented as histograms, each of which are plotted
in figures 4.1 to 4.2. Each of the plots in these figures are heat maps of the frequency of
the WEPL difference at each fraction. The darker the colour the greater the frequency of a
difference in WEPL of that magnitude. Along the x-axis is time in terms of fractions, and
along the y-axis is the change in range in millimetres. The zero line is shown so it can be
seen how the mean range error at each fraction changes with time. The plots have been
truncated down for visual purposes so outliers are not always represented here. However
all data, including outliers are included in the statistical analysis. The left hand set of the
histograms are for a beam angle of 0◦ and the right hand set of the histograms are for a beam
angle of 90◦.
To investigate how the random range error changes with both the beam angle and treatment
fraction the standard deviation of range errors for each fraction has been plotted in Figures
4.4 and 4.3 for each beam angle. The colour represents a different beam angle, with blue
being 0◦ and red being 90◦, and each plot is a single patient. Again the x-axis is time in terms
of treatment fraction, but the y-axis is the standard deviation in millimetre for that fraction
and beam angle. It can be seen there is a beam angle dependency for most of the patients
with a 0◦ beam angle resulting in random range errors with a smaller standard deviation.
4.3 Results 57
Table 4.1 Prostate data showing the target and number of fractions used. The final column
shows the percentage of fractions that did not suffer target truncation at daily imaging
Token Target Total #s Usable #s Compared to # %data usable
4 Tp 37 26 3 70.3
7 Prostate 37 13 2 35.1
21 Tp 36 12 2 33.3
34 Tp 37 26 1 70.3
36 Prostate 37 12 1 32.4
58 Tp 37 37 1 100.0
72 Prostate 37 19 2 51.4
78 Tp 37 19 4 51.4
95 Tp 37 28 2 75.7
58 Quantifying Site-specific Range Uncertainty
Fig. 4.1 Heat plots of %volume range error histograms for prostate patients 4, 7, 21, 34, and
36 (in order top to bottom) with fraction using 0 and 90 degree beam angles.
4.3 Results 59
Fig. 4.2 Heat plots of %volume range error histograms for prostate patients 58, 72, 78 and
95 (in order top to bottom) with fraction using 0 and 90 degree beam angles.
60 Quantifying Site-specific Range Uncertainty
Fig. 4.3 Plots of the standard deviation for patients 4, 7, 21, 34, and 36 (in order top to
bottom) of range error with fraction and beam angle.
4.3 Results 61
Fig. 4.4 Plots of the standard deviation for patients 58, 72, 78 and 95 (in order top to bottom)
of range error with fraction and beam angle.
62 Quantifying Site-specific Range Uncertainty
Intra-fraction
For patient 58, eight fractions were imaged before and after treatment. This has allowed
for the WEPL change within a treatment fraction to be measured for this patient. The mean
time between the start of the first image and the end of the second was 12 min (Thomas
et al., 2013). The intra-fraction results included in table 4.2 showing the patient mean and
standard deviation (σ ) in mm.
4.3.2 Head and Neck
Seven head and neck patients were analysed for five angles, 0◦, 25◦, 70◦, 110◦ and 180◦,
as shown in Figure 4.5. The high dose CTV was used as the ROI for these patients. The
high dose CTV includes the nodal region extending down the patient’s neck. Table 4.3
includes all the information about these patients including the percentage of the treatment
successfully analysed. The fractions that could not be analysed were due to part of the
region of interest being missed off the scanning image. In the head and neck cases the first
or second fraction has been used as the reference to make comparisons. Figures 4.6 to 4.8
show the heat plots of the change in WEPL histograms for each fraction and 4.9 shows the
standard deviation of WEPL change in mm at each fraction for each patient and beam angle.
4.3 Results 63
Table 4.2 Table of Intra-fraction changes in WEPL for a single prostate patient, 58, at beam
angles 0◦ and 90◦.
Angle 0◦ 90◦
Fraction Mean (mm) σ(mm) Mean (mm) σ(mm)
1 1.6 1.46 3.68 0.77
5 0.15 1.66 2.36 1.2
12 -1.42 0.57 -1.22 0.93
13 -0.13 1.72 3.84 1.48
21 -0.24 0.79 0.39 0.69
26 1.77 1.66 -5.95 0.85
27 5.6 1.5 7.3 1.41
28 1.73 2.05 1.81 1.54
Average 1.13 0.59 1.53 0.37
Fig. 4.5 Screen shot of a single slice of Patient 78 showing the high dose CTV and approxi-
mate angles used to calulate WEPL
64 Quantifying Site-specific Range Uncertainty
Table 4.3 Table of Head and Neck data showing the targets and number of fractions used.
The final column shows the percentage of fractions that did not suffer target truncation at
imaging
Token Target Total #s Usable #s Compared to # %data usable
35 CTV65 30 30 1 100.0
51 CTV65 31 21 2 96.8
52 CTV65 34 30 1 52.9
78 CTV60 30 30 1 100.0
91 CTV60 30 28 1 93.3
92 CTV65 30 30 1 100.0
102 CTV68 40 28 2 70
4.3 Results 65
Fig. 4.6 Heat plots of %volume range error histograms for all head and neck patients with
treatment fraction at 0 and 25 degree beam angles.
66 Quantifying Site-specific Range Uncertainty
Fig. 4.7 Heat plots of %volume range error histograms for all head and neck patients with
treatment fraction at 70 and 110 degree beam angles.
4.3 Results 67
Fig. 4.8 Heat plots of %volume range error histograms for all head and neck patients with
treatment fraction at 180 degree beam angle.
68 Quantifying Site-specific Range Uncertainty
Fig. 4.9 Plots of the standard deviation of range error for each patient with fraction and
beam angle.
4.3 Results 69
4.3.3 Population Data
The aim of this project was to build a better understanding of the magnitude of range uncer-
tainty from daily patient changes in proton therapy. It was also to determine a method for
this information to be used at robust optimisation and plan analysis. In treatment planning
a single composite random range uncertainty is convenient for either choosing margins, or
for using in a robust optimisation. To this end, estimates of the overall mean error, random
inter-fraction and systematic inter-fraction error for all patients have been calculated. These
have been calculated using methods similar to those set out by Greener (2003), though here
weighting has been by number of images analysed rather than number of patients for both
the overall mean and the systematic uncertainty (Σ). The systematic inter-fraction error was
calculated by plotting a histogram of the mean range change for each usable fraction across
all patients and beam angles shown in Figures 4.10 to 4.13 and taking the standard devia-
tion. Figure 4.13 shows the overall mean histogram for both the prostate data (top) and head
and neck (bottom). For all histograms, a Gaussian distribution has been fitted.
For prostate and head and neck patients, the binned data from every fraction for each beam
angle was also plotted, but as a line histogram shown in Figures 4.14 and 4.15 for each treat-
ment site respectively. The black histogram represents the average histogram of all angles,
whereas the colours represent individual beam angles. The standard deviation (σ ) were cal-
culated from all the non-binned fraction data, this means that the random uncertainty was
calculated with weighting on number of data points and not number of images. All means,
Σ’s and σ ’s for all patients for each beam angle, are shown in Tables 4.4 and 4.5 and is
representative of that treatment site.
70 Quantifying Site-specific Range Uncertainty
Fig. 4.10 The fractional means of WEPL difference for each prostate patient was plotted
as a histogram for both angles (top is 0◦ and bottom 90◦). The total mean can be taken
from these plots for each angle as well as the standard deviation which is the systematic
component of the range error for these sites.
4.3 Results 71
Fig. 4.11 The fractional means of WEPL difference for each head and neck patient was
plotted as a histogram for three angles (top is 0◦, middle is 25◦ and bottom 70◦). The total
mean can be taken from these plots for each angle as well as the standard deviation which
is the systematic component of the range error for these sites.
72 Quantifying Site-specific Range Uncertainty
Fig. 4.12 The fractional means of WEPL difference for each head and neck patient was
plotted as a histogram for three angles (top is 110◦ and bottom 180◦). The total mean can
be taken from these plots for each angle as well as the standard deviation which is the
systematic component of the range error for these sites.
4.3 Results 73
Fig. 4.13 The fractional means of WEPL difference for all prostate patients (top) and all the
head and neck (bottom) were plotted as a histogram for all angles. The total mean can be
taken from these plots as well as the standard deviation.
74 Quantifying Site-specific Range Uncertainty
Table 4.4 Table of prostate range error in mm for all patients at each beam angle.
Angle Mean (mm) Σ(mm) σ(mm)
0 -0.98 4.53 5.37
90 2.94 5.94 6.45
All 0.98 5.64 5.91
Fig. 4.14 The binned prostate data for each beam angle from each fraction for all patients
was summed and plotted as a histogram normalised to one. Each colour represents a beam
angle, and the black histogram is of all binned data.
4.3 Results 75
Table 4.5 Table of head and neck range error in mm for all patients at each beam angle
Angle Mean (mm) Σ(mm) σ(mm)
0 -1.22 3.21 4.71
25 -1.31 2.97 4.47
70 -1.91 3.20 3.93
110 -2.71 3.83 4.81
180 -1.95 3.48 5.69
All -1.82 3.39 4.72
Fig. 4.15 The binned Head and neck data for each beam angle from each fraction for all
patients was summed and plotted as a histogram normalised to one. Each colour represents
a beam angle, and the black histogram is of all binned data.
76 Quantifying Site-specific Range Uncertainty
4.4 Discussion
The method previously described in Chapter 3 has been tested using patient data to deter-
mine the range changes for a selection of prostate and head and neck patients. The range
error from inter-fraction motion has been calculated by comparing the WEPL in each pixel
within a region of interest at each fraction to the equivalent one in a fraction chosen, usu-
ally the first fraction, as a reference. Ideally, each treatment fraction would be compared to
the planning CT, as this is the reference scan the treatment will be planned to. Due to the
difference in the two modalities, kV and MV imaging, this not been carried out. However,
another solution may be to acquire an MVCT at the same time as the planning kVCT.
In the prostate data the gaps seen in the plots is where imaging data was missing for these
fractions. From the heat maps it was seen that there was a large spread in the means and
no particular trend in their values. The standard deviation was seen to be larger for 90◦, yet
there was not a distinguishable difference seen in that mean values between the two angles.
For 90◦ the standard deviation was greater than that seen for 0◦ in all except two patients,
one of which the standard deviations of both angles were comparable. The spread in means
represent the systematic component of inter-fraction error, whereas the standard deviation
represents the random component. The larger range errors seen in the 90◦ was most likely
caused by the fact that the prostate moves independently of bony anatomy and the images are
aligned to soft tissue matching. Comparing the standard deviations for individual fraction in
Figures 4.3 and 4.4 could be misleading, as the large standard deviation seen in the prostate
population data does not appear to coincide with the individual fraction data. This is due to
the large variation in the mean as seen in the heat maps. A large amount of the data could
not be analysed due to missing data within the ROI, which leads to a bias in this data set.
One patient, patient 58, had their images used to investigate intra-fraction range changes
for eight fractions. The results from are less than 1mm and as expected, comparable to the
results seen in Thomas et al. (2013).
For the head and neck data a lot more images could be analysed leading to more reliable
data. The mean range error at each fraction is seen to be closer to zero than that seen in the
prostate and there is less variability in the means with fraction. There is a noticeable shift in
the mean from zero for several patients with time which is also seen in the population results.
This shows that throughout treatment the WEPL is decreasing. This may be due to slight
weight loss and represents a risk of proton overshoot. Three patient examples are shown in
Figure 4.17, each shows an overlay of the same MVCT slice of the first and last fraction.
4.4 Discussion 77
Yellow represents pixels that match well, red represents the first fraction image and green
the final fraction image. Slight weight loss throughout treatment appears uniform across this
sample and across beam angles. There is a larger spread in ranges seen for the larger beam
angles (70◦ - 180◦). This is due to two reasons. Firstly, that the ‘beam’ is passing through a
greater volume of the patient before reaching the target, and travelling through more density
heterogeneities. Secondly, for some of the patients, the high dose CTV was on the opposite
side of the head and neck to the direction of the beam angles chosen. This has led to a
larger amount of tissue to be traversed, as well as more density heterogeneities in the proton
path. Further investigation using clinical treatment beam angles is required. Due to this
second issue a slightly smaller standard deviation for the head and neck population may be
expected. The standard deviation in the head and neck for angles 110◦ and 180◦ were larger
than those seen in the prostate data. This is most likely due to the density heterogeneities
seen in the head and neck leading to larger range uncertainties.
In several of the head and neck patients the CTV was partly in the mouth and when using
angles such as zero or 25◦ the proton path passed through the patient’s jaw and teeth. For
these two angles tooth fillings will affect the range. Figure 4.16 illustrates the problem of
tooth filling in both a kVCT image (left) and MVCT image (right). Though artefacts are
not present in the MVCT, a small setup or motion error, in any direction other than the
beam direction, will result in a large range error for a proton beam travelling through the
filling due to the high density of the filling material. Another issue is that the calibration
has not been made to include such high density material. In a clinical proton case beam
angles would be chosen to avoid passing through the teeth with fillings or the teeth may be
removed due to the risk these materials and the resulting artefacts can have on the proton
range.
78 Quantifying Site-specific Range Uncertainty
Fig. 4.16 Example of tooth filling artefact in kV CT image (left) and in MVCT image (right)
Fig. 4.17 Possible weight loss leading to negative Σ for Head and Neck. Images left to right
of patients 52, 91, 92
4.4 Discussion 79
The population results from the head and neck patient sample show range uncertainties
from inter-fraction changes are smaller than that see in the prostate data. The systematic
inter-fraction range uncertainty was largest for prostate for 90◦. The population data for
head and neck across all angles had the smallest random range uncertainties. The random
range uncertainty in the head and neck data is slightly smaller than that of the prostate. The
density heterogeneities in the head and neck prove to be as challenging as range changes due
the target moving relative to bony anatomy in the prostate. The systematic range error in
head and neck was also smaller than that of the prostate, and this was seen in the individual
heat maps. The overall means for all population data was non-zero. In the head and neck
this is believed to be potentially due to slight weight loss throughout treatment. Due to the
limited number of patients used in this study and the number of fractions that could not
be analysed for the prostate patients in particular it is expected that the values for Σ and σ
are an overestimate. As well carrying this work out with a larger patient sample, further
work to provide suitable correction strategies for determining the systematic error would
be to identify those errors unlikely to have occurred by chance and apply a correction as
suggested by Greener (2003).
These results have highlighted a limitation in this method of using of IGRT images from
conventional radiotherapy due to truncation of the region of interest due to changes in scan-
ning length. This problem has led to a lot of the fractions being unusable. Daily TomoTher-
apy IGRT allows for online correction of inter-fraction geometric uncertainties. In this work
these datasets have been used offline to quantify inter-fraction range uncertainties. In the
case of this study the data sets are used in a way they were not intended for, and it is this
reason that gives rise to their limitations for this work. For this type of data to be used
for quantifying WEPL changes set protocols on patient set-up, imaging criteria and what
constitutes as a volume of interest will be required. In a study carried out at Addenbrooke’s
by Bates et al. investigating IGRT imaging times using TomoTherapy, it was found that
for every additional 1cm of scan length resulted in a extra 8 seconds to acquire at a slice
thickness of 6mm and an extra 12 seconds to image match. Across a normal treatment day
this would add an extra 6.7 minutes per day in imaging and matching (Bates et al., 2013),
which it was deemed implementable at Addenbrooke’s for the use in IGRT studies. Another
solution would be to extrapolate the missing data after the shifts are applied to the MVCT
data and to keep this in consideration when analysing the data.
80 Quantifying Site-specific Range Uncertainty
4.5 Conclusion
This work aimed to demonstrate how the use of daily imaging data sets from previously
treated conventional radiotherapy patients could be used to quantify the random range un-
certainty in the proton radiological path-length (WEPL) within the CTV caused by inter-
fraction motion. Though limitations in the data sets exist, this method can be adjusted and
solutions have been discussed to allow it to be suitable to investigating inter- fraction range
uncertainty in proton therapy. The method has been used for two treatment sites and sev-
eral beam angles. Further work will be to carry out such calculation using clinical beam
arrangements and image guidance protocols. The same method can be applied to investigat-
ing range changes between CTV and OARs in the beam direction.
Based on these results for head and neck, assuming range uncertainty to be normally dis-
tributed, a Σ =3.4mm and a σ = 4.7mm was calculated with an overall mean of -1.82mm.
These values are represenative of the residual range uncertianty still present after daily cor-
rection has been applied. For prostate, a Σ =5.64mm and a σ = 5.91mm was calculated with
an overall mean of 0.98mm. It is deemed that there was too few prostate images analysed to
produce a non-biased determination of the population range uncertainty, however, a smaller
range uncertainty is present for an anterior beam compared to a lateral beam direction. Mea-
sured data of this kind will support more sophisticated models for future studies allowing for
limits on the actual magnitude of range uncertainty to derive more precise uncertainty mod-
els for specific tumour sites and planning protocols. These results show that inter-fraction
range changes are not insignificant and therefore, should not be excluded from either the
plan creation or plan analysis when considering robustness.
Chapter 5
Defining Robustness Protocols
The work in this chapter was carried out during a five month placement at the Proton Ther-
apy Centre at the Paul Scherrer Institute in Villigen, Switzerland. The contents of this chap-
ter, including tables and figures have been accepted for publication in Physics in Medicine
and Biology on the 8th December 2014. The permission of my co-authors has been given
to use the work here.
5.1 Introduction
Whether a plan has been created using margins or through a method of robust optimisation,
analysis of a plan’s robustness is required. Quantifying the uncertainty for given protocols
will aid in setting limits on the amount of robustness needed. As a consequence, the as-
pect of robustness should be reflected in dose reporting protocols as to better represent the
quality of the plan (Unkelbach et al., 2007). In this chapter plan robustness is retrospec-
tively analysed for a sample of clinical plans. Using this information it is demonstrated how
centre- and site-specific robustness constraints may be established for different volumes to
aid making planning decisions. 3D IMPT offers the ability to deliver highly conformal
radiotherapy to the patient, reducing integral normal tissue dose (Chang et al., 2006) and
providing greater OAR sparing through the use of multiple inhomogeneous fields. However,
without certainty in our ability to deliver these highly complex and conformal plans at each
fraction caution must be heeded when choosing beam arrangements and positioning steep
dose gradients near OARs (Albertini et al., 2010; Lomax, 2008b; McGowan et al., 2013).
As mentioned, any uncertainty might lead to severe target underdose or OAR overdose.
Several authors are proposing to include the plan robustness as an extra parameter directly
during the optimisation process (Chen et al., 2012; Fredriksson et al., 2011; Liu et al., 2012;
82 Defining Robustness Protocols
Pflugfelder, 2008; Unkelbach et al., 2007). Unfortunately when introducing the robustness
parameter into the optimisation algorithm, the nominal plan quality is often compromised
(Fredriksson et al., 2011; Unkelbach et al., 2007). Improving plan plan robustness is clearly
desirable and producing robust plans is probably the correct way to deal with the problem
of uncertainties in IMPT. However, it should not be at the cost of compromising nominal
plan quality to achieve robustness against very large and extremely unlikely uncertainties
(Unkelbach et al., 2007). A thorough understanding of the magnitude and probability of
uncertainites is required to ensure plan conformity is not being sacrificed unnecessarily. It
is extremely important that the trade-off between plan conformality and plan robustness is
thoroughly explored. Clinical experience so far has been achieved without considering the
robustness of the plan. Although the importance of including robustness during the plan
analysis is recognised, it is not yet clear how much the nominal dose distribution can be
compromised without compromising the patient’s treatment. The actual status is therefore
that the robustness analysis is not yet fully included into the clinical plan evaluation phase.
Currently no commercial planning system has a robustness evaluation module integrated
and to the best of my knowledge, the research centres that have developed their own ro-
bustness analysis tools are not yet using them clinically to make decisions. It is therefore
of paramount importance to find a way to introduce the robustness analysis into the clin-
ical plan evaluation process, without compromising the nominal dosimetric plan quality.
Specifically, site-specific robustness thresholds are needed to define a threshold between ro-
bustness and plan quality. It is important to establish the adequate level of robustness for
each volume, both the target and OAR. In doing so, this can be used on one side as an input
parameter during the robust optimisation process and on the other side as control parameter
during the plan evaluation phase. The concept of a site-specific robustness database can
indeed be used for both purposes.
5.1.1 Methods and Materials
In this chapter a method to define site specific robustness protocols is proposed based on
retrospectively analysing the robustness to both range and set-up uncertainties for a sample
of Chordoma and Chondrosarcoma plans. All plans were clinically acceptable and had been
delivered in the last five years at the Paul Scherrer Institute (PSI). All plans were calculated
using the in-house treatment planning system at PSI. The treatment planning system uses a
pencil beam model to calculate absolute dose delivered to the patient (Pedroni et al., 2005;
Scheib and Pedroni, 1992) and uses a quasi-Newton optimisation technique (Lomax, 1999).
5.1 Introduction 83
The model parameters have been adjusted to fit the measured data, which takes into account
range straggling, proton energy loss and the width of the initial energy spectrum. To take
into account density heterogeneities in the patient, CT data is converted into water equivalent
depth using a ray casting model (Schaffner et al., 1999).
To determine plan robustness the set-up and range uncertainties were modelled using the
method proposed by Albertini et al. (2011b) and later validated by Casiraghi et al. (2013).
This method will be briefly described here.
Fig. 5.1 The ebDD represents in each voxel a dose-error-bar that brackets the calculated
deviations from the nominal dose distribution. Here are the screenshots of the nominal plan
(using typical PSI beam arrangements referenced as A(4f) in this chapter) and the ebDD
calculated for ±3% HU from one patient.
Set-up Uncertainty: The effect of random set-up uncertainty on plan robustness was sim-
ulated by recalculating the nominal dose distribution on 14 isocentrically shifted CT data
sets. The shifts used were applied along the major anatomical axes and their indices in both
positive and negative directions (equivalent to the six faces and eight vertices of a cube) re-
sulting in multiple spatially shifted dose distributions associated with each shifted CT data
set. Bolsi et al. (2008) showed that the use of a remote patient positioning device reduced
systematic set-up uncertainties to negligible values (below 0.6mm), whereas random set-up
uncertainties of over 2mm can be observed, depending on tumour location and immobil-
isation device used. Using this data a three dimensional shift vector was calculated from
the standard deviation of each shift in left-right (LR), anterior-posterior (AP) and cranium-
caudal (CC) directions which represents of the radius of the spherical error-space used in
this analysis. Following the same method as Albertini et al. (2011b) a shift of magnitude
84 Defining Robustness Protocols
equivalent to a confidence interval of 85% was applied. In particular shifts of ±3.2mm,
±4.7mm were used for bite block and for head mask fixation devices respectively. The
nominal dose distribution and each recalculated dose distribution are then combined into
one error-bar dose distribution (ebDD) by storing in each voxel the difference between the
maximum (Max(Di)) and minimum (Min(Di)) value for that corresponding voxel (i) from
all the shifted plans and the nominal plan (nom). This amplitude of dose uncertainties, ∆Di,
in each voxel is calculated using formula 5.1, where Di is the dose in voxel after each shift
has been applied and nom refers to the nominal distribution. Each voxel stores a value which
represents a dose-error-bar that brackets all the possible deviations from the nominal dose
distribution that can be detected when patient shifts within an 85% confidence interval occur
(Casiraghi et al., 2013).
∆Dirandom = Max(Di+LR,D
i
−LR, ...,D
i
−CC,D
i
nom)−Min(Di+LR,Di−LR, ...,Di−CC,Dinom)
(5.1)
Range Uncertainty: Radiotherapy plans are generally produced using a single CT data set of
the patient and so any uncertainties in either the Hounsfield Units (HU), or their conversion
to relative proton stopping power, will lead to a systematic uncertainty that will propagate
throughout the treatment delivery (Lomax, 2008a). The effect of systematic range uncer-
tainty on plan robustness was simulated by recalculating the nominal dose distribution on
CT data with HUs altered by ±3%, thus to simulate an overdose or an underdose scenario
[Lomax 2008]. In this case the ebDD is calculated using formula 5.2, where DiHU refers to
the dose at voxel i after a ±HU shift has been appled.
∆Disystematic = Max(Di+HU ,D
i
−HU ,D
i
nom)−Min(Di+HU ,Di−HU ,Dinom) (5.2)
The robustness of each volume has been analysed by extracting the error-bar volume his-
togram (ebVH) . The error-bar volume histogram represents a useful and simple tool to
evaluate the quality of the treatment: the closer the histogram is to the ‘0-error’ line, the
more robust the dose is to that specific uncertainty. In particular its similarity in presen-
tation to the familiar dose volume histogram used conventionally in radiotherapy allows
for easy integration into the planning process. Indeed several metrics can be defined and
analysed as for example the maximum uncertainty, the mean uncertainty and the volume
receiving an uncertainty of 3%, 5% and 10% (Veb3%, Veb5%, Veb10%).
There is a greater clinical relevance in displaying the potential over-dosage uncertainty
5.1 Introduction 85
(highest minus nominal value) for the organ at risk and the under-dosage uncertainty (nom-
inal minus lowest value) for the target volume, instead of the more general error-bar dose
distribution discussed above, so this has been implemented when creating the database.
5.1.2 Retrospective Robustness Analysis
At PSI Chordomas and chondrosarcomas are treated to 74Gy(RBE) and 68Gy(RBE) re-
spectively. Each treatment consists of three phases delivering a given proportion of the
total prescribed dose to the target at each fraction. The first phase is generally comprised
of a 3- field uniform dose distribution (SFUD) plan used to deliver up to 34-40Gy(RBE)
of homogenous dose to the target volume. The only constraint within the SFUD plan op-
timisation is to deliver a homogeneous dose to the target volume. As a consequence, the
resulting target coverage will be robust to uncertainties, but at the cost of some compromise
of organs-at-risk (OARs) sparing (Albertini et al., 2011b). The following two phases are
comprised of 3D-Intensity Modulated Proton Therapy (IMPT) plans used to provide dose to
the target whilst prioritising OAR sparing. The second phase then will bring the dose to the
primary target up to 54Gy(RBE), sparing some critical structures. The third phase (bringing
the dose up to 68 or 74Gy(RBE)) is delivered to a reduced target volume with the IMPT
plan aiming to spare the critical structures. In all these cases, the PTV is symmetrically
grown from the CTV by 5mm. The dose constraints are mainly defined for the brainstem
(64Gy(RBE) on the surface, 53Gy(RBE) in the centre) and 60Gy(RBE) in the optical struc-
tures (Ares et al., 2009). In the IMPT plans, all the fluences of all the fields are optimized
simultaneously. Unlike SFUD treatments, IMPT can deliver a number of non-uniform fields
to produce the desired dose distribution.
The standard 3D-IMPT beam arrangement used clinically at PSI for the skull-base is a four
field beam arrangement and is illustrated in Figure 5.2. This beam arrangement consists of
two posterior oblique beams and two lateral oblique beams. The IMPT plans retrospectively
analysed in this study have been planned using these start conditions. For the purpose of this
chapter this beam arrangement is referred to as the A(4f) beam arrangement. For each IMPT
phase, the plan robustness in the CTV for an underdose scenario was analysed for both types
of uncertainty, as well as the overdose scenario in the brainstem and optical chiasm. Due
to their proximity to the target, the dose constraints defined for the brainstem and for the
optical chiasm tend to drive the optimisation outcome. The parameters measured include
the mean and maximum dose uncertainties in the volume, the volume with an uncertainty of
86 Defining Robustness Protocols
3% and 5% (Veb3% and Veb5%), and where these uncertainties occur in relation to anatomy
and field direction. From this analysis a table of robustness parameters was created to aid
planners in the selection between robustness and conformality.
5.1 Introduction 87
Fig. 5.2 The A(4f) plan consists of four non-coplanar fields, two posterior oblique and two
anterior lateral oblique fields. Clockwise from top left the beam angles and couch twists
used are (290◦, 340◦), (70◦, 20◦), (110◦, 20◦) and (250◦, 340◦), where the first angle is the
gantry orientation and the second is the table orientation.
Example Case
To illustrate how the robustness table could be clinically used in practice, one plan was
identified as an outlier. This plan was found to be un-robust in the brainstem compared to
the sample. Due to the degeneracy associated with IMPT there exist many solutions to the
given problem of ensuring target coverage whilst meeting dose and robustness constraints.
Beyond changing the optimisation itself (Chen et al., 2012; Pflugfelder, 2008; Unkelbach
et al., 2007), we can alter our starting conditions, such as number of beams and their ori-
entation to obtain solutions to satisfy both dosimetric and robustness requirements. For this
case, the plan was re-planned using four different beam arrangements. Four additional start
arrangements were investigated to determine which would be more optimal for this patient,
included B(4f), B(6f), A(6f) and C(4f) shown in Table 1 along with A(4f). Each plan was
created with the aim of keeping the dosimetric quality of the plan in the OARs and target
as similar to those achieved in the A(4f) plan so that the degeneracy of robustness to start
conditions could be investigated.
88 Defining Robustness Protocols
Table 5.1 Descriptions of each plan
Plan No. fields Beam arrangement
A(4f) 4 Two posterior oblique and two
anterior lateral oblique
B(4f) 4 Two lateral and two
posterior lateral oblique
A(6f) 6 Two posterior oblique,
two anterior lateral oblique and
two lateral fields
B(6f) 6 Two posterior oblique,
two anterior lateral oblique and
two posterior lateral obliques
C(4f) 4 Two lateral fields and two
oblique fields all coplanar
5.2 Results
5.2.1 Retrospective Robustness Analysis
Sixteen skull-base 3D IMPT plans, phase two and three, were retrospectively analysed in
terms of robustness to systematic range and random set-up uncertainties. The volumes
analysed were the brainstem, chiasm and the clinical target volume (CTV). The results are
shown in Figure 5.3, where for each box plot the central mark is the median, the edges of the
box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points
not considered outliers, and outliers are plotted individually. Outliers are considered points
with a deviation more than±2.7σ . From the box plots in Figure 5.3 it is seen that the impact
of the random set-up error appears greater for an individual fraction than that resulting from
the range uncertainty in all three volumes of interest (VOIs) for each parameter measured.
It can also be seen that the mean and standard deviation in the range error was lower for the
OARs than for the CTV; the opposite was true for the set-up uncertainty.
5.2 Results 89
Fig. 5.3 Robustness results, represented as percentage dose uncertainty for 16 skull-base
IMPT plans from 8 patients to both the range and set-up uncertainties modelled in (a) in
the brainstem, (b) in the chiasm and (c) in the CTV. On each box plot the central mark is
the median, the edges of the box are the 25th and 75th percentiles, the whiskers extend to
the most extreme data points not considered outliers, and outliers are plotted individually.
Outliers are considered points ±2.7σ .
90 Defining Robustness Protocols
From these results an example of a robustness database for 3D IMPT skull-base plans case
has been generated (Table 5.2). This table can be used to aid planners in visualising the
magnitude of effect that range and set-up uncertainties can have on patient treatment in
terms of percentage change from the nominal dose. The ranges represent the deviation from
the nominal dose in that volume caused by the uncertainty for this patient sample. The
values correspond to the 25th -75th percentiles as seen in the box plots in figure 5.3 for
an over-dosage scenario to OAR and under-dosage scenario to the CTV. The lower limits
are included to allow the planner to see where a compromise between conformality and
robustness can be made.
Table 5.2 The use of retrospective plan analysis to establish planning parameters of robust-
ness
VOI Mean Mean Veb3% Veb5% Veb5% Veb10%
range setup Range Range Setup Setup
Brainstem 1-2% 6.5-9% 2.5-17% 5% 63-85% 15-37.5%
Chiasm 1-1.5% 8-15% <5% 0% (Veb15%) (Veb17%)
<50% <4%
CTV 1.5-2.5% 5-8% 17.5-28% 2.5-11.5% 32.5-72.5% 3-25.8%
5.2.2 How to use the Robustness Database in Practice
To illustrate how the robustness database could be clinically used in practice a clinical case
is discussed. For patient X a nominal plan was generated and the robustness to set-up and
range uncertainties was calculated and compared to the robustness database table. If the
robustness of plan X is within the defined values, the nominal plan is delivered. In contrast,
if plan X has a robustness outside the database range for a given set of parameters for one
or more volumes, the nominal plan has to be improved in terms of its robustness before
delivering it to the patient. Or in the case where robustness is better than expected from the
database it must be assured that plan conformality is not being compromised.
Of the plans retrospectively analysed, one plan was shown to have a worse robustness in the
brainstem compared to the sample. In particular, this plan had a potential mean range, Veb3%
and Veb5% uncertainty in the brainstem of 3.9%, 43.5% and 15% respectively (compared to
1-2%, 2.5-17% and 5% in the database in Table 5.2). The possibility of improving the plan
5.2 Results 91
robustness, without using robust optimisation, but by changing the start conditions of the
plan has been investigated. While ensuring plan quality was maintained, the field setup was
changed to create four dosimetrically equivalent plans. For each new plan, the robustness
and conformality of the CTV and each OAR, considered during the optimisation process,
were analysed using the error-bar dose histogram and by retrieving the robustness values for
both the range and setup uncertainty from these (see Tables 5.3 and 5.4).
The error-bar dose histograms in the CTV for each plan created for this patient are shown
in Figure 5.4. The numbers in the plot key refers to the size of the area under each curve
(the smaller the area, the greater the robustness). The B(4f) plan in Figure 5.4(a) shows
worse robustness in the CTV than the other plans. A comparison between the B(4f) and
A(4f) nominal and error-bar dose distributions are shown in Figure 5.5. The B(4f) plan
in Figure 5.4(a) shows worse robustness in the CTV than the other plans. A comparison
between the B(4f) and A(4f) nominal and error-bar dose distributions are shown in Figure
5.4(c). The B(4f) plan produces an underdose uncertainty to a larger volume of the target
than the A(4f) plan. Table 5.3, shows that despite B(4f) having the worst plan robustness it
only fails one robustness constraint, the mean range uncertainty which is 2.8% for the B(4f)
plan compared to 2.1% for the A(4f) plan. The constraint set in Table 5.2 was 1.5-2.5%.
Figure 5.4(b) shows the volume histograms from the set-up uncertainty, A(4f) is the most
robust plan, though plan B(4f) is within the robustness constraints in Table 5.3. The C(4f)
plan shows to be the most robust plan to range uncertainties, though it is the least robust
plan to set-up uncertainties.
92 Defining Robustness Protocols
Fig. 5.4 Volume Histogram of (a) systematic range uncertainties and (b) random set-up
uncertainties for an underdose scenario in the CTV for five plans of different beam orienta-
tions.
5.2 Results 93
Fig. 5.5 Dose distribution for the A(4f) (top left) and B(4f) (top right) plans and ebDD of
systematic range uncertainties for the worst-case underdose scenario for the A(4f) (bottom
left) and B(4f)(bottom right) treatment plans showing the CTV in green.
94 Defining Robustness Protocols
Figure 5.7(a) shows the volume histograms of error-bars from range uncertainties in the
brainstem from all plans. The A(4f) plan is least robust to range uncertainties, as discussed
in the previous section of these results. Figure 5.6 shows the difference in ebDD between
the B(4f) and A(4f) plans indicating the greater uncertainty in the brainstem on the A(4f)
plan. Figure 5.7(c) is included to show how the plans are affected by random set-up uncer-
tainties, especially the C(4f) plan which shows a decrease in robustness when compared the
other four beam arrangements for this type of uncertainty. As seen in the A(4f) plan there
is a hot spot in the centre of the brainstem, whereas in the B(4f) the hot spot is shifted more
anteriorly due to the direction of the lateral and posterior beams. Despite this improvement
in robustness in the brainstem to range uncertainties, the B(4f) plan is less robust in the chi-
asm (Figure 5.7(c)) compared to the A(4f), whereas the B(6f) yields the greatest robustness
to range uncertainties in the chiasm of all the plans. In the case of the optic nerves again the
B(4f) plan is less robust compared to the A(4f), but this time the A(6f) plan is most robust.
These results can also be seen in Table 5.3; the chiasm fails a robustness constraint with a
V3% range uncertainty of 23%, though this uncertainty in dose has fallen to 1.4% at V5%,
matching that achieved by the A(4f) plan. Table 5.4 has been included to show the results
for four of the plans for other OARS, though constraints for these volumes have not been
established currently.
5.2 Results 95
Fig. 5.6 Dose distributions from the A(4f) (top left) and B(4f) (top right) plans and the corre-
sponding ebDD from B(4f) (bottom left) and A(4f) (bottom right) treatment plans showing
the brainstem in green and the target in yellow.
96 Defining Robustness Protocols
Fig. 5.7 Volume histogram of (a) the systematic range uncertainties, (b) the random set-up
uncertainties for an overdose scenario in the brainstem and (c) the random set-up uncertain-
ties for overdose scenario in the chiasm, for five plans of different beam orientations.
5.2 Results 97
Table 5.3 Robustness values for each parameter set in the site specific robustness database
in Table 5.2 for plans A(4f) and B(4f). Bold italic values are those that are out of tolerance
VOI Mean Mean Veb3% Veb5% Veb5% Veb10%
range setup Range Range Setup Setup
A(4f) Brainstem 3.9% 9% 43.5% 15% 84% 38%
B(4f) Brainstem 1.1% 8.7% 1.1% 1.1% 85% 33.5%
A(4f) Chiasm 1.4% 6% 2% 3.9% 1.8% 0%
B(4f) Chiasm 2.4% 8.7% 23% 1.4% 4% 0%
A(4f) CTV 2.1% 5.2% 16.5% 3% 36.5% 7%
B(4f) CTV 2.8% 5.3% 26% 8.5% 38% 10%
98 Defining Robustness Protocols
Table 5.4 This table shows robustness results from the original A(4f) plan that patient X
was treated with and those from the re-plan created for this work B(4f). Tolerances have
not been set for these OARs as they had for the brainstem and chiasm in the robustness
database, but they have been supplied here to use to further compare plans.
VOI Mean Mean Veb3% Veb5% Veb5% Veb10%
Beams range setup Range Range Setup Setup
Cochlea
B(4f) 4 6.8 64.2 32.1 68.2 5.5
A(4f) 3.8 4 84.4 24.8 20.2 1.8
B(6f) 3.7 5.8 57.8 22.9 75.2 1.8
A(6f) 3.9 3.4 75.2 24.8 10.1 1.8
Left Optic Nerve
B(4f) 3.8 5.8 36 19 47 14.2
A(4f) 3.8 5.8 36 19 47 14.2
B(6f) 4.4 4.9 23.7 8.3 38.3 14.2
A(6f) 2.1 5.8 9.5 1.2 47.5 12.5
Left Temporal Lobe
B(4f) 4.2 4.3 15.3 10.6 19.1 4.8
A(4f) 3 3.8 11.4 6.8 12.9 2
B(6f) 3.2 3.8 11.7 7.2 15.3 2.6
A(6f) 3.1 4 12.5 7.8 16.5 2.1
5.3 Discussion 99
5.3 Discussion
In this work it has been demonstrated how retrospective analysis of the robustness of clinical
plans at PSI has furthered an understanding of how random set-up and systematic range
uncertainties affect the final dose distribution. It was seen that random set-up uncertainties
resulted in a greater difference in the resulting dose distribution than that from systematic
range uncertainties. Despite this, systematic range uncertainties pose a greater concern due
to their systematic origin meaning they are present over the entire treatment leading to a
cumulative ‘shift’ of the dose distribution. In contrast, set-up uncertainties are random and
act instead to ‘blur’ the final dose distribution; therefore, the risk of reduced plan quality
associated with these uncertainties is less. This analysis assumes that systematic set-up
uncertainties are negligible. This is the case where a remote patient position device (Bolsi
et al., 2008) is used, but might not be true in centres using other set-up techniques. It
can, however, be identified and corrected with the use of image guidance techniques. The
mean values seen in the sample data show that the brainstem, chiasm and CTV appear to
be robust to the systematic range uncertainty with values of around 2% of the planned dose.
Of these volumes the target is the least robust to systematic range uncertainties. In the 3D
IMPT optimisation OAR sparing has a higher priority over target coverage. This can lead to
individual fields showing steep dose gradients within or near to the target which can result in
a greater ebDD for both range and set-up uncertainties. Through retrospectively analysing
the robustness of each VOI a site specific robustness database was established (Table 5.2). It
has been shown how the database can be used to easily identify patients that require a more
individualized plan to improve the robustness to uncertainties. In this case patient X may
have benefited from the B(4f) plan in regard to reducing range uncertainty in the brainstem.
Though for this beam arrangement there was a higher risk of range overdose in the chiasm,
this was outweighed by the improvement to the brainstem robustness and did not sacrifice
target robustness or plan conformality.
Although the defined database is specific to PSI, the idea can be used as an example for other
centres to define, for the moment, their own robustness parameters. Hopefully, with more
centres starting to introduce the robustness evaluation in the clinical process, we will be
able in the future to define an adequate level of robustness for each volume that is accepted
worldwide, in the same way as dose-volume-constraints to organs-at-risk and prescribed
dose to target volumes are generally widely accepted and applied in the different protocols.

Chapter 6
Summary
6.1 Discussions
Advances in radiotherapy, such as IMRT and IGRT, have meant that we are capable of creat-
ing and delivering highly sophisticated, and often individualised, treatment plans for many
cancer sites. These advancements are the direct result of improved technology, increased
computational power, access to published guidelines and recommendations, and also to an
extent, due to clinical experience. Theoretically, protons offer an even greater benefit for
radiotherapy patients, offering improved local tumour control and prevention or reduction
of radiation-induced side effects (Langendijk et al., 2013). As treatments progressively
move towards individualised care, balancing cost effectiveness and patient throughput is an
important factor for a health service. Within England and the NHS there is a far greater ex-
perience treating with photons in all aspects of the patient pathway. As well as prescribing
and planning itself, this also includes aspects such as carrying out patient specific quality
assurance. The greater experience and the technological advances in conventional radio-
therapy will mean for many patients protons will not necessarily be the optimal treatment
decision. Determining and prioritising indications for patient selection, to take into account
the limited capacity and high costs of proton therapy is an important challenge for NHS
England (Crellin and Burnet, 2014). Patients are currently selected based on the results
from dosimetric and observational studies, clinical trials, and models predicting reduced
side effects to ensure those treated are deemed to benefit the most (Langendijk et al., 2013).
Other challenges include establishing protocols for ensuring plan robustness and training
staff in proton therapy planning and delivery.
Optimisation in radiotherapy is a multi-parameter problem. The choice in optimisation
102 Summary
criteria and starting conditions can affect the resultant plan quality, both in its conformity
and robustness. The challenge of ensuring that plans are robust to range uncertainties in
proton therapy remains. There remains a need for established protocols for proton treatment
in the clinic for the English centres, in the same way that there already exists in conventional
radiotherapy. It has been to this end that this thesis has been carried out. Firstly, a method
of using daily IGRT to quantify changes in WEPL from inter-fraction motion has been
explored and demonstrated for several cancer sites. The quantification of uncertainties is
necessary in creating CTV-PTV (or beam specific) margins or as input into a form of robust
optimisation. Secondly, a method of robustness analysis has been described which can be
used with either PTV-based or robustly optimised plans. In this case by retrospectively
analysing plans a database of planning robustness objectives has been generated to ensure
that for each individual patient the optimal treatment is delivered.
Quantifying Random Range Uncertainty
The most significant limitations in calculating inter-fraction range uncertainty in this disser-
tation were from calibrating the TomoTherapy MVCT scanner and truncation in the imaging
data. Established proton centres such as Orsay, PSI and MD Anderson have been able to
sufficiently calibrate their TomoTherapy MVCT scanners using a stoichiometric calibration,
so these challenges would not limit other centres from using MVCT for range calculation.
At MD Anderson it was found that deviations between individual HU values and the fitted
calibration curve were smaller for the MVCT scanner, particularly within soft-tissue ranges
compared to that for the kVCT calibration. These results suggest that MVCT images may
be as accurate as kVCT images in calculating proton range (Newhauser et al., 2008), though
more research is required. If a centre were to adopt this approach of using MVCT to anal-
yse range uncertainty, calibrations would need to be carried out routinely. When analysing a
patient data set, the appropriate calibration must be applied. Site specific imaging protocols
would need to be established by widening or choosing certain slices to image on the patient
to avoid ROI truncation (Bates et al., 2013). Keeping the number of slices similar may also
be helpful, though not essential. The concomitant doses have been discussed in Chapter 3.
The concomitant dose is relatively small in comparison to the treatment, yet still requires
consideration and patient consent.
Despite these limitations the data collected in Chapter 4 from the head and neck patients was
useable for all seven patients. Though analysis of the prostate was carried out, significant
6.1 Discussions 103
numbers of fractions were missing. Missing fractions leads to an over estimate of the range
uncertainty for the sample (Greener, 2003). From the head and neck data a mean of -
1.82mm, a Σ of 3.39mm and a σ of 4.72mm was calculated for the patient sample across
all beam angles. For this method to be used clinically to produce uncertainty models to use
at either plan optimisation or evaluation a larger number of patients as well as clinical beam
angles need to be used.
The concept of deducing probability distributions of uncertainties from daily IGRT to be
used to re-optimise patient treatment is not new and has the potential to improve both TCP
and normal tissue sparing (Birkner et al., 2003). Birkner et al. (2003) re-optimised patient
treatment plans part way through treatment by including patient specific information from
portal images into the adaptive planning process. Unkelbach and Oelfke (2005) investi-
gated the same problem and used the information from daily IGRT of the first fractions to
establish probability distributions to be used in probabilistic IMRT plans. They found the
greatest limitation was the small number of images acquired to deduce the distribution. The
novelty of the work presented in this thesis is that the images have been acquired every day,
and analysed not to change that specific patient’s optimisation, but by collecting this data
for many patients with the same target sites so that probability distributions are acquired
representative of that population. This knowledge can prevent situations where plan confor-
mality may be sacrificed unnecessarily to ensure plan robustness to large or unlikely errors
(Unkelbach et al., 2007).
Part of the conclusion from the Gordon 2010 paper highlights the importance in setting the
objectives and constraints, both in terms of dose and the probability of that dose being re-
ceived by that voxel. This paper found that their robustness algorithm for producing photon
plans produced better plans than standard PTV-based plans. However, only in that optimised
plans exploited ‘slack’ in OAR doses, i.e., cases where OAR doses were below their optimi-
sation limits, to increase target coverage (Gordon et al., 2010). This is an example of when
target coverage has been prioritised over OAR sparing. Specifically, because the optimized
plans were not constrained by a predefined PTV, they were able to provide wider dosimetric
margins around the CTV, by pushing OAR doses up to, but not beyond, their optimisation
limits (Gordon et al., 2010). For such optimisation, it is imperative constraints and objec-
tives are clearly predefined prior to optimisation. The power of optimisation algorithm will
exploit any incoherence in the objective request. Unlike the planner, the algorithm will not
know what the doctor intended, just what has been asked of it.
104 Summary
Further work includes investigating how the WEPL between the distal edge of the target and
the proximal edge of an OAR in the beam direction as well as increasing patient numbers in
the samples and choosing clinical beam arangments.
Evaluating Robustness
Beyond the ability to determine the nature of uncertainties from a patient sample and then
create treatment plans robust to these uncertainties for a given percentage of treatments,
there still exists the problem in how to define what is robust enough. Due to the degen-
eracy in the inverse optimisation problem many optimal plans may be generated. The use
of a Pareto navigation surface (Chen et al., 2012) serves as an ideal way for the planner,
whether oncologist, physicist or dosimetrist, to gain experience and an ‘intuition’ for proton
planning that they may already have in conventional radiotherapy planning. When using
MCO consideration for the balance between the computationally heavy option of calculat-
ing many plans or interpolating between a small number of plans along the Pareto surface
is required, as well as an awareness of how plans on the surface are generated (Breedveld
et al., 2007). A principal advantage of probabilistic treatment planning is that it removes
the need for artificial planning volumes such as a PTV or PRV. This allows the treatment
planner to deal directly with dose distributions and metrics that incorporate the effects of
uncertainties (Gordon et al., 2010).
Choosing the specific plan to deliver relies on qualitative plan analysis, which may be sub-
ject to the planner’s experience. In chapter 5 I have sought for a quantitative method for
choosing a plan when balancing robustness and conformity, as is the crux in the multi-
criteria problem. The error-bar Dose Distribution (ebDD) was used to analyse plan robust-
ness to both systematic range and random set-up uncertainties for 16 skull-base IMPT plans
previously treated at PSI for 8 patients. A site-specific robustness database was created as
a simple solution to include robustness into the plan analysis. Such a database will also
aid planners to produce plans to meet both dosimetric and robustness criteria by establish-
ing site-specific robustness thresholds for individual volumes. Using these methods further
work can be carried out to fully explore how start conditions, and the optimisation itself,
affect both the plan robustness and conformity. This will further our ability to determine
where a trade-off can be made between these two parameters to ensure the patient receives
as optimal treatment as possible.
This quantitative method also serves as an attractive option to use in the clinic due to the
6.1 Discussions 105
similarity of the ebDVH, proposed by Albertini et al. (2011b), and robustness objectives
when compared with the familiar DVH and dose constraints used routinely in conventional
radiotherapy. Many patients have already been treated with actively scanned protons, so
using data from these patients through retrospective analysis will allow for continually up-
dating practices and therefore, continually improving practices. So far only a few centres
worldwide are performing a robustness analysis and, to the best of my knowledge, no-one
is using this information clinically. Additionally every centre is performing the plan robust-
ness analysis differently, for example by looking at the error bars, at the standard deviation,
or at the width of the DVHs. With the introduction of the robustness database as laid out in
Chapter 5 it is hoped that more centres will soon start to analyse and report plan robustness
uniformly. The ultimate goal to be able to define in the future a reference value for the plan
robustness for individual target sites and for specific OARs. In chapter 5 it was demon-
strated how, even using a small patient sample, an outlier was identified as lacking in plan
robustness when compared to the sample group. All of the plans in this study were clinically
acceptable plans and had been delivered to the patients. The clinical results at PSI achieved
for patients treated intra-cranially for chordoma and chondrosarcomas had a five year local
control of 81% for chordomas and of 94% for chondrosarcomas with only limited toxicities
(Ares et al., 2009). Plan robustness for these patients had only been retrospectively anal-
ysed; consequently the plan conformity of these plans had never been altered in favour to a
more robust plan. Without employing optimisation techniques, a plan using different start
conditions could be created that satisfied this specific patient’s need for a plan of greater in-
dividualisation. The use of a database is envisioned to be implemented alongside the use of
robust optimisation techniques to use as part of the planning workflow. The database offers
a novel tool to evaluate plan robustness as well as providing feedback that could be used to
alter robustness constraints at optimisation itself. One of the challenges of using any robust
optimisation is to find a threshold between robustness and nominal plan quality. It is hoped
that using the database to collect and analyse data from all patients treated, for many sites,
will help to find the threshold.
Further work will be to include random range uncertainty into the ebDD model and to
analyse plans that have been created using a robust optimisation as well as margin based
plans.
A final point is that there exist many optimisation algorithms in inverse planning and within
both probabilistic planning and robust optimisation themselves. These algorithms have not
106 Summary
been a topic of this work, nor has it been to investigate how exactly systematic and random
errors in proton therapy are weighted relatively at optimisation. It was seen in Chapter 5
how random range uncertainty had a greater effect on the individual fraction dose compared
to the systematic set-up uncertainty. However, over the course of the entire treatment the
systematic uncertainty would have a cumulative effect. Due to the directional nature of
range uncertainty it is imagined that the position of a voxel relative to the beam direction
should affect how the optimisation is driven.
6.2 MVCT and Proton Therapy
Practises used in conventional radiotherapy, such as IGRT and adaptive planning are also
required for proton therapy. Accurate set-up is essential in radiotherapy, and even more so
in proton therapy due to the element of range. Each centre may use other methods of IGRT,
for example PSI uses in room CT scanner on rails, performing scout scans rather than a
full CT scan. The first proton specific commercial CBCT went clinical at Penn Medicine in
September 2014 (iba: proton therapy, 2014). CBCT is becoming routine in IMRT delivery
in the UK, as it is worldwide. This development in proton therapy can be seen as delayed
and highlights the important discrepancy between what can be achieved in the research
setting for protons and the clinical one. The methods of range calculation described in
this dissertation could be applied using either MVCT or CBCT. Dose re-computation using
CBCT in conventional radiotherapy have been discussed in Murray et al. (2014); Onozato
et al. (2014). In using CBCT to perform dose re-computation for adaptive radiotherapy
these studies require pixels in the CBCT to be modified into a small number of ranges.
Using masks for HUs is required due to instability of CBCT HUs. These methods may be
suitable for conventional radiotherapy, but currently CBCT is not adequate for proton dose
re-computation or for calculating WEPL.
In conventional radiotherapy daily MVCT allows, firstly for registration to the planning
kVCT so that online patient position correction may be applied (reducing dosimetric affects
from inter-fraction changes and systematic set-up error) and secondly, the daily MVCT may
be used to assess anatomical changes and the dosimetric implications (Langen et al., 2005).
This latter application is paramount to evaluate the necessity of adapting the treatment at
some point. This will allow for compensation of dose errors due to anatomical deforma-
tions and changes (Langen et al., 2005). In this dissertation changes in range in the CTV
between fractions has been calculation from daily MVCT images. As discussed the MVCT
6.2 MVCT and Proton Therapy 107
calibration poses as a significant barrier for such practices as dose recalculation. Solutions
may include creating a planning MVCT at the same time as the planning kVCT and calcu-
lating dose on both modalities. This would be so that the MVCT doses may be used, not to
plan treatment delivery, but to carry out dose analysis on daily MVCT images. There are
two methods in conventional radiotherapy to determine the dose distribution based on daily
MVCT. One method is dose re-computation using exit dose from the treatment itself and
recalculating it on the MVCT. The other is to recalculate the planned dose distribution onto
the MVCT. Obviously, exit proton dose calculation is not an option using MVCT. However
the dose delivered can be recomputed from the proton steering files and from the gantry
itself (Meiers et al., 2014). Recalculating the recomputed distribution onto the planning
CT acts as treatment verification, but by recalculating it onto daily volumetric image would
allow the opportunity to analyse treatment delivery at each fraction and for re-planning to
take place as and when required.
Another potential benefit for using MVCT is in the case where metal implants are in or near
the target site. Metal can cause dose calculation difficulties, as accurate patient dose relies
on CT HU. Uncertainties arise both due to the high atomic number in metal causing calibra-
tion uncertainty and due to the artefacts in the planning kVCT. The sharp density interface
also affects dose homogeneity as well as there being an increased range uncertainty caused
by motion when the proton beam passes through or near metal. Patients treated with tumours
along the spinal axis may often be referred for radiotherapy after surgery. When vertebrae
have been removed titanium implants are used to stabilise the spine. Prostate or pelvic pa-
tients may be referred for radiotherapy with metal hip prosthesis and currently patients with
dental fillings may require these teeth to be removed prior to proton therapy, or beam di-
rections avoiding the teeth are chosen. Rutz et al. (2007) and Staab et al. (2011) followed
40 patients treated for extra-cranial chordoma and chondrosarcoma between 1999 and 2006
at PSI. Twelve of the thirteen local failures observed in this patient sample occurred when
stabilising implants were present. The presence of implants was a significant risk factor
(P<0.01) for local failure. Dietlicher et al. (2014) states that whether there is a clinical rea-
son for this risk, the presence of metal causes challenges in radiotherapy treatment, and even
more so in proton therapy. This group then went on to investigate the effectiveness of their
artefact correction technique using an anthromorphic phantom to reproduce calculations of
an in vivo measurement. The results from this paper suggest the poor survival for patients
with metal may not be due to the limitations of dose calculation alone, but due to the aggres-
siveness of the local tumour for this patient group. Their results also reinforce a need for
108 Summary
multiple fields in proton therapy as individual fields showed a worse gamma analysis, with
actual measured differences in dose up to -17.9% behind metal. The authors suggest there
is still a need for an improvement in artefact correction at treatment planning. Newhauser
et al. (2008) investigated using an MVCT and kVCT hybrid scan at planning to avoid the
need to correcting metal streaking and saturation artefacts. In addition to avoiding streak
artefacts, the MV-based approach yields advantages related to its greater dynamic range
(i.e. higher saturation threshold) of relative linear stopping power values than that of kV
scanning. Specifically, MVCT avoids saturation effects, enabling more accurate delineation
of the metal implant and thus theoretically enabling quantitative identification of the im-
plant material. This hybrid scan method offers a step towards reducing the subjectivity and
manual labour associated with correcting for implant-induced artefacts (Newhauser et al.,
2008). However, further comparison with streak- resistant kVCT reconstruction algorithms
(McLean, 2000; Yazdi et al., 2005) and proton CT (Sadrozinski et al., 2004) is needed
Furthermore, the use of DECT has been proposed by several groups to determine stopping
power ratios with greater accuracy (Bourque et al., 2014). The kV-MV combination makes
the DECT method more robust in resolving the effective atomic numbers for biological
tissues than the traditional kV-kV DECT method (Yang et al., 2011). Despite improvements
MVCT does produce uncertainties in stopping power ratios due to CT number uncertainties
and artefacts such as imaging noise, scatter and beam hardening effects.
6.3 Conclusion 109
6.3 Conclusion
In this dissertation the need and decisions behind establishing proton therapy in England
for NHS patients have been introduced, as have the concepts and benefits of this treatment
option. In chapter 2 the state of the current literature at the time of this work has been
presented and critically analysed. The following chapters describe the methods, results and
discussions of the subsequent work following analysis of the literature.
The importance of understanding and quantifying range uncertainty in proton therapy has
been thoroughly established and the need to quantify uncertainties to use as input into robust
optimisation has been stated. Two methods of analysis have been presented that can aid
in providing this information, namely through the use of IGRT, and by retrospective plan
analysis. Novel work has been carried out using daily MVCT data to produce probability
distributions of range uncertainty from inter-fraction motion. This knowledge is necessary
for accurate proton treatment planning.
Currently no established robustness protocols exist in proton therapy in the same way that
they do in conventional radiotherapy. Through analysing previously delivered proton plans
a database of plan robustness for different ROIs may be established. This type of database
may be used within the planning workflow as part of the plan analysis or as another method
for acquiring prior information required for probabilistic/robust planning.
The two proposed English centres, together with an additional centre recently announced at
Oxford (University of Oxford, 2014), offer the ability for the UK to make advancements in
this relatively new field of radiotherapy treatment. Further work is needed to bring IGRT and
ART in the clinical setting up to date with that which is already achieved in most state of the
art conventional radiotherapy centres. I believe there is a place for MVCT in proton therapy,
whether it is used for daily correction, dose difference analysis, DECT or for the study of
range uncertainty for different tumour sites and beam orientations. It is also essential that
methods for including uncertainties and evaluating plan robustness to these uncertainties
(without the use of a PTV) are established and become integrated into the treatment planning
workflow.

References
Albertini, F., Planning and optimizing treatment plans for actively scanned proton therapy:
evaluating and estimating the effect of uncertainties., Ph.D. thesis, 2011.
Albertini, F., E. B. Hug, and a. J. Lomax, The influence of the optimization starting condi-
tions on the robustness of intensity-modulated proton therapy plans., Physics in medicine
and biology, 55(10), 2863–78, doi:10.1088/0031-9155/55/10/005, 2010.
Albertini, F., M. Casiraghi, S. Lorentini, B. Rombi, and a. J. Lomax, Experimental ver-
ification of IMPT treatment plans in an anthropomorphic phantom in the presence of
delivery uncertainties., Physics in medicine and biology, 56(14), 4415–31, doi:10.1088/
0031-9155/56/14/012, 2011a.
Albertini, F., E. B. Hug, and a. J. Lomax, Is it necessary to plan with safety margins for
actively scanned proton therapy?, Physics in medicine and biology, 56(14), 4399–413,
doi:10.1088/0031-9155/56/14/011, 2011b.
Allen, A. M., et al., An evidence based review of proton beam therapy: The report of
ASTRO’s emerging technology committee., Radiotherapy and oncology : journal of the
European Society for Therapeutic Radiology and Oncology, 103(1), 8–11, doi:10.1016/j.
radonc.2012.02.001, 2012.
Ares, C., E. B. Hug, A. J. Lomax, A. Bolsi, B. Timmermann, H. P. Rutz, J. C. Schuller,
E. Pedroni, and G. Goitein, Effectiveness and safety of spot scanning proton radiation
therapy for chordomas and chondrosarcomas of the skull base: first long-term report.,
International journal of radiation oncology, biology, physics, 75(4), 1111–8, doi:10.1016/
j.ijrobp.2008.12.055, 2009.
Bangert, M., P. Hennig, and U. Oelfke, Analytical probabilistic modeling for radiation
therapy treatment planning., Physics in medicine and biology, 58(16), 5401–19, doi:
10.1088/0031-9155/58/16/5401, 2013.
Bates, A. M., J. E. Scaife, G. S. J. Tudor, R. Jena, M. Romanchikova, J. C. Dean, A. C. F.
Hoole, M. P. D. Simmons, and N. G. Burnet, Image guidance protocols: balancing imag-
ing parameters against scan time, The British Journal of Radiology, 86(1032), 20130,385,
doi:10.1259/bjr.20130385, pMID: 24128423, 2013.
Bert, C., and M. Durante, Motion in radiotherapy: particle therapy., Physics in medicine
and biology, 56(16), R113–R144, doi:10.1088/0031-9155/56/16/R01, 2011.
112 References
Bert, C., N. Saito, A. Schmidt, N. Chaudhri, D. Schardt, and E. Rietzel, Target motion
tracking with a scanned particle beam, Medical Physics, 34(12), 4768, doi:10.1118/1.
2815934, 2007.
Bert, C., S. O. Grözinger, and E. Rietzel, Quantification of interplay effects of scanned
particle beams and moving targets., Physics in medicine and biology, 53(9), 2253–65,
doi:10.1088/0031-9155/53/9/003, 2008.
Besemer, A., H. Paganetti, and B. Bednarz, The clinical impact of uncertainties in the mean
excitation energy of human tissues during proton therapy., Physics in medicine and biol-
ogy, 58(4), 887–902, doi:10.1088/0031-9155/58/4/887, 2013.
Birkner, M., D. Yan, M. Alber, J. Liang, and F. Nusslin, Adapting inverse planning to patient
and organ geometrical variation: algorithm and implementation, Medical Physics, 30(10),
2822, doi:10.1118/1.1610751, 2003.
Bolsi, A., A. J. Lomax, E. Pedroni, G. Goitein, and E. Hug, Experiences at the Paul Scherrer
Institute with a remote patient positioning procedure for high-throughput proton radiation
therapy., International journal of radiation oncology, biology, physics, 71(5), 1581–90,
doi:10.1016/j.ijrobp.2008.02.079, 2008.
Booth, J. T., and S. F. Zavgorodni, Set-up error & organ motion uncertainty: a review.,
Australasian physical & engineering sciences in medicine / supported by the Australasian
College of Physical Scientists in Medicine and the Australasian Association of Physical
Sciences in Medicine, 22(2), 29–47, 1999.
Bortfeld, T., Optimized planning using physical objectives and constraints., Sem Rad Oncol,
(9(1)), 20?34, 1999.
Bortfeld, T., and W. Schlegel, An analytical approximation of depth - dose distributions
for therapeutic proton beams, Physics in Medicine and Biology, 41(8), 1331–1339, doi:
10.1088/0031-9155/41/8/006, 1996.
Bourque, A. E., J.-F. Carrier, and H. Bouchard, A stoichiometric calibration method for
dual energy computed tomography., Physics in medicine and biology, 59(8), 2059–88,
doi:10.1088/0031-9155/59/8/2059, 2014.
Bragg, W. H., On the Absorption of Alpha Rays and on the Classification of Alpha Rays
from Radium, Philosophical Magazine, 8, 719–725, 1904.
Breedveld, S., P. R. M. Storchi, M. Keijzer, A. W. Heemink, and B. J. M. Heijmen, A novel
approach to multi-criteria inverse planning for IMRT., Physics in medicine and biology,
52(20), 6339–53, doi:10.1088/0031-9155/52/20/016, 2007.
Burnet, N., E. Taylor, K. Kirkby, N. Thorp, R. Mackay, and M. G., Proton beam therapy
for cancer: An important development for patients with an explicit agenda., BMJ, (344),
2488, 2012.
Burnet, N. G., et al., Practical aspects of implementation of helical tomotherapy for
intensity-modulated and image-guided radiotherapy., Clinical oncology (Royal College
of Radiologists (Great Britain)), 22(4), 294–312, doi:10.1016/j.clon.2010.02.003, 2010.
References 113
Cancer Research UK, CS_KF_CHILDHOOD, 2014a.
Cancer Research UK, Cancer statistics newsletter, 2014b.
Cancer Research UK, UK Cancer Incidence 2011 by country summary, 2014c.
Caporaso, G. J., Y.-J. Chen, and S. E. Sampayan, The Dielectric Wall Accelera-
tor, Reviews of Accelerator Science and Technology, 02(01), 253–263, doi:10.1142/
S1793626809000235, 2009.
Carabe, A., M. Moteabbed, N. Depauw, J. Schuemann, and H. Paganetti, Range uncertainty
in proton therapy due to variable biological effectiveness., Physics in medicine and biol-
ogy, 57(5), 1159–72, doi:10.1088/0031-9155/57/5/1159, 2012.
Casiraghi, M., F. Albertini, and a. J. Lomax, Advantages and limitations of the ’worst case
scenario’ approach in IMPT treatment planning., Physics in medicine and biology, 58(5),
1323–39, doi:10.1088/0031-9155/58/5/1323, 2013.
Chang, J. Y., X. Zhang, X. Wang, Y. Kang, B. Riley, S. Bilton, R. Mohan, R. Komaki, and
J. D. Cox, Significant reduction of normal tissue dose by proton radiotherapy compared
with three-dimensional conformal or intensity-modulated radiation therapy in Stage I or
Stage III non-small-cell lung cancer., International journal of radiation oncology, biol-
ogy, physics, 65(4), 1087–96, doi:10.1016/j.ijrobp.2006.01.052, 2006.
Chen, W., J. Unkelbach, T. Alexei, T. Madden, H. Kooy, B. Thomas, and D. Craft, Including
robustness in multi-criteria optimization for intensity-modulated proton therapy, Phys
Med Biol, pp. 591–608, 2012.
Chu, M., Robust intensity modulated radiation therapy treatment planning, Ph.D. thesis,
2006.
Chvetsov, A. V., and S. L. Paige, The influence of CT image noise on proton range calcu-
lation in radiotherapy planning., Physics in medicine and biology, 55(6), N141–9, doi:
10.1088/0031-9155/55/6/N01, 2010.
Cotrutz, C., M. Lahanas, C. Kappas, and D. Baltas, A multiobjective gradient-based dose
optimization algorithm for external beam conformal radiotherapy., Physics in medicine
and biology, 46(8), 2161–75, 2001.
Craft, D., T. Halabi, and T. Bortfeld, Exploration of tradeoffs in intensity-modulated radio-
therapy., Physics in medicine and biology, 50(24), 5857–68, doi:10.1088/0031-9155/50/
24/007, 2005.
Crellin, A., and N. Burnet, Proton Beam Therapy: The Context, Future Direction and Chal-
lenges Become Clearer, Clinical Oncology, 26(12), 736–738, doi:10.1016/j.clon.2014.
10.009, 2014.
De Ruysscher, D., M. Mark Lodge, B. Jones, M. Brada, A. Munro, T. Jefferson, and
M. Pijls-Johannesma, Charged particles in radiotherapy: a 5-year update of a systematic
review., Radiotherapy and oncology : journal of the European Society for Therapeutic
Radiology and Oncology, 103(1), 5–7, doi:10.1016/j.radonc.2012.01.003, 2012.
114 References
Deasy, J. O., A proton dose calculation algorithm for conformal therapy simulations based
on Molière’s theory of lateral deflections., Medical physics, 25(4), 476–83, 1998.
DeLaney, T. F., and H. M. Kooy, Proton and Charged Particle Radiotherapy, vol. 83, 87
pp., Lippincott Williams & Wilkins, 2008.
Department of Health, National Proton Addendum, 2014.
Dietlicher, I., M. Casiraghi, C. Ares, A. Bolsi, D. C. Weber, A. J. Lomax, and F. Albertini,
The effect of surgical titanium rods on proton therapy delivered for cervical bone tumors:
experimental validation using an anthropomorphic phantom., Physics in medicine and
biology, 59(23), 7181–94, doi:10.1088/0031-9155/59/23/7181, 2014.
Dosanjh, M., G. H. Corral, and L. M. M. Zentina, Development of Hadron Therapy for
Cancer Treatment in Europe, AIP Conference Proceedings, pp. 12–16, doi:10.1063/1.
2979248, 2008.
Duchateau, M., et al., The effect of tomotherapy imaging beam output instabilities on dose
calculation., Physics in medicine and biology, 55(11), N329–36, doi:10.1088/0031-9155/
55/11/N03, 2010.
Durante, M., and J. S. Loeffler, Charged particles in radiation oncology., Nature reviews.
Clinical oncology, 7(1), 37–43, doi:10.1038/nrclinonc.2009.183, 2010.
Engels, B., G. Soete, D. Verellen, and G. Storme, Conformal arc radiotherapy for prostate
cancer: increased biochemical failure in patients with distended rectum on the planning
computed tomogram despite image guidance by implanted markers., International jour-
nal of radiation oncology, biology, physics, 74(2), 388–91, doi:10.1016/j.ijrobp.2008.08.
007, 2009.
España, S., and H. Paganetti, The impact of uncertainties in the CT conversion algorithm
when predicting proton beam ranges in patients from dose and PET-activity distributions.,
Physics in medicine and biology, 55(24), 7557–71, doi:10.1088/0031-9155/55/24/011,
2010.
Fredriksson, A., A. Forsgren, and B. Hardemark, Minimax optimization for handling range
and setup uncertainties in proton therapy, Medical Physics, 38(3), 1672, doi:10.1118/1.
3556559, 2011.
Furukawa, T., T. Inaniwa, S. Sato, T. Tomitani, S. Minohara, K. Noda, and T. Kanai, Design
study of a raster scanning system for moving target irradiation in heavy-ion radiotherapy.,
Medical physics, 34(3), 1085–97, 2007.
Gammex, Reference tablel, 2009.
Glimelius, B., and A. Montelius, Proton beam therapy - do we need the randomised trials
and can we do them?, Radiotherapy and oncology : journal of the European Society for
Therapeutic Radiology and Oncology, 83(2), 105–9, doi:10.1016/j.radonc.2007.04.009,
2007.
References 115
Gordon, J. J., N. Sayah, E. Weiss, and J. V. Siebers, Coverage optimized planning: Prob-
abilistic treatment planning based on dose coverage histogram criteria, Medical Physics,
37(2), 550, doi:10.1118/1.3273063, 2010.
Grassberger, C., A. Trofimov, A. Lomax, and H. Paganetti, Variations in linear energy
transfer within clinical proton therapy fields and the potential for biological treatment
planning., International journal of radiation oncology, biology, physics, 80(5), 1559–66,
doi:10.1016/j.ijrobp.2010.10.027, 2011.
Greener, T., Geometric Uncertainties in Radiotherapy. Appendix 2c, Practical determina-
tion of systematic and random set-up errors, Σset−up and σset−up using Portal Imaging,
36-43 pp., BIR, 2003.
Grözinger, S. O., E. Rietzel, Q. Li, C. Bert, T. Haberer, and G. Kraft, Simulations to design
an online motion compensation system for scanned particle beams., Physics in medicine
and biology, 51(14), 3517–31, doi:10.1088/0031-9155/51/14/016, 2006.
Grözinger, S. O., C. Bert, T. Haberer, G. Kraft, and E. Rietzel, Motion compensation with a
scanned ion beam: a technical feasibility study., Radiation oncology (London, England),
3, 34, doi:10.1186/1748-717X-3-34, 2008.
Haberer, T., W. Becher, D. Schardt, and G. Kraft, Magnetic scanning system for heavy ion
therapy, Nuclear Instruments and Methods in Physics Research Section A: Accelerators,
Spectrometers, Detectors and Associated Equipment, 330(1-2), 296–305, doi:10.1016/
0168-9002(93)91335-K, 1993.
Harris, V., C. South, C. Cruickshank, and D. Dearnaley, 1299 poster A NATIONAL PHASE
III TRIAL OF PELVIC LYMPH NODE (LN) IMRT IN PROSTATE CANCER (PIV-
OTAL): A COMPARISON OF LN OUTLINING METHODS, Radiotherapy and Oncol-
ogy, 99, S486–S487, 2011.
Hünemohr, N., B. Krauss, C. Tremmel, B. Ackermann, O. Jäkel, and S. Greilich, Experi-
mental verification of ion stopping power prediction from dual energy CT data in tissue
surrogates., Physics in medicine and biology, 59(1), 83–96, doi:10.1088/0031-9155/59/
1/83, 2014.
iba: proton therapy, World’s First PT Specific CBCT Goes Clinical, 2014.
ICRP, The 2007 Recommendations of the International Commission on Radiological Pro-
tection, Ann. ICRP 37, ICRP Publication 103, 2–4, 2007.
ICRU 50, International Commission on Radiation Units and Measurements: Prescribing,
recording, and reporting photon beam therapy. ICRU Report no. 50, Tech. rep., ICRU,
1993.
ICRU 62, International Commission on Radiation Units and Measurements: Prescribing,
recording, and reporting photon beam therapy. ICRU Report no. 62, Tech. rep., ICRU,
1999.
ICRU 78, ICRU Report 78: Prescribing, Recording, and Reporting Proton Beam Therapy ,
Tech. rep., ICRU, 2007.
116 References
ICRU 83, The ICRU Report 83: prescribing, recording and reporting photon-beam intensity-
modulated radiation therapy (IMRT). ICRU Report no. 83, Tech. rep., ICRU, 2010.
Jette, D., and W. Chen, Creating a spread-out Bragg peak in proton beams., Physics in
medicine and biology, 56(11), N131–8, doi:10.1088/0031-9155/56/11/N01, 2011.
Jones, B., The potential clinical advantages of charged particle radiotherapy using protons
or light ions., Clinical oncology (Royal College of Radiologists (Great Britain)), 20(7),
555–63, doi:10.1016/j.clon.2008.02.012, 2008.
Khan, F. M. (Ed.), The Physics of Radiation Therapy, 4th ed., 74 pp., 2010.
Kirby, D., Radiation dosimetry of conventional and laser driven particle beams., University
of Birmingham, PhD Thesis, 2010.
Knopf, A.-C., T. S. Hong, and A. Lomax, Scanned proton radiotherapy for mobile targets-
the effectiveness of re-scanning in the context of different treatment planning approaches
and for different motion characteristics., Physics in medicine and biology, 56(22), 7257–
71, doi:10.1088/0031-9155/56/22/016, 2011.
Lambert, J., N. Suchowerska, D. R. McKenzie, and M. Jackson, Intrafractional motion
during proton beam scanning., Physics in medicine and biology, 50(20), 4853–62, doi:
10.1088/0031-9155/50/20/008, 2005.
Langen, K., N. Papanikolaou, and J. Balog, QA for helical TomoTherapy: Report of the
AAPM Task Group 148, Med Phys, (37), 4817–53, 2010.
Langen, K. M., and D. T. Jones, Organ motion and its management., International journal
of radiation oncology, biology, physics, 50(1), 265–78, 2001.
Langen, K. M., S. L. Meeks, D. O. Poole, T. H. Wagner, T. R. Willoughby, P. A. Kupelian,
K. J. Ruchala, J. Haimerl, and G. H. Olivera, The use of megavoltage CT (MVCT) im-
ages for dose recomputations., Physics in medicine and biology, 50(18), 4259–76, doi:
10.1088/0031-9155/50/18/002, 2005.
Langendijk, J. A., P. Lambin, D. De Ruysscher, J. Widder, M. Bos, and M. Verheij, Selection
of patients for radiotherapy with protons aiming at reduction of side effects: the model-
based approach., Radiotherapy and oncology : journal of the European Society for Ther-
apeutic Radiology and Oncology, 107(3), 267–73, doi:10.1016/j.radonc.2013.05.007,
2013.
Lawrence, E. O., Method and apparatus for the acceleration of ions, 1934.
Liu, W., X. Zhang, Y. Li, and R. Mohan, Robust optimization of intensity modulated proton
therapy., Medical physics, 39(2), 1079–91, doi:10.1118/1.3679340, 2012.
Lodge, M., M. Pijls-Johannesma, L. Stirk, A. J. Munro, D. De Ruysscher, and T. Jefferson,
A systematic literature review of the clinical and cost-effectiveness of hadron therapy in
cancer., Radiotherapy and oncology : journal of the European Society for Therapeutic
Radiology and Oncology, 83(2), 110–22, doi:10.1016/j.radonc.2007.04.007, 2007.
References 117
Lomax, A. J., Intensity modulation methods for proton radiotherapy, Physics in medicine
and biology, 44, 185–205, 1999.
Lomax, a. J., Intensity modulated proton therapy and its sensitivity to treatment uncertainties
2: the potential effects of inter-fraction and inter-field motions., Physics in medicine and
biology, 53(4), 1043–56, doi:10.1088/0031-9155/53/4/015, 2008a.
Lomax, a. J., Intensity modulated proton therapy and its sensitivity to treatment uncertainties
1: the potential effects of calculational uncertainties., Physics in medicine and biology,
53(4), 1027–42, doi:10.1088/0031-9155/53/4/014, 2008b.
Lomax, A. J., et al., Treatment planning and verification of proton therapy using spot scan-
ning: Initial experiences, Medical Physics, 31(11), 3150, doi:10.1118/1.1779371, 2004.
Matlab UK, The Language of Technical Computing, 2011.
Matsuura, T., et al., Apparent absence of a proton beam dose rate effect and possible dif-
ferences in RBE between Bragg peak and plateau, Medical Physics, 37(10), 5376, doi:
10.1118/1.3490086, 2010.
McGowan, S. E., C. D. Greaves, and S. Evans, An investigation into truncation artefacts
experienced in cardiac imaging using a dedicated cardiac SPECT gamma camera with
transmission attenuation correction., Nuclear medicine communications, 33(12), 1287–
91, doi:10.1097/MNM.0b013e328359db76, 2012.
McGowan, S. E., N. G. Burnet, and a. J. Lomax, Treatment planning optimisation in pro-
ton therapy., The British journal of radiology, 86(1021), 20120,288, doi:10.1259.bjr.
20120288, 2013.
McLean, D., Artefact reduction on CT images of fossils to allow 3D visualisation, 2000.
Meiers, G., A. Nanz, R. Besson, S. Safai, and A. Lomax, Logfile based dose calculations as
a quality assurance tool in scanned beam proton radiotherapy , ICTR, Geneva, 2014.
Merchant, T. E., C.-H. Hua, H. Shukla, X. Ying, S. Nill, and U. Oelfke, Proton versus
photon radiotherapy for common pediatric brain tumors: comparison of models of dose
characteristics and their relationship to cognitive function., Pediatric blood & cancer,
51(1), 110–7, doi:10.1002/pbc.21530, 2008.
Moyers, M. F., M. Sardesai, S. Sun, and D. W. Miller, Ion stopping powers and CT numbers.,
Medical dosimetry : official journal of the American Association of Medical Dosimetrists,
35(3), 179–94, doi:10.1016/j.meddos.2009.05.004, 2010.
Murray, J., D. McQuaid, A. Dunlop, F. Buettner, S. Nill, E. Hall, D. Dearnaley, and S. Gul-
liford, SU-E-J-14: A Novel Approach to Evaluate the Dosimetric Effect of Rectal Vari-
ation During Image Guided Prostate Radiotherapy, Medical Physics, 41(6), 157–157,
doi:10.1118/1.4888065, 2014.
Newhauser, W. D., A. Giebeler, K. M. Langen, D. Mirkovic, and R. Mohan, Can mega-
voltage computed tomography reduce proton range uncertainties in treatment plans for
patients with large metal implants?, Physics in medicine and biology, 53(9), 2327–44,
doi:10.1088/0031-9155/53/9/009, 2008.
118 References
Nguyen, T., A. Hoole, N. Burnet, and S. Thomas, The optimization of intensity modulated
radiotherapy in cases where the planning target volume extends into the build-up region,
Phys. Med. Biol., (54), 2511–2525, 2009.
NHS Specialised Services, NHS Specialised Services, 2011.
Nill, S., T. Bortfeld, and U. Oelfke, Inverse planning of intensity modulated proton therapy.,
Zeitschrift für medizinische Physik, 14(1), 35–40, 2004.
NRIG, International Commission on Radiation Units and Measurements: Prescribing,
recording, and reporting photon beam therapy. ICRU Report no. 50, Tech. rep., The Na-
tional Cancer Action Team, 2012.
Oelfke, U., and T. Bortfeld, Intensity modulated radiotherapy with charged particle beams:
Studies of inverse treatment planning for rotation therapy, Medical Physics, 27(6), 1246,
doi:10.1118/1.599002, 2000.
Onozato, Y., et al., Evaluation of on-board kV cone beam computed tomography-based
dose calculation with deformable image registration using Hounsfield unit modifications.,
International journal of radiation oncology, biology, physics, 89(2), 416–23, doi:10.1016/
j.ijrobp.2014.02.007, 2014.
Park, P. C., X. R. Zhu, A. K. Lee, N. Sahoo, A. D. Melancon, L. Zhang, and L. Dong, A
Beam-Specific Planning Target Volume (PTV) Design for Proton Therapy to Account for
Setup and Range Uncertainties, 2012.
Peach, K. J., et al., PAMELA Overview : design goals and principles, 2009.
Pedroni, E., S. Scheib, T. Böhringer, a. Coray, M. Grossmann, S. Lin, and a. Lo-
max, Experimental characterization and physical modelling of the dose distribution of
scanned proton pencil beams, Physics in Medicine and Biology, 50(3), 541–561, doi:
10.1088/0031-9155/50/3/011, 2005.
Pedroni, E., et al., The 200-MeV proton therapy project at the Paul Scherrer Institute: con-
ceptual design and practical realization., Medical physics, 22(1), 37–53, 1995.
Pflugfelder, D., Worst case optimization : a method to account for uncertainties in the opti-
mization of intensity, Cancer Research, 1689, doi:10.1088/0031-9155/53/6/013, 2008.
Phillips, C., What is a QALY?, 2009.
Phillips, M. H., E. Pedroni, H. Blattmann, T. Boehringer, A. Coray, and S. Scheib, Effects of
respiratory motion on dose uniformity with a charged particle scanning method., Physics
in medicine and biology, 37(1), 223–34, 1992.
Proton Therapy Today - The online magazine for proton therapy, Where to get pt?, 2014.
Ramaekers, B. L. T., M. Pijls-Johannesma, M. a. Joore, P. van den Ende, J. a. Langendijk,
P. Lambin, A. G. H. Kessels, and J. P. C. Grutters, Systematic review and meta-analysis of
radiotherapy in various head and neck cancers: comparing photons, carbon-ions and pro-
tons., Cancer treatment reviews, 37(3), 185–201, doi:10.1016/j.ctrv.2010.08.004, 2011.
References 119
Raysearch Laboratory, RayStation 4.5 news, 2014.
Rietzel, E., and C. Bert, Respiratory motion management in particle therapy, Medical
Physics, 37(2), 449, doi:10.1118/1.3250856, 2010.
Rimmer, Y. L., et al., Practical issues in the implementation of image-guided radiotherapy
for the treatment of prostate cancer within a UK department., Clinical oncology (Royal
College of Radiologists (Great Britain)), 20(1), 22–30, doi:10.1016/j.clon.2007.10.001,
2008.
Robertson, J. B., J. R. Williams, R. A. Schmidt, J. B. Little, D. F. Flynn, and H. D. Suit,
Radiobiological studies of a high-energy modulated proton beam utilizing cultured mam-
malian cells., Cancer, 35(6), 1664–77, 1975.
Rutz, H. P., et al., Extracranial chordoma: Outcome in patients treated with function-
preserving surgery followed by spot-scanning proton beam irradiation., International
journal of radiation oncology, biology, physics, 67(2), 512–20, doi:10.1016/j.ijrobp.2006.
08.052, 2007.
Sadrozinski, H.-W., et al., Toward Proton Computed Tomography, IEEE Transactions on
Nuclear Science, 51(1), 3–9, doi:10.1109/TNS.2003.823044, 2004.
Safai, S., S. Lin, and E. Pedroni, Development of an inorganic scintillating mixture for
proton beam verification dosimetry, Physics in Medicine and Biology, 49(19), 4637–
4655, doi:10.1088/0031-9155/49/19/013, 2004.
Saito, N., C. Bert, N. Chaudhri, A. Gemmel, D. Schardt, M. Durante, and E. Rietzel, Speed
and accuracy of a beam tracking system for treatment of moving targets with scanned ion
beams., Physics in medicine and biology, 54(16), 4849–62, doi:10.1088/0031-9155/54/
16/001, 2009.
Schaffner, B., and E. Pedroni, The precision of proton range calculations in proton radio-
therapy treatment planning: experimental verification of the relation between CT-HU and
proton stopping power., Physics in medicine and biology, 43(6), 1579–92, 1998.
Schaffner, B., E. Pedroni, and A. Lomax, Dose calculation models for proton treatment
planning using a dynamic beam delivery system: an attempt to include density hetero-
geneity effects in the analytical dose calculation, Physics in Medicine and Biology, 44(1),
27–41, doi:10.1088/0031-9155/44/1/004, 1999.
Schätti, A., M. Zakova, D. Meer, and A. J. Lomax, Experimental verification of motion mit-
igation of discrete proton spot scanning by re-scanning., Physics in medicine and biology,
58(23), 8555–72, doi:10.1088/0031-9155/58/23/8555, 2013.
Schätti, A., D. Meer, and A. J. Lomax, First experimental results of motion mitigation by
continuous line scanning of protons., Physics in medicine and biology, 59(19), 5707–23,
doi:10.1088/0031-9155/59/19/5707, 2014.
Scheib, S., and E. Pedroni, Dose calculation and optimization for 3D conformal voxel
scanning, Radiation and Environmental Biophysics, 31(3), 251–256, doi:10.1007/
BF01214833, 1992.
120 References
Schneider, U., E. Pedroni, and A. Lomax, The calibration of CT Hounsfield units for ra-
diotherapy treatment planning, Physics in Medicine and Biology, 41(1), 111–124, doi:
10.1088/0031-9155/41/1/009, 1996.
Schulte, R., et al., Conceptual design of a proton computed tomography system for applica-
tions in proton radiation therapy, IEEE Transactions on Nuclear Science, 51(3), 866–872,
doi:10.1109/TNS.2004.829392, 2004.
Seco, J., D. Robertson, A. Trofimov, and H. Paganetti, Breathing interplay effects during
proton beam scanning: simulation and statistical analysis., Physics in medicine and biol-
ogy, 54(14), N283–94, doi:10.1088/0031-9155/54/14/N01, 2009.
Sheehan, M., et al., Position statement on ethics, equipoise and research on charged
particle radiation therapy., Journal of medical ethics, 40(8), 572–5, doi:10.1136/
medethics-2012-101290, 2014.
Slater, J. D., C. J. Rossi, L. T. Yonemoto, D. A. Bush, B. R. Jabola, R. P. Levy, R. I. Grove,
W. Preston, and J. M. Slater, Proton therapy for prostate cancer: the initial Loma Linda
University experience., International journal of radiation oncology, biology, physics,
59(2), 348–52, doi:10.1016/j.ijrobp.2003.10.011, 2004.
Staab, A., et al., Spot-scanning-based proton therapy for extracranial chordoma., Interna-
tional journal of radiation oncology, biology, physics, 81(4), e489–96, doi:10.1016/j.
ijrobp.2011.02.018, 2011.
Stroom, J. C., H. C. de Boer, H. Huizenga, and A. G. Visser, Inclusion of geometri-
cal uncertainties in radiotherapy treatment planning by means of coverage probabil-
ity, International Journal of Radiation Oncology*Biology*Physics, 43(4), 905–919, doi:
10.1016/S0360-3016(98)00468-4, 1999.
Thieke, C., K.-H. Küfer, M. Monz, A. Scherrer, F. Alonso, U. Oelfke, P. E. Huber, J. De-
bus, and T. Bortfeld, A new concept for interactive radiotherapy planning with multi-
criteria optimization: first clinical evaluation., Radiotherapy and oncology : journal
of the European Society for Therapeutic Radiology and Oncology, 85(2), 292–8, doi:
10.1016/j.radonc.2007.06.020, 2007.
Thomas, S. J., Margins for treatment planning of proton therapy., Physics in medicine and
biology, 51(6), 1491–501, doi:10.1088/0031-9155/51/6/009, 2006.
Thomas, S. J., M. Ashburner, G. S. J. Tudor, J. Treeby, J. Dean, D. Routsis, Y. L. Rimmer,
S. G. Russell, and N. G. Burnet, Intra-fraction motion of the prostate during treatment
with helical tomotherapy., Radiotherapy and oncology : journal of the European Society
for Therapeutic Radiology and Oncology, 109(3), 482–6, doi:10.1016/j.radonc.2013.09.
011, 2013.
Towse, A., Cost effectiveness thresholds: economic and ethical issues, King?s Fund/Office
for Health Economics., London, 2002.
University of Oxford, £110m research centre in cancer medicine to be established at Oxford
University, 2014.
References 121
Unkelbach, J., and U. Oelfke, Incorporating organ movements in IMRT treatment planning
for prostate cancer: minimizing uncertainties in the inverse planning process., Medical
physics, 32(8), 2471–83, 2005.
Unkelbach, J., T. C. Y. Chan, and T. Bortfeld, Accounting for range uncertainties in the
optimization of intensity modulated proton therapy., Physics in medicine and biology,
52(10), 2755–73, doi:10.1088/0031-9155/52/10/009, 2007.
van de Water, T. a., H. P. Bijl, C. Schilstra, M. Pijls-Johannesma, and J. a. Langendijk, The
potential benefit of radiotherapy with protons in head and neck cancer with respect to
normal tissue sparing: a systematic review of literature., The oncologist, 16(3), 366–77,
doi:10.1634/theoncologist.2010-0171, 2011.
van Herk, M., Errors and margins in radiotherapy., Seminars in radiation oncology, 14(1),
52–64, doi:10.1053/j.semradonc.2003.10.003, 2004.
van Herk, M., P. Remeijer, C. Rasch, and J. V. Lebesque, The probability of correct tar-
get dosage: dose-population histograms for deriving treatment margins in radiotherapy,
International Journal of Radiation Oncology*Biology*Physics, 47(4), 1121–1135, doi:
10.1016/S0360-3016(00)00518-6, 2000.
Vees, H., G. Dipasquale, P. Nouet, T. Zilli, L. Cozzi, and R. Miralbell, Pelvic Lymph Node
Irradiation Including Pararectal Sentinel Nodes for Prostate Cancer Patients: Treatment
Optimization Comparing Intensity Modulated X-rays, Volumetric Modulated Arc Ther-
apy, and Intensity Modulated Proton Therapy., Technology in cancer research & treat-
ment, doi:10.7785/tcrt.2012.500405, 2014.
voxtox, The voxtox research programme - linking radiation dose at the voxel level with
toxicity, 2014.
Webb, S., The physical basis of IMRT and inverse planning, British Journal of Radiology,
76(910), 678–689, doi:10.1259/bjr/65676879, 2003.
www.dh.gov.uk, London, uk: Department of health; 2012, 2012.
Yadav, P., R. Tolakanahalli, Y. Rong, and B. R. Paliwal, The effect and stability of MVCT
images on adaptive TomoTherapy, Journal of Applied Clinical Medical Physics, 11(4),
2010.
Yang, M., G. Virshup, J. Clayton, X. R. Zhu, R. Mohan, and L. Dong, Does kV-MV
dual-energy computed tomography have an advantage in determining proton stopping
power ratios in patients?, Physics in medicine and biology, 56(14), 4499–515, doi:
10.1088/0031-9155/56/14/017, 2011.
Yang, M., X. R. Zhu, P. C. Park, U. Titt, R. Mohan, G. Virshup, J. E. Clayton, and L. Dong,
Comprehensive analysis of proton range uncertainties related to patient stopping-power-
ratio estimation using the stoichiometric calibration., Physics in medicine and biology,
57(13), 4095–115, doi:10.1088/0031-9155/57/13/4095, 2012.
122 References
Yazdi, M., M. Yazdia, L. Gingras, and L. Beaulieu, An adaptive approach to metal arti-
fact reduction in helical computed tomography for radiation therapy treatment planning:
experimental and clinical studies., International journal of radiation oncology, biology,
physics, 62(4), 1224–31, doi:10.1016/j.ijrobp.2005.02.052, 2005.
Zenklusen, S. M., E. Pedroni, and D. Meer, A study on repainting strategies for treating
moderately moving targets with proton pencil beam scanning at the new Gantry 2 at PSI.,
Physics in medicine and biology, 55(17), 5103–21, doi:10.1088/0031-9155/55/17/014,
2010.
Zietman, A. L., et al., Randomized trial comparing conventional-dose with high-dose con-
formal radiation therapy in early-stage adenocarcinoma of the prostate: long-term results
from proton radiation oncology group/american college of radiology 95-09., Journal of
clinical oncology : official journal of the American Society of Clinical Oncology, 28(7),
1106–11, doi:10.1200/JCO.2009.25.8475, 2010.
Appendix A
Published Work and Presentations
Presenated here is published work written during the course of this thesis.
Published review article in the Bristish Journal of Radiology. McGowan SE, Burnet NG,
Lomax AJ. (2013) Treatment planning optimisation in proton therapy The British
journal of radiology, 86(1021), 20120,288.
Oral presenation at the 20th International Conference on Medical Physics. McGowan SE,
Albertini F, Thomas SJ, Burnet NG & Lomax AJ. A retrospective study of pro-
ton treatment plan robustness in the skull base. The 20th International Conference
on Medical Physics (ICMP2013) at the Brighton Centre, Brighton, UK from 1-4th
September 2013
ePoster presentation at 2014 ESTRO Annual Meeting. McGowan SE, Albertini F, Thomas
SJ, & Burnet NG. (2014) Defining Robustness protocols: a method to include and
evaluate robustness in clinical plans.
Accepted scientific paper in Physics in Medicine and Biology Journal McGowan SE, Al-
bertini F, Thomas SJ, Lomax AJ. PMB (Accepted 2014) Defining robustness proto-
cols: a method to include and evaluate robustness in clinical plans
124 Published Work and Presentations
Fig. A.1 REVIEW ARTICLE: Treatment planning optimisation in proton therapy
125
Fig. A.2 Accepted ESTRO33 ePoster 2014
126 Published Work and Presentations
Defining robustness protocols: a method to include and evaluate robustness in clinical
plans
ABSTRACT Aim: To define a site-specific robustness protocol to be used during the clinical
plan evaluation process.
Method: Plan robustness of 16 skull base IMPT plans to systematic range and random set-up
errors has been retrospectively and systematically analyzed. This was determined by calcu-
lating the error-bar dose distribution (ebDD) for all the plans and by defining some metrics
used to define protocols aiding the plan assessment. Additionally, an example of how to
clinically use the defined robustness database is given whereby a plan with sub-optimal
brainstem robustness was identified. The advantage of using different beam arrangements
to improve the plan robustness was analysed.
Results: Using the ebDD it was found range errors had a smaller effect on dose distribution
than the corresponding set-up error, and that organs at risk were most robust to the range
errors, whereas the target was more robust to set-up errors. A database was created to aid
planners in terms of plan robustness aims in these volumes. This resulted in the defini-
tion of site specific robustness protocols. The use of robustness constraints allowed for the
identification of a specific patient that may have benefited from a treatment of greater indi-
viduality. A new beam arrangement showed to be preferential when balancing conformality
and robustness for this case.
Conclusions: The ebDD and error volume histogram proved effective in analysing plan
robustness. The process of retrospective analysis could be used to establish site specific
robustness planning protocols in proton planning. These protocols allow the planner to
determine plans that, although delivering a dosimetrically adequate dose distribution, have
resulted in sub-optimal robustness to these uncertainties. For these cases the use of different
beam start conditions may improve the plan robustness to set-up and range uncertainties.
127
A.0.1 Other Presentations
Poster presenation at Cambridge Science Festival. McGowan SE, Proton Therapy 2011
Poster presentation at the School of Clinical Medicine Research day. McGowan SE, Plan-
ning for Proton Therapy. 2012
Oral presenation at the 2014 PPRIG, NPL. McGowan SE, Albertini F, Thomas SJ, Bur-
net NG & Lomax AJ. The importance of plan robustness in proton beam therapy.
Proton Therapy Physics Workshop (NPL PPRIG) at the National Physical Labora-
tory, 12th – 13th March 2014.
Invited speaker at the Christie Hospital in Manchester McGowan SE, Treatment Plan-
ninig Uncertainites in Proton Therapy 2014
128 Published Work and Presentations
Fig. A.3 Winning popular press poster describing my research to the general public at the
School of Clinical Medicine Research day 2012
129
Fig. A.4 Copy of my certificate for presenting a popular press poster at the clinical school
Research day

Appendix B
Code
This appendix includes the main functions used to create the results in Chapter 4. In order
as presented here these include,
■ Startcons
■ First Fraction
■ Next Fraction
■ Read MVCT DICOM
■ Apply CT shifts
■ Read in structure DICOM set
■ Water equivalent path-length code
■ Results Script
A more indepth description of how principle functions work is provided below followed by
the codes. Some of these functions were created specific to either prostate or head and neck
patients. I have only included one of the versions, rather than both, as mostly the code is
identical except for directories.
Startcons
This function sets the starting conditions and is where the user can edit the parameters to
be used in the rest of the code. These parameters are split into families for programming
132 Code
ease; each family was constructed as a structure set. The four main families included Beam
Data, (Beam), Structure data, (Struct), DICOM data (Dicomstruct) and patient data, (Pat).
Through the coding process data is added and stored to these families to be called and used
later and in other functions. This way only the family structures need to be passed between
functions and not an entire lists of variables.
Beam variables include the beam angle, imaging modality and source to isocentre distance
(SID)∗.
Struct variables include the number of structures to be analysed, their names as a string, and
the number of letters in the string.
Dicomstruct variables include the pixel dimensions, number of slices in the data set and the
coordinate system.
Pat variables include the patient token, the directory names for retrieving and saving patient
data, the treatment site (i.e. prostate) and current fraction being investigated.
∗The SID has been set to seven metres in this study, this is the distance used in the proton
at CNAO, Pavia and is used simply as an example. At PSI at dose calculation the SID is
assumed infinite.
First fraction Next fraction
These are the two main scripts used to run functions. The first reads in the data required
for comparison, this is any fraction of the patient treatment that will be used to compare
the consecutive fractions to, or respective fractions depending on the investigation. This
script will first call the function to establish the start conditions. It then loops through the
folders in that patient directory to find the relevant DICOM information and Structure sets.
The next script loops through the remaining fractions and carries out all the same functions,
except the one to read in the structure set coordinates as these have already been saved in
the directory so they can just be loaded into Matlab to determine which slices on the shifted
MVCT are required for calculating the WEPL.
133
Read MVCT DICOM
This function goes to the pre-specified directory as determined in either firstfraction.m or
nextfraction.m to access the MVCT data. For some patients there may exist two sets of
MVCT with the same date label so time labels are also used when organising all the image
data. For these fractions, where two MVCT images exist, the patient was imaged prior
to and post treatment as part of a previous research study. The code allows for this and
will carry out full calculations for both imaging data sets while remaining within the same
fraction. This data will be used for investigating intra- fraction range errors by making
comparison between the prior and post treatment images.
As well reading the DICOM data and storing it as a 3D matirx, the code will also take and
store information from the DICOM header to store in Dicomstruct. This structure will hold
all the information relevant to that DICOM data needed for calculation. This information
includes, pixel size, slice thickness, number of pixels in each dimension, the coordinate
system, HU offset and the shifts that were applied to the couch in that instance. In the case
where shift data is missing, no shift will be applied.
It must be noted that DICOM data does not usually come with the shift information, this has
been retrieved directly from the TomoTherapy system to remove human error and added to
the header as part of the VoxTox study.
Apply CT shifts
For each treatment fraction on TomoTherapy the radiographer will set the patient up for
treatment and perform an MVCT scan. The radiographer will perform image registration of
the MVCT image to the planning kVCT making any of their own adjustments online and
determine the final shift required in each axis (x, yz∗, zy∗) and for the roll. These shifts
are translated to changes in the couch position and are applied. The patient will then be
treated in this position. By applying these shifts to the MVCT image we can simulate the
patient position at the time of treatment. In this study we have not included the roll shifts in
our calculation. It is also at this point that truncation is compensated for by reading in and
applying shifts to the kVCT data for prostae patients.
∗ Normal ’non-TomoTherapy’ axis – These are the axis used in script
The role of this function is to take all the information in the Dicomstruct family and use it
to apply a shift to the 3D matrix that is our patient’s image from that fraction. The shift is
134 Code
achieved by creating two coordinate systems of the same size of the matrix. The first is the
coordinate system retrieved from the DICOM header; the second is this coordinate system
with the shifts added to it in each direction. For example if say pixel (1,1,1) existed at the
coordinates (0,0,-12) and rigid coach shifts (0.5,1.5,-1) were applied, the new coordinate
system for pixel (1,1,1) would be (0.5,1.5,-13). This is done to all pixels, in pixel dimension
rather than in mm, and a 3D cubic interpolation is carried out to determine the new values
of pixel data. Instead of extrapolating the data in the z direction where it is missing, the
previous interpolated slices are repeated. In the x and y directions zeros are used to pad the
matrix out as only air existed there. These zeros are then dealt with later when all air is
assigned to the same value.
Read in structure DICOM set
A function will then read in the correct structure coordinates and store them. This function
loops through the relevant part of the structure set DICOM data to find the correct region
of interest as specified in the starting conditions. It will also match and align the structure
slices to the corresponding shifted MVCT image slices. The resolution of the MVCT and the
structure set differ in slice thickness so to determine which structure set slices are required
the inbuilt Matlab function ismember is used to pinpoint the elements of the two matrix
sets where the z coordinates are in correspondence. Only the image slices which have
structure data associated with it will be passed back to the main script. The coordinates of
the structure set will be saved to that patient directory to be used for the consecutive fractions
to save reading them in again (unless the structure set is one that has been re-drawn on each
MVCT set, it that case it is read in at each fraction). This function is only called once
when investigating the target, when investigating OAR changes for the prostate patients this
function is called at every fraction as the rectum has been redrawn on every fraction as part
of the VoxTox study which these patients are enrolled on.
It has been assumed that,
■ The target has been delineated correctly
■ the target is in the correct location
■ the target has not rotated
135
Water equivalent path-length code
In this function, the WEPL is calculated for a given beam direction, and beam parameters,
using the shifted MVCT data and the effective depth algorithm previously described. To
ensure the quality of the effective depth algorithm several verification tests have been carried
out; these included testing the WEPL is accurate by setting the source to the surface of the
CT image and the beam angle perpendicular to the pixels. This way the value of the pixel
can be checked to determine if it was as expected and that the next pixel had a value equal
to its own plus that of the previous pixel. An easy way to check the summation process is to
look at water and see that the pixel value increases by a distance of one pixel each time.
Following the calculation of the WEPL, masks are used to pull the WEPL data from pixels
that only exist inside the structure set for each fraction. Each subsequent fraction is then
compared to the first (or the one was used of the first script). All dimensions are converted
into distances and volumes and a histogram for each fraction depicting the percentage vol-
ume of the target that experiences a difference in WEPL in mm is created. A final histogram
is created representative of the entire treatment. Statitsacal parameters such as the mean and
standard deviation of each histogram is calculated, these represent the systematic and range
errors associated with the treatment, both as a whole and at each fraction for that given beam
angle.
136 Code
Start conditions
Contents
• Start conditions
• Beam data
• Struct data
• Pat data
function [ Beam, Struct, Pat ] = HNstart_cons(Beam, Pat)
%The start coonditions are set by the user here
% The beam structure contains beam information.
% The Struct struct contain info about the structure ROI.
% The Pat structure contains info about that patient.
Beam data
Beam.beam_angle=45; % define the beam angle, can define elsewhere
Beam.Type=0; % 1=orsay, 0 = tomo
Beam.source_isocentre=7000; %based on CNOA data
Beam.beam_rad=Beam.beam_angle*(pi/180);% converts to rads
Struct data
Struct.number=1; % to set how many structure sets to analyse
% Struct.volume(2).ROI=’rectum’; % define the ROI
% Struct.volume(2).let=8; %how many letters
Struct.volume(1).ROI=’ctv65’;%’CTV65’; % define the ROI
Struct.volume(1).let=3;
Pat data
Pat.token =102; % could be set else where
Pat.tok = num2str(Pat.token);
Mac=1; %set whether running on a mac or windows operating system
Windows=0;
CompU=Mac;
137
Pat.prostate=0;
Pat.headneck=1;
% creates directory paths to be used
if CompU>Windows;
Pat.direct_name = ...
([’/Users/staceymcgowan/Documents/PhD WEPL/HN/Data/’...
,Pat.tok]);
s= exist([Pat.direct_name ...
’/data_’,num2str(Beam.beam_angle)], ’dir’);
if s <1
mkdir([Pat.direct_name ’/data_’,num2str(Beam.beam_angle)]);
% creates directory if there is not one already
end
Pat.direct_nameD = ...
([’/Users/staceymcgowan/Documents/PhD WEPL/HN/Data/’...
,Pat.tok,’/data_’,num2str(Beam.beam_angle)]);
Pat.direct_nameResults = ...
([’/Users/staceymcgowan/Documents/HN/Data/’,...
Pat.tok,’/results’]);
else
Pat.direct_name = ([’R:\User Data\Stacey\Data\__Stacey\’...
,Pat.tok]);
Pat.direct_nameD = ...
([’R:\User Data\Stacey\Data\__Stacey\’,...
Pat.tok,’\data90’]);
Pat.direct_nameResults = ...
([’R:\User Data\Stacey\Data\__Stacey\’,...
Pat.tok,’\results’]);
end
Pat.Vox=1;
Pat.fract=1;
Pat.fraction =num2str(Pat.fract);
end of function
138 Code
Runs functions for the first fraction
Contents
• Runs functions for the first fraction
• read in CT data
• Apply shift to CT data
• Find slices of MVCT and structure set that match
• calculate WEPL
• load mask of structure
• pull wepl data from structure mask
• end of function
%Need to read in the Target structure from the kVCT structure set and
%determine its cordinates. Determine which z slices to use.
function [Beam, Struct, Pat] = firstfraction(Beam, Pat)
[Beam, Struct, Pat] = start_cons(Beam, Pat);
%[sort]=sort_DailyImagesFUN(Pat); %run to organise images
folders = dir(Pat.direct_name);
yy=0; %counter
i = 1;
Pat.fract=i;
Pat.fraction =num2str(Pat.fract);
Pat.foldname = strcat(num2str(i));
time_folders=dir([Pat.direct_name ’/’ Pat.foldname]);
for q = 1:length(time_folders)
if (strcmp(time_folders(q).name,...
’.’) || strcmp(time_folders(q).name, ’..’)...
|| strcmp(time_folders(q).name, ’MVCTstruct’)||...
strcmp(time_folders(q).name, ’kVCTstruct’)...
|| strcmp(time_folders(q).name, ’struct.dcm’)...
|| strcmp(time_folders(q).name,’.DS_Store’))
continue %skips files not required
end
139
read in CT data
Pat.fract_directory = ...
([Pat.direct_name ’/’ Pat.foldname...
’/’, time_folders(q).name]);
[ct_orig, Dicom_struct] = openimage_fun(Pat); %read CT
Pat.time=time_folders(q).name;
Apply shift to CT data
[shifted, Dicom_struct ]...
= applyZ_fun( Dicom_struct, Pat, ct_orig);
Find slices of MVCT and structure set that match
[Struct, b, d, ct_knownS]... %read structure set
= firstfractionStruct(Dicom_struct, Struct, Pat, yy);
%ct_knownS=ct_knownS-30;
disp(’shift applied’)
calculate WEPL
[ wepl ] = ...
Wepl_VOXfunfirstCAL( Dicom_struct, Beam, Pat, ct_knownS, yy);
load mask of structure
Target=load([Pat.direct_nameD ’/’,’1structmask.mat’]);
logic=logical(Target.thisMAT); % make mask
pull wepl data from structure mask
needed=wepl(logic);
%nnn = viewing(logic, ct_knownS, Dicom_struct)
if yy <1 % dtermines whether more than one MVCT per #
save([Pat.direct_nameD, ’/’,num2str(Pat.fract),...
’result1.mat’], ’needed’);
140 Code
yy=1;
else
save([Pat.direct_nameD, ’/’,num2str(Pat.fract),...
’result2.mat’], ’needed’);
end
end
end of function
141
Next Fraction
Contents
• Next Fraction
• Clean up slice matrix
• match up slices to Structure set
• end of function
Runs through conseuctive fractions and runs all same functions. Except reading the stuctuire
DICOm as just load the previous one. unless structure set re drawn on MVCT
function [Beam, Struct, Pat] = nextfraction(Beam, Pat)
[ Beam, Struct, Pat ] = start_cons(Beam, Pat);
%[sort]=sort_DailyImagesFUN(Pat);
folders = dir(Pat.direct_name);
no_fracts =(length(folders)-8);
for i =[2:36] %loop through fractions
yy=0;
Pat.fract=i
Pat.fraction =num2str(Pat.fract);
Pat.foldname = strcat(num2str(i));
time_folders=dir([Pat.direct_name ’/’ Pat.foldname]);
for q = 1:length(time_folders)
if (strcmp(time_folders(q).name, ’.’)...
|| strcmp(time_folders(q).name, ’..’)...
|| strcmp(time_folders(q).name, ’MVCTstruct’)||...
strcmp(time_folders(q).name, ’kVCTstruct’) ...
|| strcmp(time_folders(q).name,’.DS_Store’))
continue %skip uneeded folders
end
Pat.fract_directory = ...
([Pat.direct_name ’/’ Pat.foldname ’/’, time_folders(q).name]);
[ct_orig, Dicom_struct] = openimage_fun(Pat);
Pat.time=time_folders(q).name;
142 Code
[shifted, Dicom_struct ] = applyZ_fun( Dicom_struct, Pat, ct_orig);
disp(’shift applied’)
MVCTslicemat =...%creates matrix of MVCT slice coords
((Dicom_struct.Iz_cord:-Dicom_struct.slicethickness:...
Dicom_struct.Iz_cord-((Dicom_struct.no_files-1)...
*Dicom_struct.slicethickness)));
disp(num2str(Dicom_struct.shifts(2)))
Clean up slice matrix
for iuh= 1:Dicom_struct.no_files
if MVCTslicemat(1,iuh) <0
if MVCTslicemat(1,iuh) >(-1)
MVCTslicemat(1,iuh)=0;
end
end
end
for iuh= 1:Dicom_struct.no_files
if MVCTslicemat(1,iuh) >0
if MVCTslicemat(1,iuh) <(1)
MVCTslicemat(1,iuh)=0;
end
end
end
N= num2str(MVCTslicemat,4);
MVCTslicemat= str2num(N);
match up slices to Structure set
%need to find the coords in the MVCTstructslicemat that match the values
%in the struct mat.
%These values are the correct cordinates to be saved in the mask matrix.
load([Pat.direct_nameD ’/’,’structmat.mat’]);
load([Pat.direct_nameD, ’/1d.mat’]);
LIA=ismember(MVCTslicemat, structmat(:,oned)); % do slices match?
LIAB=ismember(structmat(:,oned), MVCTslicemat);
143
% need this to work to find the nearest values also
[c,bb]=ind2sub(size(LIA), find(LIA==1)); %where do they match?
[a,d]=ind2sub(size(LIAB), find(LIAB==1));
[w,r,numb]=size(bb);
b=fliplr(bb)
load([Pat.fract_directory ’/’,’shifted.mat’]);
%[sw,tw,yw]=size(shifted);
shifted(shifted<11)=0;
% if numb>=yw
% ct_knownS=shifted(:,:,b(1:yw));%check it is not d
% end
% if numb<yw
ct_knownS=shifted(:,:,b);
% end
%ct_knownS=ct_knownS-30;
[ wepl ] = Wepl_VOXfunfirstCAL( Dicom_struct, Beam, Pat, ct_knownS, yy);
imagesc(wepl(:,:,3));
Target=load([Pat.direct_nameD ’/’,’1structmask.mat’]);
logic=logical(Target.thisMAT);
[sw,tw,ww]=size(logic);
[sw,tw,dw]=size(wepl);
% nnn = viewing(logic, ct_knownS, Dicom_struct);
% if ww>=dw
% needed=wepl(logic(:,:,1:dw));
% end
% if ww<dw
% needed=wepl(logic(:,:,ww));
% end
needed=wepl(logic(:,:,d));
%nnn = viewing(logic, ct_knownS, Dicom_struct)%check visually
if yy <1
save([Pat.direct_nameD, ’/’,num2str(Pat.fract), ’result1.mat’], ’needed’);
yy=1;
else
144 Code
save([Pat.direct_nameD, ’/’,num2str(Pat.fract), ’result2.mat’], ’needed’);
yy=0;
end
end
end
end of function
145
Read MVCT DICOM
Contents
• Read MVCT DICOM
• read in data
• get DICOM header info
function [ct_orig, Dicom_struct] = openimage_fun(Pat)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
%directory_name = ...
%uigetdir(’/Users/staceymcgowan/dropbox/matlab/fun/’);
%directory_name = ...
%uigetdir(’C:/Documents and Settings/mcgowans/My Documents/MATLAB/’);
read in data
files = dir([Pat.fract_directory ’/*.dcm’]);
ct_orig = double(dicomread([Pat.fract_directory ’/’ files(1,1).name]));
[r,c] = size(ct_orig);
ct_orig = zeros(r,c,length(files));
Dicom_struct.no_files=length(files);
for i = 1:length(files)
filename = strcat(num2str(i),’.dcm’);
ct_orig(:,:,i) = dicomread([Pat.fract_directory ’/’ filename]);
%opens directories so user can choose ct data set, 2D or 3D,
end
get DICOM header info
infor=dicominfo([Pat.fract_directory ’/’ files(1,1).name]);
%reads dicom header
[m,n]=size(ct_orig(:,:,1));
Dicom_struct.m=m;
Dicom_struct.n=n;
146 Code
Dicom_struct.slicethickness=infor.SliceThickness;
Dicom_struct.ct_coord=infor.ImagePositionPatient(1);
Dicom_struct.pixelsize=infor.PixelSpacing(1);
Dicom_struct.huoffset=infor.RescaleIntercept; %gets the HU offset
Dicom_struct.Iz_cord=infor.ImagePositionPatient(3);
Dicom_struct.inst=infor.InstanceNumber;
%if Pat.prostate>0
Dicom_struct.shifts=infor.Private_0099_1011;
% else
% num= load(’shifts11.mat’);
% xshift=num.shiftdata(Pat.fract,1);
% yshift=num.shiftdata(Pat.fract,2);
% zshift=num.shiftdata(Pat.fract,3);
%Dicom_struct.shifts=[xshift, yshift, zshift];
if ischar(Dicom_struct.shifts)>0
Dicom_struct.shifts=[0,0,0];
end
end
End of function
147
Apply CT shifts
Contents
• The shifts
• Apply shifts
• End of function
function [ shifted, Dicom_struct ] = applyZ_fun( Dicom_struct, Pat, ct_orig)
if Pat.Vox>0
zshift=Dicom_struct.shifts(2);
xshift=Dicom_struct.shifts(1);
yshift=Dicom_struct.shifts(3);
end
if Pat.Vox<1
if Pat.headneck>0;
% [num]= xlsread(’C:/Documents and Settings/mcgowans/My Documents...
%/MATLAB/proton range study/code/HN_shifts.xls’, Pat.pt);
[num]= ...
xlsread(’/Users/staceymcgowan/dropbox/matlab/fun/HN_shifts.xls’...
,Pat.pt);
end
if Pat.prostate>0
[num]= xlsread...
(’C:/Documents and Settings/mcgowans/My Documents...
/MATLAB/proton range study/code/patient_shifts.xls’...
, Pat.pt);
%[num]= ...
xlsread(’/Users/staceymcgowan/dropbox/matlab/fun/patient_shifts.xls’...
, Pat.pt);
end
148 Code
The shifts
zshift=num(3,Pat.fract);
xshift=num(1,Pat.fract);
yshift=num(2,Pat.fract);
end
%convert from mm to pixel dimensions
dunSliceX=(xshift/Dicom_struct.pixelsize);
dunSliceY=(yshift/Dicom_struct.pixelsize);
dunSliceZ=(zshift/Dicom_struct.slicethickness);
%ct_cordz=Dicom_struct(Pat.fract).Iz_cord:-Dicom_struct.slicethickness:...
%(Dicom_struct(Pat.fract).Iz_cord-((Dicom_struct.no_files-1)...
*Dicom_struct.slicethickness));
Apply shifts
ct_cordz=1:Dicom_struct.no_files;
ct_cordx=1:Dicom_struct.n;
ct_cordy=ct_cordx;
[X, Y, Z]=meshgrid(ct_cordx, ct_cordy, ct_cordz);
%shift applied by cubic interpolation.
% 3D interp require 3D coord systems
% use meshgrid to create exsiting and new systems
Shift_cordz=ct_cordz+dunSliceZ;
Dicom_struct.zshiftmat=Shift_cordz;
Shift_cordx=ct_cordx-dunSliceX;
Shift_cordy=ct_cordy+dunSliceY;
[Xi, Yi, Zi]=meshgrid(Shift_cordx ,Shift_cordy, Shift_cordz);
% 3D interp of data in old to new coord system
% new coord system is old one plus the shifts
149
shifted=interp3(X, Y, Z, ct_orig, Xi, Yi, Zi, ’cubic’);
% no extrapolating use, so buffer with zeroes in x and y
% assume no pt trunction
% For missing data in z, copy last exsiting slice
%So...
if zshift == 0
else
if zshift <0
first =0;
%find the cordinate of the first value >1
if first < 1;
for I=1:Dicom_struct.no_files
if Shift_cordz(I) >1
first=I;
break
end
end
end
%jt=abs(floor(Shift_cordz(1)));
j=first-1;
w=ones(1,j).*(j+1) ;
shifted(:,:,1:j)=shifted(:,:,w);
end
if zshift >0
jt=abs(ceil(Shift_cordz(Dicom_struct.no_files)));
j=abs(Dicom_struct.no_files-jt);
q=Dicom_struct.no_files-(j-1):Dicom_struct.no_files;
w=ones(1,j).*Dicom_struct.no_files-j;
shifted(:,:,q)=shifted(:,:,w);
end
end
150 Code
shifted(isnan(shifted)) = 0; %get rid of NaNs
\subsection{If prostate read in kV and apply shifts and make masks}
if Beam.beam_angle>0
MVCT=shifted;
[Dicom_struct, Pat, kVCT]=kvctreadin(Pat, Dicom_struct);
[Dicom_struct, Pat, kvCT_shifted]=kvctshift(Pat, Dicom_struct,kVCT);
width = Dicom_struct.kvm;
height = Dicom_struct.kvn;
radius = 193/Dicom_struct.kvpixelsize;
centerW = width/2;
centerH = height/2;
[W,H] = meshgrid(1:width,1:height);
mask = sqrt((W-centerW).^2 + (H-centerH).^2) < radius;
maskop=mask.*10;
maskop(maskop<10)=1;
maskop(maskop>1)=0;
mask_mat=ones(Dicom_struct.kvn,Dicom_struct.kvm,Dicom_struct.kvno_files);
for j=1:Dicom_struct.kvno_files
mask_mat(:,:,j)=maskop;
end
kv_mask=mask_mat.*kvCT_shifted;
kv_mask(kv_mask(:,:,:)<100)=0;
kv_mask(kv_mask(:,:,:)>0)=924;
Dicom_struct.kv_mask=kv_mask;
%if Beam.beam_angle>0
[Dicom_struct, Pat, MVCT] = painted(Dicom_struct, Pat, MVCT);
%end
one=Dicom_struct.need*2;
x=one+1:Dicom_struct.total;
y=Dicom_struct.need:(Dicom_struct.total-Dicom_struct.need-1);
151
shifted=MVCT(y,x,:);
end
clear ct_orig
% s= exist([Pat.direct_name ’\’ Pat.tok ’\’ Pat.fraction], ’dir’);
% if s <1
% mkdir(Pat.direct_name, Pat.tok);
% mkdir([Pat.direct_name ’/’ Pat.tok], Pat.fraction);
% end
if Pat.prostate>0
save([Pat.fract_directory ’/’ ,’shifted.mat’], ’shifted’);
else
j=Pat.foldname;
save([Pat.direct_name ’/’ j ’/’ ,’shifted.mat’], ’shifted’);
end
end
End of function
152 Code
Read in structure DICOM set
Contents
• Read in structure DICOM set
• Get coordinates out of header
• Create z cordinate matix of MVCT & struture set
• Storing the slices required
• end of function
%function loops through structure DICOM header
% pulls out coordinates for pre defined ROI
% stores this data in 3D matrix
% the slices in this martrix that match the MVCT slices are correlated
Get coordinates out of header
article graphicx color
function [Struct, b, d, ct_knownS, test]...
= HNfirstfractionStruct(Dicom_struct, Struct, Pat)
metadata = dicominfo([Pat.direct_name ’/kVCTstruct.dcm’]);
%struct if HN
for j = 1:length(struct2cell(metadata.ROIContourSequence))
% loop through number of structures
%j=48; %Loop modified from Nick Early’s Loop
z = eval([’metadata.StructureSetROISequence.Item_’ num2str(j) ’.ROIName’]);
if strncmpi(z,Struct.volume(1).ROI,Struct.volume(1).let);
% if j==2;
disp(z);
disp(j);
g = eval([’length(struct2cell(metadata.ROIContourSequence.Item_’...
num2str(j) ’.ContourSequence))’]); % finds number of slices
test.Storage(1).stores=zeros(Dicom_struct.m, Dicom_struct.n, g);
%test.Storage(1).DEstores=zeros(Dicom_struct.m, ...
Dicom_struct.n, Dicom_struct.no_files);
%stores=zeros(100,3,g);
153
%storey=zeros(g,1);
for k = 1:g % loop through all these slices
%k=27; %slice of interest number
s = [’metadata.ROIContourSequence.Item_’ num2str(j)...
’.ContourSequence.Item_’ num2str(k) ’.ContourData’];%get cords
CD = eval(s);
CD = [CD; CD(1:3)]; %makes co-ords end where they start
st = -1 + length(CD)/3;%co-ords are all jumbled up (xyzxyz...)
d = []; %this code splits them into 3 arrays d,e,f
e = [];
f = [];
for i = 0:st
d(i+1) = CD(1+(3*i),1);
Struct.cord(Pat.fract).volume(1).xpos(i+1)=d(i+1);
e(i+1) = CD(2+(3*i),1);
Struct.cord(Pat.fract).volume(1).ypos(i+1)=e(i+1);
f(i+1) = CD(3+(3*i),1);
Struct.cord(Pat.fract).volume(1).zpos(i+1)=f(i+1);
end
cord1=d’;
cord2=e’;
cord3=f’;
if k==1
Pat.first=cord3(1);
% Pat.first=3;
end
if k==2
Pat.structslice= abs(Pat.first-cord3(1));
%Pat.structslice=1;
end
A=[cord1, cord2, cord3];
pos_x = (cord1-Dicom_struct.ct_coord+Dicom_struct.pixelsize)...
/Dicom_struct.pixelsize;
pos_y = (cord2-Dicom_struct.ct_coord+Dicom_struct.pixelsize)...
154 Code
/Dicom_struct.pixelsize;
%pos_y = 512-pos_ys;
pos_z=k;
%((cord3+Dicom_struct.Iz_cord+Pat.structslice)/Pat.structslice)
roi_mask = poly2mask(pos_x, pos_y, 512, 512);
test.Storage(1).stores(:,:,k)=roi_mask;
end
end
end
\subsection*{Storing the slices required}
MVCTslicemat =(((Dicom_struct.Iz_cord:-Dicom_struct.slicethickness...
:Dicom_struct.Iz_cord-((Dicom_struct.no_files-1)*...
Dicom_struct.slicethickness)))-Pat.eee);%took-1 out
[hjk,mvctslice]=size(MVCTslicemat);
%Dicom_struct.no_files=mvctslice;
for iuh= 1:Dicom_struct.no_files
if MVCTslicemat(1,iuh) <0
if MVCTslicemat(1,iuh) >(-1)
MVCTslicemat(1,iuh)=0;
end
end
end
for iuh= 1:Dicom_struct.no_files
if MVCTslicemat(1,iuh) >0
if MVCTslicemat(1,iuh) <(1)
MVCTslicemat(1,iuh)=0;
end
end
end%+Dicom_struct.shifts(1,3)
N= num2str(MVCTslicemat,4);
MVCTslicemat= str2num(N);
structmat=(Pat.first:Pat.structslice:Pat.first+(Pat.structslice*g-1));
155
%MVCTslicemat(MVCTslicemat<Pat.end)=NaN; % use to avoid slice (shoulders)
%structmat(structmat<Pat.end)=NaN;
LIA=ismember(MVCTslicemat, structmat); %finds matches
LIAB=ismember(structmat, MVCTslicemat);
[c,bb]=ind2sub(size(LIA), find(LIA==1));
[a,dd]=ind2sub(size(LIAB), find(LIAB==1));
[r,numb]=size(bb);
%d=fliplr(dd)
%b=fliplr(bb);
d=dd;
thisMAT=zeros(Dicom_struct.m, Dicom_struct.n, numb);
thisMAT=test.Storage(1).stores(:,:,d);
%disp(d);
oned=d%+facto;
save([Pat.direct_nameD, ’/’,num2str(Pat.fract), ’d.mat’], ’oned’);
j=Pat.fraction;
load([Pat.direct_name ’/’ j ,’/shifted.mat’]);
shifted(shifted<11)=0;
b=fliplr(bb);
%b=bb;
ct_knownS=shifted(:,:,b);%check it is not d
save([Pat.direct_nameD ’/’ , ’structmat.mat’], ’structmat’);
save([Pat.direct_nameD ’/’ ,num2str(1),’structmask.mat’], ’thisMAT’);
nnn = viewing(thisMAT, ct_knownS, Dicom_struct)
%clear test;
clear roi_mask
end
156 Code
This LaTeX was auto-generated from an M-file by MATLAB.
Water equivalent path-length code
Contents
• Water equivalent path-length code
• Set up beam conditions
• WEPL loop
• end of function
%Wepl_VOXfun is a function that requires beam, Dicom and patient data to
%then determine the wepl path legnth of a patient Data set. The
%calibration details for converting HU to relative proton stopping
%power for MVCT are included, kVCT could also be added with new cal.
%This code can be used for photons by adding a new calinration curve.
Set up beam conditions
function [ wepl ] = Wepl_VOXfunfirstCAL...
( Dicom_struct, Beam, Pat, ct_knownS,yy)
%Modified from Dr Simon Thomas’ code
xiso=(Dicom_struct.n/2+1);
ziso=(Dicom_struct.m/2+1);
sid=(Beam.source_isocentre/(Dicom_struct.pixelsize));
%sets up centre of data and isocentre
xf=xiso +(sid*sin(Beam.beam_rad));
%finds x coord of focus (centre of rotation is at centre of CT)
zf=ziso-(sid*cos(Beam.beam_rad));
%finds z coord of focus
xxp=1:Dicom_struct.n;
xp=repmat(xxp, 1, Dicom_struct.n);
%each element is the x coord of the mxn matrix and is in 1 row
nB=1:Dicom_struct.m;
Brep=repmat(nB,Dicom_struct.n,1);
zp=reshape(Brep,1,((Dicom_struct.m)*(Dicom_struct.n)));
157
%each element is the z coord of the mxn matrix and is in 1 row
dx=xf-xp;
%distance from focus to point of interest in x
dz=zf-zp;
%distance from focus to point of interest in y
steps=max((max(abs(dx),abs(dz))));
%number of intervals for later
rl=sqrt(dx.*dx + dz.*dz);
%distnce from point of interest to source
dr=rl/steps;
deltax=dx/steps;
deltaz=dz/steps;
x=xp;
z=zp;
%so can keep original values of xp and zp (zp,xp is coords)
rpath=0;
clear Brep
[a,b,no_slices]=size(ct_knownS);
v=ones((Dicom_struct.m+2),1)*(1);
t=ones(1,Dicom_struct.n)*(1);
C=zeros(Dicom_struct.n+2, Dicom_struct.m, no_slices);
G=zeros(Dicom_struct.n+2, Dicom_struct.m+2, no_slices);
for tt = 1:no_slices
C(:,:,tt)=[t;ct_knownS(:,:,tt);t];
G(:,:,tt)=[v C(:,:,tt) v];
end
ct=permute(G, [3 1 2]);
%reorders pt dtata in the form of slicesx261244
Gdim=max(size(ct));
Gsize=Gdim*ones(1,2);
WEPL loop
%Loop structure from Dr Simon Thomas
for i=1:steps
158 Code
%loop to calculate wtaer equivlent path to the point of interest
iz=round(z);
iz(iz<1)=1; iz(iz>Gdim)=1;
ix=round(x);
ix(ix<1)=1; ix(ix>Gdim)=1;
use 1 as flag both sizes
if max(iz)==1, break, end
if max(ix)==1, break, end
hu=double(ct(:,sub2ind(Gsize,iz,ix)))+Dicom_struct.huoffset;
%gets hu value from ct
%hu=ctimages(:,sub2ind(ctsize,iz,ix));
hu(hu<=(-900))=-1000;
%for air so it does contributes to the water equivelent path
fhu=double(hu);
if Beam.Type<1 %new cal
rhop=0.9652+(0.0009914*fhu);
rhop(rhop<0)=0; %trap negative densities
end
end
if Beam.Type>0 %orsay cal
rhop=1.0153+(fhu/1000); % A
rhop2=1.0152+(fhu.^3*4e-08)+(fhu.^2*6e-06)+(0.0006*fhu);
rhop3=1.0356+(fhu/1111.11);
rhop(rhop<0.812)=0;
rhop(rhop>0.812)=rhop2(rhop>0.812);
rhop(rhop>1.065)=rhop3(rhop>1.065);
end
rpath=rpath+rhop;
%sums the previous relative electron densities and
%will get multiplied by distance dr later
x=x+deltax;
z=z+deltaz;
end
a=ones(no_slices,1);
dr2=a*dr; %
159
rpath=rpath.*dr2;
%water equivelent path length for all slices
wepl=permute(rpath, [2 1]);
wepl=reshape(wepl,Dicom_struct.n, Dicom_struct.m, no_slices);
wepl=permute(wepl, [2 1 3]);
%puts matirx into conventional order of (XxY)xZ
if yy <1
save([Pat.direct_nameD ’/’, num2str(Pat.fract),’wepl1’], ’wepl’);
yy=1;
else
save([Pat.direct_nameD ’/’, num2str(Pat.fract),’wepl2’], ’wepl’);
end
end of function
Results code Fraction Data
Contents
• Stats for fraction bit!
• Plot bit fraction bit
• Cumlative
• Entire Treatment bit
• codes
function [Beam] = ResultsOut( Beam)
count=0;
statsdatastore=zeros(5, 12);
for beam_angle=[0,90];
count=count+1;
Beam.beam_angle=beam_angle;
Pat.token=21;
bnbn=36;
Dicom_struct.pixelsize=0.7647;
Dicom_struct.slicethickness=6;
160 Code
Dicom_struct.no_files=bnbn;
storingHIST=zeros(Dicom_struct.no_files+2, 58);
storingCUML=zeros(Dicom_struct.no_files+2, 58);
[ Beam, Struct, Pat ] = start_cons(Beam, Pat);
% if Beam.beam_angle > 0
% sheet=Pat.token;
% else
% sheet=Pat.token+1;
% end
%folders = dir(Pat.direct_nameD);
no_fracts = bnbn;
getfirst=load([Pat.direct_nameD ’/’ ,num2str(1),’result1.mat’]);
[resSize,tyh]=size(getfirst.needed);
gotn=zeros(resSize,1);
goth=zeros(resSize,1);
got_again=zeros(resSize,no_fracts-1);
counter=0;
hgf=[17 36 15 18 11 19 3 10 27 35 12 32];
for i = [17 36 15 18 11 19 3 10 27 35 12 32]
%2:3,7,10:12,16,17,23,31,33,34 to
%1?5,6,8,9,15,18:20,24:30,36 to4? -pt 78
counter=counter+1;
getnext=load([Pat.direct_nameD ’/’ num2str(i),’result1.mat’]);
neededfirst=getfirst.needed.*Dicom_struct.pixelsize;
needednext=getnext.needed.*Dicom_struct.pixelsize;
got=(needednext-neededfirst);
gotn=gotn+got;
Stats for fraction bit!
medata=mean(got);
got_again(:,i-1)=got;
mi=min(got);
mx=max(got);
161
st=std(got);
varian=var(got);
storingHIST(i,52)=medata;
%storing(i,103)=mod;
storingHIST(i,53)=varian;
storingHIST(i,54)=st;
storingHIST(i,55)=mi;
storingHIST(i,56)=mx;
%v=linspace(mi,mx,100);
Plot bit fraction bit
v=-50:2:50;
[numbe,a]=size(got);
%h_vol=(numbe*(Dicom_struct.pixelsize*Dicom_struct.pixelsize...
%*Dicom_struct.slicethickness));
whole_vol=(numbe*(Dicom_struct.pixelsize*Dicom_struct.pixelsize...
*Dicom_struct.slicethickness))/1000;%cm cubed
h_vol=(Dicom_struct.pixelsize*Dicom_struct.pixelsize*...
Dicom_struct.slicethickness);%mm cubed
totalvol=h_vol*numbe;%mm cubed
totalvoll=repmat(totalvol,1,51);
h=got;
hh=hist(h, v);
normed=((hh.*h_vol)./totalvoll).*100;
storingHIST(1,1:51)=v;
storingHIST(i,1:51)=normed;
handle(i)=figure(i);
plot(v, normed);
xlabel(’Change in WEPL (mm)’);
ylabel(’%Volume’);
%title([’Range Volume Histogram: Interfraction WEPL changes in fraction’,...
% num2str(i), ’ - pt ’,num2str(Pat.token), ’ prostate ’, ...
%num2str(Beam.beam_angle), ’ coplanar beam’]);
saveas(handle(i),[Pat.direct_nameD, ’/’ , num2str(Beam.beam_angle)...
162 Code
,Pat.tok,’_’, num2str(i),’_’, ’Hist.png’], ’png’);
Cumulative
totalvoll=repmat(totalvol,1,51);
cuh=abs(got);
goth=goth+cuh;
cvg=0:0.5:25;
cuhh=(hist(cuh, cvg));
cuhc=cumsum(cuhh);
cumnormedd=(((cuhc.*h_vol)./totalvoll).*100);
cumnormed=(cumnormedd*(-1))+100;
%cumnormed(cumnormed>100)=100;
%cumnormed(cumnormed<0.5)=0;
handle(i+40)=figure(i+40);
plot(cvg, cumnormed);
xlabel(’Change in WEPL (mm)’);
ylabel(’%Volume’);
%title([’Cumulative Range Volume Histogram: Interfraction WEPL ...
%changes in fraction’,...
% num2str(i), ’ - pt ’,num2str(Pat.token), ’ prostate ’...
%, num2str(Beam.beam_angle), ’ coplanar beam’]);
saveas(handle(i+40),[Pat.direct_nameD, ’/’ , num2str(Beam.beam_angle...
),Pat.tok,’_’, num2str(i),’_’, ’Cuml.fig’], ’png’);
storingCUML(1,1:51)=cvg;
storingCUML(i, 1:51)=cumnormed;
ggett=[2,3,5,10,20]; % vol getting range mm
pervols=interp1(cvg, cumnormed, ggett);
storingCUML(i,52:56)=pervols;
end
163
Entire Treatment bit
v=-50:2:50;
disp(counter)
gotT=gotn./(counter-1);
%gotT=sum(got_again,2);
[numbe,a]=size(gotT);
%h_vol=(numbe*(Dicom_struct.pixelsize*Dicom_struct.pixelsize*...
%Dicom_struct.slicethickness));
whole_vol=(numbe*(Dicom_struct.pixelsize*Dicom_struct.pixelsize...
*Dicom_struct.slicethickness))/1000;%cm cubed
h_vol=(Dicom_struct.pixelsize*Dicom_struct.pixelsize...
*Dicom_struct.slicethickness);%mm cubed
totalvol=h_vol*numbe;%mm cubed
totalvoll=repmat(totalvol,1,51);
h=gotT;
hh=hist(h, v);
normed=((hh.*h_vol)./totalvoll).*100;
% storingHIST((i+1),:)=sum(storingHIST(2:no_fracts+1,:),1);
% storingHIST((i+2),:)=storingHIST(i+1,:)./no_fracts;
handle(31)= figure(31);
plot(v, normed);
% stdevT=std(gotn/36)
% VarT=var(gotn/36)
xlabel(’Change in WEPL (mm)’, ’fontsize’, 12);
ylabel(’%Volume’, ’fontsize’, 12);
%title([’Range Volume Histogram: Interfraction WEPL changes in pt ’...
%,num2str(Pat.token), ’ sacrum, ’, num2str(Beam.beam_angle)...
%, ’degree beam’], ’fontsize’, 12);
saveas(handle(31),[Pat.direct_nameD, ’/’ , num2str(Beam.beam_angle)...
,Pat.tok,’_’, ’THist.png’], ’png’);
medata=sum(gotn)./(36*resSize);
handle(32)=figure(32);
164 Code
scatter(hgf, storingHIST(hgf, 54));
lsline
xlabel(’fraction’, ’fontsize’, 12);
ylabel(’stdev’, ’fontsize’, 12);
%title([’Pt ’, num2str(Pat.token), ’,’ num2str(Beam.beam_angle),
%’degree beam - stdev of change in range in CTV
%relative to first fraction’], ’fontsize’, 12);
saveas(handle(32),[Pat.direct_nameD, ’/’ ,...
num2str(Beam.beam_angle),Pat.tok,’_’, ’StDev.png’], ’png’);
%mod=mode(gotn)
mi=min(storingHIST(2:no_fracts,55));
mxsys=max(storingHIST(2:no_fracts,52));
misys=min(storingHIST(2:no_fracts,52));
mx=max(storingHIST(2:no_fracts,54));
meanT=mean(gotT);
stdev_thistimeB=sqrt(sum(storingHIST(i,53)));
stdec_thistimeC=std(gotT);
stdev_thistimeD=mean(storingHIST(i,54));
storingHIST(i+2,52)=meanT;
storingHIST(i+2,53)=sum(storingHIST(i,53));
storingHIST(i+2,54)=stdev_thistimeB;
storingHIST(i+2,55)=misys;
storingHIST(i+2,56)=mxsys;
storingHIST(i+2,57)=stdec_thistimeC;
storingHIST(i+2,58)=(sum(storingHIST(i,53)))/(sqrt(counter-1));
cuh=abs(gotT);
cuhh=(hist(cuh, cvg));
cuhc=cumsum(cuhh);
cumnormedd=(((cuhc.*h_vol)./totalvoll).*100);
cumnormed=(cumnormedd*(-1))+100;
handle(33)=figure(33);
plot(cvg, cumnormed);
xlabel(’Change in WEPL (mm)’, ’fontsize’, 12);
165
ylabel(’%Volume’, ’fontsize’, 12);
%title([’Cumulative Range Volume Histogram: Interfraction
%WEPL changes in pt ’,num2str(Pat.token), ’ Sacrum, ’,
%num2str(Beam.beam_angle), ’degree beam’], ’fontsize’, 12);
saveas(handle(33),[Pat.direct_nameD, ’/’ , ...
num2str(Beam.beam_angle),Pat.tok,’_’, ’TCuml.png’], ’png’);
pervolsT=interp1(cvg, cumnormed, ggett);
storingCUML(i+2,52:58)=[pervolsT,0,0];
%storingHIST(1,52:58)=[’mean’,’var’,’std’,’min’,’max’,’C’,’SE’];
% figure(i+41)
% plot(cv, storingCUML(i+2,1:51));
%storingHIST(1, 52:56)=([’Mean’, ’std’, ’Max’, ’Min’, ’hmm’]);
%storingCUML(1, 52:56)=([’V3%’, ’V5%’, ’V7%’, ’V10%’, ’V50%’]);
storingone=[storingCUML];
storing=storingHIST;
statsdata=[storingHIST(i+2,52:58) , pervolsT];
statsdatastore(count, :)=statsdata;
%save([Pat.direct_nameD, ’/’ , num2str(Beam.beam_angle), ’stats.mat’],
%’statsdata’);
% for bb=[2:32];
% savefig([Pat.direct_nameD, ’/’ , num2str(Beam.beam_angle),
%Pat.tok, num2str(bb), handle(bb), ’Hist.fig’]);
% end
% for bb=[41:72];
% savefig([Pat.direct_nameD, ’/’ , num2str(Beam.beam_angle),
%Pat.tok, num2str(i), handle(bb), ’cuml.fig’]);
% end
% savefig([Pat.direct_nameD, ’/’ , num2str(Beam.beam_angle),Pat.tok,
%num2str(i), handle(73), ’stdev.fig’]);
%savefig([Pat.direct_nameD, ’/’
%, num2str(Beam.beam_angle), handle(2), ’figcuml.fig’]);
titleh = [’hist’,num2str(Pat.token),’_’,num2str(Beam.beam_angle)];
%titlel = [’HN2cuml’,num2str(Pat.token),’_’,num2str(Beam.beam_angle)];
166 Code
dlmwrite(titleh, storing, ’delimiter’, ’\t’);
xlswrite([titleh,’.xlsx’], storing);
%xlswrite([titlel,’.xlsx’], storingone);
%close(handle)
%save([Pat.direct_nameResults ’/’,Pat.tok, ’
%NcalHist_data’,num2str(Beam.beam_angle)], ’storing’);
end of function
article graphicx color
function [Beam] = ResultsOut_all( Beam)
got_all=zeros(7,2);
anothercounter=0;
for beam_angle=[0, 25,70,110,180];
anothercounter=anothercounter+1;
Beam.beam_angle=beam_angle;
counter=0;
for token= [35,51,52,78,91,92,102]
counter=counter+1;
Dicom_struct.pixelsize=0.7647;
Pat.token=token;
[ Beam, Struct, Pat ] = HNstart_cons(Beam, Pat);
co=1;
fractions=[2:30];
if token==51
fractions=[3:6,8,10,15:20,22:25,27:31];
co=2;
end
if token==102
co=2;
fractions=[3:5,7:9,11,12,18:21,23:25,27:35,37,39:40];
end
getfirst=load([Pat.direct_nameD ’/’ ,num2str(co),’result1.mat’]);
[resSize,tyh]=size(getfirst.needed);
167
got_again=zeros(resSize,37);
if token==52
fractions=[4:18,20:34];
end
if token==91
fractions=[2,4:6,8:30];
end
for i = fractions
getnext=load([Pat.direct_nameD ’/’ num2str(i),’result1.mat’]);
neededfirst=getfirst.needed.*Dicom_struct.pixelsize;
needednext=getnext.needed.*Dicom_struct.pixelsize;
got=(needednext-neededfirst);
got_again(:,i-1)=got;
end
Store.token(counter).squareit=sum(sum(got_again.^2));
Store.token(counter).summer=sum(sum(got_again));
Store.token(counter).numall=length(fractions)*resSize;
end
Tl=[Store.token(1).squareit,Store.token(2).squareit,...
Store.token(3).squareit,Store.token(4).squareit,...
Store.token(5).squareit,Store.token(6).squareit,...
Store.token(7).squareit];
Sl=[Store.token(1).summer,Store.token(2).summer,...
Store.token(3).summer,Store.token(4).summer,...
Store.token(5).summer,Store.token(6).summer,...
Store.token(7).summer];
Nl=[Store.token(1).numall,Store.token(2).numall,...
Store.token(3).numall,Store.token(4).numall,...
168 Code
Store.token(5).numall,Store.token(6).numall,...
Store.token(7).numall];
N=sum(Nl);
S=sum(Sl);
T=sum(Tl);
mean = S/N;
sd=sqrt((T/N)-((S/N)^2));
got_all(anothercounter,:)=[mean, sd];
end
end
end of function