Computer-assisted delineation of lung tumor regions in treatment planning CT images with PET/CT image sets based on an optimum contour selection method.

To assist radiation oncologists in the delineation of tumor regions during treatment planning for lung cancer, we have proposed an automated contouring algorithm based on an optimum contour selection (OCS) method for treatment planning computed tomography (CT) images with positron emission tomography (PET)/CT images. The basic concept of the OCS is to select a global optimum object contour based on multiple active delineations with a level set method around tumors. First, the PET images were registered to the planning CT images by using affine transformation matrices. The initial gross tumor volume (GTV) of each lung tumor was identified by thresholding the PET image at a certain standardized uptake value, and then each initial GTV location was corrected in the region of interest of the planning CT image. Finally, the contours of final GTV regions were determined in the planning CT images by using the OCS. The proposed method was evaluated by testing six cases with a Dice similarity coefficient (DSC), which denoted the degree of region similarity between the GTVs contoured by radiation oncologists and the proposed method. The average three-dimensional DSC for the six cases was 0.78 by the proposed method, but only 0.34 by a conventional method based on a simple level set method. The proposed method may be helpful for treatment planners in contouring the GTV regions.

INTRODUCTION

The primary purpose of radiation therapy is to administer the highest possible dose to tumors while minimizing the dose to normal tissues and organs at risk (OARs), such as the spinal cord, esophagus and heart. The stereotactic body radiotherapy (SBRT) technique has recently been developed to meet this goal of maximal and minimal doses to the tumor and normal tissue, respectively. In performing radiotherapy, the gross tumor volume (GTV), which is defined by the gross demonstrable extent and location of a malignant growth [1], should be accurately determined, because the prescribed dose distribution in radiation treatment planning is determined individually for each tumor region. Therefore, GTV regions have been manually delineated by radiation oncologists using treatment planning computed tomography (CT). However, the subjective manual contouring of a tumor region is tedious and time-consuming, and its reproducibility is relatively low, which could cause inter- and intra-variability of tumor regions [2-4]. Consequently, it would be not easy for many radiation oncologists to accurately delineate the GTV regions with high reproducibility and low variability.

Computer-assisted delineation methods can assist radiation oncologists in overcoming these limitations of manual contouring [3, 5–7]. In several computer-assisted methods, both morphological and functional images (such as those from CT and 18F-fluorodeoxyglucose positron emission tomography (FDG-PET)) have been employed for improving lung tumor delineation. CT is the essential modality for radiation treatment planning, because a CT image contains anatomical information and it is employed to obtain the electronic density map that is used for dose calculation. On the other hand, PET is a functional imaging modality that delineates the active tumor region, and thus it has higher sensitivity and specificity for detecting tumor tissue than CT images. El Naqa et al. proposed a multivalued level set method (LSM) for GTV extraction by using active contours for radiotherapy treatment planning [5]. Day et al. reported a region-growing method for tumor volume segmentation on FDG-PET images for rectal and anal cancer patients. In their method, the variation in the GTV was 9%, versus 37% using a simple thresholding method [3]. Strassmann et al. developed an atlas-based method for semi-automatically extracting tumor regions for head and neck tumors, and this method cut planning time approximately in half [6]. Ballangan et al. reported a region-based active contour approach for improving the accuracy of tumor delineation by using both CT and PET images for cases in which tumors abut or involve the chest walls or mediastinum [7].

This study focuses on selection of a global optimum object contour from among multiple possible delineations around a tumor. For that purpose, we have proposed an automated contouring algorithm based on an optimum contour selection (OCS) method using a LSM for treatment planning CT images combined with PET/CT images. The OCS method retrospectively determines a global optimum objective contour from among multiple active delineations around a tumor. In addition, the PET images were employed for determination of initial GTV regions to be used in the OCS method.

MATERIALS AND METHODS

