Active contouring and 3D model deformable registration of radiotherapy planning and cone-beam computed tomography images

Introduction

In order to treat a patient’s disease via safe and effective modern radiotherapy techniques, radiation therapy (RT) procedures must include various stages: patient data acquisition, treatment planning, treatment simulation, patient setup, and target localization before and during treatment. High accuracy is required to align the patient with respect to the radiation beam during treatment. There are several types of image techniques for the visualization of soft tissue structures. The most extensively used of these techniques is image-guided radiation therapy (IGRT), which includes image guidance procedures for target localization before and during treatment. These procedures use imaging technology to identify and correct problems arising from inter- and intra-fractional variations in patient setup and anatomy, including planning target volumes (PTVs), the shapes of the treatment target, and the surrounding normal tissues, especially that of organs at risk (OARs).

When easily correctable problems are recognized, such as dislocation, the couch can be shifted to match the original region. There are, however, some problems that cannot be solved directly, such as when a tumor changes its shape, size, or relative location to OARs, especially when said tumor is outside the treatment field. In such cases, re-alignment of the contour and treatment volume is the solution. The general procedure consists in arranging a computed tomography (CT) simulation and then contouring and planning again. Recreating a comprehensive treatment plan is complicated as well as time consuming and generally requires at least 4 to 6 days. The new plan cannot be immediately used on board when errors are detected, and the patient has to either be still treated using the old plan at incorrect locations or the treatment has to be temporarily halted, which decrease its radiobiological effects. When the new plan is available, the tumor might change its size, shape or location again.

Ideally, we want change to contour and plan immediately when an error is found. This would require recontouring straight away and comparing the two isodose plans to enable professionals to make setup corrections or modify the treatment parameters to minimize the variations between the planned and the actual treatment. Such procedures fall into a category called image-guided adaptive radiation therapy (IGART). For such techniques, online rapid recontouring would be very convenient. However, there are some limitations: (I) the suitability of the software used for online contouring; (II) the image resolution of cone-beam computed tomography (CBCT) is not very high; (III) the time and personnel required for contouring.

We designed a novel system that combines: (I) active contouring, (II) 3D model registration to increase the signal-to-noise ratio of the CBCT images; (III) demons deformable image registration (DIR) to transfer the treatment contour from simulated CT images to the corresponding CBCT images. Based on current clinical practices, CBCT is the most commonly used IGRT technique during radiation therapy. Although CBCT shares many similarities with traditional (fan beam) CT, there are important differences, particularly for reconstruction. CBCT has the advantages of having a lower cost, a shorter scanning time (it can be completed in less than one rotation), and lower radiation doses (1). The wider collimation in CBCT leads to increased scatter radiation and degradation of image quality, as evidenced by artifacts and decreased contrast-to-noise ratio. The time required for image reconstruction is longer for CBCT (1 minute) compared with that for traditional CT (real time) because of the computationally demanding cone-beam reconstruction algorithms used.

There are several research groups trying to improve the CBCT technique. For the reconstruction algorithm, there analytic 3D reconstruction algorithms (2,3) and statistical iterative reconstruction algorithms (4,5) have been developed. As for the problem of scattered doses, new methods have also been tried (6). To the best of our knowledge, there are no techniques or systems developed to combine the remodel CBCT images with auto-contouring to match the traditional CT images. In actual practice, transformations between CT and CBCT images are associated with deformations. The first step consists in increasing the resolution and performing contrast enhancements. The second step consists in matching structures with the same isocenter and then perform an auto-contour.

In the present study, we established the 3D reconstruction algorithm and demonstrated the results as a basis for further research.

Methods

Loading the CT and CBCT images

Different medical institutions and linear accelerators for radiotherapy employ different secure data systems and unique recognition procedures. CBCT and traditional CT use different file formats (CBCT data uses the extension *.his, whereas CT data uses *.dcm). They are both based on the Digital Imaging and Communications in Medicine (DICOM) standard. We need to develop a new system for loading the database.

For loading the image database, we wrote a program in MATLAB language to read and convert CBCT images into cross-sectional (tomographic) images similar to those obtained via CT.

Increasing resolution and contrast enhancement

Previously reported methods for increasing contrast are histogram equalization (HE) and wavelet decomposition and reconstruction. They are used for enhancing medical images, but they decrease saturation at the same time. Adaptive histogram equalization (AHE) (7,8) was developed to solve this problem.

