Abstract
Image-guided radiotherapy (IGRT), adaptive radiotherapy (ART), and online reoptimization rely on accurate mapping of the radiation beam isocenter(s) from planning to treatment space. This mapping involves rigid and/or nonrigid registration of planning (pCT) and intratreatment (tCT) CT images. The purpose of this study was to retrospectively compare a fully automatic approach, including a non-rigid step, against a user-directed rigid method implemented in a clinical IGRT protocol for prostate cancer. Isocenters resulting from automatic and clinical mappings were compared to reference isocenters carefully determined in each tCT. Comparison was based on displacements from the reference isocenters and prostate dose-volume histograms (DVHs). Ten patients with a total of 243 tCTs were investigated. Fully automatic registration was found to be as accurate as the clinical protocol but more precise for all patients. The average of the unsigned and offsets and the standard deviations (σ) of the signed offsets computed over all images were (avg. ± σ (mm)): 1.1 ± 1.4, 1.8 ± 2.3, 2.5 ± 3.5 for the clinical protocol and 0.6 ± 0.8, 1.1 ± 1.5 and 1.1 ± 1.4 for the automatic method. No failures or outliers from automatic mapping were observed, while 8 outliers occurred for the clinical protocol.
1. Introduction
Image-guided radiotherapy (IGRT) [1], off-line adaptive radiotherapy (ART) [2], and online reoptimization [3] involve pretreatment imaging, taken here to be CT imaging. A procedure held in common by all three methods is registration of the planning (pCT) and treatment (tCT) images to map the planned isocenter to treatment space. Accuracy and precision of this step are important for delivering an accumulated dose distribution that closely matches the treatment plan. Mapping methods involve at least rigid registration. Ideally a nonrigid step would be included to account for differences in organ shape between planning and treatment times (Figure 1). The composite of the rigid, and possibly nonrigid, matrices is then used to map the planned isocenter to the tCT.
(a)
(b)
Quantitative evaluation of isocenter mapping methods is muddled by the lack of gold standards [4]. In the absence of standards, this work retrospectively compared a fully automatic method against a user-assisted procedure used in a clinical protocol for IGRT for prostate cancer. The study focused on tCTs from a conventional diagnostic scanner due to the availability of clinical data. However, the approach applies to kilovoltage cone-beam CT (CBCT) images as illustrated later.
2. Methods
2.1. Clinical Protocol
The clinical protocol (Table 1) was practiced over the period 2005–2008 at the University of North Carolina (UNC) as part of the routine workflow during evaluation of a CT-on-rails system [5] (Primatom, Siemens Medical Solutions, Concord, CA). The Primatom system consists of a conventional CT scanner, a set of rails between the scanner and the linac, and a dual-purpose imaging/treatment couch that moves along the rails. Patients under the IGRT protocol were imaged before treatment, and the couch was rolled from the CT scanner into the image-guided treatment position. The fixed geometry between the scanner and linac allows cross calibration of the patient laser alignment systems, and the rails facilitate moving the couch from the scanner into treatment position with minimal time delay or mechanical disturbance of patient geometry.
All patients underwent simulation using a Philips AcQSim CT system [6] (Philips Healthcare, Andover, MA). The planning isocenter was localized during simulation, and the anterior and lateral laser crosshairs were tattooed according to standard practice. Treatment planning and calculation of hypothetically delivered dose were accomplished with PlanUNC [7, 8]. (PLanUNC is a set of modular software tools for external beam treatment planning and dose calculation development at UNC.) During planning the dosimetrist placed the computer crosshairs at the planned isocenter and then marked reference points in the pCT at the intersections of the computer crosshairs with the skin (Figure 2). These reference points defined the geometrically correct positions of the anterior and lateral laser crosshairs based on patient geometry at planning time. Before acquisition of each tCT, the patient was positioned by aligning the tattoos on the patient’s skin with the CT laser crosshairs. After alignment, steel BBs were taped to the skin at the center of each laser’s crosshairs (AP, R&L lat). Assuming accurate laser calibration, this placement allows the treatment isocenter for the initial patient position to be inferred from the imaged BBs. Immediately after imaging, the tCT was imported to the planning system and the physicist defined the BB centers via point and click on a computer display screen. The planning system then automatically registered the tCT and pCT by matching the BBs in the tCT with the ideal locations marked in the pCT. The algorithm minimized (distance between corresponding points)2. The final step was interactive rigid registration of the prostate ROIs using the planning contours for reference. Table tilt around the lateral axis and rotation around the craniocaudal axis were not allowed in the final two steps because rotational errors were not handled during patient setup. Since the prostate was not contoured in the tCT, the final step relied on human judgment to match the pCT intensities and contours with the tCT intensity patterns in the prostate ROI. Table displacements were computed by comparing the coordinates of the mapped isocenter after BB registration with the coordinates after interactive prostate-based registration. The displacements were used to reposition the treatment table. Assuming accurate registration and repositioning, and no changes in patient anatomy between tCT imaging and treatment times, this procedure registered the mapped planned isocenter with the treatment isocenter. In typical practice of IGRT, displacements are implemented only when they exceed a predefined threshold. In the UNC protocol, the threshold was 3 mm along a given axis. However, for comparison purposes this study assumed that the displacements were applied without error regardless of magnitude.
(a)
(b)
2.2. Automatic Mapping Method
Automapping (Table 1) was performed using a beta version of ARTSuite (Morphormics, Inc. Chapel Hill, NC) installed at UNC. (Effective July 16, 2012, Morphormics, Inc. became a wholly owned subsidiary of Accuray, Inc., Sunnyvale, CA.) ARTSuite is a DICOM-RT compliant software system developed to support IGRT, ART, and online reoptimization with tools for autosegmentation of pCTs and tCTs, nonrigid mapping between treatment and planning spaces, and dose accumulation and analysis.
Automapping (Figure 3) involved three major steps: multiscale rigid registration of the previously segmented pCT to the tCT, model-based autosegmentation of the prostate in the tCT [9, 10], and nonrigid mapping of the planning isocenter to the tCT using a transformation that included the rigid registration matrix and a nonstatistical, nonrigid diffeomorphism determined from correspondences between the pCT and tCT prostate models.
(a)
(b)
(c)
(d)
(e)
(f)
2.3. Patient Data
Ten patients with a total of 243 tCTs were investigated. All pCTs were acquired with contrast in the bladder, and one patient had contrast in the bladder and rectum. The pixel dimensions in the axial plane were , and the slice thickness was 3 mm for pCTs and tCTs.
All patients were treated with step-and-shoot IMRT using an anterior and six oblique fields to ~75 Gy prescribed to a point or isodose curve. Treatment typically included a boost starting ~50 Gy to reach the final dose. Margins of 5 mm and 3 mm were applied to the planning target volume for the initial and boost plans, respectively.
2.3.1. Model Fitting to Prostate Contours
This study assumed that the manually drawn prostate contours in the pCT were true at treatment time. To facilitate model-based segmentation of the tCTs as described below, a prostate model was fit to each set of planning contours, using an approach based on prostate shape statistics described by Merck et al. [11]. This process yielded a custom prostate atlas for each patient that was used during registration and segmentation of tCTs. Contours were represented as many short line segments joined together. This representation caused small differences between a contour and a model at sharp vertices where two line segments joined. These differences were less than 0.2 mm per contour on average and were caused by smoothness constraints that forced the model to have continuous curvature.
2.3.2. Multiscale Rigid Registration
The overall rigid registration approach has similarities with that of Court and Dong [12] and Smitsmans et al. [13], discussed later in Section 3. To minimize compute time, registration progresses from coarse to fine scale in three steps. The output of each step serves as a prior for the following step. Step aligns the skin bounding boxes and registers the images by sliding the image data in the bounding box in the pCT along the box for the tCT. The registration algorithm computes the voxel count per slice in a predefined intensity window and aligns the images by finding the best match between graphs of voxel count versus slice position for each image pair. A bone window performs best but can fail when a “bright” contrast medium is used for the pCT. In such cases, Step is ignored and the algorithm starts over at Step , which optimizes global mutual information (MI) using a gradient descent approach [14]. To avoid convergence to a distant optimum, the algorithm is run multiple times with different starting points. The result with the best score over all runs is selected as the output. Step is similar to Step but with MI computed over an ROI defined by the atlas prostate model. As in the clinical procedure, rotations during automatic registration were not allowed explicitly. However, the segmentation step treated rotation as a nonrigid deformation in images with adequate intensity information. In Figure 4, for example, the rigidly registered prostate model overlaps pubic bone anteriorly, causing the model to deform in a manner that avoids overlap with both the bone and the gas bubble in the rectum.
2.3.3. Model-Based Autosegmentation
The model used to represent the prostate consists of a chain of so-called medial atoms (Figure 5). A full collection of atoms for an organ is called a medial-representation (m-rep) [9]. The chain configuration is well suited for objects that are more or less tube-like, including objects with closed ends. Each prostate atom has a hub and 16 spokes radiating to the organ surface. Additional hubs and spokes can be interpolated as needed. The skeletal framework serves as an organ-relative coordinate system with a formalism for converting back and forth between model and image coordinates [9, 15]. After deformation of the starting model in a target image, corresponding positions are defined by pairing points in the starting and deformed models that have the same model-relative coordinates (Figure 5). The nonrigid transformation matrix can be specified in image coordinates in terms of a standard displacement vector field [16], where each vector originates on a voxel in the reference image and terminates on its postdeformation position in the target image. In contrast to voxel-scale deformations, the model approach is statistical at the scale of an organ but nonstatistical at voxel scale, eliminating small-scale artifacts [16].
(a)
(b)
(c)
(d)
Segmentation of the prostate in the tCT is necessary to determine corresponding points in the pCT and tCT. The algorithm transfers the atlas model to the tCT using the rigid registration matrix, and autosegmentation proceeds in a statistical framework based on Bayes’ theorem [9]. A conjugate gradient algorithm seeks to find the optimal model such that where is the currently deformed model in the trained shape space [11], is the target image intensity pattern relative to is the probability of (geometric typicality), and is the probability of given (image match [17]).
2.3.4. Mapping the Isocenter
The rationale for using a model to map the isocenter stems from several considerations: during planning the isocenter is positioned relative to the prostate; a point in an image can be more accurately found by relying on regional image features that are correlated spatially with the point rather than using local information near the point itself [18]; and the trainable models used in this study provide a means for determining both the rigid and nonrigid components of the mapping transformation. The mapping step is straightforward and involves labeling the point in the tCT that has the same prostate-relative coordinates as the planned isocenter using (2): where models in pCT and tCT, respectively, and are corresponding positions in - and -relative coordinates, VAL = value of scalar, for example, label, dose, or intensity, in the pCT at -relative position , and mapped to -relative position .
2.4. Reference tCT Isocenters
The clinical and automapping protocols yielded two independent sets of tCT isocenters. To compare these sets, a reference isocenter was determined for each tCT by repeating the automapping procedure with human supervision. The main task was to reposition and/or edit the autosegmented prostate as necessary to achieve the best match with the tCT image data while generally preserving the global shape and volume defined by the planning contours. Every rigid registration and prostate segmentation were evaluated and edited based on human judgment, after which the tCT isocenter from (2) was accepted without modification. This procedure was performed without the pressure of clinical time constraints over ~8 months by two physicists and a dosimetrist working as a team. In general, the dosimetrist and one of the physicists made the initial pass, and the results from that pass were reviewed at a later time by the second physicist. The team met about once a week to discuss and review ongoing progress. Results from the dosimetrist/physicists team were periodically evaluated by one or two radiation oncologists based on the criterion of clinical reasonableness; that is, given the planning contours as truth, would the radiation oncologist judge the location and shape of the prostate in the tCT to be clinically reasonable? This criterion eliminated a potential source of interobserver bias and variability and effectively served the need for clinical standards. Only a few tCTs per patient were reviewed because of the large number of cases and the fact that the position and shape of the prostate in a given tCT are expected to be strongly correlated with other tCTs for the same patient.
2.5. Calculation of Hypothetical Delivered Dose Distributions
Treatment dose distributions for each tCT and each mapped isocenter were computed using PlanUNC assuming all beams were delivered as planned for each of the three isocenters (Figure 6). This was accomplished by importing each tCT to PlanUNC, registering the ensemble of planned beams to each of the treatment isocenters in turn, and calculating the delivered dose assuming all beams were delivered as planned. After dose calculation the DICOM-RT files for the planned and treatment dose distributions were imported to ARTSuite and dose distributions were mapped from tCTs to the pCT using (2) in the reverse direction. Except for a rind ~2–5 mm thick around the prostate, dose to interstitial tissues was not mapped since (2) applies only to modeled tissues. Dose to the rind was mapped by extrapolating the skeletal framework a small distance beyond the prostate surface. The mapped treatment doses were summed over all tCTs and all isocenters and resampled to the grid of the planning dose to simplify comparison of planned and treatment dose distributions. Figure 6(a) shows the planned isodose curves for the 10, 25, 50, 85, 90, and 95% levels on the pCT for the initial (nonboost) portion. The proximal two slices (6 mm length) of the seminal vesicles were included in the PTV. Figure 6(b) shows the same isodose levels from the cumulative dose for the clinical protocol over 15 fractions. Error dose distributions (Figure 6(c)) were computed by subtracting the scaled planned dose distribution from the summed treatment dose distribution. The scaling factor was determined by the number of fractions contributing to the summed treatment dose. Isodose curves are shown for −100, −50, and −25 cGy. Figure 6(d) shows the same isodose levels as in Figure 6(a) computed from one of the 15 tCTs. The frequency distribution in Figure 7 illustrates the expected underdosing of the prostate due to differences between prostate shape and location at planning and treatment times.
(a)
(b)
(c)
(d)
For dose accumulation purposes, the spatial accuracy required for nonrigid registration of a point in the tCT with its corresponding point in the pCT depends on the tolerable dose error and the steepness of dose gradients near the points. Assuming that the spacing between calculation points in the dose grid is matched to the dose gradients [19, 20], correspondence errors should be small compared to the grid spacing. The clinical planning grid spacing was 5 mm in this study. Figure 8 shows color-coded maps of differences between the positions of two points, one point for each method, resulting from a single point in the tCT. The differences were ~1 mm along the and axes and ~2 mm along the axis for the prostate ROI. The larger difference along the axis is attributed to the 3 mm slice thickness. These findings support the use of m-reps for dose accumulation for this study.
(a)
(b)
(c)
3. Results
3.1. Distance Metrics
Frequency histograms (Figure 9) were computed from the signed differences (Figure 10) between the reference isocenters and the clinical and automatic isocenters for all 243 tCTs. The bin width for and was 0.5 mm. To maintain comparable counts per bin, the width for was chosen to be 1.0 mm. Summary statistics for all ten patients are given in Table 2. Figure 9 shows that the distributions for clinical and automatic isocenters are centered near the reference values in approximately Gaussian fashion, supporting the utility of the reference values.
(a)
(b)
(c)
(a)
(b)
Comparison of the frequency histograms for clinical and automatic protocols shows that automatic mapping is robust, as accurate as the clinical protocol along all three axes, and more precise than the clinical protocol, where accuracy is the average of the unsigned compared to the reference values, and precision is the spread (standard deviation) for each axis. Even though the averages for the automatic method are smaller (closer to the references) than the clinical protocol, greater accuracy is not claimed because the references are not golden. On the other hand, the standard deviation is characteristic of the registration method and independent of the reference values. The ANOVA -test is < for all three axes, demonstrating that human and automatic variances are significantly different. As seen from Table 2, these observations apply for each individual patient. Moreover, the clinical protocol yielded 8 outliers, defined here as differences > ( = standard deviation of clinical protocol), with the largest being almost . In comparison, the largest difference for the automatic approach is slightly less than along the axis for patient 10, and for this case is small (≈1 mm).
Table 2 also gives results for Court and Dong [12] and Smitsmans et al. [13]. Both studies evaluated automatic localization of the prostate in tCTs via multiscale rigid registration with pCTs. The values given in Table 2 for these studies are differences in prostate position, as opposed to isocenter position, between automatic methods and manually prepared references. Court and Dong reported results for two patients: patient A had 22 tCTs and was considered less challenging than patient B, who had 21 tCTs. Smitsmans et al. looked at a collection of 19 patients with 8–13 tCTs each. The results in Table 2 are for 91% of the tCTs. The remaining 9% were considered outliers and quantitative results were not reported. The automatic rigid registration method presented in this paper appears qualitatively to be comparable in performance to both Court and Dong and Smitsmans et al. No direct quantitative comparisons are possible however due to differences in study designs.
3.2. Prostate DVHs
3.2.1. Initial Portion of Treatment Regimen
Prostate DVHs for doses accumulated using the clinical, automatic, and reference isocenters for two typical patients are compared against the scaled planned DVHs in Figure 11. These cases illustrate (i) small differences among the three isocenters, (ii) degradation of the shoulder region, and (iii) a decrease in delivered versus planned dose of ~100–200 cGy scaled to the full nonboost portion. These general findings were consistent across all ten patients, but the severity of shoulder degradation was patient specific. Automatic mapping performed better in the shoulder region (Figure 11(a)) than the clinical protocol in about half the cases and as well in the other half (Figure 11(b)). The absence of more significant differences is attributed to margins (5 mm) that were relatively insensitive to isocenter mapping variations on the order of a few mm, a finding expected for properly designed margins.
(a)
(b)
Figure 7 shows the differential DVH for the error dose for the patient in Figure 11(a). Error isodose curves in Figure 6(c) show that the underdosed region occurs at the apex. The superior shift of the delivered dose distribution displayed on the tCT in Figure 6(d) was also present on other tCTs, explaining the origin of the under dose. Underdosing at the base was not observed for this patient because the planning target volume was enlarged superiorly to include the proximal SVs, providing extra protection at the base. However, underdosing was observed at both the base and apex for all patients whose SVs were not included in the target volume.
3.2.2. Boost Portion of Treatment Regimen
TCTs for boost fractions were available for only two patients. DVHs in Figure 12 were computed from tCTs acquired for four of twelve boost fractions. The DVHs for both patients show the same general trends observed for the initial treatment portion (Figure 11). However, the clear separation between the descending portions of the clinical and automatic isocenters in Figure 12 suggests that automatic mapping may reduce the overall prostate under dose compared to the clinical protocol. If true, this finding would not be surprising since the boost margin (3 mm) would be expected to be more sensitive to image registration errors. However, further study is needed to test this finding.
4. Discussion and Conclusions
This work presents a formalism for mapping the planned isocenter to a tCT based on correspondence properties of a deformable prostate organ model that is used for registration of the pCT and tCT and for segmentation of the prostate in the tCT. The fully automatic mapping algorithm is as accurate as the clinical protocol but more precise. The algorithm had no failures or outliers for the tCTs studied. Better precision can be explained in terms of the robust properties of the algorithm and the absence of intra- and interuser variabilities. Moreover, human registrations were made under the pressure of clinical time constraints that can hasten decisions and lead to suboptimal results.
In dosimetric comparisons for the prostate, automatic mapping showed less degradation in the shoulder region of DVHs for 50% of patients in this study. In the other 50%, degradation was no worse than the clinical protocol. The absence of large differences in dose-volume metrics is attributed to prostate margins that were relatively insensitive to variations in isocenter mapping on the order of a few mm.
The clinical significance of the observed dosimetric improvements was not addressed but appears to be modest for the patients in this study. This conclusion however depends on the prostate margin and dose fractionation scheme as suggested by Figure 12. In particular, dosimetric improvements might be significant for less forgiving forms of treatment delivery such as stereotactic body radiotherapy. Also the finding that accurate pCT and tCT image registration does not fully compensate for geometric variability supports conclusions of other studies [2, 3] that full compensation for patient-specific geometric changes requires off-line adaptive planning or online reoptimization.
The overall conclusion is that the automatic algorithm robustly maps the planned isocenter to a position close to the correct location in a tCT and thus is well suited to augment human judgment in the clinical setting. Furthermore, the algorithm offers the potential for reducing registration outliers. During the workflow for mapping the isocenter, all of the essential image processing steps for calculation of delivered dose and mapping the delivered dose to planning space are performed.
Conflict of Interests
Edward Chaney is an employee of Morphormics, Inc., and Principal Investigator of NIH Grant no. R44 CA141941 that supported the research reported in this paper. Mark Foskey is an employee of Morphormics, Inc. Andrea Kress and An Qin were partially supported under subcontracts funded by NIH Grant no. R44 CA141941. The remaining authors declare no conflict of interests.
Acknowledgment
The research reported in this paper was supported by the National Cancer Institute of the National Institutes of Health under Award no. R44 CA141941. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.