Modeling geometrical uncertainties for radiotherapy plan optimization without margins

(1)

Modeling

geometrical uncertainties

for radiotherapy

plan optimization

without margins

Eka Budiarto

(2)

MODELING GEOMETRICAL UNCERTAINTIES

FOR RADIOTHERAPY PLAN OPTIMIZATION

WITHOUT MARGINS

(3)

(4)

Modeling geometrical uncertainties for

radiotherapy plan optimization

without margins

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus Prof. ir. K. C. A. M. Luyben, voorzitter van het College voor Promoties,

in het openbaar te verdedigen

op dinsdag 6 januari 2015 om 15:00 uur

door

Eka BUDIARTO

Master of Science in Industrial Mathematics van Technische Universit¨at Kaiserslautern, Duitsland

(5)

Prof. dr. ir. A. W. Heemink Prof. dr. B. J. M. Heijmen Copromotor: Dr. ir. M. Keijzer

Samenstelling promotiecommissie: Rector Magnificus, voorzitter

Prof. dr. ir. A. W. Heemink, Technische Universiteit Delft, promotor Prof. dr. B. J. M. Heijmen, Erasmus Universiteit Rotterdam, promotor Dr. ir. M. Keijzer, Technische Universiteit Delft, copromotor Prof. dr. ir. D. den Hertog, Tilburg University

Prof. dr. M. B. van Herk Nederlands Kanker Instituut

Prof. Dr. rer. nat. U. Oelfke The Institute of Cancer Research, Engeland Prof. dr. ir. G. Jongbloed, Technische Universiteit Delft

Dr. ir. P. R. M. Storchi heeft als begeleider in belangrijke mate aan de totstandko-ming van het proefschrift bijgedragen. (art. 16.4)

Dit onderzoek kwam tot stand met steun van Erasmus Medical Center Cancer Insti-tute in Rotterdam

Eka Budiarto,

Modeling geometrical uncertainties for radiotherapy plan optimization without margins, Ph.D. Thesis, Delft University of Technology,

with summary in Dutch.

Keywords: radiotherapy, geometrical uncertainties, population-based model, princi-pal component analysis, treatment-plan optimization

ISBN: 9789461087898 Copyright:

c

2014 by Eka Budiarto c

2011, 2014 by Institute of Physics and Engineering in Medicine (chapters 3 and 4) All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted in any form or by any means without the written permission of the copyright owner.

Typeset by the author with LA_{TEX 2εdocumentation system}

Cover design by Bo Logiantara

(6)

For my family thank you for your continuous support

(7)

(8)

Summary

Radiotherapy is one of the methods used to treat cancer. One common ap-proach for radiotherapy is exposing the patient to external beams of high-energy X-ray photons. Using the intensity-modulated radiation therapy (IMRT) technique, the fluence (radiation energy per unit area) of the radiation beams can be modulated to optimally shape the 3D high-dose volume to the tumor shape. The term radiotherapy will henceforward refer to this technique.

Before radiotherapy can be administered, a treatment plan must be gen-erated to accomplish high tumor-dose and maximum sparing of surrounding healthy tissues. This planning is usually based on a computer-tomography (CT) scan showing the tumor and healthy organs-of-interest, called the planning CT scan. The treatment itself is normally given in several fractions over several days (e.g. about 40 for prostate cancer), where for each fraction a proportional part of the dose is given.

Due to movements and deformations, the geometrical positions of the organs in the body at the time of the treatment sessions might not be the same as in the planning CT scan. These geometrical uncertainties can make radiation beams miss the target (the tumor site, usually called the clinical target volume, CTV), making the actually delivered dose in the CTV less than what is prescribed. So, eradication of the tumor may fail. To prevent this, the irradiated area is usually enlarged with a margin, following a margin recipe. Although using a margin recipe can prevent underdosage in the CTV, unfortunately it can also potentially harm the healthy organs-at-risk (OAR) around the CTV. The larger the margin, the higher the probability of damaging the OAR.

The aim of this thesis is to develop a method to calculate the optimized flu-ence profiles of the radiation beams, while taking into account the geometrical uncertainties. The irradiated area will also be enlarged as with the margin-recipe, but now more locally adapted to the movements and deformations of the organs, sparing the OAR as much as possible, while delivering a high enough dose to the CTV. In this thesis, the movements and deformations are first mod-eled for a whole patient population. Then a useful method is described with which in this uncertain geometry the expected radiation dose and its corre-sponding variance can be calculated. Finally the correcorre-sponding radiation plan is optimized.

(9)

The first step is to model the movements and deformations of the organs. For this, a method called principal component analysis (PCA) has been used, which extracts the dominant modes of movements of the voxels (the volume discretization units) of the organs, as well as estimating their probabilities. Unfortunately PCA can usually not be done directly to a new patient because there is generally just one CT scan (or at most a few) available for this patient, which gives not enough data on movements and deformations. Therefore in this thesis a scheme for combining the movement data for the organs of interest from a population of patients has been developed, and tested for a prostate cancer site. PCA has then been applied to the resulting database, and it turns out that the dominant modes have plausible physical explanations. Further verification has been obtained by carrying out an error analysis, which shows that the dominant modes are also shared by prostates which are not included in the database.

The movements and deformations have been used to calculate the expected value (mean) and the variance of the dose. To do this, an integration over the probability space is needed, which is computationally expensive. If the mean and the variance of the dose are going to be used in the objective functions and constraints of the fluence optimization, they have to be calculated for every optimization iteration. To shift the burden of probability integrations away from the optimization iterations, dynamic dose deposition matrices have been developed. Multiplication of these matrices with the fluence profile vectors gives the expected value and variance of the voxel doses.

The dynamic dose deposition matrices also help to calculate the derivatives of the expected value and the variance of the dose with respect to the fluencies. These derivatives are usually needed in optimizers, such as Erasmus-iCycle (an optimization suite developed in Erasmus MC Cancer Institute in Rotterdam, The Netherlands), which is used in this thesis. Numerical derivations without dynamic dose deposition matrices, for example using the Monte Carlo method, are computationally expensive, especially since the derivatives must be calcu-lated for each optimization iteration as well.

Once the dynamic dose deposition matrices are calculated, inclusion of the expected value and the variance of the doses in the fluence optimization algo-rithm (called here dynamic optimization) is computationally as costly as the usual fluence optimization (called static optimization). In this thesis, dynamic optimization is done by substituting the dose (in static optimization) with the expected value of the dose in the objective functions and constraints. The variance is included in the form of the average of the voxel variances. The preprocessing effort to calculate the average of the variances is much less than to calculate the variances at every voxel separately.

In Erasmus-iCycle, the formulation of the optimization criteria is done using a wish-list. There the objective functions (costlets) are ranked according to their priorities, and each costlet has its own goal. For example for a prostate CTV, a costlet for dynamic optimization can be to maximize the mean dose

(10)

Summary

(expected dose), with the goal of 78 Gy. In the wish-list, some hard constraints are also prescribed, e.g. the minimum mean dose in the prostate CTV can be 74.1 Gy (95% of 78 Gy). Using an optimization method called 2pc, the costlets are optimized with respect to the constraints in two steps, so that a Pareto-optimal solution is obtained, where no improvement of a costlet can be made without deteriorating the others.

To evaluate the results of the fluence optimization, a new evaluation tool called dose probability volume histogram (DPVH) is introduced in this thesis to complement the conventional dose volume histogram (DVH). While in dynamic optimization the DVH depicts the volume percentage of the organ that receives a certain expected dose, the DPVH shows the probability that the delivered dose in the organ actually fulfills the optimization criteria.