This study was performed with approval by the Institutional Review Board of the Kyushu University Hospital. We employed a computer with four cores of a 3.2-GHz central processing unit (CPU) (Intel i7-3930) and a 16-GB memory on a linux system (Ubuntsu) on a Windows 7 operating system. The overall scheme of the proposed method for automated contouring of a GTV is shown in Fig. 1. First, the PET image was registered with the planning CT image through a diagnostic CT image of the PET/CT dataset [8, 9]. Initial GTV regions were obtained by thresholding the PET image at 80% of the maximum standard uptake value (SUVmax) within a rectangular volume of interest (VOI), which had the same geometric position as the VOI in the planning CT image. Each initial GTV location was corrected in the VOI. Finally, the GTV region was segmented using the OCS method.

Fig. 1.

An overall scheme of the proposed method for automated contouring of a GTV based on an optimum contour selection (OCS) method.

Clinical cases

Datasets consisting of planning CT and PET/CT images of six lung cancer patients (mean age: 74 years; range: 65–86 years; females: 3; males: 3; mean effective diameter of GTV: 23.8 mm; range: 17.7–29.4 mm) who had received SBRT were selected for this study. The patient characteristics are summarized in Table 1. All tumors were solid type cancers. Planning CT images of the patients were acquired with breath-holding at the end of expiration to delineate the tumor region and to calculate the dose distribution for each patient with a radiation treatment-planning (RTP) system, whereas PET/CT image sets were obtained to help the treatment planners delineate the tumor region from a functional point of view.

Table 1.

Summary of patient characteristics

Patient No.	Gender	Age (years)	GTV sizea (mm)	Tumor location	SUVmaxb	Time differencec (days)	Tumor CT imaging characteristics
1	female	71	17.7	RULd	8.43	20	homogeneous irregular
2	female	81	25.8	RUL	4.43	2	homogeneous irregular pleural indentation
3	male	65	25.3	RUL	6.79	20	inhomogeneous irregular
4	female	67	24.2	RUL	12.2	34	homogeneous irregular vascular
5	male	75	20.2	LULe	8.74	5	inhomogeneous irregular adjacent pleural
6	male	86	29.4	LUL	9.68	20	cavity irregular

aEffective diameter. bMaximum standardized uptake value. cTime difference between the planning CT scan and PET/CT scans. dRight upper lobe. eLeft upper lobe.

Planning CT images were acquired from a four-slice CT scanner (Mx 8000; Philips, Amsterdam, The Netherlands) with dimensions of 512 × 512 pixels, an in-plane pixel size of 0.977 mm, and a slice thickness of 2 mm. The original planning CT images were converted to isotropic images with matrix sizes of 512 × 512 × 234–313 and an isovoxel size of 0.977 mm using a cubic interpolation method.

PET imaging was performed on an integrated PET/CT scanner (Discovery STE; General Electric Medical Systems, Milwaukee, WI). Patients were scanned with their arms down and breathing free 60 min after the FDG injection. The PET data were acquired in the 3D mode and reconstructed with corrections for attenuation, scatter, decay, random, and dead time using a 3D ordered subsets-expectation maximization algorithm (VUE Point Plus; GE Healthcare), resulting in a 128 × 128 matrix with an in-plane pixel size of 5.47 mm and a slice thickness of 3.27 mm. Meanwhile, diagnostic CT images were acquired from a 16-slice CT scanner (a 512 × 512 matrix, an in-plane pixel size of 0.977 mm, and a slice thickness of 5 mm) in the PET/CT system.

The matrix size and isovoxel size of the PET image and diagnostic CT image in a PET/CT dataset were conformed to those (512 × 512 × 234–313 and 0.977 mm) of the planning CT image for registration of the PET image to the planning CT image through the diagnostic CT image. After having changed the voxel sizes, the matrix sizes of the PET and diagnostic CT images were 716 × 716 × 793–920 and 512 × 512 × 234–313, respectively. Because the isotropic PET images were larger than the isotropic diagnostic CT images, the PET images used for this study were cropped from the isotropic PET image by use of a template-matching technique with normalized mutual information [10], in which the diagnostic CT images were used as templates. In this template matching, image positions of the PET and diagnostic CT images were used as initial points.

