PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Med Image Anal. Author manuscript; available in PMC 2010 October 1.
Published in final edited form as:
PMCID: PMC2771756
NIHMSID: NIHMS143644

Constrained non-rigid registration for use in image-guided adaptive radiotherapy

Abstract

A constrained non-rigid registration (CNRR) algorithm for use in prostate image-guided adaptive radiotherapy is presented in a coherent mathematical framework. The registration algorithm is based on a global rigid transformation combined with a series of local injective non-rigid multi-resolution cubic B-spline Free Form Deformation (FFD) transformations. The control points of the FFD are used to non-rigidly constrain the transformation to the prostate, rectum, and bladder. As well, the control points are used to rigidly constrain the transformation to the estimated position of the pelvis, left femur, and right femur. The algorithm was tested with both 3D conformal radiotherapy (3DCRT) and intensity-modulated radiotherapy (IMRT) dose plan data sets. The 3DCRT dose plan set consisted of 10 fan-beam CT (FBCT) treatment-day images acquired from four different patients. The IMRT dose plan set consisted of 32 cone-beam CT (CBCT) treatment-day images acquired from 4 different patients. The CNRR was tested with different combinations of anatomical constraints and each test significantly outperformed both rigid and non-rigid registration at aligning constrained bones and critical organs. The CNRR results were used to adapt the dose plans to account for patient positioning errors as well as inter-day bone motion and intrinsic organ deformation. Each adapted dose plan improved performance by lowering radiation distribution to the rectum and bladder while increasing or maintaining radiation distribution to the prostate.

Keywords: Non-rigid image registration, B-spline free form deformation, Rigid registration, Image-guided radiotherapy, Adaptive prostate radiotherapy

1. Introduction

Prostate cancer is the most commonly diagnosed cancer among men in the United States. The tumor is localized to the prostate for 90% of the cases, making external beam radiotherapy (EBRT) a viable, non-invasive treatment option for the majority of patients (Potosky et al., 2000). For this reason, EBRT has become the primary treatment modality for prostate cancer (Potosky et al., 2000).

The latest advancement of image-guided adaptive radiotherapy (IGART) integrates an in-room cone-beam CT (CBCT) imager with radiotherapy linear accelerators for treatment-day imaging. The treatment-day CBCT images enable the radiotherapist to rigidly adjust the position of the patient to best align the patient to the position used during planning. Meijer et al. showed that the prostate planning margin can be reduced from 11 mm to 8 mm when aligning the patient using rigid registration to align bony anatomy (Meijer et al., 2008).

Treatment-day images allow changes in tumor position, size, and shape that take place over the course of radiotherapy to be measured and accounted for in order to boost geometric accuracy and precision of radiation delivery. Increased precision and accuracy are expected to augment tumor control, reduce incidence and severity of unintentional toxicity, and facilitate development of more efficient treatment methods than currently available (Dawson and Sharpe, 2006).

Very little work has been done to develop an on-line method to non-rigidly update dose plans, despite the push in the literature (Li et al., 2008; Frank et al., 2008). Foskey et al. developed a non-rigid registration (NRR) technique which first shrinks bowel gas prior to registration, allowing for a diffeomorphism between the domain of the two images (Foskey et al., 2005). Foskey used the resulting non-rigid transformation to segment the organs from the treatment-day CT in order to determine the dosimetric effect of a given plan in the presence of non-rigid tissue motion, not to update the dose plan to account for such motion.

In this paper, we present a constrained non-rigid registration (CNRR) algorithm combined with a bone segmentation algorithm designed to meet the demands of IGART. Our algorithm segments the pelvis, left femur, and right femur, registers the planning-day fan-beam CT (FBCT) image to the treatment-day image, then updates the IMRT dose plan using the estimated transformation. A hierarchical composite transformation model was developed to capture the global and local motion of the patient anatomy. The global motion from the planning-day image to the treatment-day image is modeled using a rigid transformation, while the unpredictable rigid bone motion and internal organ deformation is captured using a series of non-rigid cubic B-spline Free Form Deformation (FFD) transformations. The composite transformation from the CNRR algorithm maps the original dose plan into the treatment-day domain, effectively adapting the dose plan to the treatment-day position of the organs and bones. An updated dose plan that accurately accounts for bone and organ motion will allow for tighter planning margins and greater radiation dose, while maintaining or lowering bladder and rectum toxicity levels (Dawson and Sharpe, 2006; Klein et al., 2005; Li et al., 2008; Frank et al., 2008).

Our previous work in (Greene et al., 2008a) was only on FBCT-FBCT matching and 3DCRT plans and did not incorporate automatic bone segmentation and matching; the current work addresses the more difficult FBCT-CBCT matching problem in addition to FBCT-FBCT matching, works with both forward-looking IMRT and 3DCRT dose plans and describes a robust algorithmic approach that includes both automated segmentation of key bone structure and automated, constrained non-rigid mapping.

2. Methods