Dynamic optimization using the dynamic deposition matrices has been tested for a simple cubic geometry and a prostate case. The results show that the margin-recipe solution prescribes a larger irradiated area for the same DVH in the CTV compared to dynamic optimization using expected values of the doses. Consequently using expected doses in the dynamic optimization damages the OAR less than the margin-recipe solution.

Adding costlets on the average value of the variances improves the DVH and DPVH for the CTV. Unfortunately our results are not yet conclusive for the OAR in this case. It has further been shown though, that the local movements are even more taken into account, when the variances are included.

(11)

(12)

Samenvatting

Radiotherapie is een van de methoden om kanker te behandelen, en dan wordt de pati¨ent meestal bestraald met r¨ontgenstraling met hoog-energetische foto-nen. Als daarbij intensity-modulated radiation therapy (IMRT) wordt toege-past, dan kan de fluentie (stralingsenergie per eenheid oppervlakte) van de stralingsbundels zo geregeld worden, dat het deel van het weefsel dat een hoge dosis krijgt zo goed mogelijk overlapt met het tumorgebied. Deze techniek zal hier verder met de term radiotherapie aangeduid worden.

Voordat de behandeling kan beginnen, moet er een behandelplan opgesteld worden zodat de tumor een hoge dosis krijgt en tegelijkertijd de omliggende gezonde weefsels maximaal gespaard worden. Het plan wordt meestal gemaakt op basis van een computertomografie-scan (CT-scan), waarop de tumor en de nabije gezonde organen aangegeven zijn. Deze CT-scan heet de plannings-CT. De behandeling wordt daarna meestal in een aantal fracties gegeven, verdeeld over meerdere dagen. Radiotherapie voor prostaatkanker kan bijvoorbeeld uit 40 fracties bestaan, waarbij in elke fractie een proportioneel deel van de dosis wordt gegeven.

Door verschuivingen en vervormingen zal de geometrie van de organen in het lichaam tijdens de bestralingen niet altijd hetzelfde zijn als die van de plannings-CT. Door deze geometrische onzekerheden kan het gebeuren dat sommige stralingsbundels het doelvolume (het tumorgebied, oftewel het cli-nical target volume, CTV) missen, zodat dat uiteindelijk een lagere dosis krijgt dan was voorgeschreven. Hierdoor zou het uitschakelen van de tumor kunnen mislukken. Om dit te voorkomen, wordt het bestraalde gebied meestal ver-groot met een marge, die bepaald wordt met een marge-recept. Een marge kan voorkomen dat het CTV een te lage dosis krijgt, maar jammer genoeg kun-nen daardoor ook gezonde orgakun-nen rond het CTV, de organs-at-risk (OAR) beschadigd worden. Hoe groter de marge, des te hoger de kans op schade aan de OAR.

Het doel van het onderzoek in dit proefschrift is om een methode te ontwik-kelen voor het berekenen van de optimale stralingsbundels, die met name re-kening houdt met de geometrische onzekerheden. Het bestraalde gebied wordt vergroot, net zoals met een marge, maar nu op een manier die lokaal beter af-gestemd is op de bewegingen en vervormingen van de organen, zodat de OAR

(13)

zo veel mogelijk gespaard kunnen blijven, terwijl het CTV een dosis krijgt die hoog genoeg is. In dit proefschrift worden eerst de locale bewegingen en vervormingen gemodelleerd voor een hele pati¨entenpopulatie. Dan wordt een handige methode beschreven waarmee in deze onzekere geometrie de verwachte stralingsdosis en de bijbehorende variantie kunnen worden berekend. Tenslotte wordt hiermee het bestralingsplan geoptimaliseerd.

Eerst worden de bewegingen en vervormingen van de organen gemodelleerd. Hiervoor is de principale-componentenanalyse (PCA) gebruikt. Het volume van de patiënt wordt gediscretiseerd in voxels, en met de PCA wordt geschat wat de belangrijkste bewegingsmodes zijn van de voxels en hoe waarschijnlijk die modes zijn. Helaas kan bij een nieuwe patiënt de PCA meestal niet uitgevoerd worden, omdat er meestal maar een paar (of zelfs maar één) CT-scans beschik-baar zijn voor die patiënt, en dat geeft veel te weinig informatie voor een PCA. Daarom is in dit onderzoek een manier ontworpen om de bewegingen van or-ganen van een populatie van patiënten te kunnen vergelijken en te gebruiken. Deze methode is getest op een reeks CT-scans van prostaten. Met de PCA werden zo de belangrijkste bewegingsmodes van prostaten bepaald, en deze kunnen verklaard worden als reële bewegingen en vervormingen. De resultaten zijn verder bevestigd door een foutenanalyse: de belangrijkste modes werden ook gevonden in de bewegingen van prostaten die niet in de database waren opgenomen.

De bewegingen en vervormingen worden gebruikt om de verwachtingswaarde en de variantie van de stralingsdosis in de prostaat te berekenen. Hiervoor moet over de waarschijnlijkheidsruimte worden ge¨ıntegreerd, wat veel rekentijd kost. Als de verwachtingswaarde en de variantie van de dosis gebruikt worden in de doelfuncties en de voorwaarden bij het optimaliseren van de fluentie, dan moe-ten ze in elke iteratiestap van de optimalisatie berekend worden. Om de last van deze integraties uit de optimalisatie te halen, zijn dynamische dosisdepositie-matrices ge¨ıntroduceerd. De verwachtingswaarde en de variantie van de dosis kunnen nu snel berekend worden door deze matrices te vermenigvuldigen met de fluentie-vector.

Bovendien zijn de dynamische dosisdepositie-matrices handig om van de verwachtingswaarde en de variantie van de dosis, de afgeleiden naar de flu-enties te berekenen. Deze afgeleiden zijn voor de meeste optimalisatiepro-gramma’s nodig, ook voor het in dit onderzoek gebruikte Erasmus-iCycle (een optimalisatie-omgeving ontwikkeld in het Erasmus MC Kanker Instituut in Rotterdam, Nederland). Zonder dynamische dosisdepositie-matrices is het nu-meriek berekenen van afgeleiden, bijvoorbeeld met de Monte-Carlomethode, heel rekenintensief, vooral omdat het ook in elke optimalisatie-iteratiestap moet gebeuren.

Als de dynamische dosisdepositie-matrices eenmaal berekend zijn, dan kost het optimaliseren van de fluentie waarbij de verwachtingswaarde en de variantie van de dosis gebruikt worden (hier de dynamische optimalisatie genoemd) niet meer rekentijd dan de gebruikelijke fluentie-optimalisatie (de statische

(14)

optima-Samenvatting

lisatie). In de doelfuncties en de voorwaarden van de dynamische optimalisatie wordt in dit onderzoek de dosis (van de statische optimalisatie) vervangen door de verwachtingswaarde van de dosis. Van de variantie van de dosis wordt de gemiddelde waarde over de voxels gebruikt. Het vooraf berekenen van de ge-middelde waarde van de variantie kost veel minder rekentijd dan het berekenen van de varianties in alle voxels apart.

