Automated high-dose rate brachytherapy treatment planning for a single-channel vaginal cylinder applicator

In the conventional manual planning approach at our institution, CT images are imported to the BrachyVision 13.6 (Varian Medical System, Palo Alto, CA) treatment planning system (TPS) from which the planner creates a treatment plan. To automate this process, we developed the infrastructure illustrated in figure 1, which consists of three main components. The first component is a data server that stores CT image files and the plan files that are generated based on them. The second component of the system is a web-based user interface that allows users to browse the patient list and select a patient, as well as to specify treatment parameters such as the prescription dose, treatment depth, treatment length, and cylinder diameter. The third component is an in-house developed MATLAB-based computation module that performs treatment planning.

Zoom In Zoom Out Reset image size

After the CT scan of the patient is exported to the data server, a treatment planner uses the web-browser interface to select the patient and to input the patient-specific treatment parameters. The AutoBrachy software receives treatment parameters from the web server, reads CT image files from the data server, and then automatically performs all relevant treatment planning calculations without additional human guidance. The plan and structure files are then imported into the TPS to verify plan quality before being approved by a physician.

Within the AutoBrachy system, there are four major steps to develop a treatment plan for an HDR treatment with a cylinder applicator. First, a template is generated based on the input planning parameters; the template contains all of the geometrical components of a treatment plan, e.g. applicator volume, relevant points, dwell positions, etc. Second, the cylinder, based on the given parameters, is oriented and positioned to the CT image set using rigid registration techniques. All the geometrical components are transferred to the CT volume by matching the applicator volume in the template to that in the CT image set. Third, the dwell times are determined via a quadratic optimization method. Fourth, the automatically planned geometric and dosimetric data are transferred to the clinical TPS via a DICOM-RT interface. These four steps are described in detail in the subsequent sections.

We built a web-based user interface, as shown in figure 2. The user can access this webpage through a standard browser to perform treatment planning. There are two afterloaders at our institution. The user first selects an afterloader on the top left of the page. Treatment parameters including prescription dose, treatment depth, length, and cylinder diameter are entered in the lower left region. On the right side of the interface, the user selects the patient.

Zoom In Zoom Out Reset image size

Per our institutional protocol following American Brachytherapy Society guidelines for HDR brachytherapy treatment with vaginal cylinder applicators (Small et al 2012), a treatment plan contains the following geometry components: a clinical target volume (CTV), relevant reference dose points (apex point, dose calculation point, and vaginal surface point), cylinder central channel, and dwell positions. The CTV is defined as the expansion volume from the cylinder applicator for a specified treatment depth and length. The apex point is defined as the tip of the cylinder applicator. The location of the dose calculation point is chosen to be on the CTV outer surface half way from the apex point to the total treatment length. The vaginal surface point is on the cylinder applicator surface at the same level as the dose calculation point.

Given that all geometrical components are defined according to their geometric relations with respect to the cylinder applicator, a template is generated containing the applicator and all other components. Specifically (1) A 3D template model of the cylinder is created with the principal axis along the z direction (purple surface in figure 3(A)) based on the user-specified cylinder diameter; (2) The CTV is generated by expanding the cylinder volume along  +Z and radially from the central axis according to the user-specified treatment depth and length (yellow surface in figure 3(A)) and then subtracting the original cylinder volume; (3) Locations of the dose calculation, apex and vaginal surface points are determined according to their geometric relations with respect to the cylinder applicator (figure 3(A)); (4) The applicator channel is drawn along the central axis of the cylinder, as shown by the green line in figure 3(B). In our practice, segmented cylinder applicators are used. The central channel starts 3.43 mm from the tip of the cylinder (Varian Medical Systems 2013); and (5) Starting from 4 mm away from the applicator channel tip, consecutive dwell positions are placed to cover the targeted treatment length. Successive dwell positions are separated by 5 mm according to the Varian VS2000 source used in our clinic (red bars in figure 3(B), and red dots in figure 3(C)).

Zoom In Zoom Out Reset image size

To place the cylinder template in the patient CT image set, a rigid registration process is automatically performed. After that, all other geometric components that were generated in the previous section are transferred from the template to the CT volume. The cylinder applicator template in figure 3 should be oriented and translated within the 3 dimensions to match with the actual cylinder applicator in the patient CT images. In principle, this can be achieved solely by calling a rigid registration function with the template and full CT image set. However, we found that simply calling the rigid registration function on the whole CT volume is slow due to the iterative nature of the algorithm. In addition, without appropriate initialization of the template position and orientation, it is difficult to achieve accurate registration results. To overcome these challenges, an automatic procedure to perform the template-matching step was developed. The main idea is to first perform a rough initialization of the template in a smaller sub-region of the CT volume and then call the registration function to fine tune the results; these methods are described below.