The SUV was employed for identification of initial GTV regions, and it was calculated as a ratio of the radioactivity concentration of tissue at one timepoint to the injected dose of radioactivity concentration at that timepoint, divided by the body weight [11]: where C represents the radioactivity concentration in kBq/ml obtained from the pixel value in the PET image multiplied by a cross calibration factor, D is the injected dose of 18-FDG administered in MBq (decay corrected), and W is the body weight of a patient in kilograms.

The treatment-planning CT images were obtained from digital imaging and communications in medicine (DICOM) files with the personal information removed. In accordance with the definition of JCOG0403, GTV contours were determined based on a consensus between two experienced radiation oncologists with reference to the PET/CT images using an RTP system (Varian, Eclipse Ver. 6.5).

Two-step registration of PET image to planning CT image

The CT scan has the advantage of providing high-resolution images with anatomical detail, but it is not good at giving functional information about tumors. In contrast, PET images reflect the functional processes in biological bodies. To utilize both functional and anatomical information in the same coordinate system, the PET images were registered with the planning CT images.

Figure 2 illustrates a two-step registration of a PET image to a planning CT image. In the first step, a diagnostic CT image of PET/CT was registered with a planning CT image by using an affine transformation, which was calculated by manually selecting nine feature points. These points are the apexes of the left and right lungs, the most anterior points of the left and right lungs, the most posterior points of the left and right lungs, the most lateral points of the left and right lungs, and the most superior point of the diaphragm from treatment-planning CT images and diagnostic CT images of PET/CT image sets. Twelve elements in the affine transformation matrix were calculated using a least squares method based on a singular value decomposition. In the second step, the PET images were registered to the treatment-planning CT images by using the same transformation matrix, because the PET image was scanned using the same coordinate system as used for the diagnostic CT image of PET/CT. However, the two-step registration was performed in terms of the left and right lungs, rather than with reference to a lung tumor. Therefore, an initial lung tumor identified on the PET image should be locally aligned with the planning CT image, as described in a later section.

Fig. 2.

An illustration of a two-step registration of a PET/CT image set to a planning CT image.

Identification of initial GTV regions

Initial GTV regions were determined by thresholding the PET images by a certain percentage of the SUVmax, because brighter tumor regions in PET images indicate regions where tumor cells may be active.

First of all, a rectangular VOI, which was slightly larger than a circumscribed parallelepiped of a tumor, was determined in the planning CT image by manually selecting the minimum and maximum coordinates of the circumscribed parallelepiped, and calculating its widths in the x, y and z directions. The PET image was thresholded at a certain percentage of the SUVmax within the VOI to identify the initial GTV region. We can assume that the region of interest (ROI) size and location would not affect the segmentation results, because the initial GTV regions were determined based on the tumor SUV, which was independent of the ROI.

In a preliminary study, the suitable threshold percentages of the SUVmax for segmentation of initial tumor regions on the PET images were determined for all cases by changing the percentage of SUVmax from 30 to 80% as a threshold value in the VOI. Figure 3 shows an example that compares the tumor regions (white regions) segmented by different threshold values and the GTV contours (red lines) determined by radiation oncologists for Case 6. The comparison for all cases revealed that regions segmented at 80% of the SUVmax were suitable as the initial GTV regions for the proposed method, because the initial regions should be inside the desired contours in the proposed method. Figure 4 illustrates the initial regions (red regions) in a VOI. All initial regions were located inside the tumor region in this case, and they should be as large as possible to prevent potential cavities due to inhomogeneous inner regions of lung tumors, and to reduce the calculation time.

Fig. 3.

A comparison between the tumor regions (white regions) segmented by different threshold values and a GTV contour (red line) determined by radiation oncologists for Case 6.

Fig. 4.

An illustration of initial regions (red regions) in a VOI.