In Erasmus-iCycle worden de optimalisatiecriteria in een wish-list gezet. De doelfuncties daarin (costlets) staan op volgorde van prioriteit, en elke costlet heeft een eigen streefwaarde. In een dynamische optimalisatie kan een costlet voor het prostaat-CTV bijvoorbeeld het maximaliseren van de verwachte dosis zijn, met een streefwaarde van 78 Gy. In de wish-list worden ook een aan-tal harde voorwaarden gesteld, bijvoorbeeld dat de verwachtingswaarde van de dosis in het CTV minimaal 74.1 Gy is (oftewel 95% van 78 Gy). Met de 2pc-optimalisatiemethode worden de costlets binnen de voorwaarden geopti-maliseerd in twee stappen. De gevonden oplossing is Pareto-optimaal, wat wil zeggen dat geen costlet meer kan worden verbeterd zonder andere costlets te verslechteren.

Om de resultaten van de fluentie-optimalisatie te kunnen beoordelen, wordt in dit proefschrift een nieuw evaluatie-instrument ge¨ıntroduceerd, het dose pro-bability volume histogram (DPVH), wat de informatie van het gebruikelijke dose volume histogram (DVH) compleet maakt. In een dynamische optimali-satie wordt in een DVH het percentage van het orgaanvolume getoond dat een bepaalde verwachte dosis krijgt. Een DPVH is een histogram van de kans dat de uiteindelijke dosis in de orgaanvoxels voldoet aan de optimaliseringsvoor-waarden.

De dynamische optimalisatie met dynamische dosisdepositie-matrices is ge-test op een eenvoudige kubusgeometrie en op een prostaatcasus. De resultaten laten zien dat de statische marge-oplossing een groter bestraald gebied geeft, vergeleken met de dynamische optimalisatie met de verwachtingswaarde van de dosis, terwijl het DVH in het CTV hetzelfde blijft. Daarom zal het bestra-lingsplan van de dynamische optimalisatie minder schade aan de OAR geven dan dat van de statische optimalisatie.

Het toevoegen van een doelfunctie op de gemiddelde waarde van de varian-ties verbetert het DVH en het DPVH voor het CTV nog verder. Helaas kunnen er nog geen conclusies worden getrokken voor de OAR in deze casus. Als de varianties worden meegenomen in de optimalisaties, dan is in onze resultaten wel te zien dat er nog meer rekening wordt gehouden met de locale bewegingen.

(15)

(16)

CHAPTER

1

Introduction

Abstract – This chapter is an introduction to cancer and radiotherapy. The treatment procedures, the concept of intensity modulated radiotherapy (IMRT), as well as the evaluation of the treatments using both physical and biological parameters are described. Some common terms which are used in this thesis are explained. The main problem of geometrical uncertainties in ra-diotherapy optimization is then described, followed by a review of current work on this problem. The chapter continues with a description of the approaches used in this thesis, and the way the thesis is organized.

1.1 Cancer

Cancer has become a leading cause of death. The World Health Organization (WHO) estimated that in 2008, 7.6 million people died of cancer: 13% of all deaths worldwide at that time [WHO 2013]. In 2012, the International Agency for Research on Cancer (IARC) reported that the death number has increased to 8.2 million [IARC 2012].

There are many types of cancer. According to the WHO [2014], the main types of cancer, in descending order of number of deaths, are lung, liver, stom-ach, colorectal, breast, and oesophageal cancers. Although not in the previous list, the American Cancer Society [2014] reported that prostate cancer, which is the test case in this thesis, is the second most frequently diagnosed cancer in men.

The cause of cancer is a malfunction of genes that control cell growth and division. According to the American Cancer Society [2014], this malfunction can be caused by external factors (e.g. tobacco, infectious organisms, chemicals,

(21)

and radiation) or internal factors (e.g. inherited mutations, hormones, immune conditions, or metabolism). The resulting abnormal cells can spread, not only in one organ to enlarge the tumor (the collection of these abnormal cells), but also to other organs. If the abnormal cells disrupt critical functions of the body, death can occur.

For most cancers, a staging system called TNM is used to describe the extent or spread of cancer at the time of the diagnosis [American Cancer Society 2014]. The TNM system assesses the cancer based on the extent of the primary tumor (T), the presence of cancer cells in the regional lymph node (N), and the presence of metastases (M), which is the spread of the cancer cells in distant organs. Based on these TNM categories, the cancer can be divided into stage 0, I, II, III, or IV (with IV the most advanced disease). The choice of therapy (surgery, radiation therapy, chemotherapy, hormone therapy, immune therapy, or targeted therapy) and the prognosis of the cancer are based on this staging.

1.2 Radiotherapy for cancer treatment

In this thesis, the focus is on cancer treatment using radiation (radiotherapy). There are two ways to deliver the radiation. It can be delivered by implanting (temporarily) a radioactive source inside the body: brachytherapy. Another way is to administer the radiation by using external beams of high-energy photons (X-rays), electrons, or protons, usually referred to as external beam radiation therapy (EBRT). In this thesis only EBRT with high-energy X-rays will be dealt with, and henceforward the term radiotherapy will refer to this modality of radiation.

Using X-ray beams for radiation, high-energy photons will ionize atoms in the cells. The electrons released by this ionization are the ones that actually damage the cells. The amount of damage is proportional to the amount of energy per unit mass which is absorbed in the irradiated tissue, usually called a dose and presented in the unit of Gray (1 Gy = 1 J/kg). The calculation of the dose is explained further in section 2.4. All irradiated tissues, including healthy tissues, may be damaged. The goal of radiotherapy is to kill the cancer cells, while sparing the healthy tissues as much as possible.

1.2.1 Treatment procedures

When a cancer patient is planned to have radiotherapy, a computer tomography (CT) scan of the to-be-irradiated area is made. The organs of interest are delineated: both the tumor (target) and the surrounding organs at risk (OAR). A treatment plan is then made to specify the settings of the treatment device, with the competing goals to deliver a high enough dose to target, and at the same time to spare the OAR as much as possible.

(22)

1.2 Radiotherapy for cancer treatment

Usually the whole treatment is given in several (daily) sessions (called frac-tions), e.g. 40. The dose at each fraction is scaled down from the total dose: the more fractions, the smaller the dose at each fraction. The purpose of this fractionation is to give the healthy cells time to recover.

1.2.2 Intensity-modulated radiation therapy (IMRT)

Intensity-modulated radiation therapy (IMRT) is one of the great developments in radiotherapy (for a historical review see [Webb 2003, Bortfeld 2006]). In IMRT, the fluence (radiation energy per unit area) of the radiation beams are modulated to achieve good radiation dose conformality to the target shape and optimally avoid dose delivery to the surrounding healthy tissues. Although there are other methods to do IMRT, for example using tomotherapy or robotic arms to deliver radiation (see for a review [Ahnesj¨o et al. 2006]), in this thesis the focus is on the delivery with a multi-leaf collimator (MLC). In this method, the radiation beam is led through a gantry which can rotate around the patient. The couch on which the patient lies can also rotate horizontally. By changing the rotation angles of the gantry φ and the couch θ, the patient can be irradiated from any direction (figure 1.1).

φ

θ

radiation beam

couch gantry

Figure 1.1 – Illustration of rotation angles for the gantry φ and the couch θ. The arrows start from 0◦ and point to the positive direction. The gantry angle φ ranges from 0◦to 360◦in clockwise direction, while the couch angle θ ranges from −90◦(to the left) to 90◦(to the right).

The collimator (MLC) on the gantry consists of several leaves made from tungsten (a highly radiation-absorbent material). There are typically around

(23)