De-noising is another practical approach. Some de-noising methods, such as spatial filtering, Fourier-based spectral filtering, mean filtering (9), and Gaussian filtering, can successfully decrease noise but they make the margins of the image blurry. Median filtering preserves the margins but can overcut the image.

Image registration

Where is the true location of the patient in real time? Unlike simulated CT images, the patient’s location in real CT images might have multi-directional errors in rotation, yaw, and/or pitch. The relationship between locations in two images can usually be determined via rigid-body transformations and by matching the isocenters of the two images (10). Auto-registration has been used in some algorithms, such as the Bone (uses the chamfer matching algorithm), GreyValue (11), and Demons algorithms (12,13). We tried these three methods and compared their advantages and disadvantages.

Auto-contouring

Active Contour model (Snakes)

In order to extract the region of interest in the image, an image segmentation technique is required. Two types of methods are available: (I) identifying pixels and then carrying out a reconstruction; (II) contour the outline of the item of interest. The active contour model, also called the Snakes model (14), is usually used for this purpose. This algorithm can clearly determine the edges of soft tissues. Because they are based on different formulas, there are several subtypes of the Snake model. One of them is the Chan-Vese model, which is based on the Mumford-Shah functional and can be used to solve energy minimization problems.

The core equation of the Snake model is shown in Equation [1].

F(c1,c2,C)=∫Ω1=ω(u0(x,y)−c1)2dxdy+∫Ω2=Ω−ω(u0(x,y)−c2)2dxdy+v|C|[1]

Equation [1]: Core equation of the Snake model.

C indicates the contour. C1 is the average radiodensity (CT number, Hounsfield units) inside the contour and C2 is that outside the contour. Ω means area; Ω1 is the area inside the contour and Ω2 is that outside the contour. µ₀ is the origin of the image and ω is the target contour.

Equation [1] can be simplified as shown below in Equation [2]:

F1(C)+F2(C)=∫inside(C)|uo−c1|2dx+∫outside(C)|uo−c2|2dx[2]

Equation [2]: Simplified core equation of the Snake model.

F1 corresponds to the construction of contour C and F2 corresponds to the dilation of contour C. The equation is balanced when the contour matches the margins of the item of interest.

Actual practice

Using formal origin contour

Figure 1 shows an example situation. Figure 1A shows what happens when the construction force is higher than the dilation force (F1 > F2), resulting in contour construction. In Figure 1B, this situation is reversed as F1 < F2, resulting in contour dilation.

Figure 1 Example situation of the active contour model.

Figure 1C shows an unbalanced situation, with the lower left region exhibiting contour dilation and the upper right region showing contour construction. Figure 1D indicates that the equation is balanced, resulting in no changes to the contour.

It is not easy to change the location and size of a round origin contour. On the other hand, rectangular origin contours do not usually match the target shape and are time consuming.

Using reformed origin contour

The first slice can be used to reform the contour from the formal origin contour. Then, the results (reformed origin contour) can be transferred to the next slice directly. In this way, the shape of the template contour can be more similar to the target than round or rectangular shapes and achieve faster results.

In actual practice, using reformed origin contour then remodeling as the sequence of steps presented in Figure 2.

Figure 2 The sequence of steps for using reformed origin contour.

Remodeling the reformed origin contour

Between slices, the image might undergo construction or dilation. The algorithm addresses certain problems during image construction in the next slide if the previous slice image is used as a template. These problems are, for example, when the margins cannot be found or when the contour is not correctly constructed. An addition algorithm was employed to detect changes in area between slices. If the area of the next slice is smaller, the contour is directly constructed on the template image. This addition results in a smoother and more ideal remodeling procedure.

Extraction of the main body

The resulting image from the previous step is recognized, remodeled, and filled up. The filled-up image is compared and multiplied with the original image, allowing us to extract the main body (Figure 3).

Figure 3 The result of extraction of the main body.

3D model registration

The patient’s location cannot always remain the same in all directions. This also happens when matching CBCT images and simulated CT images. Therefore, we need 3D model registration to registration in giant view.

For dislocation between the images as in Figure 4, we used the iterative closest point (ICP) method (15) to reconstruct 3D models of the images to make a general comparison. The ICP method converts the images into a cloud of points and calculate the distances between the points to then construct 3D models, as shown in Figure 5. The 3D models are then compared in terms of offset angle and size to carry out any necessary corrections.

Figure 4 Dislocations between CT and CBCT images. CBCT, cone-beam computed tomography.

Figure 5 3D models constructed from a cloud of points identified in the images.

