Towards the elimination of Monte Carlo statistical fluctuation from dose volume histograms for ...

27 Pages·1999·0.59 MB·English

Phys.Med.Biol.45(2000)131–157.PrintedintheUK PII:S0031-9155(00)06802-0
Towards the elimination of Monte Carlo statistical fluctuation
from dose volume histograms for radiotherapy treatment
planning
JosepSempau†‡andAlexFBielajew†
†DepartmentofNuclearEngineeringandRadiologicalSciences,TheUniversityofMichigan,
AnnArbor,Michigan,USA
‡InstitutdeTe`cniquesEnerge`tiques,UniversitatPolite`cnicadeCatalunya,Diagonal647,
08028Barcelona,Spain
Received9August1999
Abstract. TheMonteCarlocalculationofdoseforradiotherapytreatmentplanningpurposes
introduces unavoidable statistical noise into the prediction of dose in a given volume element
(voxel).Whenthedosesinthesevoxelsaresummedtoproducedosevolumehistograms(DVHs),
thisnoisetranslatesintoabroadeningofdifferentialDVHsandcorrespondinglyflatterDVHs.
Abruteforceapproachwouldentailcalculatingdoseforlongperiodsoftime—enoughtoensure
thattheDVHshadconverged.Inthispaperweintroduceanapproachfordeconvolvingthestatistical
noisefromDVHs,therebyobtainingestimatesforconvergedDVHsobtainedabout100timesfaster
thanthebruteforceapproachdescribedabove.Therearetwoimportantimplicationsofthiswork:
(a)decisionsbaseduponDVHsmaybemademuchmoreeconomicallyusingthenewapproachand
(b)inversetreatmentplanningoroptimizationmethodsmayemployMonteCarlodosecalculations
atallstagesoftheiterativeproceduresincetheprohibitivecostofMonteCarlocalculationsatthe
intermediatecalculationstepscanbepracticallyeliminated.
(Somefiguresinthisarticleappearincolourintheelectronicversion;seewww.iop.org)
1. Introduction
Withtheadventofaffordablehigh-speedcomputersandthedevelopmentoffasterMonteCarlo
based algorithms (Neuenschwander and Born 1992, Neuenschwander et al 1995, Keall and
Hoban1996a,b,Kawrakowetal1996,Sempauetal2000),theMonteCarlo(MC)methodis
becominganimportanttoolformedicalphysicistsperformingphotonandelectrontreatment
planningforradiotherapycancertreatment(Bielajew1994,1995,Mohan1997).
Theinformationprovidedbyatreatmentplanningsystemcanberepresentedgraphically
inavarietyofdifferentforms,forexampleisodoselinesanddosevolumehistograms(DVHs).
GivenatargetvolumeV withinapatient’sbodyandcertainirradiationconditions,theDVH,
say p.x/, is defined such that p.x/dx equals the fraction of V with a dose in the interval
[x;x+dx]. Thus,theDVHdefinestheprobabilitydensityfunction(PDF)offindingadosex
inV. TheintegralDVH
Z
x
F.x/D dx0p.x0/ (1)
0
is also an important representation, since F.x / − F.x / represents the fraction of V that
2 1
receivesadosebetweenx andx .
1 2
0031-9155/00/010131+27$30.00 ©2000IOPPublishingLtd 131
132 JSempauandAFBielajew
DVHshavebeenusedasatesttocomparetheaccuracyandperformanceofMCcodes
(Hartmann-Siantaretal1997). Tothisend,thepatient’sgeometry,includingthetargetvolume,
is divided into voxels as provided by CT or MRI scans. Voxels are normally shaped as
rectangular boxes with a square base with characteristic lengths of the order of millimetres.
As the particle transport simulation is carried out, the dose in each voxel is scored in an
accumulator. Whenasufficientnumberofparticlehistorieshasbeenprocessed,thecalculation
stops and each voxel is assigned to a dose bin, to which it contributes with a count. The
histogram of counts in the different bins corresponds to the discretization of the underlying
continuousPDF,p.x/. Wenotethatthediscretizationisperformedbothinspace(voxels)and
indose(bins).
The remainder of the paper is organized as follows. In section 2 we give a general
discussionofthestatisticalnoiseinthecalculateddoseinavoxel. Insection2.1weprovide
a general proof that the variance of the dose is proportional to the dose and argue that
the distribution of dose should be distributed according to a Gaussian. We demonstrate in
section 2.2 that calculated DVHs and true (or converged) DVHs, those produced after an
infinitenumberofparticlehistorieshavebeensimulated,arerelatedbyaFredholmequation
of the first kind, an integral equation that yields the calculated (finite number of particle
histories) DVH after integration over the true DVH and a ‘response function’, the Gaussian
derivedinsection2.1. Insection2.3wetestourGaussianresponsemodelanddemonstrateits
usefulnessandinsection2.4wefindpracticallowerlimitsofthenumberofhistoriesforwhich
themodelcanbeemployedaswellasadeterminationofthenumberofhistoriesrequiredto
ensureconvergenceofthecalculatedDVHs. Insection3wedevelopadeconvolutionmethod
thatallowsustoestimatetheconvergedDVHfromanapproximateDVHobtainedwithorder-
of-magnitudefewerparticlehistories. Insection4weshowhowthenewmethodperformson
avarietyoftestcasesbeforetheconclusions.
2. MonteCarlogeneratedstatisticalnoise
2.1. Dosevarianceinavoxel
Letussuppose,asisusuallythecase,thatMCsimulationsarecarriedoutwiththetransport
ofelectronsandpositronsturnedon,eitherbecausetheyaretheprimaryparticlesorbecause
theyareproducedwhenphotonsinteractwithmatter. Inanycase,thisimpliesthatallenergy
depositionsareactuallymadebychargedparticles. Forthesakeofbrevity,wewillreferonlyto
electronshereafter,sincepositrontransportinvolvessimilarmethodsexceptforthepossibility
ofannihilation,whichdoesnotimplyanydirectenergydeposition.
VariousofthecurrentlyavailableMCcodesrelyon‘mixed’schemestosimulateelectron
energy losses (Nelson et al 1985, Baro´ et al 1995, Sempau et al 2000). The basic idea
is to separate inelastic interactions into two categories, depending on whether the energy
lost is above (hard) or below (soft) some threshold, E . Hard events are simulated in a
cut
detailed way, that is, they are processed one at a time. Soft events, on the contrary, are
simulated in a condensed way, which means that many soft interactions are accumulated in
a single fictitious event. Delta rays are generated only by hard events and therefore energy
losses larger than E involve the generation of secondary particles but no deposition on
cut
the spot. This motivates the definition of the restricted stopping power, S , or mean energy
r
lost per unit step length for energy losses below E . For the majority of applications,
cut
including radiotherapy, it is plausible to set E equal to the absorption energy, E ,
cut abs
below which electrons are considered to be absorbed and their kinetic energy deposited
locally.
Dosevolumehistogramsforradiotherapytreatmentplanning 133
MC calculations provide an exact solution to the Boltzmann transport equation except
forthefluctuationsinherentintherandomnatureofthemethod. Thesefluctuationsmaybe
made arbitrarily small by simulating a sufficiently large number of particle histories. Let x
be the dose deposited per primary particle in a given voxel when the number of simulated
particlehistoriesisinfinite,thatis,whentherearenostatisticalfluctuations,andlet1V and
(cid:26) bethevoxelvolumeandaveragedensityrespectively. Aftercompletionofafinitenumber
ofhistories,N,theaveragedosex definedasthequotientbetweentheMCsimulatedenergy
depositionperhistory,1E=N,andthevoxelmass,(cid:26)1V,isanunbiasedestimatorofx.
A voxel will be said to have been ‘hit’ when an electron, either primary or secondary,
traverses it. Hereafter, n and (cid:15) stand for the random variables representing the number of
hits produced in a given voxel and the energy deposited by a particle that hits this voxel,
respectively. As1Eandconsequentlyx arerandomvariables,theyhaveanassociatedmean
andvariance. Toinvestigatethesequantities,considerxexpressedintermsofnand(cid:15),namely
P
1 1E 1 n (cid:15)
x (cid:17) D iD1 i (2)
N (cid:26)1V N (cid:26)1V
where(cid:15) istheenergydepositedbytheithhit.
i
Forvoxelsizesoftheorderofsomemillimetres,theprobabilityofobservingamultiple
hit on a single voxel is small compared with the probability of the voxel receiving a single
hit,P. Hence,nfollowsapproximatelyabinomialdistributionwithN trials. Moreover,for
broadbeams,suchasthosefoundinradiotherapyapplications,P issmallforthemajorityof
voxelsandNP isratherlarge. Theseconsiderationsleadtotheconclusionthatthebinomial
distributionofncanbeapproximatedbyaGaussianwithvariance,(cid:27)2.n/,givenby
(cid:27)2.n/Dhni: (3)
Ifnisdistributedaccordingtoabinomialdistribution,itisshownintheappendixthat
(cid:28) (cid:29)
Xn
(cid:15) Dhnih(cid:15)i: (4)
i
iD1
Thatistosay,theaveragevalueoftheenergydepositedinavoxelisgivenbytheproductof
theaveragenumberofhitsmultipliedbytheaverageenergydepositionperhitonthevoxel.
Thisallowsequation(2)tobewrittenas
hni h(cid:15)i
hxiD (cid:17)x (5)
N (cid:26)1V
wherewehaveredefinedxtomakethenotationalittlelesscumbersome. Intheappendixitis
alsoproventhat
(cid:18) (cid:19)
Xn
(cid:27)2 (cid:15) Dhni(cid:27)2.(cid:15)/+h(cid:15)i2(cid:27)2.n/ (6)
i
iD1
indicating that the variance of the Monte Carlo estimated energy deposition in the voxel
originatesfromthetwodistributedquantities,nand(cid:15).
Substitutingfromequation(3)resultsin
(cid:18) (cid:19)
Xn
(cid:27)2 (cid:15) Dhni.(cid:27)2.(cid:15)/+h(cid:15)i2/Dhnih(cid:15)2i (7)
i
iD1
whichleadstoanexpressionfortheestimatedvarianceofthedoseinthevoxel
hnih(cid:15)2i x h(cid:15)2i
(cid:27)2.x/D D : (8)
N2(cid:26)21V2 N1V (cid:26)h(cid:15)i
134 JSempauandAFBielajew
3
)
g
m/2 2
c
V
e
M

(
r S/r 1
0
0 4 8 12 16 20
Energy (MeV)
Figure1. Restrictedmassstoppingpowerforwater(fullcurve),bone(shortdashedcurve)and
air(longdashedcurve). Ecut wassetto200keV.Datahavebeenobtainedfromthepreparation
programassociatedwiththeDPMcode(Sempauetal2000).
We note that only energy losses due to soft collisions are involved in the distribution of (cid:15).
Theenergystragglingduetothesesoftcollisions,whichaccountsforadispersionin(cid:15) when
electronswithidenticalinitialconditionstravelagivenpathlength,canbeneglectedforlengths
equal to 1 mm or larger. This ‘continuous slowing down approximation’ (CSDA) has been
adoptedsuccessfully(RogersandBielajew1990)formanymedicalphysicsapplicationsby
theEGS4MonteCarlocode(Nelsonetal1985),foranypathlength. Massrestrictedstopping
powers(S =(cid:26))formaterialsrelevanttoradiotherapyapplicationshaveaverysmallvariation
r
withenergy,exceptatlowenergies,andarefairlyindependentofthematerialaswell. Thisis
seeninfigure1. ItisplausibletoapproximateS =(cid:26) byaconstantvalueforallvoxelsinthe
r
estimateofthevarianceofthedosetoavoxel.
Takingtheseapproximationsintoaccount
(cid:18) (cid:19)
S
(cid:15) D(cid:26) r ‘ (9)
(cid:26)
where‘standsforthepathlengthtraversedthroughavoxel. Consequently,equation(8)can
berewrittenas
(cid:18) (cid:19)
x S h‘2i .S =(cid:26)/(cid:17)
(cid:27)2.x/D r Dx r (10)
N1V (cid:26) h‘i N1V2=3
wherethedimensionlessparameter
1 h‘2i
(cid:17)(cid:17) (11)
1V1=3 h‘i
can,inprinciple,changefromvoxeltovoxeldependingonthevoxelshapeandonthedegree
ofisotropyoftheradiationfieldatthevoxellocation.
Forallpracticalradiotherapysituations,however,(cid:17)hasaveryrestrictedrangeofvariation
and can be considered as approximately constant inside the target volume. To justify this
assertionitwillbeconsideredthat,duetothevoxelsmallness,electronstravelalongalmost
straightlinesandthatparticlefluxishomogeneousinsideagivenvoxel. Thus,h‘2i,h‘iand
then (cid:17) can be readily obtained by means of an MC numerical integration for any radiation
Dosevolumehistogramsforradiotherapytreatmentplanning 135
field. Usingthisapproach,wefindformonodirectionalbeamsandcubicvoxels,that(cid:17)spans
theinterval[0:86;1:00]asthebeamdirectionchanges. Asanyelectronradiationfieldcanbe
thoughtofasacontinuoussuperpositionofmonodirectionalbeams,(cid:17)isboundedbythesame
intervalforanyfield.
Inpractice,collimatedradiationcomingoutofanacceleratorheadwillenterthepatient’s
geometry nearly parallel to a voxel face, say z D 0 where the z-axis is directed along the
patient’s body. This includes the cases in which the direction is normal to x D 0 or y D 0
(whichgives(cid:17)D1),oritgoesalongthediagonal(1,1,0)((cid:17)D0:94). Aselectronspenetrate
deeperintothepatient,theelectronangulardistributionbecomesprogressivelymoreisotropic
as the kinetic energy decreases and the electrons scatter. For a completely isotropic field,
(cid:17)D0:90. Itisinterestingtonotethat,intheisotropiccase,themeanvalueof‘equals4V=S
foranyregionofvolumeV limitedbyaconvexsurface†ofareaS. Whentheregionisacube
withunitside,h‘iD4V=S D2=3.
Hence, (cid:17) is likely to have values close to 1 in the region where the radiation enters the
patient’sbodyandprogressivelydecreaseuntilnearlyreachingtheisotropicvalue0.90atthe
end of the electron range. Taking into account that target volumes cover only a fraction of
the region where radiation is transported, it is clear that this range of variation will usually
be further limited. Moreover, the decrease in (cid:17) may somewhat compensate for the increase
in S =(cid:26) at low energies, keeping the product even less sensitive to variations across voxels.
r
Alltheseconsiderationsjustifytheplausibilityoftheassumptionthat(cid:17) maybetreatedasa
constant.
Takingthisassumptionintoaccount,equation(10)canberecastas
C
(cid:27)2.x/Dx (12)
N
wheretheparameter
.S =(cid:26)/(cid:17)
C (cid:17) r (13)
1V2=3
isaconstantforagiventarget. Thisequationdemonstratesthattheestimateofthevariance
in dose in a voxel is approximately proportional to the dose itself, inversely proportional to
thenumberofhistoriesandinverselyproportionaltothevoxelsurface. Fromthecentrallimit
theorem(Lindeberg1922,Feller1967)itcanbeconcludedthat,whenthenumberofsimulated
historiesislargeenough,x hasaGaussiandistribution. Forthesakeofnotationalsimplicity,
thequantity(cid:27).x/willbedenotedby(cid:27).x/hereafter.
Itisimportanttonoticethatwhenelectrontransportisturnedoffsothatphotoninteractions
produceenergydepositionsonthespot,equation(9)lossesitsvalidityandtheexpressionfor
(cid:27).x/ is no longer applicable. This is due to the fact that, in this case, (cid:27).(cid:15)/ comprises the
fluctuationsof‘and thefluctuationsoftheenergydepositedbyphotonsalongagivenpath.
Thesecondcontributioncanbemadenegligibleforphotonsaswellbyusinginteractionforcing
orsomeequivalentmethodandthevalidityoftheformeranalysiscanberestored. Astudyof
theinfluenceofvariancereductiontechniquesonthepresentmodelwillnotbepursuedherein.
It suffices to say that equation (12) seems unlikely to hold in general when non-analogue
methods are applied unless all particles entering any voxel carry similar statistical weights,
whichwillthenappearasamultiplicativefactorinequation(12)modifyingthevalueofC.
† Everyelementofsurfaceareamustbeabletobejoinedtoanyothersurfaceelementbyanuninterruptedstraight
line.
136 JSempauandAFBielajew
0.05 0.05
N=105 N=106
0.04 0.04
V) V)
e 0.03 e 0.03
k k
g/ g/
DF ( 0.02 DF ( 0.02
P P
0.01 0.01
0 0
0 40 80 120 160 0 40 80 120 160
dose (keV/g) dose (keV/g)
0.05 0.05
N=107 N=108
0.04 0.04
V) V)
e 0.03 e 0.03
k k
g/ g/
DF ( 0.02 DF ( 0.02
P P
0.01 0.01
0 0
0 40 80 120 160 0 40 80 120 160
dose (keV/g) dose (keV/g)
Figure2. MCsimulatedDVHsfor10MeVelectronsimpingingnormallyonawaterphantom
formedby128(cid:2)128(cid:2)1281mm3 cubicvoxels. Theelectronfieldonthewatersurfacewas
collimatedtoa5(cid:2)5cm2square.Thetargetvolumewasdefinedasafield-centredparallelepiped
of4(cid:2)4(cid:2)2cm3,withtheshortestdistanceinthedirectionoftheincidentbeamandlocated
between2and4cmunderthewatersurface.Thedosebinwidthwasroughly1%ofthefullscale.
DifferentcasescorrespondtodifferentnumbersofhistoriesN.
2.2. DVHconvolution
The statistical ‘noise’ produced by the fluctuation of x about its mean can have a profound
effect on the DVH. An ideal calculation of the dose in each voxel would reproduce its true
value, say y, and the corresponding true DVH, which will be represented by p. The MC
simulation represents such an ideal only when the number of histories, N, tends to infinity.
ForfiniteN,thecalculateddosex fluctuatesaboutitsmeany withthevariance(cid:27)2 givenby
equation(12). Consequently,theMCcalculatedDVH,pQ,maydifferfromp.
Letusdefinearesponsefunctiong.xIy/suchthatg.xIy/dx givestheprobabilitythat
a voxel with a true dose y yields a measure in the interval [x;x +dx] for a given value of
N. TheMCcalculatedfractionofvoxelswithadoseinthissameintervalispQ.x/dx. This
fraction includes contributions from voxels that have a true dose y not in the interval under
consideration, but that have been associated with it by virtue of the statistical fluctuations.
Withthisinterpretationofg,pQ andparerelatedby
Z
pQ.x/D dyp.y/g.xIy/ (14)
D
withDrepresentingthedomainwherep.y/ 6D 0. Consideringg andpQ asknownquantities,
theformerisaFredholmequationofthefirstkind(ArfkenandWeber1995)intheunknownp.
Dosevolumehistogramsforradiotherapytreatmentplanning 137
N=105
0.02
V)
e
k
g/
(
DF 0.01
P
0
0 40 80 120 160
dose (keV/g)
Figure3. DVHsobtainedwithDPM(fullcurve)andPenelope(dashedcurve)for105histories.
Source,targetandphantomarethesameasinfigure2.
Alternatively,g.xIy/canbeviewedastheMCcalculatedpQobtainedwhentheunderlying
pcorrespondstoadeltafunctioncentredaty,thatis,whenthedosesdepositedinthedifferent
voxelsofthetargetvolumehaveexactlythesameexpectedvaluey. Accordingtoequation(12),
thisimpliesthatallthevoxelshavenearlythesamedepositeddosePDF,aGaussiancentred
aty withavariancegivenbythatequation. Therefore,fromthispointofviewitisclearthat
g.xIy/isalsorepresentedbythissameGaussian,namely
s
(cid:18) (cid:19) (cid:18) (cid:19)
1 .x−y/2 N N.x−y/2
g.xIy/D p exp − D exp − : (15)
(cid:27).y/ 2(cid:25) 2(cid:27)2.y/ 2(cid:25)Cy 2yC
AsN tendstoinfinity,(cid:27) tendstozero,g.xIy/collapsesintoadeltafunction,(cid:14).x−y/,and
pQ convergestop,whichwillbecalledtheconvergedDVHhereafter. Forvaluesofy much
largerthanC=N,g.xIy/canbeapproximatedbyafunctionof.x−y/onlyandtheFredholm
equation (14) becomes a familiar convolution integral. Although in practical situations this
limitmaynotbeattained,thetermconvolutionwillberetainedtodesignatetheprocessthat
transformspintopQ.
The effect of this convolution is to make pQ wider than p. Since the variance decreases
withincreasingN,thiswideningisexpectedtobesignificantatrelativelysmallvaluesofN.
Infigure2acomparisonofDVHsobtainedforthesameconditionsbutwithdifferentnumber
ofhistoriespresented. EvenforthesamevalueofN,differentmechanismsofsimulatingsoft
energylossescanproducedifferentvaluesofCandthusyieldapparentlydifferentDVHs. For
instance,inhomogeneousmedia,electronstepscanbetakenindependentlyofvoxelboundaries
andsoftenergylossesproducedalongagivensteplengthcanbedepositedatsomepointchosen
randomlyontheelectronpath. Thismeansthatalltheenergylostinastepisdepositedinonly
one voxel, irrespective of the number of them that have been traversed. (This is the current
schemeadoptedinthePenelopecode.) Alternatively,theparticlecanbestoppedeverytimeit
crossesavoxelboundary,inwhichcasesomeenergydepositionisobtainedeverytimeavoxel
ishit. (ThisisthecurrentschemeadoptedintheDPMcode.) Givenasufficientnumberof
histories,bothmethodsproducethecorrectmeanvalue,buttheirvariancesaredifferent. The
Penelope scheme produces a larger variance. Our analysis strictly applies only to the DPM
approach. However, when DVH convergence is employed to compare the reliability of MC
codesthesubtletiesmustbeproperlytakenintoaccount.
138 JSempauandAFBielajew
Figure3presentsanexampleofthisdanger. TwodifferentMCcodes,DPMandPenelope,
are used to calculate a DVH in water. These codes have shown very good agreement in
calculationsofdepthandlateraldoseinwaterphantoms. Surprisinglyenough,thetwoDVHs
present significant differences for the same number of histories. The reason is that, in this
particularcase,Penelope’senergy-depositionmodelislessefficientthanDPM’s.
Itisimportanttorealizethattheeffectoftheconvolutionisessentiallydifferentfromthe
customary MC uncertainty that appears in the form of error bars, e.g. in the calculation of
depthdoseprofiles. Inthiscase,uncertaintiesproduceonlyfluctuationsoftheMCcalculated
datapointsinthedirectionoftheordinateaxisandthereforetheycanbedenotedwithvertical
errorbars. Incontrast,theconvolutionofDVHsisaprocessinwhichtheuncertaintyofthe
doseinavoxelaffectstheabscissavaluethatwillbeassignedtoit. Hence,thisuncertaintyis
horizontal.
The convolution of DVHs is a phenomenon somewhat analogous to the response of a
scintillationdetector,wherethespectrumofthedepositedenergyisconvolvedwitharesponse
function, analogous to the action of the integral in equation (14), to produce the spectrum
actually measured. In this case, the fluctuations in the light yield and the response of the
photoelectrictubeareresponsibleforthenearlyGaussianshapeoftheresponsefunction,which
also has a width roughly proportional to the deposited energy. In contrast to the scintillator
case,however,thewidthoftheresponsefunctioncanbechangedintheprocessofcalculating
aDVHbychangingthenumberofhistoriessimulatedorthevoxeldimensions.
2.3. Testingthemodel
Inthissectionwepresentsomeresultsintendedtoclarifythemeaningofthemodelintroduced
aboveandverifyitsaccuracy. AllofthemhavebeenobtainedwithanewMCcodecalledDPM
(Sempauetal2000),anacronymforDosimetryPlanningMethod. DPMisbasedonthephysics
contained in EGS4 plus some newer developments intended to improve electron multiple
scattering and thus speed up the simulations considerably. It is optimized for radiotherapy
applications,andinthetestsperformedsofarithasbeenfoundtoagreetowithabout3%with
otherwell-knownstandardcodes. However,theinaccuraciesthatDPMoranyotherMCcode
may have are irrelevant here, for what is being investigated refers to the basic nature of the
MCmethod.
Inordertoverifytheaccuracyoftheapproximationsleadingtoequation(12),thevariance
ofthedepositeddose,x,hasbeencalculatedusingMCmethodsforeachvoxel. Astheabscissa,
ordoseaxis,x,oftheDVHisdividedintoanumberofequispacedbins,ameanvariancefor
each bin can be obtained by averaging this quantity over all voxels falling within the bin.
Fittingastraightlinetothedatapointsrepresenting(cid:27)2versusxyieldsavalueoftheconstant
C andthethevalidityofthelinearfitcanbejudged.
Thevarianceofeachindividualvoxel,(cid:27)2,canbeestimatedusing
(cid:20)P (cid:18)P (cid:19)(cid:21)
1 N x2 N x 2
(cid:27)2 ’ iD1 i − iD1 i (16)
N N N
wherex representsthedosescoredinthevoxelduringthesimulationoftheithhistoryand
i
N isthenumberofhistoriessimulated. Inprinciple,tocarryoutthesumsinthenumerators,
threedifferentcounters,sayxtmp,xscoreandx2scoremustbedefined. xtmpaccumulates
theindividualenergydepositionsproducedbyalltheparticlesofonlyPoneshower. Oncethe
shower is completed, its value is addePd to xscore, which contains xi, and its square is
added to x2score, which represents x2. Then xtmp is set to zero before a new shower
i
starts. Althoughthissimpleprocedurerepresentsaninsignificantcomputationalcostforthe
Dosevolumehistogramsforradiotherapytreatmentplanning 139
0.20
0.16
2
g)
V/
e 0.12
k
(
e
c
n 0.08
a
ari
v
0.04
0.00
30 40 50 60 70 80 90 100
dose (keV/g)

Figure4. Thevarianceindoseforthesamearrangementasinfigure2andN D 108. Dots
represent‘experimental’MCdataandthefulllineisalinearfitofthetype(cid:27)2 D xC=N. The
valueofthefittingparameterCandthecorrelationcoefficientare181:6MeVg−1andrD0:9998
respectively.
majorityofMCcalculations,thisisnottrueforradiotherapyapplicationssincethenumberof
countersisthreeforeachvoxel,andthetotalnumbermayamounttomanymillions. Accessing
thecomputermemoryaftereachhistorycanslowdowntheexecutionofthecoderendering
thecalculationinfeasible.
Toovercomethisdifficulty,weemployedthefollowingmethod†. Wedefineacounter,say
ihist,thatidentifiesthecurrentshowerhistory. Thenoveltyofthemethodconsistsinusing
an integer label, saylast, for each voxel in addition to the three counters described above.
lastidentifiesthelastshowerthatcontributedtothevoxelbystoringitsihistvalue. Every
timethereisadeposition(cid:15)inagivenvoxel,thecurrentihistiscomparedwithlast. Ifthey
match, thedepositionisaddedtothetemporarycounterxtmp. Iftheydiffer, xtmpisadded
toxscore,itssquareisaddedtox2score,thetemporarycounterxtmpissetto(cid:15) andlast
isupdatedtothecurrentihist. Thiscountingmechanismensuresthateverytimeashower
contributestoavoxel,thecontributionfromanyformershowerhasbeenstoredcorrectly. The
simulationtimeisnotseverelyaffected,becausethenewmethodinvolvesonlyanadditional
ifeverytimeadepositionismade,plusanadditionalarraytostoretheintegerlabellastfor
eachvoxel.
Acalculationofvarianceispresentedasafunctionofthedosex infigure4. Ourmodel,
equation(12),predictsalinearvariation,ingoodagreementwithwhatitisobserved.
SincetheGaussianshapeofgisrelatedtothevalidityofequation(12),whichstatesthat
allvoxelswiththesamedosehavethesamevariance,atestofthemodelconsistsincalculating
thedispersionofthisvariance. SupposethatK ofthevoxelsinthetargetvolumehappento
beinacertaindosebin. Consideringtheirvariances(cid:27)2withk D1:::K asasampleofsome
k
population,thequantity
(cid:20)P (cid:18)P (cid:19)(cid:21)
K (cid:27)4 (cid:27)2 2
s2.f(cid:27)2 g/D k − k (17)
kD1:::K K −1 K K
isameasureofthedispersioninthevariance.
† WeareindebtedtoProfessorFrancescSalvatoftheUniversityofBarcelonaforsupplyinguswiththeideaforthis
method.
140 JSempauandAFBielajew
10
8
105
)
%
( 6
2
s)/
2 4 106
s(
s
2 107
108
0
0 40 80 120 160
dose (keV/g)
Figure5.Relativedispersionofthevariancefordifferentnumberofhistories,N.
3
)
g
2/ 2
m
c
V
e
M
(
se 1
o
d
target
0
0 1 2 3 4 5 6
depth (cm)
Figure6. Theaxialdepth-doseprofileforthesameproblemconditionsasinfigure2. Thetwo
verticaldashedlinesindicatethelocationofthetargetvolumeandthehorizontaldashedlinesets
thelimitthatdividesthetargetareaintworegions,onewithaone-to-onecorrespondencebetween
depthanddose(belowthehorizontaldashedline)andtheotherwithatwo-to-onecorrespondence
(abovethehorizontaldashedline).
Figure 5 presents the dispersion of the variance as a function of the dose. The small
relative dispersion seen validates, at least for this particular arrangement, the premise upon
whichequation(15)wasderived.
NotethatasN increasestherelativedispersionstopsdecreasingindicatingthattheDVH
hasconvergedandhence,thefluctuationsbecometoosmalltomovevoxelstoneighbouring
dose bins. The bump observed at high doses in the cases 107 and 108 for a dose around
85 keV g−1 can also be understood as a consequence of the model. Figure 6 shows the
depth dose curve for the problem under consideration. For doses below a certain limit, the
correspondencebetweendoseanddepthissuchthattoeachdosetherecorrespondsonlyone
depth. Abovethislimit,however,everydosecorrespondstotwopossibledepths,onedeeper
thanthepointwherethemaximumoccursandtheotheratashallowerdepth. Thesetworegions