40-80 leaves on each side of the (radiation) treatment field. By moving these leaves, the shape of the radiation field can be changed, for example to match the shape of the target projected in that beam angle (beam’s-eye-view, BEV). In IMRT, the collimator is also used to change the fluence profile of the radiation beam. There are two main ways to create the fluence profile using the collimator [Bortfeld 2006, Ahnesj¨o et al. 2006]. The first is the dynamic mode, where the collimator leaves move when the beam is still on. The second is the step-and-shoot method, where the radiation is given in multiple differ-ently shaped beams (segments) which are created by moving the collimator leaves when the beam is off. The fluence profile is realized by superposition of the radiation fields of these segments. The main criteria for the choice of segmentation are time efficiency and the amount of the unavoidable leakage radiation. Usually the smoother the fluence profile is, the easier the realization will be.

In practice, the desired dose distribution is obtained from a superposition of radiation from several beam directions (figure 1.2). Using several beams makes it easier to achieve better conformity of the dose to the target and/or better sparing of the healthy tissues.

Target

Figure 1.2 – Illustration of IMRT with several beams. The schematic shows an axial cut of the patient’s body, with the target in the middle. The four beams used are shown with the corresponding fluences and directions.

The treatment plan must specify the settings for each beam direction. Plan-ning usually starts with the calculation of the fluence settings, and then on how this fluence profile can be realized by sequencing the collimator-leaves. There is also an approach to directly optimize a series of realizable MLC settings (apertures) with respect to the many criteria for the dose distribution (also called column generation [e.g. Shepard et al. 2002, Romeijn et al. 2005].)

The calculation of the treatment plan (the fluence profile in this thesis) is sometimes called inverse planning, because it starts from the end-product

(24)

1.2 Radiotherapy for cancer treatment

(the required dose) and ends with input (the fluence profile). Although the relation between the dose and the fluence can be stated in a linear form (see section 2.4), it is not straightforward to find the right fluence profile for a given dose distribution. One of the problems is that a fluence can not be negative. Therefore, inverse planning is usually done with an optimization algorithm [Webb 2003]. One of the challenges is how to formulate the correct objective functions and constraints [Ahnesj¨o et al. 2006], besides of course finding a suitable fast and robust optimization algorithm.

In this thesis, we focus on the fluence optimization for a number of beams with already fixed angles, so beam-angle optimization (see e.g. [Breedveld et al. 2012]) is not considered here. In fluence optimization, settings for the smoothing of fluence profiles are such that the resulting profiles are realiz-able. We do not do the extra segmentation optimization (sometimes called collimator-leaf sequencing).

1.2.3 Evaluation of radiotherapy treatment plans

Usually the quality of a radiotherapy plan is evaluated by its dose volume histogram (DVH). This is a graph which shows the volume percentage of the organ that receives a dose equal or lower than a certain value. The ideal is to have 100 % of the target volume receive the prescribed dose (Dprescribed), and

0% of the OAR volume receive any dose (see figure 1.3). The percentage of the

D_prescribed 0

100

Volume (%)

Dose (Gy) ideal DVH for OAR

ideal DVH for target

Figure 1.3 – Ideal dose volume histogram (DVH) for target (full line) and organs at risk (dashed line). Ideally, 100% of target volume should receive the prescribed dose, and 0% of the organs at risk (OAR) should be irradiated.

target volume that receives a dose higher than Dprescribedis ideally 0% to avoid

overdosage. If the percentage of the target volume that receives Dprescribed is

less than 100%, then there is underdosage. Of course in reality, the ideal case is never achieved. For prostate cases, a 7% higher dose (maximum dose

(25)

constraint) and a 5% lower dose (minimum dose constraint) than the prescribed in the target volume are usually tolerated.

Some researchers have developed models for the biological response of liv-ing cells to radiation. The models usually assume that the cell survival after radiation follows binomial or Poisson statistics. Based on the death or sur-vival of the cells, the responses of the organs of interest are thus quantified, usually assuming that statistically all cells and patients respond identically. Examples of such quantities are the tumor control probability (TCP) and the normal tissue complication probability (NTCP), which quantify the effects of radiation dose to the tumor and healthy tissues respectively (see for a review [Schultheiss 1994]). The concept of TCP and NTCP is further developed into equivalent uniform dose (EUD), which for a given dose distribution estimates the uniform dose with the same biological response [Niemierko 1997]. The EUD can thus be used as a tool to compare and evaluate the quality of differ-ent treatmdiffer-ent plans.

Although the concept of using biological response-functions (e.g. EUD, TCP, NTCP) to evaluate radiotherapy treatments is powerful, their use in flu-ence optimization is not yet popular. One problem is the highly non-linear forms of the functions in terms of the fluence, which complicates optimization. The other is the validity of the models and the determination of their parame-ters. In this thesis, we will use only the measures based on the resulting physical dose distribution, such as DVH, to evaluate the radiotherapy treatment plans.

1.3 Geometrical uncertainties in radiotherapy

In the body the organs are moving and deforming. Some of these movements and deformations are involuntary, and therefore can not be stopped during the radiotherapy treatment. Due to these, the positions and shapes of the organs are not known for sure during the treatment. These uncertainties are usually called the geometrical uncertainties.

Geometrical uncertainties can reduce the quality of a radiotherapy treat-ment. IMRT technique has made it possible to deliver a radiation dose dis-tribution which conforms well to the desired dose disdis-tribution, even in quite complicated cases. Unfortunately, geometrical uncertainties make the actual treatment geometry not exactly the same compared to the CT scan which has been used to plan the treatment (planning CT scan). This can make the radiation beam miss the target and accidentally hit the healthy organs. Con-sequently the target can potentially receive less dose than required and the cancer might not be successfully eradicated. The main target to be consid-ered here is the clinical target volume (CTV), which is defined by International Commission on Radiation Units and Measurements (ICRU) [2010] as a volume of tissue that contains the gross demonstrable extent and location of the tumor

(26)

1.3 Geometrical uncertainties in radiotherapy

(called gross tumor volume, GTV) and/or subclinical malignant disease with a certain probability of occurrence considered relevant for therapy. The notion of subclinical malignant disease includes the microscopic tumor spread at the boundary of the primary-tumor GTV (thus outside of what can be observed, palpated, or visualized in a particular imaging modality), the possible regional infiltration into lymph nodes, and the potential metastatic involvement of other organs, despite their normal appearance on clinical and radiological examina-tions. If some healthy tissues (organs at risk, OAR) receive damaging dose, new health problems can arise, which will reduce the patient’s quality of life. It can even be fatal if the accidentally irradiated healthy tissue is a critical organ, such as the spinal cord.

The geometrical uncertainties are usually divided into systematic errors and random errors [van Herk 2004]. A systematic error occurs because the organ shape depicted in the planning CT scan is not the true mean shape of the organ. It is an unavoidable error because any CT scan can just capture a momentary shape of the organ. To get the true mean shape, a lot of CT scans should be taken for planning, which usually is not feasible in the practice. Setup errors, which introduce misalignments of the target with respect to the planning CT scan, will of course worsen the situation. Besides human errors, misalignments can also occur because patients are usually not stressed when CT scans are taken, but do feel stressed when the treatments are given. Stress can induce some involuntary muscle contractions, that can lead to changes in the position or shape of the organs. Systematic errors can also be caused by delineation errors, either due to the limited resolution of the CT scans or different human interpretations of the same scan.