Matching the most similar images

We used the matchTemplate method and the square difference matching method (method = CV_TM_SQDIFF) (16) to match the most similar images. We chose the least difference criteria to match the images. The procedure is presented in Figure 6.

Figure 6 Procedure for matching the most similar images.

DIR

We need an improved system and a method for precisely transforming CT images into CBCT images so as to enable the real-time treatment of the information comprising the treatment plan and CBCT image data.

How can CT image data be transferred to CBCT image data? We used the DIR method, which is based on the Demons algorithm. The latter was developed by Thirion in 1998 and is used to calculate the gradience and fluence of every spot in images for comparison purposes. It then indicates shifts in the shapes and locations of the image.

The algorithm consists in the three following steps and Figure 7:

Figure 7 The deformable image registration method based on the Demons algorithm.

• Calculate momentum of every spot in the images according to the following equation:

  $${\text{dr}}^{(k+1)}=\frac{(I_m^{(k)}-I_S)\nabla I_S}{{(I_m^{(k)}-I_S)}^2+{\left|\nabla I_S\right|}^2}+\frac{(I_m^{(k)}-I_S)\nabla I_m^{(k)}}{{(I_m^{(k)}-I_S)}^2+{\left|\nabla I_m^{(k)}\right|}^2}$$[3]

Where dr = (dx, dy, dz) indicates the location of each spot;

• Use a Gaussian equation to smooth the calculated dr;
• The smoothed dr is added to the vector field v(x) and then Im is renewed.

These three steps are repeated until the equation is balanced.

Results

The proposed novel system integrates data loading, image enhancement, reforming, auto-contouring, and image matching to make the auto-registration between CT and CBCT images possible.

After increasing resolution, contrast enhancement and image registration, the auto-contouring was developed by the Snake model.

After solving several practical problems, the system was put into actual practice, as shown in Figure 8.

Figure 8 The result of auto-registration between CT and CBCT images. CBCT, cone-beam computed tomography.

Reformed CT image was made using planning CT image and is the result of comparing and reforming. The shape was matched to that of CBCT image to achieve the clearer resolution and higher quality offered by CT images.

Discussion

Radiotherapy had widely used on-board CBCT scanning as image guidance procedure for target localization before and during treatment. This image technique has the nature advantages of having a lower cost, a shorter scanning time and lower radiation doses (17). It also has some limitation for the use in the practice, for example scatter radiation, degradation of image quality, different image database from traditional CT and lacking technique to match and modify these two kinds of images.

Some researchers used learning-based method to effectively capture the relationship between the planning CT and CBCT. It reduce scatter artifacts improving CBCT image quality to a level close to planning CT but not directly fit with the planning images (18). There was a study using MV CBCT images generated by proprietary 3D reconstruction software based on the FDK algorithm for megavoltage treatment planning (19). It is feasible after phantom evaluation but the planning only based by low quality CBCT image. Most of the papers declared as adaptive radiotherapy and used offline CBCT data for check the dose quality control with the CTV-to-PTV margin and fraction dose recalculations. Based on my knowledge, there was no acceptable technique directly transfers online CBCT to be modified, registration and auto-contouring (20).

In order to change contour and plan immediately when an error is found, we require recontouring straight away to minimize the variations between the planned and the actual treatment.

We develop the new system for loading the CBCT and traditional CT database. Increasing resolution and contrast enhancement of CBCT image. Auto-contouring use the Active Contour model (Snakes) then remodeling the reformed origin contour for the actual practice. We used the ICP method (15) to reconstruct 3D models and the square difference matching method to matching the most similar images. Adaptive RT aiming to modify RT target volumes during RT course is based on registration of CBCT and planning CT images. This new technique can automatically online rapid re-contouring and match the CBCT and traditional CT images, saving the time and cost.

The dose calculation is based on an accuracy Hounsfield Unit (HU) and electron density. We modified the images from the planning CT (traditional CT). The HU data remained to depend on the original CT images even we changed their shapes. We had checked HU values after modifying the images. The HU values were similar to the planning CT except when the volume changed significantly for the reasons such as marked changes in body weight. In that case, we arranged re-simulation and used the HU values of the new planning CT for calculation.

The next step of the IGART is the re-planning. The re-planning procedures had been developed by some company, for example the new Varian Linear Accelerator “Halcyon”. The combined of them might make the IGART in the real practice.

This isolated system can also be linking with the other radiation treatment procedure. It might make adaptive radiotherapy with proton or heavy ion technique possible.