2.4.1. Region of interest extraction.

The first step is the extraction of a region of interest (ROI) that fully encompasses the cylinder applicator from the original CT images. Extracting the ROI reduces the problem size and helps increase the computational efficiency of the subsequent processing steps. A typical patient CT image in a sagittal view is shown in figure 4(A). The applicator central channel is clearly visible as a line of highlighted voxels because a radio opaque wire was placed in the applicator at CT simulation stage. To automate this segmentation, the voxels with the highest intensity are identified via a threshold of 2000 hounsfield units (HU); with a properly calibrated CT simulation, this threshold limits the results largely to metallic or radiopaque objects such as the guide wire and markers. Detailed discussion regarding this issue will be presented later. Typical results are presented in figure 4(B) where the result contains a group of voxels lying on a straight line corresponding to the central applicator channel on the CT image, as well as other voxel clusters indicating radiopaque markers.

Zoom In Zoom Out Reset image size

To distinguish the true central channel voxels automatically from other high-HU objects, spectral cluster analysis is performed (von Luxburg 2007). The number of eigenvalues equaling zero for the Laplacian matrix based on the distances between voxels of interest provides an estimate of the number of unique cluster groups. The voxels of interest are then clustered into this number of groups. The result of this step is shown as dots with different colors in figure 4(C). To automatically identify which cluster is the central channel guide wire, each cluster is given a score $s=N{{\left(\frac{d}{a}\right)}^{2}}$, where $N$ is number of voxels in the group, $d$ is maximum separation of points within a cluster group, and $a$ the shortest distance from the cluster to the center of the CT simulation. The definition of this score $s$ is empirically chosen, but we found it performed well. The cluster with the highest score is identified as the central channel.

After the central channel has been detected, the minimum bounding box for the voxels containing the central channel is computed. The maximum cylinder radius is 17.5 mm and the angle of the cylinder applicator central axis from the patient superior-inferior (SI) direction is generally less than 45°. Hence, expanding the bounding box by 30 mm in patient lateral directions yields  >3.5 mm padding radially for the largest cylinder for the ROI. This ensures the cylinder is fully encompassed within the ROI. Since it is critical that the whole length of the cylinder is fully encompassed for the subsequent automatic planning sections, in the SI direction the bounding box is expanded by 50 mm in both directions. Using a threshold of 150 HU, we render the boundary of bones, additional metal markers and the whole applicator inside the newly created ROI, as shown in figure 4(D).

2.4.2. Determination of cylinder orientation.

2.4.3. Determination of cylinder translation.

The next step is to perform a translation on the rotated cylinder template to accurately match it with the true cylinder location. The template is first approximately translated by matching the tip location of the applicator channel in the template with the first (most superior) voxel on the cylinder central channel in the CT image. Then a rigid registration function is launched to finely tune the template location by matching edges detected on the template image and the segmented CT image (as shown in figure 4(D)). Finally, another 1D translational match is performed along the cylinder orientation direction to match the cylinder tip in the template with that in the CT image. The tip of the cylinder in the CT image is identified as the location where the CT number gradient reaches its maximum value along this axis. Theoretically, due to the finite separation between the applicator channel tip and the cylinder tip (3.43 mm for our applicator type) there should be two drops in CT number along the axis: one from the applicator channel to the cylinder and one from the cylinder to the tissue. The location of the second drop is the correct cylinder tip position. However, due to the finite CT image resolution along the SI direction (typically 2 mm), the two drops are blurred and we can only observe one distinct CT number drop. Matching the highest gradient point with the applicator tip was found to be acceptable in our tested cases.

Each of these steps together yields a well-matched cylinder template with the cylinder applicator in the CT image, as shown in figure 5(A). Meanwhile, the translation vector, $T$, and the rotation matrix, $R$, that are determined in this process define a transformation from the template to the CT images coordinates. By applying this transformation, all other geometric components in the template are transferred to the CT images (figure 5(B)).

Zoom In Zoom Out Reset image size