A CNRR algorithm for use in updating prostate radiotherapy dose plans has been developed. The CNRR algorithm has four steps:

  1. Segmentation. The bones are segmented using an automated fast level set segmentation algorithm. The organs are hand segmented from both the planning-day and treatment-day images prior to the start of the registration process.
  2. Global rigid registration. A global rigid registration is used to capture the gross motion between the planning-day image I0 and the treatment-day image Id.
  3. Iterative non-rigid registration. The images are initially down-sampled to a low resolution (8 mm). At each iteration, planning-day organs {prostate, rectum, bladder} are independently registered to treatment-day organs using a non-rigid transform. Planning-day bones {pelvis, left femur, right femur} are independently registered to treatment-day bones using a rigid transform. The constrained registration results are then used to constrain the objects in the CNRR to their individually estimated transformations, generating T0d, the composite transformation from the planning-day image I0 to the treatment-day image Id. At the start of the next iteration, the images are sampled at a higher resolution, updated to account for the registration results from the previous iteration, and the registration process is repeated.
  4. Update dose plan. Once the registration process is complete, the composite transformation T0d is used to update the planning-day dose plan P0 to account for articulated rigid bone motion and organ deformation, creating Pd, the adapted dose plan.

A detailed description of the algorithm is provided in the following sections.

2.1. Segmentation

The soft tissue organs were hand segmented in the planning-day FBCT images and treatment-day images by a radiation oncologist. The planning-day and treatment-day segmented images are denoted as S0obj and Sdobj, respectively, where obj [set membership] {prostate, rectum, bladder, bones}. The segmented images are all binary. The bones were segmented using a fast level set method (Shi and Karl, 2005). This technique avoids directly solving the PDEs connected with the level set method resulting in a reduction of computation time. In this method, contours are defined by two sets of points that are located on the inner and outer edges of the boundary. The speed function consists of an evolution speed term derived from image data and a regularization speed term that defines the boundary smoothness. The implementation is briefly outlined.

Let an image IRN be defined on a uniformly sampled grid. Any point x = (x1,x2, . . . , xN) on the grid is denoted by its coordinates. An object boundary represented implicitly by the zero level set of a function ϕ, lying on the grid, is defined uniquely by two lists of points Łin and Łout that are formally defined as:

Łin={xϕ(x)<0andyN(x)suchthatϕ(y)>0}Łout={xϕ(x)>0andyN(x)suchthatϕ(y)<0}

where N(x) is a discrete neighborhood of x and is defined as N(x)={yI;yx=1}xI. The boundary interface lies. The boundary interface lies between Łin and Łout. Computation time is minimized by limiting the values of ϕ to a discrete set of values: 3 if x is exterior to Łout and is not enclosed by the boundary, 1 if x is on Łout, –1 if x is on Łin and –3 if x is interior to Łin and is enclosed by the boundary. This is schematically illustrated in Fig. 2.

Fig. 2
A sample segmentation figure. The level set function ϕ assumes discrete values of [–3,–1,0,1,3] depending upon the voxel's location.

The interface evolution is driven by the velocity field v in the level set equation: ϕt+vϕ=0. If the value of v(x) of a point x [set membership] Łout is positive, the boundary moves forward and x is transferred to Łin and ϕ(x) is set to –1. Any gaps in Łout are reconnected by including points in N(x) that are equal to 3 and reset to the value of Łout, i.e. 1. The evolution of Łin is handled in a converse manner. At this stage, both the inward and outward boundaries are smoothed with an appropriate Gaussian filter. Once the curve evolutions are completed, the final boundary of the bony structure ϕ(x) is returned.

2.2. Global rigid registration

The global motion from the planning-day image I0 to the treatment-day image Id is modeled using a rigid transformation T0dr,global. To estimate T0dr,global, we use a normalized cross-correlation (NCC) match metric CNCC.

Given that the rigid registration (RR) is an inter-modality registration between FBCT and CBCT, we can not use a direct comparison of image intensities, i.e., sum of squared differences, as a similarity measure. Although the image intensities for both imaging modalities are in CT numbers, their intensities differ by a linear transformation and Gaussian noise due to X-ray scatter during CBCT acquisition. An alternative voxel-based similarity measure is normalized cross-correlation (NCC) which has been shown to accurately and robustly align inter-modality images (Hill et al., 2001; Holden et al., 2000). Furthermore, NCC has successfully been used to register FBCT to CBCT for use in prostate IGRT (Yang et al., 2007). As well, NCC has been recommended over MI for registering FBCT to CBCT because it is computationally more efficient (Thilmann et al., 2005).

CNCC(A,B)=1Ni=0N(A(i)A)(B(i)B)σAσB

where i indexes over the voxels. A, B denote the average image intensity and σA, σB represent the intensity standard deviation. If both images are aligned, CNCC is maximized.

A conjugate gradient method (Press et al., 1992) is used to estimate the rigid transformation T0dr,global that maps I0 to Id such that –CNCC is minimized

T^0dr,global=argminT0dr,globalCNCC(I0(T0dr,global),Id).
(1)

2.3. Resample and map images

After the global rigid registration, the images are sampled to a low resolution at the start of the first non-rigid registration iteration, then sampled at a higher resolution at the start of each additional iteration.

Once the images have been resampled they are transformed to account for the registration results from the previous iteration.

I0=I0(T0d)Sobj0=Sobj0(T0d)