The random error, on the other hand, quantifies the deviations from the mean shape, i.e. the movements and deformations of the organs. For the prostate case for example, a momentary contraction or expansion of the rectum and or the bladder can push or compress the prostate that lies between these organs. Soft tissues such as seminal vesicles can also move easily. Momentary stress can also cause organ deformation as mentioned above.

All the causes of geometrical uncertainties can happen randomly. This means that both systematic and random errors have statistical properties.

There are some efforts to minimize these geometrical uncertainties. Fixing the patient on the treatment table using a body cast is one example. Nowadays, fiducial markers (that can be clearly seen in a CT scan) are implanted in the organ, so that alignment of organs does not only depend on bone structures, but also on these markers, reducing the errors. Nevertheless there are always some uncertainties remaining.

To fully handle these stochastic geometrical errors, they must be incorpo-rated into the inverse planning calculation of the radiation dose. The planning CT scan can not be relied on to give the actual position and shape of the organs:

(27)

there will certainly be some movements and deformations. On-line adaptive image-guided radiotherapy is the ultimate solution. Unfortunately doing on-line imaging of the target while the therapy is being done and calculating the required beam settings in real time (directly on-line during the treatment) are still not practical or too expensive. Consequently off-line adaptation is still widely used. There are several ways proposed to do this, which are surveyed in the following section.

1.4 Survey of methods to deal with geometrical

uncertainties in radiotherapy

There are roughly three schools of thought on how to incorporate geometrical uncertainties in the inverse planning calculations:

• Enlarge the target volume with a margin, computed using a margin recipe • Calculate the dose for a number of scenarios, which are simulated possible

treatment geometries

• Use the probability distribution function of the random geometrical un-certainties directly in the objective functions and constraints of the opti-mization algorithm

These approaches are summarized in the following subsections. 1.4.1 Margin calculation using a margin recipe

One way to incorporate geometrical uncertainties in the planning process is to enlarge the CTV with a margin. The enlarged volume, usually called the planning target volume (PTV), is then used as the target in the fluence-profile optimization.

An example of the margin method is the one proposed by Stroom et al. [1999]. This method is based on the coverage probability, which is the probabil-ity of a voxel (discretization volume) in the treatment space to be occupied by the CTV. The calculation of this coverage probability is based on rigid body movements of the organs (no deformations are considered). Systematic and random errors are taken into account, which are both assumed normally dis-tributed. The relation between the margin and the variation of the systematic and random errors is called the margin recipe. Using the coverage probability, Stroom et al. estimate how large the margin should be so that the probability of a large percentage of the CTV to receive a high (physical) dose is high enough. They prescribe the margin recipe as 2Σ + 0.7σ, where Σ and σ are the standard deviations of the systematic and random errors respectively. These parameters are obtained from measurements from the patient population. Both Σ and σ

(28)

1.4 Survey of methods to deal with geometrical uncertainties in radiotherapy

can be different for different directions, but they are not prescribed differently for each voxel.

Instead of using the coverage probability, van Herk et al. [2000] compute directly the probability distribution of cumulative dose over a population of patients to derive a margin recipe. The movements are assumed to be only translations with no deformations (rigid body). Again both systematic and random errors are considered, and they are both assumed to be normally dis-tributed. The margin needed to have a minimum dose of 95% of the prescribed dose for 90% of the patients for a 3D object using this method is prescribed as 2.5Σ + 1.64(σ − σp), where Σ, σ, and σp are the standard deviations of

the systematic errors, random errors, and the penumbra respectively. Further σ2_{= σ}2

m+σs2+σ2p, where σmand σsare the standard deviations due to random

organ movements and random setup errors respectively. If σp is taken to be

3.2 mm, then the margin recipe can be simplified to 2.5Σ + 0.7σ, where here σ =pσ2

m+ σ2s.

A modification on the margin recipe is made by van Herk et al. [2002] by evaluating the effects of the systematic (rigid body) displacements with biological indices, such as the average tumor control probability (TCP) and the average equivalent uniform dose (EUD). With these modifications, the margin recipe is 2.5Σ + 0.7σ − 3 mm. There are several other margin recipes, which are summarized in [van Herk 2004].

The margin recipes are widely used in practice because they are simple and usually give good-enough dose distributions. Nevertheless, the method does not yet handle deformations (i.e. non-rigid movements) of the organs. It also does not yet take into account the correlation of movements or deformations of different organs, which can cause conflicting margin requirements, for example for the CTV and the neighboring OAR.

1.4.2 Weighted scenarios

Instead of using a fixed margin recipe to enlarge the irradiated volume, some re-searchers use scenarios to take into account the geometrical uncertainties. The scenarios sample the transformed shape of the patient’s organs. The trans-formation is based on (simulated) movements and detrans-formations of the organs. The scenarios are often weighted with the probabilities of the corresponding movements and deformations. Fluence profiles can be computed for each sce-nario, and methods exist to evaluate or compute the to-be-used fluence profile from these scenario-based fluence profiles [e.g Mageras et al. 1996, McShan et al. 2006]. The scenarios can also be used to compute the objective functions and constraints which will be used to compute the fluence profile [e.g Birkner et al. 2003, Yang et al. 2005, Chu et al. 2005]. By lumping together the geome-tries of the organs from the scenarios, an enlarged PTV area can be determined [Yan et al. 2000].

(29)

Mageras et al. [1996] propose a method to evaluate a fluence profile based on scenarios. First, for the patient at hand (the study patient ), an optimized fluence profile is computed without considering movements and deformations, called the nominal fluence profile. Then, by comparing the planning CT scans and the CT scans taken during treatments of a group of previous patients, called the reference patients, the movements and deformations of each refer-ence patient are computed. Scenarios are made by adjusting the movements and deformations of each reference patient to the study patient, by assum-ing uniform expansion or contraction between the organs of the two patients. For each scenario, the dose distribution for the adjusted study-patient’s organ is computed using the nominal fluence profile. These scenario-based dose-distributions are then ranked with respect to certain quantities, for example the dose at a specific point in the organ, the tumor control probability (TCP), the normal tissue complication probability (NTCP), etc. By showing the nom-inal dose volume histogram (DVH), together with those corresponding to the 10th _{and 90}th _{percentiles of the rank list, a so called confidence-limit DVH is}

created with 10% and 90% confidence limits respectively. This confidence-limit DVH is then used to evaluate the nominal fluence profile.

Also in the multiple instance geometry approximation (MIGA) method from McShan et al. [2006], scenarios are used to compute the to-be-used fluence profile. Here scenarios are made from replicas of the geometry of the patient which are subjected to pre-defined movement-deformation transformations. In general, non-linear deformations can also be accommodated, but McShan et al. only use rigid-body transformations in their publication. The optimized fluence profile for each beamlet is then computed for each scenario. The to-be-used fluence profile is the sum of all the scenario-based fluence profiles, weighted with the corresponding probabilities of the scenarios. McShan et al. assume this probability distribution to be a uniform over 7 geometries (two samples for translations in three directions + original).

The method from Birkner et al. [2003] differs from MIGA in the fluence com-putation part. No scenario-based optimized fluences are computed. Instead, the to-be-used fluence profile is directly computed in an optimization, where for each iteration a fluence profile is evaluated. Given a fluence profle, dose at voxel ν (D(ν)) can be computed. D(ν) changes due to geometrical uncertain-ties, which is simulated by some scenarios, taken from the first M treatment fractions of the same patient. These M geometry instances are given equal probabilities. The mean dose µ(ν) at voxel ν from these geometrical scenar-ios is used as the input for the evaluation function g(µ(ν)), for example TCP, NTCP, or maximum dose. The function g is then used as the objective function to find the optimal fluence profile.