Correction of initial GTV location in a ROI of a planning CT image

To register PET/CT images and planning CT images, Cai et al. performed a validation study of CT and PET lung image registration based on a chamfer-matching method, with registration errors of 2–3 mm in the transverse plane, 3–4 mm in the longitudinal direction, and about 1.5 degrees in rotation [12]. Mattes et al. proposed a free-form deformation method, with overall errors ranging from 0 to 6 mm [13]. Figure 5 illustrates the discrepancies in location and shape of tumor regions between the CT and PET images after the two-step registration mentioned above. These discrepancies might be caused by respiratory motion, inhomogeneous radioactivity concentrations in the GTV regions, or the large time intervals between the PET/CT and planning CT scans.

Fig. 5.

Illustrations of discrepancies in location and shape of tumor regions between the CT and PET images after the two-step registration.

To correct the tumor location described above, the centroid of an initial GTV region obtained from a PET image was aligned with the center of the ROI in a planning CT image on a slice-by-slice basis. The ROI was the same as the transverse plane of the VOI used for the identification of the initial GTV regions. In this process, only the ROI region (not the entire PET image) was moved to adjust the tumor location. Figure 6 illustrates the correction of an initial GTV location in an ROI of a planning CT image by alignment of the initial region obtained from a PET image with the center of the ROI in each slice. The location discrepancy between the initial GTV region and the lung tumor was successfully corrected after the alignment.

Fig. 6.

An illustration of the correction of an initial GTV location in a ROI of a planning CT image by alignment of the initial region obtained from a PET image with the center of the ROI on a slice-by-slice basis.

Segmentation of GTV regions using an OCS method

The final GTVs were segmented by applying the OCS method to the initial region as determined on the PET images. The basic concept of the OCS method is to retrospectively select a global optimum object contour from among multiple active delineations with an LSM around tumors. In the OCS method, the LSM [14] is employed for searching for the optimum object contour in the relationship between the average speed function value on an evolving curve and the evolution time.

In the first step, the GTV contour and the speed function value obtained by the LSM were recorded at each evolution time from the initial GTV region until the evolution time reached 10 000 or the evolving curve reached the edge of the ROI in the planning CT image. The level set function was updated from the initial GTV region contour by using the following discrete partial differential equation: where n is the evolution number, t is the evolution time, is the evolution time interval, and is the speed function. The evolution time is the time of the contour deformation in updating the discrete partial differential equation. The zero level set of , which corresponds to the contour of the segmented region, moves according to the speed function in the 3D level set function. The zero level set function, i.e. the evolving curve, moves according to the following speed function : where is the function of the edge indicator, is the Gaussian function, is the planning CT image to be processed, ⊗ denotes convolution, v is a constant, and is the curvature. The edge indicator function and speed would be small around the edge, whereas the functions and would be large in relatively homogeneous regions.

In the second step, the GTV contour was determined from the optimum contour derived using the LSM by searching for the minimum point in the relationship between the evolution time and the average speed function value on an evolving curve, based on the steepest descent method (SDM). To avoid local minimum traps, the average speed function was smoothed by a median filter, and the smoothed function was resampled by a larger interval than the original one, before applying the SDM.

Evaluation of the proposed method

The performance of the proposed method was evaluated using a Dice similarity coefficient (DSC) [15]. The DSC denotes the degree of region similarity between the GTV gold standard region and the GTV region obtained by the proposed method. The DSC was calculated by the following equation: where T is the GTV gold standard region determined by two radiation oncologists, D is the GTV region contoured by the proposed method, n(T) is the number of pixels in the region T, n(D) is the number of pixels in the region D, and n(T∩D) is the number of logical AND pixels between T and D. The DSC ranges from 0 (no overlap between T and D) to 1 (T and D are identical). A DSC value of greater than 0.7, which denotes a good agreement between T and D, has been generally accepted in the medical field [16-18].