The dwell time of a treatment plan is determined by solving an optimization problem (Xing et al 1999). The planning objective of this cylinder treatment is to cover the CTV surface uniformly with the prescription dose $p$ (Small et al 2012). As such, a set of $m$ anchor points on the surface of the CTV are sampled first (blue stars in figure 3(C)), and a dose vector, $b$, has the dose at each anchor point set to the prescription dose. Due to rotational symmetry, it is only necessary to place the points on one side of the CTV surface. The separation of these points is 5 mm along the cylindrical portion of the CTV, and then at 15, 30, 45, and 90° for the hemispherical portion of the CTV. The dose deposition coefficient matrix $D$ is computed following the task group report #43 of the American Association of Physicists in Medicine (Rivard et al 2004). Its element ${{D}_{ij}}$ represents the contribution from the ith dwell position to the jth anchor point for a unit dwell time. The vector of non-negative dwell times that is to be determined in this optimization problem is denoted by $x$. A quadratic programming approach is employed to determine the dwell time of each dwell position by solving an optimization problem as:

$$\begin{eqnarray}&&{{\min}_{x}}\,\left|Dx-b\right|^{2},~\text{s}.\text{t}.\,Dx-b\geqslant 0,\,x\geqslant 0\end{eqnarray} \tag{ 1 }$$

where $|.{{|}^{2}}$ denotes the standard L-2 norm. This optimization problem is solved using the interior-point algorithm (Bazaraa et al 2006).

The plan and structure information including all structures, points of interest, channel positions, dwell positions and dwell times are written into standard DICOM-RT files. No private tags are used to ensure compatibility with multiple treatment planning systems. This completes the automatic treatment planning process. These files are imported into FDA approved TPS. After a physician draws other organ contours for treatment planning per our planning protocol, a physicist inspects the plan and makes changes to the dwell times if the dose to OARs is above acceptable limits.

To test our system, nine patients previously treated at our institution were selected. These patient cases had different treatment prescriptions, lengths, depths, and cylinder diameters and hence were expected to represent a large patient population. The prescription doses ranged from 5.25 Gy to 6 Gy, treatment lengths ranged from 2 cm to 5 cm, the treatment depths were either 0 mm or 5 mm, and cylinder diameters were 2.5 cm or 3 cm. For each case, we generated a treatment plan using the developed module. Each case also had a clinical plan that was generated manually by physicists covering clinical HDR service. Specifically, the physicist first selected a cylinder plan template from our clinical TPS according to the applicator size, prescribed treatment length and depth. The location and orientation of the applicator were adjusted to match that in the CT images. The template contained predefined dwell times, which were scaled according to the prescription to yield the final dwell times in the plan.

The automatically generated and manually generated plans were compared. To evaluate geometric accuracy, the identified applicator tip positions and central channel orientations in the automated plans were compared with the corresponding manual plans. To evaluate the dosimetric quality, comparisons between automated and manual planning were performed in terms of the minimum dose received by at least 5% (D5) and 95% (D95) of the CTV surface, as well as the dose homogeneity index (DHI), defined as $1-(D5-D95)/p$, where p is the prescribed dose. A DHI closer to unity indicates a better dose homogeneity on the CTV surface. A dose surface histogram (DSH) was also generated for the outer CTV surfaces in both automated and manual plans. We would like to remark that the use of DSH was for the purpose of examining the homogeneous coverage of the outer surface of the CTV, which was the goal of treatment planning in this case per American Brachytherapy Society guideline(Small et al 2012).

Table 1 shows the 3D distances between the automatically and manually digitized tip positions and the angle between the automatically and manually digitized central channels. The mean 3D distance between tips was found to be 0.75  ±  0.39 mm (typical voxel dimensions are $1.17~\times ~1.17~\times ~2.0~$ mm³). The mean angle between central channels was found to be 0.69  ±  0.43°. These results demonstrate that the geometric difference between the two planning methods is small.

Figure 7 shows D5, D95 and DHI of the two planning methods for each case with the prefix 'm' and 'a' for the manual and automatic plans, respectively. In seven out of the nine cases (Cases 2–5, 7–9), the automated plans achieved DHI greater than the manually generated plans. In Case 1, the DHI for the manual plan was slightly higher, however, this plan underdosed the CTV surface compared to the automated plan, as indicated in the DSH shown in figure 8(A). In Case 6, the DHI values differ by 0.001 for the manual and automatic planning. A sample case (Case 2) is illustrated in figure 8(B). For Cases 6 and 8, the manual plan overdosed the CTV surface. Case 8 is illustrated in figure 8(C).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The nine patients chosen for this study were meant to represent a large range of parameter space to test the automatic planning module. The patient data have not been sorted according to any pattern. We used these patient cases as a proof of concept to show that the module works from end-to-end over a range of applicator diameters, orientations, CTV depths, treatment lengths, and prescription doses. However, as an evaluation on the module for clinical usage, this number of patients is probably insufficient. We are currently in the process of performing a test retrospectively on all cylinder applicator patients treated in our institution for the previous 2.5 years (~75 patients) for a comprehensive test of the module.