Yang et al. [2005] use an objective function that maximizes the number of geometrical scenarios that meet a set of (biological) dose constraints, both for the CTV and OAR. The geometrical scenarios are made of replicas of the

(30)

geometry of the patient, which is subjected to some movement or deformation transformations, just like in MIGA. The to-be-minimized objective function is the sum of costlets (an objective function for one organ) evaluated for different scenarios, and weighted with their corresponding probabilities. The costlet can only have value 0 or 1, and it is only 0 if the constraints for the target (target dose larger than prescribed dose) and OAR (OAR dose lower than the tolerated dose) are fulfilled for the corresponding scenarios. Yang et al. assign the same probabilities for all scenarios.

In the work of Chu et al. [2005], the scenarios are used to incorporate the geometrical uncertainties into the constraints and objective function of the fluence-profile optimization. Chu et al. use rigid body translations of the origi-nal geometry as scenarios, where all the replicas are given the same probability weights and the original geometry a higher one. For each scenario, the dose for a given fluence profile (the evaluated fluence profile for each optimization iteration) is computed. The mean dose µ and the variance σ of the dose are computed by weighing the scenario-based dose with the corresponding probabil-ity of the scenarios. The constraints are made in the form that the probabilprobabil-ity of the accumulated dose (after all scenarios) in the CTV (OAR) is less (more) than a certain limit m, should not be more than a certain small threshold δ. To compute this probability, Chu et al. assume that the accumulated dose after all the scenarios is normally distributed. By using the distribution of probabil-ity P of standard normal random variable Z, the constraints are transformed into algebraic equations involving m, µ, σ, and z1−δ, which is chosen so that

P (Z > z1−δ) = δ. Dose volume histogram constraints can also be imposed by

restricting the sum of excess dose, i.e. dose that exceeds a certain a thresh-old. The seriousness of violating these constraints is measured by pre-defined penalty weights. The to-be-minimized objective function is the weighted sum of these constraint violations.

A different use of scenarios is proposed by Yan et al. [2000]. Here scenarios are used to compute the margin to make a PTV out of the CTV. The PTV will be then used as the new target like in margin-recipe method (section 1.4.1). Yan et al. use CT scans from the first several treatment fractions as scenarios. First k of these CT scans are aligned with respect to bony structures. A bounding volume P T V0(k) is then constructed as the convex hull (a kind of

union) that contains all the k CT scans. The bounding box P T V0(k) is then

used as the target volume for the dose calculation. The minimal number of initial CT scans k∗ is chosen so that if the same dose distribution (computed using bounding box P T V0(k∗)) is kept for the remainder of the treatment (the

CTV moves according to treatment fractions k + 1, k + 2, etc.), the maximal dose reduction in the CTV due to movements and deformations (i.e. due to CTV elements that lie outside the bounding box P T V0(k∗)) should be less

than or equal to ∆% of the prescribed dose. To take care of the rigid body motion of the bony structures with respect to the treatment room coordinates,

(31)

a margin is given around the P T V0(k∗). The effects of daily setup errors

are additionally convoluted with the planned dose distribution which uses the P T V0(k∗) + margin as the target, creating the overall dose distribution. The

optimal margin is chosen as the minimal margin that gives at most δ% overall dose reduction in the P T V0(k∗). ∆ and δ are prescribed a priori to be 2% and

1% respectively. The final PTV is the P T V0(k∗) + optimal margin.

1.4.3 Direct use of the probability distribution function

The methods reviewed in this section directly use the probability distribu-tion funcdistribu-tion of the geometrical uncertainties in the calculadistribu-tion of the objec-tive function and constraints in the fluence profile optimization. The prob-ability distribution is usually applied to compute the expected values (over time/treatment-fractions), and in some cases, the variance of some parame-ters. It is common to assume normal probability distributions, although there are some researchers who try to use histograms from measurements to estimate the probability distribution.

Witte et al. [2007] propose the use of the mean value (over time) of the tumor control probability (TCP) and of the normal tissue complication proba-bility (NTCP) in the objective function. The calculation of these mean values uses the direct integration of the probability distribution functions of the geo-metrical uncertainties. Both the systematic and the random errors are assumed to be normally distributed, with zero means and given standard deviations. The resulting integrals are too complicated, so they used simplifications: the random error is only assumed to blur the biologically effective dose (a part of TCP and NTCP formulations), whereas the systematic error is only assumed to shift (not rotate) the clonogenic cell density (the area where the tumor cells are most likely be cloned). The biological dose is assumed not to depend on the systematic error, and similarly the clonogenic cell density is not influenced by the random error. The other biological parameters (e.g. radiosensitivity parameters) are also assumed to be normally distributed.

In [Baum et al. 2006], the mean value of a costlet (local objective func-tion) is computed by multiplying the value of the costlet of an organ with the coverage probability of the voxels of the corresponding organ. The total objective function is the sum of all the mean values of these costlets (here in the form of equivalent-uniform-dose, EUD). Baum et al. expand the concept of coverage probability from Stroom et al. [1999] (subsection 1.4.1). It has now two parts: setup errors and organ motion or deformation errors. Both of them are stochastic variables. Organ motion/deformation is evaluated on a patient-based coordinate system (PCS). The probability PP CS(x) for a voxel

at position x in the PCS to be covered by the CTV or the OAR due to organ motion/deformation is computed by summing the number of (preliminary) CT scans where this voxel is covered by the CTV or the OAR, divided by the total

(32)

number of (preliminary) CT scans. On the other hand, the setup errors de-scribe the movement of the PCS with respect to the treatment room coordinate system, and are assumed to have a normal distribution. The mean and variance for the setup errors are taken from measurements in the patient population. The final coverage probability Pcoverage is then computed as a convolution of

both these errors:

Pcoverage(x) =

Z

Psetup(y)PP CS(x − y)dy (1.1)

where Psetupand PP CSare probabilities of setup errors and organ motion errors

respectively.

For a prostate case, Unkelbach & Oelfke [2005a] utilize the probability distributions of the geometrical uncertainties in the calculation of the mean value µiand the variance σi2of the dose at voxel i. The probability distributions

are assumed to be normal, and there is no distinction between systematic and random errors. The objective function for the fluence profile optimization is in the form of an integral over the whole volume of interest V :

E = X i∈CT V αCT V (µi− DCT V)2+ 1 Nf rac σ2_i +X n " X i∈OARn αn(µi− Dnmax) 2 + # (1.2) where DCT V is the prescribed dose at the CTV, and Dnmax is the maximum

dose at the n-th OAR. The term (µi− Dmaxn )+= µi− Dnmaxif µi> Dnmaxand

zero otherwise. The objective function E is minimized, and the weights αCT V

and αn are given.

In their other works, Unkelbach & Oelfke [2004, 2005b] use only the first two terms of equation 1.2 in the objective function of the fluence profile opti-mization (only minimize the mean and the variance of the dose in the CTV). In [Unkelbach & Oelfke 2004], the systematic and random error probabilities are separated for a simple circular CTV with no OAR, where the total prob-ability is the product of the probabilities of both these errors. For this same CTV geometry, Unkelbach & Oelfke [2005b] have used Bayesian inference to incorporate the probability that the modeled mean and variance are true given the data, which is then used as an additional probability to compute the mean and the variance of the dose.