The GTV regions were obtained from the DICOM for radiation therapy (DICOM-RT) files. The isotropic GTV region was used as the gold standard, which was produced by use of a shape-based interpolation [19] for matching with the isotropic planning CT image with an isovoxel size of 0.977 mm.

The DSC for each case was calculated in two different ways. First, the average value of the DSCs was calculated based on a 2D slice, yielding the 2D average DSC. Second, the DSC was calculated from a 3D planning CT image to obtain the 3D DSC.

To evaluate the usefulness of the proposed method, a conventional method was compared with the proposed method. In the identification of an initial GTV region step, the centerpoint of each ROI was used as the initial region without using PET images. In the segmentation of GTV regions using an OCS method step, a fast LSM [20] was used instead of the OCS method. In the fast LSM, the threshold for the value of average speed function was manually optimized based on all six cases, and was 0.015 for this study. In contrast, the proposed method adaptively determines the threshold value depending on the GTV regions, which is one of the advantages of the proposed method. Those two methods were fully automatic except for determination of the VOI around the GTV regions.

RESULTS

Figure 7a illustrates the GTV contours, which were delineated a number of times by the proposed method on the planning CT image for Case 6 at evolution times of 500, 2000, 4000 and 6000. Figure 7b shows the relationship in the LSM between the evolution time and the average speed function value on an evolving curve. In Case 6, the average speed function as a function of evolution time t converges to the global minimum between evolution times of 2000 and 4000. At an evolution time of 500, the segmented GTV region is relatively small and corresponds to a larger as a result of the low contrast inside the segmented region. For evolution times of 2000 and 4000, the GTV contours seem to be well extracted and approach the minimum value of due to the tumor edge. At an evolution time of 6000, the estimated GTV contour was too large, and it had a greater value due to the low contrast outside the tumor region. Therefore, the optimum contour can be determined by detecting the minimum point in the relationship between the average speed function and the evolution time.

Fig. 7.

Illustrations of (a) contours on a planning CT image at evolution times of 500, 2000, 4000 and 6000, and (b) the relationship in the LSM between the evolution time and the average speed function on a moving front line.

The resultant GTVs were visually investigated and categorized in terms of the tumor CT imaging characteristics [21], which were determined by a radiologist (HY) with more than 5 years of experience. Figure 8 compares the segmentation results of the proposed method and the conventional method for six types of lung tumors classified according to tumor CT imaging characteristics. For Case 1, the proposed method provided a GTV contour that is closer to the GTV gold standard on the irregular tumor margin. The evolving curve stopped at almost the same position as the GTV gold standard when the proposed method was used, but it was terminated before reaching the GTV gold standard when the conventional method was used, due to a fixed threshold value for the speed function. For Case 2, the proposed method better segmented a homogeneous and irregular tumor with pleural indentation with a higher DSC of 0.89. For the inhomogeneous and irregular tumor region in Case 3, the proposed method showed a higher DSC of 0.84 with the GTV gold standard, whereas the conventional method overestimated the GTV because of the unclear edge of the tumor. For Case 4, the results of both the proposed method and conventional method were comparable for homogeneous and irregular tumors adjacent to vascular tissues. For Case 5 in Fig. 8, neither method was able to correctly extract reasonable GTVs, because the tumor contacted a chest wall, and there appeared to be no boundary between the tumor and the chest wall. To extract this type of GTV, we will need to develop a method for estimating the chest wall boundary in the future. For Case 6, the proposed method extracted only a part of the GTV region with cavities, whereas the evolving curve in the conventional method was terminated inside a very small region of the GTV. Consequently, we need a method for more accurately segmenting initial regions, so that they cover the cavity regions to avoid the effect of low speed values inside tumors.

Fig. 8.

A comparison between results of the proposed method and conventional method in terms of tumor CT imaging characteristics.

Table 2 shows the DSCs between the GTV gold standard region and the GTV segmented by the proposed method and the conventional method for six cases. The average 3D DSC between the GTV gold standards contoured by radiation oncologists and the GTV regions obtained by the proposed method was 0.78 for six cases. This result was higher than the average DSC of 0.34 obtained by the conventional method.