While the CT value is related to the material and should thus exist in a well-defined range of values for a given material, CT scanner calibration inaccuracy and partial volume averaging effects may lead to wider ranges of values. To make our algorithm robust, we studied the distribution of CT numbers in a number of patients. A typical histogram of the extracted ROI is shown in figure 9. A peak in the distribution related to the metallic objects was observed at 3000 HU, but to account for the partial volume effects and calibration variation, we found 2000 HU was a good choice of threshold to separate the metal from other structures. Hence this number was used in our algorithm for the ROI extraction (section 2.4.1). Similarly, the choice of 150 HU, used to discriminate the applicator from surrounding soft tissue in the image registration after the ROI was selected (section 2.4.3), was chosen because it consistently represented a local minimum in the distribution above the expected CT value range for tissues.

Zoom In Zoom Out Reset image size

The cylinder applicator planning module was written for CT images. However, it can be potentially expanded to work with other planning image modalities, such as MRI. When using a different image modality, the most critical steps are to identify the applicator channel and the cylinder edge. This may be achieved by analyzing the image intensity distribution and determining appropriate thresholds. With new thresholds, the applicator channel can be determined and the rest of the algorithms will perform in the same manner as they did for CT images.

In our automatic planning module for cylinder applicators, the dwell times are determined using an iterative quadratic optimization approach. However, this is not strictly necessary. Since there are a finite number of combinations of treatment parameters (e.g. treatment length, depth, and cylinder applicator diameter), it is possible to enumerate them, and to add the optimal dwell times for each of them to the template. Nonetheless, the automatic planning module is a testbed for future, more complex automatic planning modules, so we chose to solve the optimization problem each time, because these techniques will be required for the other types of applicators.

Solving the optimization problem was not found to be a significant computational burden in the current study. For the cylinder applicator case, the problem size is small. This is because of the simple planning objective. Typically used constraints in more complex cases, e.g. organ dose-volume constraints, are not needed either. In addition, the rotational symmetry of the applicator also allowed using a small number of anchor points.

We have observed improvements in DHIs for the CTV surfaces in automatically generated plans (figure 7, bottom) even though the manually generated plans were also clinically acceptable and previously delivered to patients. The difference between manually and automatically generated plans can be ascribed to the plan quality variation in the templates used in the manual planning process. In fact, when a template was generated, the physicist manually tuned the prescription iso-dose line to match with the CTV surface. Because of the sharp dose gradient, a slight misalignment between the two led to large difference in dose. When planning a patient case, the manual plan was essentially a scaled version of the template, while the automatically generated plan strictly enforced dose level at the CTV surface. Hence, dosimetric differences could be observed when comparing DSH.

Likewise, both the manually and automatically created plans were acceptable in terms of geometrical accuracy, although differences between the two were seen. There was a case with relatively large discrepancies between apex point and central channel orientation, case 7 (table 1). After re-examining this case, we found that the guide wire was slightly off the center of the cylinder (<1 mm). The manual plan digitized the actual guide wire, whereas the automatic plan used the center of the cylinder as the channel position. Therefore, a slight difference was observed. However, in terms of channel position, both the automatic and manual plans were of sufficient quality to be acceptable for clinical treatment.

As well as extending the methods for additional planning image modalities and adding modules for more complex applicators, there exists the potential to integrate these methods directly into TPS. The current AutoBrachy system is a standalone program that communicates with the clinical TPS through the DICOM-RT interface. Since no non-standard fields are used with AutoBrachy, the plans and structure files should be compatible with multiple TPS suites. While keeping AutoBrachy as a separate module is clinically feasible, it introduces an additional step of manually importing plan data to the TPS. It is desirable to incorporate the automatic planning module within the TPS. At present, we are in the process of developing interfaces with the BrachyVision TPS via the Eclipse API. Upon completion, this will allow a user to call the automatic planning module from BrachyVision and the results will be available inside the TPS immediately. It is expected that this step will further improve the workflow, making the AutoBrachy module more suitable for clinical applications.

To the best of our knowledge, this is the first study describing a CT-based automated treatment planning for HDR brachytherapy. Our study has demonstrated that auto-planning for HDR brachytherapy is feasible and clinically attractive. There is potential to expand our AutoBrachy platform to accommodate more complex cases utilizing other brachytherapy applicators, e.g. tandem-and-ovoid (T&O) and tandem-and-ring (T&R) treatments.