Conclusions

This novel established technique consisting of active contouring with 3D model and deformable image registration. This technique would enable the on-line radiation treatment planning for adaptive radiotherapy.

Acknowledgments

We thank the staffs in graduate institute of the National Taiwan University of Science and Technology for the management and data collection for this study.

Funding: This research is supported by research program of National Taiwan University of Science and Technology and Mackay Memorial Hospital under grant number “MMH-TW-10601 and 10302”.

Footnote

Conflicts of Interest: YJC serves as an Editor-in-Chief of Therapeutic Radiology and Oncology. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Ning R, Chen B, Yu R, et al. Flat panel detector-based cone-beam volume CT angiography imaging: system evaluation. IEEE Trans Med Imaging 2000;19:949-63. [Crossref] [PubMed]
2. Johnson RH, Hu H, Haworth ST, et al. Feldkamp and circle-and-line cone-beam reconstruction for 3D micro-CT of vascular networks. Phys Med Biol 1998;43:929-40. [Crossref] [PubMed]
3. Turbell H. Three-dimensional image reconstruction in circular and helical computed tomography. Linköping: Linköpings Universitet, 1999.
4. Grass M, Köhler T, Proksa R. 3D cone-beam CT reconstruction for circular trajectories. Phys Med Biol 2000;45:329-47. [Crossref] [PubMed]
5. Fessler JA, Hero AO. Space-alternating generalized EM algorithms for penalized maximum-likelihood image reconstruction. Available online: https://web.eecs.umich.edu/~fessler/papers/lists/files/tr/94,286,sage.pdf
6. Huang K, Zhang D, Li M. Scatter detection and correction method of FPD-based DR/CT imaging systems. Available online: https://ieeexplore.ieee.org/document/5274066
7. Laine AF, Schuler S, Fan J, et al. Mammographic feature enhancement by multiscale analysis. IEEE Trans Med Imaging 1994;13:725-40. [Crossref] [PubMed]
8. Stark JA. Adaptive image contrast enhancement using generalizations of histogram equalization. IEEE Trans Image Process 2000;9:889-96. [Crossref] [PubMed]
9. Zhong J, Ning R, Conover D. Image denoising based on multiscale singularity detection for cone beam CT breast imaging. IEEE Trans Med Imaging 2004;23:696-703. [Crossref] [PubMed]
10. Hill DL, Batchelor PG, Holden M, et al. Medical image registration. Phys Med Biol 2001;46:R1-45. [Crossref] [PubMed]
11. Smitsmans MH, De Bois J, Sonke JJ, et al. Automatic prostate localization on cone-beam CT scans for high precision image-guided radiotherapy. Int J Radiat Oncol Biol Phys 2005;63:975-84. [Crossref] [PubMed]
12. Lu W, Chen M, Chen Q, et al. Adaptive fractionation therapy: I. Basic concept and strategy. Phys Med Biol 2008;53:5495-511. [Crossref] [PubMed]
13. Godley A, Ahunbay E, Peng C, et al. Automated registration of large deformations for adaptive radiation therapy of prostate cancer. Med Phys 2009;36:1433-41. [Crossref] [PubMed]
14. Kroon DJ. Multimodality non-rigid demon algorithm image registration. Robust Non-rigid Point Matching 2008;14:120-6.
15. Feldmar J, Declerck J, Malandain G, et al. Extension of the ICP algorithm to nonrigid intensity-based registration of 3D volumes. Comput Vis Image Underst 1997;66:193-206. [Crossref]
16. Bradski G, Kaehler A. Learning OpenCV: Computer vision with the OpenCV library. Sebastopol, CA: O'Reilly Media, Inc., 2008.
17. White LR. Mol Imaging Biol 2012;14:397-8; author reply 399-401.
18. Yang X, Lei Y, Higgins KA, et al. A Leaning-Based Method to Improve Cone Beam CT Image Quality for Adaptive Radiation Therapy. Int J Radiat Oncol Biol Phys 2017;99:S224. [Crossref]
19. Thomas TH, Devakumar D, Purnima S, et al. The adaptation of megavoltage cone beam CT for use in standard radiotherapy treatment planning. Phys Med Biol 2009;54:2067-77. [Crossref] [PubMed]
20. Posiewnik M, Piotrowski T. A review of cone-beam CT applications for adaptive radiotherapy of prostate cancer. Phys Med 2019;59:13-21. [Crossref] [PubMed]