Table 2.

DSCs between GTV gold standard regions and GTV regions segmented by the proposed method and the conventional method for six cases

Method	DSC	Case number						Mean
		1	2	3	4	5	6
Conventional method	2Da	0.30	0.38	0.54	0.49	0.46	0.44	0.44
	3Db	0.40	0.12	0.52	0.08	0.40	0.52	0.34
Proposed method	2D	0.74	0.63	0.80	0.83	0.78	0.65	0.74
	3D	0.81	0.71	0.80	0.83	0.81	0.71	0.78

aMean for 2D DSCs in all slices with GTV contours between the two methods. b3D DSC for whole GTV volumes between the two methods.

DISCUSSION

To evaluate the usefulness of initial regions obtained from PET images, a method without using PET images was compared with the proposed method. In this method, the centerpoint of each ROI was used instead of using the initial regions obtained from PET images, but the other steps of the method were exactly the same as the proposed method. This method differed from the conventional method, which was used for comparison in Table 2. As a result, the mean 3D DSC of this method was 0.56, which was much lower than the 0.78 obtained by the proposed method. Therefore, we believe that the initial regions on the PET images could be useful for segmentation of lung tumors in the proposed method.

In this paragraph, three studies [5, 22, 23] on automatic delineation methods for lung cancer radiation therapy were compared with the proposed method, because they were evaluated using the DSC. El Naqa et al. [5] developed a segmentation method based on CT, PET and magnetic resonance (MR) images by applying a multivalued LSM. The average DSC was 0.90 for a phantom study, not patient studies. In general, MR imaging is not applied for lung cancer patients. Wanet et al. [22] validated a gradient-based segmentation method for performing GTV delineation on FDG-PET of non-small-cell lung cancer patients using surgical specimens in comparison with threshold-based approaches and CT. The average DSC was 0.58 to 0.62 on CT images and 0.62 to 0.68 on PET images. Markel et al. [23] developed an approach based on a combined decision tree with k-nearest neighbors to classify voxels based on their first-, second-, and higher-order statistics as well as structural and Tamura features. The approach was able to segment tumors with an average DSC of 0.607 for 31 lung cancer patients. On the other hand, the average DSC obtained by the proposed method was 0.78 for six patients, which was lower than for the phantom study, but higher than for the patient study.

The mean calculation time of the proposed method was 16 min 27 s for the six cases using the general personal computer. This was less than the mean spending time of the conventional method, which was 27 min 54 s, because the proposed method could terminate the evolution of moving contours when the tumor edges were blurred, but the conventional method could not, as is shown for Case 3 in Fig. 8. Therefore, the proposed method was more efficient than the conventional method. However, the computation time can be reduced with the use of a general-purpose graphical processing unit (GP-GPU) for clinical practice.

Six cases with various types of tumor CT imaging characteristics were selected for this study, but the proposed method should be applied to many types of tumors in the future so as to improve the segmentation accuracy, especially for cavity-type tumors and tumors abutting the chest wall. In addition, independent databases obtained from different institutions and planning CT scanners should be used in order to improve the robustness of the proposed method.

CONCLUSION

We proposed an automated contouring method based on the OCS method for treatment-planning CT images with PET/CT images. In the proposed method, initial GTV regions were determined from 80% of SUVmax on the PET images, and then global optimum GTV contours from among multiple possible delineations around lung tumors were selected. The proposed method was more accurate and efficient than the conventional method, and thus may be useful for assisting radiation oncologists in contouring lung GTV regions.

FUNDING

Funding to support this study and the payment of the Open Access publication charges for this article was provided by the Okawa Foundation for Information and Telecommunications, the Ministry of Education, Culture, Sports, Science, and Technology (MEXT), and by Grants-in-Aid for Scientific Research on Innovative Areas (24103707) and Scientific Research (c) (22611011), 2012.