Calculating tumor trajectory and dose-of-the-day using cone-beam CT projections

1.7 MB·English

by  Bernard L. Jones

Calculating Tumor Trajectory and Dose-of-the-day Using Cone-Beam CT
Projections
Bernard L. Jones,1,a) David Westerly,1 and Moyed M. Miften1
Department of Radiation Oncology, University of Colorado School of Medicine
Purpose: Cone-beam CT (CBCT) projection images provide anatomical data in real-time over several
respiratory cycles, forming a comprehensive picture of tumor movement. We developed and validated a
method which uses these projections to determine the trajectory of and dose to highly mobile tumors during
each fraction of treatment.
Methods: CBCT images of a respiration phantom were acquired, the trajectory of which mimicked a
lung tumor with high amplitude (up to 2.5 cm) and hysteresis. A template-matching algorithm was used
5 to identify the location of a steel BB in each CBCT projection, and a Gaussian probability density function
1 for the absolute BB position was calculated which best fit the observed trajectory of the BB in the imager
0
geometry. Twomodificationsofthetrajectoryreconstructionwereinvestigated: first,usingrespiratoryphase
2
information to refine the trajectory estimation (Phase), and second, using the Monte Carlo (MC) method to

n sample the estimated Gaussian tumor position distribution. The accuracies of the proposed methods were
a
evaluated by comparing the known and calculated BB trajectories in phantom-simulated clinical scenarios
J
using abdominal tumor volumes.
4
Results: With all methods, the mean position of the BB was determined with accuracy better than 0.1
1
mm, androot-mean-square (RMS) trajectory errorsaveraged3.8±1.1%of themarkeramplitude. Dosimetric

] calculations using Phase methods were more accurate, with mean absolute error less than 0.5%, and with
h error less than 1% in the highest-noise trajectory. MC-based trajectories prevent the over-estimation of dose,
p
but when viewed in an absolute sense, add a small amount of dosimetric error (<0.1%).
-
d Conclusions: Marker trajectory and target dose-of-the-day were accurately calculated using CBCT pro-
e jections. Thistechniqueprovidesamethodtoevaluatehighly-mobiletumorsusingordinaryCBCTdata,and
m could facilitate better strategies to mitigate or compensate for motion during SBRT.
. This manuscript was submitted to Medical Physics
s
c
i
s I. INTRODUCTION using surrogates), and to treat only when the tumor lies
y
within some specified region (i.e. gating)6,7. More recent
h
p As radiotherapy becomes more conformal, the nega- advancesinvolvemodifyingorre-directingthetreatment
[ tive effects of tumor motion become more pronounced1. beam in real-time in order to compensate for intrafrac-

1 Highly targeted techniques, such as Stereotactic Body tion motion8,9.
v Radiotherapy (SBRT), rely on the establishment of a More accurate methods to account for tumor motion
6 very steep dose gradient around the tumor in order to are based on knowing the instantaneous tumor location
0 reduce normal tissue effects2. Motion blurs the bound- over time (i.e. the tumor trajectory). Several meth-
4 ary between tumor and normal tissue, making this dose ods exist to determine this. If the tumor is implanted
3
gradient harder to establish. Interfraction motion can withradiopaquefiducialmarkers,thepositioncanbetri-
0
largelybeaddressedbyimage-guidedtherapytechniques angulated using stereoscopic imaging from two or more
.
1 such as cone-beam CT (CBCT)3; however, intrafraction view angles10-13. Similarly, the tumor can be implanted
0
motion presents a more complicated problem. Specifi- with electromagnetic transponders7,14, or can be moni-
5
cally, motion makes it more difficult to ensure that the tored in the treatment suite using MRI images15. One
1
: targetreceivesthefullprescribeddose,andalsoaddsun- major drawback of these methods is that they rely on
v certainty when calculating the dose-of-the-day based on technologiesortreatmentsnotcommonlyavailabletothe
i
X pre-treatment imaging. majorityofradiotherapycenters,orareonlyavailablefor
A common approach to account for regular intrafrac- use in certain tumor sites.
r
a tionmotionistouse4-DimensionalCT(4DCT)imaging
CBCT is a widely used pre-treatment imaging modal-
to characterize the full range of tumor motion, and to
ity which captures a relatively low-quality 3D image of
treat a target volume which encompasses this range4,5.
the patient anatomy in the treatment position3. Typi-
This guarantees full coverage of the tumor but increases
cal clinical CBCT systems capture roughly 650 projec-
dose to normal tissue in the tumor vicinity. Another ap-
tionimages,whichareusedtoreconstructthe3Dimage.
proachistomonitorthetumorposition(eitherdirectlyor
These projections are typically acquired over a 30-120
second timescale, and the reconstructed images exhibit
significant blurring due to respiratory motion. However,
the projection images themselves provide a comprehen-
a)Author to whom correspondence should be addressed: sive picture of the anatomical movement, as they typ-
[email protected] ically observe the motion of several respiratory cycles.
2
For tumors implanted with fiducial markers, these pro- II. METHODS
jections can reveal the motion of the target volume dur-
ingtreatment. Unfortunately, theclinicalutilityofthese
To quantify the motion of the target, the location of a
projections is often limited as these data are incomplete
marker in 3D phantom geometry is calculated by exam-
with respect to the tumor position, since the process of
ining the position of that marker in a series of orbiting
projectioncondensesthelocationofanobjectin3Dspace
2Dprojectionimages. Ourmethodbuildsontheworkof
(x,y,z) to only two dimensions (u,v) in the imager ge-
Poulsenet al22 andproceedsasfollows. First,weacquire
ometry. Calculation of the tumor trajectory from these
CBCTprojectionsofaphantominwhichthe position of
projections would allow for determination of tumor mo-
a small fiducial marker is visible. Next, the position of
tionanddose-of-the-dayestimationfortumorsimplanted
thefiducialmarkerintheseprojectionsisdeterminedau-
with fiducial markers, which is often the case for tumors
tomatically using a template-matching algorithm. Then,
of the lung, liver, or pancreas16-18. Additionally, this
the acquired images are binned into different phases of
information could be used in the context of Adaptive
the respiratory cycle. An estimate of the likely distri-
Radiotherapy19-21.
bution of marker positions in each phase is built using
a Gaussian approximation of uncertainty22,25. Finally,
this Gaussian distribution is used to estimate the unre-
The goal of this work was to calculate the tumor tra- solved component of position, either by calculating the
jectory using CBCT projections, and to validate the ac- mostlikely position orbyusingMonteCarlomethodsto
curacy of this trajectory in calculating the tumor dose. sample this distribution.
Mathematically, determining the instantaneous 3D po-
sition of the tumor using the 2D projections acquired
during CBCT is an under-defined problem; in each pro-
jection we have two measured data points (u,v position A. Definition of geometry
within the imager geometry) yet wish to calculate three
quantities (x,y,z position of the tumor). The projection Let (x,y,z) represent an orthogonal basis in 3D space,
measurement confines the position to lie along a line be- withtheorigincorrespondingtothetreatmentandimag-
tweentheimagingsourceandthepointofdetection, but ing isocenter. In this work, we consider the patient to
cannotbyitselfdeterminethedepthofthepositionalong lieinthehead-firstsupineposition,andchoose(x,y,z) to
this line. Methods have been developed to estimate the correspondtotheanterior-posterior(AP)orvertical,left-
tumor position using orbiting CBCT projections, mostly right (LR) or lateral, and superior-inferior (SI) or longi-
in the context of tracking/targeting mobile tumors in tudinaldirections,respectively. ThekVsource/imageris
real-time22-24. These methods choose the tumor location oriented orthogonal to the z-axis, and orbits the patient
by maximizing the expectation value of tumor position inthex/y plane(rotationaboutthez-axis). θrepresents
(which maximizes the probability of hitting the tumor). the angle of the kV imaging source with respect to the x
However,thismaynotbeappropriatefordosimetriceval- axis. Let(u,v) describetheaxesofthe2Dprojectionim-
uation, since the most likely tumor location is generally age, where u is parallel to the x/y plane and v is parallel
in the center of the planning target volume (PTV) and to the z-axis. Note that this coordinate system matches
mayleadtoanover-estimationofthetumordose. More- those typically described in the radiotherapy literature
over, tumors which undergo large respiratory motion of- (e.g. the seminal work of Feldkamp et al26); however,
ten exhibit a cyclical trajectory. Previous methods are particularchoiceoftherelativedirectionsofx, y, z, u, v,
built around the assumption that tumor position is dis- andθ mayaltertherelationshipand/orsigndependence
tributedrandomlyaboutthedailymeanlocation,leading of these quantities.
to a Gaussian distribution of random position errors25.
However, this assumption is broken in the presence of
non-random respiratory motion, and a Gaussian distri-
B. Projection operator
bution of these positions will contain relatively large un-
certainty. For cyclical motion, the true uncertainty is
much smaller, since there is a large correlation between LetPdefinetheprojectionofaposition(x,y,z)intothe
position and respiratory phase. (u,v) plane at an angle θ; in other words, P represents
theacquisitionofaCBCTprojection. SDDisthesource-
to-detector distance, p is the image pixel size, SAD is
the source-to-axis distance, and o ,o are the projected
In this work, we extended these methods for more ac- u v
locationsoftherotationalisocenterintheu,v coordinates
curate calculation of tumor dosimetry in highly-mobile
(the so-called piercing-point or principal point). In Eqs
tumors by incorporating information regarding the res-
1-3,u andv aregiveninunitsofpixelsinimagingspace.
piratory phase. Using a motion phantom, we validated
theaccuracyofthesemethods,andquantifiedtheuncer-
tainty in using CBCT projections to calculate dose-of-
the-day in mobile tumors. (u,v)=P(x,y,z,θ) (1)
3
(cid:115)
detB
K = e−(p(cid:126)+eˆµ−r(cid:126)0)TB(p(cid:126)+eˆµ−r(cid:126)0)/2 (7)
SDD xsin(θ)+ycos(θ) (2π)3
u= +o (2)
p SAD−[xcos(θ)−ysin(θ)] u
g(t)=Ke−(t−mu)2/(2σ)2 (8)
SDD z
v = +o (3)
p SAD−[xcos(θ)−ysin(θ)] v Here, eˆ describes a unit vector pointing from the ob-
served marker position to the imaging source, and p(cid:126) de-
Consider the following clinical scenario. A patient is
scribes the point on the imager at which the marker was
implantedwithseveralradiopaquefiducialmarkersinor-
detected (converted to patient geometry).
der to assist in the delineation and localization of this
Usingtheobservedmarkerpositions(u ,v )ineachim-
i i
tumorfortreatment. ACBCTisacquiredinordertovi-
age i, optimization was performed to estimate the mean
sualizethepatientanatomyonthedayoftreatment,and
positionandcovariancematrixwhichdescribestheprob-
involves the acquisition of several hundred images which
abilistic marker positions. By integrating Eq. 8, it can
contain the projection of that anatomy into 2D. The ac-
be shown that the probability of finding the marker at
quisitionofthisCBCTscanisrepresentedbyperforming
a given detector position is linearly proportional to the
P(θ) for values of θ from 0◦ to 360◦, and for each fidu-
product K ·σ. The optimization finds a mean position
cial marker, each image contains a measurement of the
and covariance matrix which maximizes the total prob-
marker position u and v .
θ θ ability of finding the markers in positions observed in
each image i. The optimization was performed using the
Nelder-Mead simplex search method27.
C. Marker position
(cid:88)
WeusedtheformalismofPoulsenetal22toconstructa [r(cid:126),B]=argmin −log(K σ ) (9)
0 i i
3DGaussianprobabilitydensityfunction(PDF)describ-
r(cid:126)0,B i
ing the marker position. A full derivation of this Gaus-
sian formalism is given in Ref 22, and for convenience of
thereader,weusethesamenotationwhenpossible. This
method uses maximum-likelihood optimization to calcu- D. Respiratory phase
late a PDF which best describes the observed motion of
a marker in projection images. Then, using this PDF,
Respiratory phase was determined using the SI dis-
the unresolved component of marker position is deter-
placement of the fiducial markers. The SI position is
mined by computing the expectation value of this PDF
well-determined in the CBCT projections for locations
alongthelinebetweenthex-raysourceandthedetection
near the rotational isocenter, since the only error derives
point within the image (Fig 1).
frommagnificationeffectsrelatedto(x,y)errors(Eq. 3),
A 3D Gaussian PDF is defined by 9 quantities: mean
which are typically on the order of 1% or less. This in-
(cid:126)
location (r0)=(x0,y0,z0), variance (σx,σy,σz), and co- ternalsurrogateinformationprecludedtheneedtodefine
variance (σxy,σyz,σzx). The variance and covariance to- respiratoryphasebasedonsomeothersurrogate,suchas
gether form the covariance matrix A. From this, the 3D chest markers or diaphragm movement.
(cid:126) First, the marker trajectory in the v direction was
PDF for marker position (r)=(x ,y ,z ) was given by:
0 0 0
smoothed using a five-sample moving average, and dis-
crete velocity was calculated by taking the difference be-
(cid:115)
tween smoothed positions in sequential projections. It
detB
G(x,y,z)= e−((cid:126)r−r(cid:126)0)TB((cid:126)r−r(cid:126)0)/2 (4) is assumed that these marker positions represent tumor
(2π)3
trajectoryinahead-firstsupineorientation, withthe+v
axis pointing in the superior direction. End-exhale (v )
Here, B denotes the inverse covariance matrix B = ee
wasidentifiedfromprojectionswherethemarkervelocity
−1
A , and T denotes the transpose operator. The func- switched from positive to negative; likewise, end-inhale
tional form of Eq. 4 along the line from the imaging (v ) was identified where velocity went from negative to
ei
source to detection point (defined here as g) was de- positive.
scribed by the following: Projections were binned into one of 10 respiratory
phases based on the velocity and location of the marker
(cid:113) relative to the end-exhale and end-inhale positions.
σ = eˆTBeˆ (5) Phases 1-5 (inhale) consisted of projections with nega-
tivevelocity,groupedintouniformlyspacedbinsbetween
subsequentv andv coordinates. Phases6-10(exhale)
ee ei
µ= (r(cid:126)0−p(cid:126))TBeˆ (6) consisted of projections with positive velocity, grouped
into linearly spaced bins between subsequent v and v
eˆTBeˆ ei ee
4
2DvProjection 2DvProjection
3DvGaussian DistributionsvofvTumorv
Distributionvof Positionvforveachv
TumorvPosition RespiratoryvPhase
MeasuredvPosition MeasuredvPosition
(u,v) (u,v)
ExpectedvPosition ExpectedvPosition
(PhasevResolved)
FIG.1. Determinationoftumorpositionusing3DGaussianprobability. A)Usingall2Dprojections,a3DGaussiandistribution
isfoundwhichbestfitstheobservedmarkertrajectory. Foreachmeasuredposition(i.e. eachprojectionimage),thepositionis
determined by finding the expectation of this Gaussian distribution along the line between the x-ray source and the detection
point. B) In the phase-resolved method, a 3D distribution is independently calculated for each phase of the respiratory cycle.
coordinates. Projectionsatthebeginning(orend)ofthe tween 0 and 1. The position of the marker along the ray
motion trace (i.e. those that did not contain measure- was chosen to be the point T which satisfied CDF(T)=
ments of the full cycle) used v (or v ) coordinates of ξ, and this solution was found numerically.
ei ee
the nearest full phase. To increase the accuracy of tracking for targets with
high respiratory motion, the trajectory of the marker
was calculated using a phase-resolved Gaussian approxi-
E. Position determination mation (Method 3, Phase). All projections were binned
into one of 10 respiratory phases, and Eq. 9 (optimiza-
Four methods were used to determine the 3D marker tion of Gaussian distribution) was applied to data from
trajectory from the projection images. These meth- each phase independently. This resulted in a set of 10
odswereimplementedusingin-housesoftwaredeveloped phase-specific mean positions r(cid:126) and inverse covariance
p
with Matlab (The Mathworks, Natick MA). In Method matrices B , or in other words, 10 distinct distributions
p
1 (Base), position estimation followed the method of describing the position of the marker during each phase
Poulsen et al22. Data from all projections was used dur- of respiration (Fig 1B). In Method 4 (PhaseMC), both
ing optimization to calculate a 3D Gaussian distribution thePhaseandMCmethodswereappliedsimultaneously.
which described the probability of finding the marker in
a given 3D location (Eq. 9). This equation was parame-
terized along the detection ray (the line between the ra-
F. Motion phantom
diation source and the imager detection point), and the
markerpositionalongthislinewasgivenbyt=µ,where
In order to validate the proposed methods, a phantom
µ equals the expectation value of the distance along this
was used to simulate respiratory tumor motion (Quasar
line (Eqs. 5-8).
Phantom, Modus Medical Devices, Inc.). This phantom
In Method 2 (MC), Monte Carlo methods were used
contains a 1-mm-diameter steel BB, the trajectory of
to select the marker position in 3D space in place of the
whichisshowninFig2forvarioussuperior-inferior(Sup-
expectation value µ. First, the detection ray was dis-
Inf) amplitudes of motion. The trajectory is similar to
cretized into 10,000 points, spaced linearly between the
that of a lung tumor trajectory with high amplitude and
source and detector. The PDF of the 1D Gaussian (Eq.
hysteresis11. ACBCTscanofthisphantomwasacquired
8) was integrated across the discretized ray to form a
using the imaging capabilities of a Varian Truebeam ac-
normalized cumulative probability distribution function.
celerator (Varian Medical Systems, Palo Alto, CA). 656
projection images were acquired using a half-fan trajec-
(cid:80)i=Ig(t ) tory (θ = −180◦ → 180◦). The phantom was set to a
CDF(t )= i=1 i (10)
i (cid:80)i=imaxg(t ) rate of 11 breaths-per-minute, with a maximum range
i=1 i in the SI direction of 2.5 cm (resulting in maximum AP
Next, a random number ξ was uniformly chosen be- andLRrangesof1.2and0.6cm, respectively). Multiple
5
Amplitude
5 mm
10 mm
15 mm
20 mm
25 mm
0
-0.5
)
m
c
(

st -1 1.5
o
P 1
-
nt 0.5
A
-1.5
0
-0.5 Sup-Inf (cm)
-1
-2
2.5 2 1.5 -1.5
Left-Right (cm)
FIG. 2. Trajectory of Motion Phantom. The motion of the 1-mm-diameter steel BB embedded within the motion phantom
is shown relative to the isocenter of the CBCT imaging system. Trajectories are shown for amplitudes of motion in the Sup-
Inf direction ranging from 5 25 mm. For clarity, the trajectories are also shown projected onto the XY, XZ, and YZ axes.
Trajectories are shown relative to a standard head-first supine orientation (i.e. the primary axis of motion in the Sup-Inf
direction).
projection image datasets were acquired with trajectory G. Dosimetric model
amplitudes ranging in 1 mm increments from 1 mm to
2.5 cm.
The goal of this study was to evaluate the accuracy
In each image, the location of the BB was determined of dosimetric calculations based on reconstructed trajec-
automatically using a template matching algorithm28. A tories. To that end, the dosimetry of the reconstructed
discretized 2D Gaussian kernel was used as a template trajectories was compared to the dosimetry of the true
image, with width equal to the radius of the BB. In each trajectory. To evaluate the dosimetry of the true and
image, the normalized cross correlation was computed reconstructed phantom motion trajectories, a phantom-
between the projection image and the template29. The simulatedclinicalscenariowasconsidered. Aclinicaltar-
locationoftheBBintheimagewasgivenbythelocation get volume (CTV) was assumed to undergo the motion
where the normalized cross correlation was maximized. trajectories under analysis. Treatment was prescribed
In cases where the template matching algorithm failed, to a planning target volume (PTV), which included the
the position was determined manually. CTV plus an isotropic margin. To simplify the model,
UsingthepositionsoftheBBintheprojectionimages, andtoexaggeratedifferencesbetweendifferentmethods,
thetrajectoryoftheBBin3Dspacewascalculatedusing the dose within the PTV was assumed to equal the pre-
methods 1-4 (Sect 2.5). The accuracy of these methods scription,withzerodoseoutsidethePTV.ThePTV(and
wascalculatedbycomparingthecalculatedtrajectoryto resulting dose distribution) was stationary over time,
the true trajectory (shown in Fig 2). To simulate more while the CTV was translated within the PTV based on
difficult and unpredictable clinical scenarios, additional either the true phantom trajectory or the reconstructed
(simulated)trajectorieswerecreated. Gaussiannoisewas trajectories. For each of the 656 points in the trajectory
added to the true marker trajectory by adding normally (corresponding to the 656 projection images), the posi-
distributed random values to the x,y, and z locations at tion of the CTV relative to its center was given by the
each point in time. Simulated trajectories were created instantaneous trajectory position relative to the mean
by varying the standard deviation of this noise from 0 position. Basedon this position, the dosewascalculated
to 1 mm in 0.1 mm increments between datasets. These for each voxel in the CTV (i.e. based on the binary dose
altereddatasetswereprojectedinto2Dimagespace(Eqs. model, or in other words assigning 1/656th of the pre-
1-3), and Methods 1-4 were used to calculate the marker scription dose for voxels inside the PTV, and zero for
trajectory as described previously. those outside). Independent simulations were performed
6
2.5 10
A- B- C-
h1 Base
10
m- MC
m 2 m 8 d Phase
ErrorvO1.5 vev3m 6 bserve10h2 PhaseMC
oryv Abo yvO
MSvTraject0.51 FvErrorsv 24 Frequenc10h3
R
0 0 10h4
0 5 10 15 20 25 0 5 10 15 20 25 0 0.2 0.4 0.6 0.8 1
TrajectoryvAmplitudevOmm- TrajectoryvAmplitudevOmm- TrajectoryvErrorvOcm-
FIG. 3. Errors in calculating trajectory using orbiting views. Results were calculated with no corrections (Base), with phase
corrections (Phase), with Monte Carlo sampling (MC) and with both phase and Monte Carlo (PhaseMC). Results are shown
for phantom motion amplitudes in the SI direction ranging from 1 mm to 25 mm. A) Root-mean-squared (RMS) errors in
the trajectory reconstruction using different methods. Phase and PhaseMC methods reduce RMS errors by half. B) Relative
fraction of trajectory errors greater than 3 mm. Phase methods reduce the fraction of large errors, and Monte Carlo methods
increasethepercentageoflargeerrorsslightly. C)Histogramoftrajectoryerrorsinallprojectionsacrossallacquireddatasets.
Phase methods demonstrate increased accuracy.
for CTV-to-PTV margins ranging from 0 mm to 15 mm, trajectories. Themeandoseerrorwascalculatedforeach
and were deliberately set less than the amplitude of mo- volume as the mean dose difference between calculated
tion so as to allow for differences between the proposed dose and reference dose. The mean absolute dose er-
methods to become apparent. The goal of our analysis rorwascalculatedasthemeanofthemagnitudeofthese
was to quantify the accuracy of dosimetric calculations dosedifferencesacrossallvoxelofthevolume. Meandose
using CBCT-based trajectory reconstruction. To that error gives a sense of the dosimetric trends of each tra-
end, our binary dose model serves as a conservative esti- jectory reconstruction method, and can reveal whether
mateoftheinaccuraciesofthisapproach,asanyerrorsin a given method underestimates or overestimates dose.
position determination will be reflected in the idealized Meanabsolutedoseerrordenotestheoverallaccuracyof
nature of the target dosimetry simulation. a method in determining accurate dose in each voxel (by
preventing positive and negative errors from cancelling).
In addition to these metrics, dose-volume histograms of
H. Validation metrics each volume were computed for each method, and com-
pared to the reference.
We analyzed the accuracy of the proposed methods
(Base, MC, Phase, PhaseMC) by comparing these tra-
III. RESULTS
jectories to the true phantom motion trajectory. Each
CBCT projection dataset (656 images) constituted 656
measurements of the marker position in phantom geom- A. Accuracy of trajectory reconstruction
etry. For each measurement, the error was given as the
magnitude of the 3D vector difference between the cal- The accuracy of the trajectory calculations are shown
culated position and the true position. The root-mean- in Fig. 3; results were calculated without respiratory
squared trajectory error was given for each trajectory as phaseorMonteCarlocorrections(Base),withphasecor-
the square root of the mean of the squared errors for all rectionsonly(Phase),withMonteCarlocorrectionsonly
measurements in that trajectory. Additionally, the fre- (MC),andwithbothphaseandMonteCarlocorrections
quency of large errors was computed by calculating the (PhaseMC). Trajectory amplitude ranged from 0 to 2.5
fraction of measurements with errors exceeding 3 mm. cm.
Theaccuracyofdosimetriccalculationswascomputed Withallmethods,themeanpositionofthemarkerwas
using our phantom-simulated clinical scenario. Since the determined with accuracy better than 0.1 mm. In Fig.
dosimetric results depend on the shape of the target, 32 3a, itcanbeseenthatthePhaseandPhaseMCmethods
abdominal/lungtumorvolumeswereusedtoevaluatethe reduced the root-mean-squared (RMS) trajectory errors
accuracy of the proposed methods (mean volume, 25 ± by 50% compared to the Base and MC methods, with
23 cm3; range 2.0 - 80 cm3). The reference dose was cal- an average RMS error of 3.8 ± 1.1% of the trajectory
culated by performing this analysis using the exact tra- amplitude. In Fig 3b, the Phase and PhaseMC meth-
jectory, and the calculated dose used the reconstructed ods decreased the fraction of large errors (defined here
7
xF 1 AF 1 CF 1 EF
a
M
4P 0.5 0.5 0.5
orPN
Err 0 0 0
eP
s
o
D -0.5 -0.5 -0.5
nP
a
e
M -1 -1 -1
0 5 10 15 20 25 0 5 10 15 0 0.2 0.4 0.6 0.8 1
xF
Ma 2 BF 2 DF 2 FF
4P Base
orPN 1.5 MPhCase 1.5 1.5
Err PhaseMC
eP 1 1 1
s
o
D
s.P 0.5 0.5 0.5
b
A
nP
ea 0 0 0
M 0 5 10 15 20 25 0 5 10 15 0 0.2 0.4 0.6 0.8 1
MotionPAmplitudePNmmF PTVPMarginPNmmF TumorPPositionPNoisePNmmF
FIG.4. Meanandmeanabsolutedoseerrorinthereconstructedtrajectories. Errorswerecalculatedusing32clinically-drawn
abdominal/lungCTVs. Motionamplituderangedfrom0to2.5cm,andresultswerecalculatedwitharangeofisotropicmargins
and for trajectories with added noise. Dosimetric calculations using phase-resolved methods were more accurate, with mean
absolute error less than 0.5% with zero margin, and with error less than 1% in the high-noise calculation. As noise increases,
expectation-based methods generate a trajectory which over-estimates the dose to the CTV. Monte Carlo-based trajectories
under-estimatedosebyanaverageof0.2%ofprescriptiondose,evenastrajectorynoiseincreases. InA)andB),thenoiseis0
mm and the margin is 4 mm. In C) and D), the amplitude is 2.5 mm and the noise is 0 mm. In E) and F), the amplitude is
1.2 mm and the margin is 4 mm.
as trajectory error > 3 mm) by an average of 80%. MC tory estimation methods. This is likely due to the lower
methods increased both the RMS trajectory error and RMS trajectory error demonstrated by these methods.
percentage of large errors slightly. Fig. 3c shows his- As gaussian trajectory noise increased, the expectation-
tograms of trajectory errors in all projections across all based methods Base and Phase tended to overestimate
acquired datasets (amplitudes ranging from 1 mm to 25 dose, while Monte Carlo methods MC and PhaseMC
mmin1mmincrements). Trajectoryerrorwaslessthan tendedtounderestimateit. WhileMonteCarlomethods
1mm in 92% of projections (Phase and PhaseMC) and prevented the overestimation of dose, they also slightly
82% of projections (Base). increased the overall absolute error (<0.1%).
Fig. 5 shows the accuracy of each method in de-
terming dosimetric parameters D90 and D95 for vary-
B. Dosimetric Accuracy ing amplitudes, PTV margin sizes, and trajectory noise.
One can see that the error of Phase methods was lower
Fig. 4 shows the distribution of mean dose errors thantheerrorinmethodswhichignorerespiratoryphase.
and mean absolute dose errors using trajectories recon- Expectation-basedmethodsagainoverestimatedthedose
structedwiththeBase,MC,Phase,andPhaseMCmeth- (relativetothetruevalue)andMonte-Carlomethodsun-
ods. Results are shown across variations in trajectory derestimated it, but the Base and Phase methods had a
amplitude, PTV margin size, and amount of added tra- lower absolute error. The greatest accuracy was seen in
jectory noise. thePhasemethod,whichwasabletodetermineD90and
D95 of the CTV to within 0.5% relative to prescription.
Doseerrorsincreasedsharplyasthemotionamplitude
increased,althoughinthePhaseandPhaseMCmethods,
absolute error never exceeded 0.5% of the prescription
dose. Ingeneral,doseerrorsdecreasedasthemarginsize IV. DISCUSSION
increased, mainly due to an overall reduction in the con-
tribution of motion seen by having more margin around In this work, we have developed and tested methods
the CTV. However, it can also be seen that mean dose to evaluate dose to a moving tumor using data from a
error was lower when using Phase and PhaseMC trajec- standard clinical cone-beam CT. In particular, we have
8
2 2 2
eh Ah Ch Eh
m
u 1 1 1
ol
V
N. 0 0 0
0.F
9-1 -1 -1
D
n.
or.i-2 -2 -2
Err
-3 -3 -3
0 5 10 15 20 25 0 5 10 15 0 0.2 0.4 0.6 0.8 1
Base
eh 2 Bh 2 MC Dh 2 Fh
m Phase
u 1 1 1
ol PhaseMC
V
N. 0 0 0
5.F
9-1 -1 -1
D
n.
or.i-2 -2 -2
Err
-3 -3 -3
0 5 10 15 20 25 0 5 10 15 0 0.2 0.4 0.6 0.8 1
Motion.Amplitude.Fmmh PTV.Margin.Fmmh Tumor.Position.Noise.Fmmh
FIG. 5. Accuracy of calculating clinical dose metrics D90 and D95 from reconstructed trajectories. Errors are calculated in
32 clinically-drawn CTVs. Motion amplitude ranged from 0 to 2.5 cm, and results were calculated with a range of isotropic
marginsandfortrajectorieswithaddednoise. Phase-resolvedmethodsincreasedtheaccuracyofthesedosemetriccalculations.
Astrajectorynoiseincreased,theBaseandPhasemethodstendedtooverestimateD90andD95,whileMCmethodsprevented
overestimation of dose. In A) and B), the noise is 0 mm and the margin is 4 mm. In C) and D), the amplitude is 2.5 mm and
the noise is 0 mm. In E) and F), the amplitude is 1.2 mm and the margin is 4 mm.
tailoredthesemethodstotheproblemofcalculatingdose limitations; it cannot be used for tumors which undergo
tohighlymobiletumorsimplantedwithfiducialmarkers, non-cyclical motion, such as prostate tumors. It is also
such as lesions in the lung, pancreas, or liver16-18. The difficult to determine the respiratory phase in real-time.
trajectory and subsequent dose received by the target
Given the prevalence of CBCT, this method may pro-
volume was determined with minimal error. Errors in
vide a useful technique to clinically estimate the dose-
mean position did not exceed 0.1 mm, and the RMS tra-
of-the-day received by a mobile tumor, since over half of
jectory error averaged 3.8% of the total range of motion.
radiation oncology clinics perform some cone-beam CT
Using the Phase method, the mean dose and mean ab-
imaging for localization30. Our analysis is directed to-
solute dose were calculated with over 99% accuracy, and
wards tumors implanted with fiducial markers, but also
D90/D95 were determined to within 0.5% of the true
holds in cases without fiducials as long as the tumor lo-
value.
cation can be discerned in each projection (as has been
By incorporating information related to the respira- demonstratedinlungtumors31,32). Intheory,ourmethod
tory cycle, the accuracy of tumor trajectory reconstruc- allows for very accurate trajectory/dose reconstruction,
tion and subsequent dosimetric analysis was improved. sincedataareobtainedusingamuchmorecomprehensive
The 3D Gaussian PDF utilized in this study stems from motion profile than other techniques. A CBCT projec-
the work of Papiez and Langer25, where it was assumed tion dataset constitutes upwards of 600 samples of tu-
that any systematic position errors are corrected with mor position, whereas 4D CT typically uses only 10 res-
dailyimageguidance,andthatonlyrandomfluctuations piratory phases, and has been shown to not accurately
remain. Thusthecentrallimittheoremdictatesthatthe represent daily intrafraction motion of the abdomen33,34.
distribution of the remaining random errors is Gaussian. Asradiotherapybecomesmoreconformal, andhypofrac-
However,thisassumptionisbrokenfortumorswhichun- tionationbecomesmorecommonplace,itbecomescrucial
dergo regular respiratory motion, and the distribution to understand and account for the effects of tumor mo-
of the trajectory about the mean position is no longer tion on dosimetry. This method could be used to design
well definied as Gaussian. By creating a unique Gaus- patient-specific, optimized margins for treatment, or to
sian distribution for each portion of the respiratory cy- verifytheadequacyofthemarginsused. Itcouldalsobe
cle, this effect is minimized, and the accuracy of trajec- used to decide between motion management strategies
tory reconstruction is improved. This method does have on a patient-specific basis (such as breath-hold, gating,
9
or free-breathing). It should be said that while we have tion ray. Yet one expects this determination to exhibit
validatedtheaccuracyofourmethodinasimpleclinical a Gaussian error profile, as the tumor position under-
scenario, true dose-of-the-day reconstruction would re- goes random fluctuations in addition to the cyclic respi-
quireamoresophisticateddosimetricapproach. Bycon- ratory motion. MC sampling addresses this by sampling
structing a simplified clinical scenario which maximizes a greater portion of the tails of this distribution (rela-
the impact of tumor motion, our results indicate that tive to using the expectation value), thereby generating
tumor trajectories can be derived from CBCT projec- a distribution of positions with a more realistic propor-
tion datasets with sufficient accuracy to determine true tionoflargedeviations. OnemayexpecttheMCmethod
dose. An additional benefit of this method is that the to increase trajectory error; however, under certain con-
processing time is short. Template-matching to deter- ditions one would expect a more complete sampling of
mineseedlocationscanbeperformedinreal-time,andit the error distribution to result in a more accurate calcu-
takes roughly 2 seconds to calculate the trajectory once lation of derived metrics (such as target dosimetry). As
the projected fiducial locations are known. However, to thetrajectorynoiseincreases, theBaseandPhasemeth-
build the 3D PDF, one must know the entire set of pro- ods to overestimate the true tumor dose, as the expecta-
jectedlocations;inotherwords,thePDFcannotbebuilt tion value of the distribution likely lies within the PTV.
until after CBCT acquisition is completed. Our method MCsamplingslightlyincreaseserrorinanabsolutesense
could be applied to generate an accurate estimate of tu- (Fig. 4),butisaconservativeestimateofD90/D95target
mordosimetryshortlyafterthecompletionoftheCBCT dose metrics in noisy trajectories (Fig. 5). MC methods
acquisition. maybeusefulforcalculatingdosetotheprescriptionvol-
ume,andcouldconstituteaworst-casescenarioestimate
The conditions of this study were chosen to test the
of tumor dose. However, for organs-at-risk and other
boundaries of dose reconstruction in highly mobile tu-
volumes,onewouldlikelyuseexpectation-basedposition
mors. Ranging up to 2.5 cm, the amplitude of motion
estimation (i.e. Base or Phase), which yields a more ac-
in the SI direction corresponded to the displacement of
curate result.
very highly mobile lower lobe lung tumors, and exceeds
the range of liver and other abdominal tumors11,35. Hys-
teresis was also exaggerated, and motion in the AP and
LR directions (up to 1.2 and 0.6 cm) exceeded that seen
in clinical cases11,35. In addition to the large range of
V. CONCLUSIONS
motion, the effects of that motion on the dose were ex-
aggerated by assuming perfect dose fall-off outside the
PTV. Our methods were also tested over a wide range Wehavedevelopedmethodsforaccuratelydetermining
of tumor volumes, and with up to 1 mm of Gaussian the trajectory of highly-mobile abdominal tumors using
noise added to each point in the trajectory. Even un- cone-beam CT projections. These methods have been
der these circumstances, dose metrics of the tumor were applied to reconstruct the 3D trajectory of a respiratory
determinedaccurately. OnemaynotethattheRMStra- motion phantom, and to calculate the dose received un-
jectory error reported in this work (1-2 mm) exceeds the der that trajectory. Results were calculated without res-
mean error in the similar work of Poulsen et al36. This piratory phase or Monte Carlo corrections (Base), with
canbeexplainedbythelargeamplitudeofmotioninthe phasecorrectionsonly(Phase),withMonteCarlocorrec-
current study, and agrees well with the maximum RMS tions only (MC), and with both phase and Monte Carlo
error (> 2 mm) reported for lung tumors in that work. corrections (PhaseMC). In all methods, errors in mean
Additionally, these results give insight as to the scale of positiondidnotexceed0.1mm, andtheRMStrajectory
motion which causes the greatest dosimetric inaccuracy error was only 3.8% of the total range of motion. Using
foreachmethod. InFig. 4C,errorsincreasefortheBase phase-resolved methods Phase and PhaseMC, the mean
andMCmethodsasthemaginsizeapproachestherange dose and mean absolute dose were calculated with over
of 5-8 mm, while the Phase and PhaseMC see the great- 99% accuracy, and D90/D95 were determined to within
esterrorfor4-5mmmargins. Sinceourdosimetricmodel 0.5% relative to the prescription dose. These methods
is sensitive to motion equal to or greater than the PTV canbeappliedtoestimatethedose-of-the-dayformobile
margin, andthefrequencyofmotionerrorsdecreasesex- tumors, and can aid in the selection of motion manage-
ponentially for larger errors (Fig 3C), these results sug- ment strategies.
gest that motion on the order of 7 mm has the grestest
effectfortheBaseandMCmethods, andontheorderof
4 mm for the Phase and PhaseMC methods.
Itisclearfromtheseresultsthatincorporationofrespi-
VI. ACKNOWLEDGEMENTS
ratory phase information is beneficial when reconstruct-
ing the tumor trajectory of highly mobile tumors. How-
ever, the role of MC sampling is mixed. In the Base The authors would like to thank Hirotatsu Armstrong
scenario, the position of the tumor is determined as the for several helpful discussions and comments regarding
expectation value of the Gaussian PDF along the detec- this work.
10
REFERENCES tional Journal of Radiation Oncology* Biology* Physics 82,
1665-1673 (2012).
14. A. P. Shah, P. A. Kupelian, B. J. Waghorn, T. R.
1. P. J. Keall, G. S. Mageras, J. M. Balter, R. S. Emery,
Willoughby, J. M. Rineer, R. R. Maon, M. A. Vollenweider
K.M.Forster,S.B.Jiang,J.M.Kapatoes,D.A.Low,M.J.
andS.L.Meeks,”Real-TimeTumorTrackingintheLungUs-
Murphy and B. R. Murray, ”The management of respiratory
inganElectromagneticTrackingSystem,”InternationalJour-
motion in radiation oncology report of AAPM Task Group
nal of Radiation Oncology* Biology* Physics (2013).
76,” Medical Physics 33, 3874 (2006).
15. J. Yun, K. Wachowicz, M. Mackenzie, S. Rathee, D.
2. S.H.Benedict,K.M.Yenice,D.Followill,J.M.Galvin,
Robinson and B. Fallone, ”First demonstration of intrafrac-
W.Hinson,B.Kavanagh,P.Keall,M.Lovelock,S.Meeksand
tional tumor-tracked irradiation using 2D phantom MR im-
L. Papiez, ”Stereotactic body radiation therapy: the report
ages on a prototype linac-MR,” Medical Physics 40, 051718
ofAAPMTaskGroup101,”MedicalPhysics37,4078(2010).
(2013).
3. D. A. Jaffray, J. H. Siewerdsen, J. W. Wong and A.
16. W.G.Park,B.M.Yan,D.Schellenberg,J.Kim,D.T.
A. Martinez, ”Flat-panel cone-beam computed tomography
Chang,A.Koong,C.PatalanoandJ.VanDam,”EUS-guided
for image-guided radiation therapy,” Int J Rad Oncol* Biol*
gold fiducial insertion for image-guided radiation therapy of
Phys 53, 1337-1349 (2002).
pancreatic cancer: 50 successful cases without fluoroscopy,”
4. T. Pan, T.-Y. Lee, E. Rietzel and G. T. Chen, ”4D-
Gastrointestinal endoscopy 71, 513-518 (2010).
CT imaging of a volume influenced by respiratory motion on
17. W.Wunderink,A.MndezRomero,W.DeKruijf,H.De
multi-slice CT,” Medical Physics 31, 333 (2004).
Boer,P.LevendagandB.Heijmen,”Reductionofrespiratory
5. H. H. Liu, P. Balter, T. Tutt, B. Choi, J. Zhang, C.
livertumormotionbyabdominalcompressioninstereotactic
Wang, M. Chi, D. Luo, T. Pan and S. Hunjan, ”Assessing
bodyframe,analyzedbytrackingfiducialmarkersimplanted
respiration-inducedtumormotionandinternaltargetvolume
in liver,” International Journal of Radiation Oncology* Biol-
using four-dimensional computed tomography for radiother-
ogy* Physics 71, 907-915 (2008).
apy of lung cancer,” International Journal of Radiation On-
18. M. Imura, K. Yamazaki, H. Shirato, R. Onimaru, M.
cology* Biology* Physics 68, 531-540 (2007).
Fujino, S. Shimizu, T. Harada, S. Ogura, H. Dosaka-Akita
6. C. R. Ramsey, D. Scaperoth, D. Arwood and A. L.
and K. Miyasaka, ”Insertion and fixation of fiducial markers
Oliver, ”Clinical efficacy of respiratory gated conformal radi-
for setup and tracking of lung tumors in radiotherapy,” In-
ation therapy,” Medical Dosimetry 24, 115-119 (1999).
ternationalJournalofRadiationOncology*Biology*Physics
7. P. Kupelian, T. Willoughby, A. Mahadevan, T. Djemil,
63, 1442-1447 (2005).
G. Weinstein, S. Jani, C. Enke, T. Solberg, N. Flores and D.
19. C.Wu,R.Jeraj,G.H.OliveraandT.R.Mackie,”Re-
Liu,”Multi-institutionalclinicalexperiencewiththeCalypso
optimization in adaptive radiotherapy,” Physics in Medicine
System in localization and continuous, real-time monitoring
and Biology 47, 3181 (2002).
of the prostate gland during external radiotherapy,” Interna-
20. G. D. Hugo, D. Yan and J. Liang, ”Population and
tional Journal of Radiation Oncology* Biology* Physics 67,
patient-specific target margins for 4D adaptive radiotherapy
1088-1098 (2007).
to account for intra-and inter-fraction variation in lung tu-
8. M. J. Murphy, ”Tracking moving organs in real time,” mour position,” Physics in Medicine and Biology 52, 257
Seminars in Radiation Oncology 14, 91-100 (2004). (2007).
9. P. J. Keall, H. Cattell, D. Pokhrel, S. Dieterich, K. 21. C. Ozhasoglu and M. J. Murphy, ”Issues in respira-
H. Wong, M. J. Murphy, S. S. Vedam, K. Wijesooriya and tory motion compensation during external-beam radiother-
R. Mohan, ”Geometric accuracy of a real-time target track- apy,”InternationalJournalofRadiationOncology*Biology*
ing system with dynamic multileaf collimator tracking sys- Physics 52, 1389-1399 (2002).
tem,”InternationalJournalofRadiationOncology*Biology* 22. P. R. Poulsen, B. Cho and P. J. Keall, ”Real-time
Physics 65, 1579-1584 (2006). prostatetrajectoryestimationwithasingleimagerinarcra-
10. H. Shirato, K. Suzuki, G. C. Sharp, K. Fujita, R. On- diotherapy: A simulation study,” Physics in medicine and
imaru, M. Fujino, N. Kato, Y. Osaka, R. Kinoshita and H. biology 54, 4019 (2009).
Taguchi, ”Speed and amplitude of lung tumor motion pre- 23. P. R. Poulsen, B. Cho, K. Langen, P. Kupelian and
cisely detected in four-dimensional setup and in real-time P. J. Keall, ”Three-dimensional prostate position estimation
tumor-tracking radiotherapy,” International Journal of Ra- with a single x-ray imager utilizing the spatial probability
diation Oncology* Biology* Physics 64, 1229-1236 (2006). density,” Physics in medicine and biology 53, 4331 (2008).
11. Y.Seppenwoolde,H.Shirato,K.Kitamura,S.Shimizu, 24. R. Li, B. P. Fahimian and L. Xing, ”A Bayesian ap-
M. van Herk, J. V. Lebesque and K. Miyasaka, ”Precise and proach to real-time 3D tumor localization via monoscopic x-
real-time measurement of 3D tumor motion in lung due to ray imaging during treatment delivery,” Medical physics 38,
breathingandheartbeat,measuredduringradiotherapy,”In- 4205 (2011).
ternationalJournalofRadiationOncology*Biology*Physics 25. L. Papiez and M. Langer, ”On probabilistically de-
53, 822-834 (2002). fined margins in radiation therapy,” Physics in Medicine and
12. B. Cho, P. R. Poulsen, A. Sloutsky, A. Sawant and P. Biology 51, 3921 (2006).
J. Keall, ”First demonstration of combined kV/MV image- 26. L.A.Feldkamp,L.C.DavisandJ.W.Kress,”Practical
guided real-time dynamic multileaf-collimator target track- cone-beam algorithm,” J Opt Sci Am A 1, 612-619 (1984).
ing,” International Journal of Radiation Oncology* Biology* 27. J. C. Lagarias, J. A. Reeds, M. H. Wright and P. E.
Physics 74, 859-867 (2009). Wright,”ConvergencepropertiesoftheNelder–Meadsimplex
13. K. Malinowski, T. J. McAvoy, R. George, S. Diet- methodinlowdimensions,”SIAMJournalonoptimization9,
richandW.D.DSouza,”Incidenceofchangesinrespiration- 112-147 (1998).
induced tumor motion and its relationship with respiratory 28. J. Lewis, ”Fast template matching,” Vision interface
surrogates during individual treatment fractions,” Interna- 95, 15-19 (1995).