Chan et al. [2006] use the histograms of voxel displacements (from several measurements) to estimate upper and lower bounds of the probabilities of voxel positions. A nominal probability is chosen from the average of the measure-ments. In the objective function, the mean dose is computed by allowing the

(33)

voxels to move, computing the dose for each resulting position, and then sum-ming the dose for each movement weighted with the nominal probability of that movement. For the constraints however, the probability is taken from the upper or the lower bound, such that the worst cases are considered.

The fluence optimization presented in [Moore et al. 2012] incorporates the random and systematic errors differently. Here two stages of optimization are used. In the first stage the initial radiotherapy plan (fluence field assignments) is computed. This stage starts with an assumption of a uniform fluence field, which is then optimized by taking into account the random and systematic errors. First the fluence is convolved with the assumed normally distributed probability of the movements caused by the random errors. Thus the fluence is blurred with the random errors. Next the geometrical uncertainties due to the systematic errors are sampled iteratively in the optimization process, doubling the number of samples for each iteration. The systematic errors are used to shift the dose distribution. After the initial plan is obtained, the weights of each objective function is adjusted by using an iterative binary search until the satisfactory solution is found.

A different incorporation method for the random and systematic errors is also used by Bohoslavsky et al. [2013], who use both the physical dose and the EUD (equivalent uniform dose) in the objective functions. The random errors are used to blur the dose distribution, and it is done by convolving the dose distribution with the probability distribution of the movements, which is assumed to be normal with possibly different standard deviations in different directions. The organs are assumed to be rigid (no deformation), and rotational random errors are treated as localized translations. The blurred dose and the blurred EUD (computed by substituting the dose in the EUD formulation with the blurred dose) are then rotated and shifted according to the systematic errors. Substitutions of the blurred, shifted, and rotated dose and EUD into the original original objective function create the blurred, shifted, and rotated objective function. The probability pj of the j-th sample of systematic error is

used to weigh the corresponding contribution to the probabilistic k-th objective function Fk

P (different objective functions are prescribed for different regions

of interest): F_Pk = n X j=1 pjFk ∗ j (1.3) where Fk∗

j is the blurred (due to the random errors) and displaced (according to

the j-th systematic-error sample) k-th objective function. P is the probability corresponding to the desired confidence limit for the systematic errors, and the number n is chosen such that the cumulative probability Pn=

h Pn j=1pj i ≤ P . Fk

(34)

1.5 Approach used for this thesis

In this thesis the focus is on the fluence optimization for IMRT, while taking into account the geometrical uncertainties. The radiation beam angles are as-sumed to be given, and no separate optimization is done here for the realization of the resulting fluence profile in the collimator sequencing. In the fluence op-timization, the fluence profiles are already smoothed such that the optimized fluence profile is realizable.

Although promising, the use of biological parameters such as EUD, TCP, and NTCP, are not yet clinically popular. Optimization using these parameters are more difficult due to the non-linearity and complexity of the functions. This is also the reason why in this thesis the main focus is on using physical dose parameters.

As described in section 1.4, many methods were developed to deal with geometrical uncertainties, although clinically, the margin recipe as explained in section 1.4.1 is still the most popular approach. The method developed in this thesis mainly follows and expands the method of Unkelbach & Oelfke (section 1.4.3).

The focus of this thesis is on accounting for random geometrical errors in treatment planning. Systematic errors, although reduced by the application a of body cast, external markers, or implanted fiducial markers, are generally not zero. In this thesis they have been ignored. Further research is needed to also account for these errors.

To deal with random errors, the following steps are used:

1. Make a (statistical) model of the movements and deformations, using principal component analysis, aided by a shape transformation to com-bine information from a population of patients.

2. Develop a strategy to use the model for computations of relevant vari-ables (mean dose and variance of the dose), using dynamic dose deposition matrices. The goal of using these matrices is to be able to use the exist-ing fluence optimization algorithms, now with geometrical uncertainties included.

3. Incorporate the model in the fluence optimization, in this case using the 2-phase -constraint (2pc) and primal-dual interior point methods. The results are compared to conventional margin recipe solutions.

1.6 Structure of the thesis

The main content of this thesis is a collection of papers, which are made ac-cording to the steps explained in section 1.5:

(35)

1. A population-based model to describe geometrical uncertainties in radio-therapy: applied to prostate cases, published in Phys. Med. Biol. 56 (2011) 1045-1061

2. Computation of mean and variance of the radiotherapy dose for PCA-modeled random shape and position variations of the target, published in Phys. Med. Biol. 59 (2014) 289-310

3. Radiotherapy fluence optimization using principal component analysis to model the geometrical uncertainties, submitted to Phys. Med. Biol. To make the thesis self-contained, a summary of the tools used in the papers is presented in chapter 2. The presentations are more in-depth than the ones in the papers, but they are not meant to be exhaustive. Interested readers are referred to the references for more detailed information.

At the end of the thesis, general conclusions and recommendations for future work are presented.

(36)

CHAPTER

2

Tools

Abstract – In this chapter several tools for the main results of the thesis (chapters 3, 4, and 5) are described. The first tool is the shape transformation to create the database for the statistical analysis, followed by the principal component analysis (PCA) which is the statistical tool used to analyze the database of organ shapes. The third tool is Erasmus-iCycle, the software suite used in this thesis for dose calculation and fluence optimization. The methods used in Erasmus-iCycle for these two purposes are the fourth and the fifth tools described in this chapter.

2.1 Shape transformation

Geometrical uncertainties are modeled in this thesis with a principal compo-nent analysis (PCA), which is a kind of statistical analysis (see section 2.2). Unfortunately, implementing this idea is not straightforward in practice. The number of CT-scans available for one patient is usually very small, not enough to do a statistical analysis like PCA. For a new patient, usually even only one CT-scan is available for the radiotherapy planning. The CT-scans from other patients cannot be directly used as additional data, because the shape of the organs involved are not compatible: they come from different persons with different heights, weights, ages, etc.

One way out of this problem is to transform the organ shapes from the CT-scans of other patients so that they can be compared with the data of the patient at hand (see section 3.2.4). The tool used in this thesis for this purpose is the thin-plate-spline (TPS) robust-point-matching (RPM) method of Chui & Rangarajan [2003], which has been further developed for radiotherapy applications by V´asquez Osorio et al. [2009].

(37)

There are two main tasks in matching organ shapes from two CT-scans: the correspondence and the transformation. The first task is to determine which points in one scan correspond to the ones in the other scan, even when the shapes have different numbers of discretization points, which happens in practice. The second task is to find the correct equations that transform a point in one scan to the corresponding point in the other scan. Solving both problems simultaneously is very difficult, if not impossible. TPS-RPM solves these two problems iteratively in alternation. The solution of one problem becomes the initial condition of the other one.

For TPS-RPM, two sets of points are defined: the reference set and the deforming set. The reference set consists of coordinates of the points on the surface of the target shape, and is denoted in this section as points xi ∈ R3,

i = 1, · · · , Nx. The deforming set points lie on the surface of the original shape,

and are denoted by points uj∈ R3, j = 1, · · · , Nu. The goal is to transform the

original shape to match the target shape. The number of points in reference set Nxis not always the same as the number of points in the deforming set Nu.

The points in both sets are not always uniformly distributed. This irreg-ular distribution is not desirable, because correspondence between irregirreg-ularly spaced points can be problematic and error prone. Therefore, before the trans-formation is made between shapes, some uniformly distributed control points are created first. The control points on the target and original shape surfaces then replace the reference and deforming sets respectively. The number of control points should be kept small to reduce the computational effort.