where T0d is the composite non-rigid transformation that maps the planning-day image I0 into the treatment-day image Id. In the case of the first iteration, the images are transformed using the global rigid transformation T^0dr,global.

2.4. Independent organ registration

The segmented organ images are independently registered using a non-rigid transform based on a multi-resolution cubic B-Spline FFD. A FFD based transform was chosen because a FFD is locally controllable due to the underlying mesh of control points which are used to manipulate the image (Rogers, 2001).

A cubic B-spline FFD is defined by designating the image volume as Ω={(x,y,z)0x<X,0y<Y,0z<Z}. Φ denotes a (l + 3) × (m + 3) × (n + 3) mesh of control points ϕi,j,k with uniform spacing δx, δy, δz = δ. The parameter domain of the image is defined as ϴ={(u,v,w)0ul,0vm,0wn}. The FFD can then be written as a 3D tensor product of 1D cubic B-splines

T0d(x,y,z)=i=03j=03k=03Bi(uu)Bj(vv)Bk(ww)ϕ(u+i)(v+j)(w+k)

where T0→d maps I0 to Id, and u=xδx, v=yδy, w=zδz, and, Bi represents the ith basis function of the cubic B-spline.

B0(u)=(1u)36B2(u)=(3u3+3u2+3u+1)6B1(u)=(3u36u2+4)6B3(u)=u36.

Cubic B-splines have limited support, such that changing one control point only affects the transformation in the local neighborhood of the manipulated control point (Rogers, 2001; Greene et al., 2007).

Pseudovertices are used to ensure that the control volume being deformed by the control point mesh spans the entire image volume. The pseudovertices are the last control points along each axis and do not lie within the image domain. In order to guarantee that the control volume of the control point mesh, when deformed, spans the entire image volume, the last three control points, two contained within the image domain and one outside the image domain (pseudovertex), must not be manipulated (Rogers, 2001; Greene et al., 2007).

The degree of non-rigid motion captured depends on the spacing of the control points in the mesh Φ. In order to create an algorithm with an adequate degree of non-rigid deformation, a hierarchical multi-resolution approach is implemented in which the control point resolution is increased in a coarse to fine fashion at each iteration (Rueckert et al., 1999; Greene et al., 2007; Greene et al., 2008a,b).

Let Φ1, Φ2, . . . , ΦN denote a hierarchy of control point meshes, each improving in resolution over the previous mesh. T1, T2, . . . , TN are the transformations associated with each control point mesh. The composition of all transformations T(Ω)=TNT2T1(Ω), is a map that defines the deformation of image volume Ω (Rueckert et al., 1999; Greene et al., 2007; Greene et al., 2008a).

Given that the registration results from the previous iteration are used to update the current segmented organ, the composite transformation T0dobj which maps the planning-day organ to the treatment-day organ is

T0dobj(x,y,z)=T^0dobjnT0d(x,y,z)

where obj [set membership] {prostate, rectum, bladder}, and T0d is the constrained composite mapping from the previous iteration, or T^0dr,global in the case of the first iteration.

In order to determine T^0dobjn, we developed an objective function based on a sum of squared differences intensity matching term CSSD and a transformation penalty term Cinjective.

The sum of squared differences intensity match metric CSSD is written as

CSSD(A,B)=1Ni=0N[A(i)B(i)]2

where i indexes over the voxels. Although the planning and treatment-day images are inter-modal, the object segmentations are binary which allows a direct comparison of image intensities to be used as a registration metric (Hill et al., 2001).

In the following we will develop a cost function to ensure the FFD transformation is both smooth and injective. A FFD transformation is a tensor product of differentiable 1D cubic B-splines, which gives a FFD transformation an inherent smoothness property because polynomials are always continuously differentiable. Choi and Lee (1999) presented sufficient conditions for the local injectivity of a 3D uniform cubic B-spline function. An injective function is a function which has a unique output for each unique input. Given

T:UVS:VU

then T is injective if S(T(u)) = u for every u in U. The injective conditions presented by Choi and Lee (1999) are used in this algorithm to guarantee the injectivity of the FFD transformation.

Without loss of generality, we will assume that the control points are arranged on a lattice with unit spacing. Let Δϕijk = (Δxijk, Δyijk, Δzijk) be the displacement of control point ϕijk. Let δx = max|Δxijk|, δy = max|Δyijk|, δz = max|Δzijk|.

Property. A FFD based on cubic B-splines is locally injective over all the domain if δx<1K, δy<1K, and δz<1K, where K ≈ 2.48.

Choi and Lee (1999) have determined a value of K ≈ 2.48, thus the maximum displacement of control points given by the bound 1K is approximately 0.40. This means that the maximum displacement of control points is determined by the spacing of control points in the lattice. For example, a lattice with 20 mm control point spacing will have a maximum control point displacement of 20 mm × 0.40 = 8 mm. Note, this condition applies only to individual transformations. Larger effective displacements are possible when transformations are composed, such as in our multiresolution formulation.

To ensure the transformation is injective, each control point is restricted to move within a sphere of radius r < R, where R ≈ 0.40 × δ and δ is the uniform control point spacing (Choi and Lee, 1999; Rueckert et al., 2006; Greene et al., 2007; Greene et al., 2008a,b). With this, the injective penalty term is posed as

Cinjective(ΔΦ)=kf(ΔΦk)

where k indexes over the control points and

f(ΔΦk)={0ifΔΦk<0.40×δΔΦkotherwise}

At each iteration n, T^0dobjn is estimated by using a conjugate gradient method to minimize the following objective function with respect to T0dobjn

T^0dobjn=argminT0dobjn[CSSD(Sdobj,S0obj(T0dobjn(Φobjn)))+Cinjective(ΔΦobjn)]
(2)

The organ control point mesh Φobj–n is used to constrain obj in the CNRR.

2.5. Bone registration

The segmented bone images are registered at each iteration using a rigid transformation to estimate the composite transformation which maps S0bone into Sdbone. Given that the registration results from the previous iteration are used to update the current segmented bone images, the composite transformation T0dr,bone is written as

T0dr,bone(x,y,z)=T^0dr,bonenT0d(x,y,z)

where T^0dr,bonen is the nth rigid update to the composite transformation.

To estimate T^0dr,bonen, we minimize CSSD with respect to T0dr,bonen using the conjugate gradient method

T^0dr,bonen=argminT0dr,bonenCSSD(Sdbone,S0bone(T0dr,bonen))
(3)

The rigid transformation T^0dr,bonen is applied to an associated control point mesh Φbone–n, which is used to constrain the bones in the CNRR.

2.6. Constrained non-rigid registration

The constrained non-rigid registration algorithm uses the same FFD transform developed in Section 2.4 to model the constrained local motion which maps I0 into Id. The rigid and local motion from I0 to Id are modeled using a composite transformation T0d

T0d(x,y,z)=T^0dnT^0d1T^0dr,global(x,y,z)

where T^0dn is the nth local motion update to the composite transformation. T0d is updated at each iteration such that

T0d(x,y,z)=T^0dnT0d(x,y,z)

where T0d is set to T^0dr,global prior to the first iteration. To estimate T^0dn, we developed a cost function based on CNCC, Cinjective, and an object matching metric Cobject.

The object matching metric Cobject aligns the organs and bones by forcing the control points that lie within the segmented objects to their position estimated by the individual object registration.

Cobject(ΔΦn)=objkobj[ΔΦknΔΦknobj]2

where k indexes over the control points, obj indexes over the constraint objects (prostate, rectum, bladder, pelvis, left femur, right femur).

At each iteration, T^0dn is estimated by using a conjugate gradient method to minimize the following objective function with respect to T0dn

T^0dn=argminT0dn[Cobject(ΔΦn)CNCC(Id,I0(T0dn(Φn)))+Cinjective(ΔΦn)]
(4)

2.7. Adapting and assessing the dose plan

Once the registration process has completed, the composite transformation T0d is used to map the original dose plan into the treatment-day domain, effectively adapting the dose plan to account for patient positioning errors as well as inter-day bone and intrinsic organ motion.

Pd=P0(T0d)

where P0 and Pd represent the unadapted planning-day dose plan and the adapted treatment-day dose plan, respectively.

The adapted dose can be used directly as the dosimetric goal in dose plan optimization for the treatment-day anatomy, however, this approach is not yet practical for real-time treatment adaptation. The time needed for dose plan optimization is still relatively long and quality assurance checks of the new plan are required before it can be delivered to the patient. Alternatively, a library of pre-approved dose plans can be established based on the variations of tissue deformation observed from a similar population of patients. On a given treatment day, the adapted dose plan and anatomy can be used to select an appropriate dose plan from the library for treatment delivery (after the approval of a radiation oncologist). This approach, while not perfect, has the potential for implementing adaptive radiotherapy in real-time.

The dose distributions generated by the adapted dose plan Pd in the tumor and organs at risk are analyzed by computing tumor control probability (TCP) and normal tissue complication probability (NTCP). TCP is a measure of the effectiveness of the administered therapy in destroying the malignant cells within the target volume, and NTCP is an estimate of the risk to the neighboring normal organs. NTCP is computed for both the rectum and bladder. These two parameters are used to predict the success of a dose plan (Lyman, 1985; Niemierko and Goitein, 1993; Jackson et al., 1995; Chelikani et al., 2006). An optimal dose plan maximizes TCP while minimizing both rectal and bladder NTCP.

To compute NTCP, we employed the widely used Lyman model (Lyman, 1985). To compute TCP, we used the Niemierko model developed by Niemierko and Goitein (1993), which is a mechanistic cell kill model based on Poisson statistics (Levegrun et al., 2001). For more details on NTCP and TCP, the reader is referred to the literature (Lyman, 1985; Niemierko and Goitein, 1993).

3. FBCT and 3DCRT results

We used 10 treatment-day FBCT images acquired from four different patients undergoing prostate EBRT at Yale-New Haven Hospital. Two patients each had three treatment-day FBCT images, and an additional two patients each had two treatment-day FBCT images. Each of the four patients analyzed had an associated planning-day 3DCRT dose plan and planning-day FBCT image. Please refer to Fig. 1 for an example of a planning-day 3DCRT dose plan.The prostate, rectum, and bladder were hand segmented by a qualified clinician. The pelvis, left femur, and right femur were segmented using a semiautomatic thresholding algorithm available in BioImage Suite (Papademetris et al., 2008). The automated fast level set segmentation algorithm incorporated into the CNRR algorithm was not used because the bones had been previously segmented.