To create the control points, a starting point is chosen on the triangulated surface of the shape. This point becomes the center of a sphere with radius r. All points of the set that lie inside this sphere are replaced by the centroid of these points. This procedure is repeated for the remaining points on the surface, until in every sphere there is only one point. The bigger the radius r, the fewer the control points are. In this thesis the radius r is chosen such that there are around 600-700 control points for every shape.

The correspondence matrix is ideally a binary matrix: 0 for not correspond-ing, and 1 for corresponding. However, integer computation makes the opti-mization method involved in the calculation difficult. Therefore in TPS-RPM, a soft-assign approach is used, where the correspondence matrix starts with continuous values between 0 and 1, making partial matches between points in reference and deforming sets possible. For example, one point in the deforming set can then have a correspondence with a cluster of points in the reference set. An element of this correspondence matrix mij describes the correspondence

between reference point xi and deforming point uj. To make sure that every

(38)

2.1 Shape transformation

on mij:

Nx

X

i=1

mij = 1 for all j and

Nu

X

j=1

mij = 1 for all i. (2.1)

The correspondence matrix slowly evolves into a binary (or nearly binary) one using a so-called simulated annealing process. Here a variable T called temperature is introduced. The annealing process starts with a large value for T . This means the correspondence between points is loose (fuzzy), i.e. the correspondences between points ui and vj are not influenced very much by the

distance between the points: mij = 1 T exp −(xi− ϕ(uj)) 2 2T (2.2) where ϕ is the transformation function. The lower the temperature T , the more restricted the correspondence to the points in close proximity. Special arrangements are made for outliers, which can be found in [Chui & Rangarajan 2003].

TPS-RPM tries to find a transformation of the form:

ϕ(uj) = Auj+ Wφ(uj) (2.3)

which transforms the deforming point uj to match the reference point xi.

Ma-trix A is a linear transformation maMa-trix, whereas φ is the thin-plate-spline kernel which is weighted by the entries of matrix W. With the second term, non-rigid transformations (allowing shape deformations) are also accommo-dated. The transformation function ϕ is obtained by solving the minimization problem of the energy function E:

min ϕ E(ϕ) = minϕ   Nv X j=1 kyj− ϕ(uj)k2+ αT kΨϕk2   (2.4) where yj = Nx X i=1 mijxi (2.5)

is the estimated position of the point-set in the reference set (target shape) that corresponds to the deforming point uj. The other parameters are stiffness

parameter α and smoothing operator Ψ applied to ϕ (see for details [Chui & Rangarajan 2003, V´asquez Osorio et al. 2009]).

The higher the temperature T , the smoother the energy function E. In this thesis, the initial temperature is set to be the maximal distance between points belonging to the same organ. Decreasing T slowly will help the minimization

(39)

to avoid the traps of local minima. The temperature decreases over time with rate ra:

Tk+1= raTk (2.6)

where k is the iteration number. ra is chosen to be 0.93 in this thesis.

The iteration starts with a large stiffness parameter α, restricting the trans-formation to mostly affine (preserving collinearity, i.e. all points lying on a line initially still lie on a line after transformation, and proportions on lines, e.g. the midpoint of a line remains the midpoint after transformation), and therefore avoids a too wild mapping. This parameter α is decreased with the same annealing rate ra as temperature T , allowing gradually more flexibility.

The initial α is chosen to be 2 in this thesis.

Once the transformation function is computed with equation 2.4 using the initial correspondence values, temperature T is decreased according to equation 2.6. The new temperature T and transformation function ϕ is used to update the correspondence elements mij according to equation 2.2. To conform mij

to the constraints in equation 2.1, a normalization is done. The iteration is repeated until the final temperature Tfinalis reached (Tfinal= 5 in this thesis).

Tfinal is chosen to be the desired mean square distance between the nearest

neighbors of the transformed points [V´asquez Osorio et al. 2009].

2.2 Principal component analysis

Using the transformation described in the previous section, the shapes of the organ-of-interest taken from the CT scans of different patients can be made compatible. In chapter 3, a scheme to extract and combine the information about the movements and deformations of that organ from a population of patients is described. As a result of the scheme, we will have a set of vectors that show the deviations from the mean shape of the organ. Principal component analysis (PCA) is then used to model the movements and deformations of the organ by extracting the dominant modes from the deviation vectors.

2.2.1 Derivation of the principal components and their probability distributions

The shape vector p ∈ R3M _{is the collection of 3D-coordinates of all M}

dis-cretization points of the observed organ. Due to geometrical uncertainties, p is a random variable. The covariance matrix Ξ of p is defined as:

Ξ = h(p − p0)(p − p0)Ti (2.7)

(40)

2.2 Principal component analysis

The probability distribution P of shape vector p is usually assumed to be normal, and therefore can be written as:

P (p) = 1 (2π)3M2 |Ξ|1/2 exp −1 2(p − p0) T_Ξ−1_{(p − p} 0) . (2.8) Here |Ξ| is the determinant of covariance matrix Ξ.

In practice, covariance matrix Ξ is not known. It can be estimated though by the sample covariance matrix C, defined as:

C = 1 S − 1 S X i=1 (pi− p0) · (pi− p0)T. (2.9)

where pi is the i-th sample of shape vector p and S is the total number of

samples. The mean shape p0 is also estimated using the samples, i.e. p0 = 1

S

PS

i=1pi. From now on in this section we will take C ≈ Ξ.

The approximation of the covariance matrix Ξ by the sample covariance matrix C confines the observed probability space into a subspace of the original space. The dimension of the observed probability space is now S instead of 3M . With the principal component analysis (PCA) one tries to find an even smaller subspace which still contains the important features. This is done by using a transformation for random variable vector (p − p0) so that only a few

(let’s say L < S) transformed variables are needed to represent the (sample) covariance matrix C. This is a further reduction of the probability space, now to a subspace with dimension L.

The new transformed variables are called principal components:

ck= qTk(p − p0), k = 1, · · · , L (2.10)

where qk is a transformation vector. This transformation should fulfill several

requirements:

1. The variances of each ck (V ar(ck)) should be maximized (maximizing

the important features contained in the subspace).

2. Each principal component ck should be uncorrelated to the other

princi-pal components ci, i.e. Cov(ck, ci) = Cov(ci, ck) = 0; i 6= k (cross

covari-ances are zero).

3. Each transformation vector qk is normalized, i.e. qTkqk= 1.

It turns out that the suitable transformation vectors qkare the orthonormal

Modeling geometrical uncertainties for radiotherapy plan optimization without margins

Modeling

geometrical uncertainties

for radiotherapy

plan optimization

without margins

Eka Budiarto

MODELING GEOMETRICAL UNCERTAINTIES

FOR RADIOTHERAPY PLAN OPTIMIZATION

WITHOUT MARGINS

Modeling geometrical uncertainties for

radiotherapy plan optimization

without margins

PROEFSCHRIFT

Eka BUDIARTO

Summary

Samenvatting

Contents

CHAPTER

1

Introduction

1.1

Cancer

1.2

Radiotherapy for cancer treatment

φ

θ

1.3

Geometrical uncertainties in radiotherapy

1.4

Survey of methods to deal with geometrical

uncertainties in radiotherapy

1.5

Approach used for this thesis

1.6

Structure of the thesis

CHAPTER

2

Tools

2.1

Shape transformation

2.2

Principal component analysis