Fig. 1
This figure shows an example of a IMRT dose plan and a 3DCRT dose plan created for the same patient undergoing EBRT for prostate cancer at Yale-New Haven Hospital. The dose plans are overlaid on the planning-day FBCT image. The IMRT dose plan is on the ...

3.1. Registration models

The FBCT images had a spatial resolution of 1.17 mm × 1.25 mm × 1.17 mm. The initial rigid registration was performed at full resolution, while the non-rigid registration was performed using a multi-resolution framework where the first iteration was performed at a low image resolution 8 mm3 and low control resolution 40mmcp. The image resolution and control point resolution were each increased in subsequent iterations to 4 mm3 and 20mmcp by the third non-rigid iteration. Details of the registration parameters are listed in Table 1.

Table 1
CNRR algorithm parameters for 3DCRT data set.

The CNRR algorithm was tested using six different combinations of organ and bone constraints to determine if constrained organs and bones improved the registration of the unconstrained organs and bones. In order to compare the CNRR algorithm results to those achieved using standard registration algorithms, both a single-resolution rigid registration (RR) and multi-resolution non-rigid registration (NRR) were tested. The NRR test was performed using the same framework as the CNRR algorithm without organ or bone constraints. The six different constraint scenarios used to test the CNRR algorithm and both standard registration algorithms used are listed in Table 2.

Table 2
This table lists the eight different registration tests used to register the patient data. The RR and NRR tests were performed to assess how well the CNRR algorithm performs in comparison to standard registration algorithms. The CNRR algorithm was tested ...

3.2. Registration results

At the end of each registration iteration, the global composite transformation T0d was used to map the full resolution planning-day FBCT image into the treatment-day space. The quality of the registration was then assessed by computing the percent overlap of the mapped planning-day organs and bones and the treatment-day organs and bones.

The average registration and adapted dose plan results are presented in Fig. 3.

Fig. 3
(a) Registration results are averaged over 10 treatment-day FBCT images acquired from four different patients. The average object overlap increase is plotted for all 8 registration tests. (b) 3DCRT dose plan results from 10 different treatment days across ...

From prior experience with our CNRR algorithm (Greene et al., 2008a), we discovered that including the bones as rigid constraints does not impact the transformation in the region of the organs and further carry the organs along. However, for completeness, we included the pelvis, left femur, and right femur as rigid constraints. From Fig. 3, it is clear that including the bones as constraints in registration test CNRR-3 did not improve the average prostate, rectum, or bladder overlap when compared to NRR results.

The best prostate registration results occurred when the prostate was used as the only constraint or both the prostate and rectum were constrained. The maximum average prostate overlap was achieved for CNRR-2, however, there was no statistical difference between the prostate overlap results for all registration tests which constrained the prostate.

The greatest rectum overlap occurred for CNRR-2, CNRR-4, and CNRR-6, when the rectum was held as a constraint.

CNRR-5 used the prostate and bladder as constraints in addition to the left femur, right femur and the pelvis. The goal of CNRR-5 was to determine if only one organ was left out, would the overlap of the unconstrained organ improve when compared to RR and NRR results. From the results in Fig. 3a it is clear that constraining the organs and bones in the local area surrounding the rectum did not significantly improve rectum registration results when compared to RR and NRR results. Rectum overlap was not significantly improved by CNRR tests which did not use the rectum as a constraint when compared to RR and NRR results. The local control domain of the control points was not large enough to influence the registration of any surrounding unconstrained organs.

On average, bladder overlap improved by 10.20% for CNRR-4, CNRR-5, and CNRR-6 which all used the bladder as a constraint. Leaving the bladder out as a constraint and constraining the registration to the prostate and rectum during CNRR-2 did not significantly improve bladder overlap when compared to RR and NRR results.

3.3. Dose plan results

After each registration iteration, the global composite transformation T0d was used to map the 3DCRT dose plan into the treatment-day space. Improvements in the dose plan are quantified by comparing the unadapted treatment-day tumor control probability (TCTP) and normal tissue complication probability (NTCP) values to the adapted treatment-day values.

The 3DCRT dose plan results associated with each registration test are presented in Fig. 3b. The average adapted dose plan values (TCP, rNTCP, bNTCP) are normalized by dividing by the unadapted treatment-day value. Fig. 3d plots the normalized average dose plan values for each iteration of CNRR-6. In addition, the organ overlap results from CNRR-6 are plotted for each registration iteration in Fig. 3c to show the improvement in organ overlap associated with each iteration.

The normalized average adapted TCP value did not significantly change for any registration tests performed, despite an average 19.07% improvement in prostate overlap for registration tests which constrained the prostate. The planning margins are statistically chosen such that there is a low probability that the prostate will move outside of the planning margins throughout the course of treatment.

Rectal NTCP (rNTCP) improved for all registration tests performed and approached a steady state normalized average value of 0.60 as rectum overlap approached 80%. As rectal overlap improved beyond 80% in subsequent iterations, rNTCP did not further decrease.

Similarly, bladder NTCP (bNTCP) approached a steady state average normalized value of 0.70 as bladder overlap approached 92%, which represents a 30% improvement in bNTCP. Further improvement in bladder registration did not further decrease bNTCP. For registration tests which did not approach a 92% bladder overlap, no significant change in bNTCP was seen.

The data indicates that the greatest benefit from adapting the dose plan is decreasing the radiation dosage delivered to the rectum and bladder. The steady state phenomena for bNTCP and rNTCP, and the lack of improvement in TCP occurs as a result of the 3DRT planning margins (Chen et al., 2000, 2007). The planning margins are statistically chosen to account for prostate, rectum and bladder motion to ensure adequate dosage to the prostate and minimize exposure to the organs at risk throughout the course of treatment. This explains why a 19.07% increase in prostate overlap does not significantly increase TCP. As long as the prostate does not move outside the planning margins, there is no room to improve TCP. Similarly for the rectum and bladder, once the registration carries the organs outside of the planning margins, further decreases in NTCP can not be achieved.

4. CBCT and IMRT results

4.1. Patient data

We tested our CNRR algorithm on 32 CBCT treatment-day images acquired from four different patients undergoing prostate IGRT at Yale-New Haven Hospital. Each patient data set consisted of a planning-day FBCT, eight weekly treatment-day CBCTs, and a planning-day IMRT dose plan. Please refer to Fig. 1 for an example of a planning-day IMRT dose plan.

4.2. Registration models

The planning-day and treatment-day images were down-sampled from 1.25 × 1.25 × 2.0 mm3 and 2.0 × 2.0 × 5.0 mm3, respectively, to 4.0 × 4.0 × 4.0 mm3 for the initial rigid registration. Three CNRR iterations were performed at 8.0 mm3, 6.0 mm3, and 4.0 mm3 with a control point spacing of 35mmcp, 30mmcp, and 25mmcp, respectively. Please see Table 3 for further registration parameter details.

Table 3
This table lists the CNRR algorithm parameters used for each registration iteration.

The CNRR algorithm was tested using six different combinations of organ and bone constraints to determine if constrained organs and bones improved the registration of the unconstrained organs and bones. In order to compare the CNRR algorithm results to those achieved using standard registration algorithms, both a single-resolution RR and multi-resolution NRR were tested. The NRR test was performed using the same framework as the CNRR algorithm without organ or bone constraints. The six different constraint scenarios used to test the CNRR algorithm and both standard registration algorithms used are listed in Table 4.

Table 4
This table lists the eight different registration tests used to register the patient data. The RR and NRR tests were performed to assess how well the CNRR algorithm performs in comparison to standard registration algorithms. The CNRR algorithm was tested ...

4.3. Registration results

At the end of each iteration, the global composite transformation T0d was used to map the original full resolution planning-day FBCT image into the treatment-day space. The quality of the registration was then assessed by computing the percent overlap of the mapped planning-day organs and bones and the treatment-day organs and bones.

The planning-day and treatment-day images were highly aligned due to manual patient alignment prior to CBCT imaging and dose delivery, so the RR test did not change average bone or organ alignment.

As before, for completeness, we used the bone as a constraint (CNRR-10) on two patient data sets (16 treatment-days). The CNRR-10 test improved bone overlap when compared to RR results (p < 0.001, paired t-test) and NRR results (p < 0.001, paired t-test). However, from Fig. 4b, it is clear that including bone as a constraint did not improve the average prostate, rectum, or bladder overlap when compared to NRR results from the same two patient data sets. CNRR-9 results tested across the same two patient data sets were included in Fig. 4b to show average prostate, rectum, and bladder results when they are directly constrained.

Fig. 4
(a) Average registration results from 32 different treatment days across four different patients. The change in percent organ overlap is plotted for five different registration tests. (b) Average registration results with bone from 16 different treatment ...

The registration results presented in Fig. 4a were tested across four patient data sets for a total of 32 different treatment-days. Because the prostate was constrained for all four CNRR tests, all CNRR tests significantly increased prostate overlap when compared to NRR results. CNRR-8 increased the average prostate overlap the most when compared to NRR results (p < 0.001, paired t-test), taking the average prostate overlap from 65.01% to 91.57%.

Average initial rectum overlap was very low (3.55%) due to organ deformation, so large improvements in overlap were not expected. CNRR-2 and CNRR-7 increased rectum overlap the most when compared to NRR results, increasing rectum overlap to 8.58% and 8.62%, respectively. The CNRR-7 test significantly outperformed the NRR algorithm at aligning the rectum (p < 0.0191, paired t-test), as did CNRR-2 (p < 0.0195, paired t-test). CNRR-1, which did not constrain the rectum, did not improve rectum alignment.

Average initial bladder overlap was high (94.25%), so there was not much room for improvement. The greatest average increase in bladder overlap occurred for CNRR-8 which increased bladder overlap by 2.71%, however, there was no significant statistical improvement when compared to NRR results which increased bladder overlap by 2.09%. CNRR-8 combined all three organs into a single constraint image and increased average bladder overlap from 94.25% to 96.80%.

4.4. Dose plan results

After each iteration, the global composite transformation T0d was used to map the planning-day IMRT dose plan into the treatment-day space. Please refer to Fig. 5 for an example of an adapted IMRT dose plan overlaid on a treatment-day CBCT image. Improvements in the dose plan are quantified by comparing the unadapted treatment-day TCP and NTCP values to the adapted treatment-day TCP and NTCP values.

Fig. 5
This figure shows an adapted IMRT dose plan overlaid on the associated treatment-day CBCT image. The IMRT dose plan was warped using the CNRR-3 estimated displacement field. The contour at the top outlines the prostate and the contour below it outlines ...

Dose plan results from CNRR-7 are presented in Fig. 4d. The average adapted dose plan values (TCP, rNTCP, bNTCP) are normalized by dividing by the unadapted treatment-day value and are plotted for each iteration of the CNRR algorithm. In addition, the organ overlap results from CNRR-7 are plotted in Fig. 4c to show the improvement in organ overlap associated with each iteration. TCP is maximized to 1.036 after the second iteration. Rectum NTCP (rNTCP) reached a minimum steady state value of 0.787 by the third iteration. Bladder NTCP (bNTCP) reached a minimal value of 0.781 by the second iteration, and did not significantly change for further iterations. Organ registration is improved at each iteration, yet the dose plan values do not improve at each iteration once the dose plan values approach steady state. The plateau in NTCP achieved is again due to the relatively wide planning margins used to account for organ motion. The planning margins also help explain why an average increase on 26.56% in prostate alignment only increases TCP by 3.6%. Once the organs are contained within the treatment margins, further registration accuracy does not improve TCP or NTCP. The statistical placement of the planning margins explains why such a small change in average organ overlap caused a large decrease in NTCP measures. A 5.05% increase in rectum overlap decreased rNTCP by 21.3% and a 2.27% increase in bladder overlap decreased bNTCP by 21.9%.

5. Summary

Our proposed CNRR algorithm allows accurate registration by incorporating organ and bone constraints in a hierarchical composite non-rigid B-spline transformation (Greene et al., 2008a). The approach includes automated bone segmentation and allows radiotherapy dose plans to be carried forward. The CNRR algorithm is described in a coherent mathematical framework that easily permits incorporation of new constraints into the processing.

Testing of our approach on both 3DCRT plans with fan-beam CT data and IMRT plans with cone-beam CT data showed that all combinations of constraints resulted in improvement in organ overlap compared to unconstrained rigid and multiresolution non-rigid alignment. While bone constraints were not found to be helpful in aligning the organs, good organ overlap was achieved when prostate, bladder and rectum were included as constraints. Prescribed planning margins, however, limit the corresponding improvement in TCP, as the prostate is ordinarily contained within the relatively wide planning margins. Rectum and bladder NTCP, on the other hand, improve due to the constraints of these organs in the registration process. The resulting adapted dose can be used as a dosimetric goal for dose plan optimization or to facilitate plan selection from a library of pre-approved plans.

The ultimate goal was to develop an accurate and robust algorithm such that the dose plan needed to deliver the CNRR mapped dose distribution on the treatment day can be directly determined in real-time from the tissue deformation and the initial planning-day dose plan approved during treatment planning.

References

  • Chelikani S, Purushothaman K, Knisely J, Chen Z, Nath R, Bansal R, Duncan J. A gradient feature weighted minimax algorithm for registration of multiple portal images to 3DCT volumes in prostate radiotherapy. Int. J. Radiat. Oncol. Biol. Phys. 2006;65(2):535–547. [PMC free article] [PubMed]
  • Chen Z, Song H, Deng J, Yue N, Knisely J, Peschel R, Nath R. What margin should be used for prostate IMRT in absence of daily soft-tissue imaging guidance? Med. Phys. 2000:31.
  • Chen Z, Deng J, Nath R, Peschel R. A two-stage imrt treatment protocol provides more robust radiotherapy for prostate cancer in presence of inter- and intra-fractional organ motions. Int. J. Radiat. Oncol. Biol. Phys. 2007;69:S713–S714.
  • Choi Y, Lee S. Local injectivity conditions of 2D and 3D uniform cubic b-spline functions. Proc. Pacific Graph. 1999:302–311.
  • Dawson LA, Sharpe MB. Image-guided radiotherapy: rationale, benefits, and limitations. Lancet Oncol. 2006;7:848–858. [PubMed]
  • Foskey M, Davis B, Goyal L, Chang S, Chaney E, Strehl N, Tomei S, Rosenman J, Joshi S. Large deformation three-dimensional image registration in image-guided radiation therapy. Phys. Med. Biol. 2005;50:5869–5892. [PubMed]
  • Frank SJ, Dong L, Kudchadker RJ, Crevoisier RD, Lee AK, Cheung R, Choi S, O'Daniel J, Tucker SL, Want H, Kuban DA. Quantification of prostate and seminal vesicle interfraction variation during imrt. Int. J. Radiat. Oncol. Biol. Phys. 2008;71(3):813–820. [PubMed]
  • Greene WH, Chelikani S, Papademetris X, Knisely JPS, Duncan J. A constrained non-rigid registration algorithm for application in prostate radiotherapy. ISBI. 2007:740–743. [PMC free article] [PubMed]
  • Greene WH, Chelikani S, Purushothaman K, Chen Z, Knisely JPS, Staib LH, Papademetris X, Duncan J. A constrained non-rigid registration algorithm for use in image-guided radiotherapy. MICCAI. 2008a;5241:780–788. [PMC free article] [PubMed]
  • Greene WH, Chelikani S, Papademetris X, Staib LH, Knisely JPS, Duncan J. Tracking organ overlap for a constrained non-rigid registration algorithm. ISBI. 2008b:1159–1162. [PMC free article] [PubMed]
  • Hill DLG, Batchelor PG, Holden M, Hawkes DJ. Topical review: medical image registration. Phys. Med. Biol. 2001;46:R1–R45. [PubMed]
  • Holden M, Hill DLG, Denton ERE, Jarosz JM, Cox TCS, Rohlfing T, Goodey J, Hawkes DJ. Voxel similarity measures for 3-D serial MR brain image registration. IEEE Trans. Med. Imag. 2000;19(2):94–102. [PubMed]
  • Jackson A, Haken RKT, Robertson JM, Kessler ML, Kutcher GJ, Lawrence TS. Analysis of clinical complication data for radiation hepatitis using a parallel architecture model. Int. J. Radiat. Oncol. Biol. Phys. 1995;31(4):883–891. [PubMed]
  • Klein EE, Drzymala RE, Purdy JA, Michalski J. Errors in radiation oncology: a study in pathways and dosimetric impact. J. Appl. Clin. Med. Phys. 2005;62:1517–1524. [PubMed]
  • Levegrun S, Jackson A, Zelefsky MJ, Skwarchuk MW, Venkatraman ES, Schlegel W, Fuks Z, Leibel SA, Ling CC. Fitting tumor control probability models to biopsy outcome after three-dimensional conformal radiation therapy of prostate cancer: pitfalls in deducing radiogiologic parameters for tumors from clinical data. Int. J. Radiat. Oncol. Biol. Phys. 2001;51(4):1064–1080. [PubMed]
  • Li HS, Chetty IJ, Enke CA, Foster RD, Willoughby TR, Kupellian PA, Solberg TD. Dosimetric consequences of intrafraction prostate motion. Int. J. Radiat. Oncol. Biol. Phys. 2008;71(3):801–812. [PubMed]
  • Lyman JT. Complication probability as assessed from dose-volume histograms. Radiat. Res. 1985;104:S13–S19. [PubMed]
  • Meijer GJ, Klerk J, Bzdusek K, Berg HA, Janssen R, Kaus MR, Rodrigus P, Toorn P. What CTV-to-PTV margins should be applied for prostate irradiation? four-dimensional quantitative assessment using model-based deformable image registration techniques. Int. J. Radiat. Oncol. Biol. Phys. 2008;72(5):1416–1425. [PubMed]
  • Niemierko A, Goitein M. Implementation of a model for estimating tumor control probability for an inhomogeneously irradiated tumor. Radiat. Oncol. 1993;29:140–147. [PubMed]
  • Papademetris X, Jackowski M, Rajeevan N, Constable R, Staib L. BioImage Suite: An Integrated Medical Image Analysis Suite, Section of Bioimaging Sciences. Dept. of Diagnostic Radiology, Yale School of Medicine; 2008. < http://www.bioimagesuite.org>. [PMC free article] [PubMed]
  • Potosky AL, Legler J, Albertsen PC, Stanford JL, Gilliland FD, Hamilton AS, Eley JW, Stephenson RA, Harlan LC. Health outcomes after prostatectomy or radiotherapy for prostate cancer: results from the prostate cancer outcomes study. J. Nat. Cancer Inst. 2000;92(19):1582–1592. [PubMed]
  • Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in C. Cambridge University Press; Cambridge, UK: 1992.
  • Rogers D. An Introduction to NURBS: with Historical Perspective. first ed. Morgan Kaufmann Publishers; San Francisco: 2001.
  • Rueckert D, Sonoda L, Denton E, Rankin S, Hayes C, Leach MO, Hill D, Hawkes DJ. Comparison and evaluation of rigid and non-rigid registration of breast MR images. SPIE. 1999;3661:78–88. [PubMed]
  • Rueckert D, Aljabar P, Heckemann R, Hajnal JV, Hammers A. Diffeomorphic registration using b-splines. MICCAI. 2006:702–709. [PubMed]
  • Shi Y, Karl WC. A fast level set method without solving PDEs. ICASSP. 2005:97–100.
  • Thilmann C, Nill S, Tucking T, Hoss A, Hesse B, Dietrich L, Bendl R, Rhein B, Haring P, Thieke C, Oelfke U, Debus J, Huber P. Correction of patient positioning errors based on in-line cone beam CTs: clinical implementation and first experiences. Radiat. Oncol. 2005;63(2):s550–s551. [PMC free article] [PubMed]
  • Varian Eclipse Treatment Planning System . version 8.1. Varian Medical Systems; Palo Alto, CA: 2007. < http://www.varian.com>.
  • Yang Y, Schreibmann E, Li T, Wang C, Xing L. Evaluation of on-board kV cone beam C (CBCT)-based dose calculation. Phys. Med. Biol. 2007;52:685–705. [PubMed]