REZART BELSHI, PH.D., DOMINIQUE PONTVERT, M.D., JEAN-CLAUDE ROSENWALD,
PH.D. AND GENEVIEVE GABORIAUD, M.Sc.
Department of Radiotherapy, Institute Curie. Paris, France

Purpose: A method is provided for the automatic calculation of the Clinical Target Volume (CTV) by automatic three dimensional (3D) expansion of the Gross Tumor Volume (GTV), keeping a constant margin M in all directions and taking into account anatomic obstacles.

Methods and Materials: Our model uses a description of the GTV from contours (polygons) defined in a series of parallel slices obtained from Computed Tomography (CT) or Magnetic Resonance Imaging (MRI). Each slice is considered sequentially, including those slices located apart from the GTV at a distance smaller than M. In the current slice a two-dimensional (2D) expansion is performed by transforming each vertex of the polygon into a circle with a radius equal to M, and each segment into a rectangle with a height equal to 2M. A cartesian millimetric grid is then "projected" onto the slice and a specific value is assigned at each point depending if the point is internal to the 2D expansion. The influence in the current slice of any slice located at a distance D z smaller than M is taken into account by applying a 2D expansion using a margin

Additional contours representative of various "barriers" stopping the expansion process can also be defined.

Results: The method bas been applied to cylindrical and spherical structures and bas proven to be successful, provided that the slice thickness is small enough. For usual slice thicknesse3 and margins, it gives a slight over-estimation of the additional volume (around 5%) due to the choice that the calculated target volume would not be less than the expected volume. It bas been shown that for a spherical volume, a 2D expansion performed slice by slice leads to a volume Up to 80% smaller than that obtained by 3D expansion.

Conclusions: This tool, which mimics the tumor cell spreading process, bas been integrated in our treatment planning software and used clinically for conformal radiotherapy of brain and prostatic tumors. It bas been found to be extremely useful, not only saving time but also allowing a precise determination of the CTV which would be impossible to do manually. © 1997 Elsevier Science Inc.

Determination of the target volume is an essential step in planning conformal radiotherapy. The precision of its delineation directly influences the quality of treatment (2, 5, 6.8, 9.13, 14). According to International Commission on Radiation Units and Measurements (ICRU) (7), one can distinguish the Clinical Target Volume (CTV), which is obtained by adding to the Gross Tumor Volume (GTV) a margin taking into account the spread of the infraclinic disease. and the Planning Target Volume (PTV), which encompasses the CTV and includes any uncertainty due to patient motion or beam setup. In this article we will only address the problem of CTV determination.

Delineation of the GTV is performed by the radiotherapist and/or radiologist. usually on axial CT or NIRI scans. To provide consistency between the conformal treatments of a series of patients, it is necessary to define the CTV using a constant margin around the GTV. In principle, this margin must be the same in all directions, including the longitudinal one. This proves to be difficult or even impossible to do manually. At the most. one can use a constant lateral margin around the each GTV outline drawn in each cross-section. However, such an approach is time consuming and does not account for the 3D characteristics of the GTV. In particular, it does not provide any CTV description in slices longitudinally adjacent to the GTV where the risk of microscopic spread is the same as laterally.

We have, therefore. developed an automatic 3D calculation model that expands any volume described as a series of polygons in parallel planes into another volume, described in the same way, and obtained by adding a constant margin in all directions around the initial volume. In addition, this model is intended to mimic the cell spreading process and takes into account some predefined 'obstacles" that confine the expansion within given limits. A detailed description of our automatic volume expansion algorithm is given in the following sections. The practical use of this model in conformal radiotherapy of brain tumors is also discussed.

The external surface. anatomical data, and target volume were determined from CT and/or MRI images. A complete volume definition required the acquisition of thin adjacent axial slices. A thickness of 2 to 5 mm was generally sufficient. However, to correctly delineate small cerebral lesions, the zone of interest had to be explored by 1 mm slices (4, 10, 11, 14). It was recommended to make sure that there was no gap between the slices. If it is not the case, the axial positions and the thickness were preprocessed in order to fill in the gaps, but the results can become inaccurate.

The structures were delineated on the image display slice by slice, either manually or by using density threshold techniques (3.4). They were represented as polygons consisting of a maximum of 500 points (cumulative limit for all contours associated with each slice). For the expansion process only one polygon (the first one) was accepted for a given structure. The polygons defined adjacent prisms with a thickness equal to that of the CT slice. Each structure was composed of a series of prisms (12). The various structures were distinguished by allocating a code, a name, a color, and a density.

The objective of our model was to create a new structure from an existing structure. with a predefined margin M and ensuring regular expansion in three dimensions.

Application of the margin around the structure

This step began with de termination of the smallest parallelepiped P having its sides perpendicular to the principal directions and containing the expanded volume. It was easily obtained by adding the margin M to each of the six sides of the parallelepiped. which includes the original structure. The expansion was then performed slice by slice within this parallelepiped P

Our model has in common with models previously described for automatic generation of beam apertures (1) the use of a cartesian millimetric grid (pixel array) held in memory. For a given current slice S_C, the polygon which represented the contours of the initial structure. was 'projected" onto the grid, point by point and segment by segment. By applying the margin M, each point was transformed into a circle with a radius M and each segment became a rectangle with a height equal to 2M. The area of the circles and rectangles. covering the grid, formed a cloud of points corresponding to the 2D expansion of the polygon delimiting the structure in S_C (Fig. 1). A matrix POINT_INTERNAL (i, j) was associated with the grid and consisted of values either equal to 1 or 0, depending on whether the pixel (i, j) belonged to the expanded polygon.

However, the expansion in S_C was also influenced by the overlying and underlying slices passing through the structure and situated at a distance smaller than M. Let us assume provisionally that the slices are infinitely thin and consider one of them S located at a distance D z smaller than M from S_C (Fig. 2). In slice S, the expansion process can be represented by spheres and cylinders of radius M replacing, respectively, the 2D circles and rectangles in S_C referred to above. Its influencing on slice S_C can be thought of as the intersection of these spheres and cylinders with the current slice.

This results in the polygon in S influencing the expansion in S_C according to a margin m obtained from the relationship:

m being equal to M for the current slice (case treated above) and equal to 0 for D z ³ M.

To "record" this influence of S on S_c, all points included in circles or rectangles built with a margin m were 'back-projected" onto the grid of the current slice and contributed to the cloud of all points internal to the expanded structure as represented by the matrix POINT_INTERNAL (i, j). At the beginning of the whole expansion process, this matrix was initialized with a 0 value (= "false"). For each slice S a value of 1 (= "true") was then assigned to the points belonging to the cloud, this process being repeated for all slices S within the parallelepiped P.

The last step consisted of plotting the final contour which envelopes the cloud of points POINT_INTERNAL (i, j) = 1, representative of the expanded volume in the current slice. The matricial data (cloud of points) had to be transformed into vectorial data (contour). To do this, a first borderline point (i, j) of the cloud was initially selected by systematic search from an initial external point. The adjacent borderline point was then found by searching in the four main directions around this element (i, j), starting from the point opposite to the incoming direction and always tuning Clockwise. The search stopped at the adjacent point (k. 1) where POINT_INTERNAL (k,l) = l. This process. repeated until closure, allowed us to obtain a step contour of the final volume in the current slice. A reduction algorithm was then applied to minimize the number of points of the contour. The first phase of this algorithm proceeded with triplet of points (1, 2, 3), taking out point 2 if the angle (l® 2 l® 3) was smaller than a fraction of degree. A second pass lead to further reduction, when the algebraic distance between point 2 and segment 1-3 was smaller than a preset limit (1 mm). The sign of this limit was such that the resulting contour enclosed an area, always larger than the cloud area (Fig. 3). This "maximized" contour was attributed to the expanded volume and characterized by its code, name, and color. In addition, for dose calculation purposes. a specific density, obtained from the average density of the enclosed original pixels, was eventually associated.

The same algorithm was applied to each slice S_C within the parallelepiped P.

In the previous section, we assumed that each slice had a 0 thickness. However, this was not the case, and the exact meaning of D z has to be clarified. If we consider two slices of longitudinal position, z_i and z_j and thickness, t_i and t_j, we can define D z either as the distance |Z_i - Z_j| between the central planes of the two slices. or as the distance { |Z_i - Z_j| - 0.5 (t_i + t_j) } between their planar sides. If we further assume that a structure seen on a slice is attributed to its whole thickness, the second solution is better because. contrary to the first one, it pro-vides at the top and at the bottom of the initial structure a margin really equal to M (Fig. 4). Actually, this solution leads to the minimum "staircase" volume that completely encompasses the 'true" smoothed expanded volume. As illustrated in Fig. 4, it can be shown that the resulting overestimation is negligible as soon as the slice thickness t is small enough and that it is much smaller than the volume underestimation corresponding to the other solution. using the same slice thickness t. The optimal choice of t will be discussed later.

In practice, it is well known that the expansion is not necessarily the same in all directions. In particular, it respects anatomical limits such as the external surface of the body or, for the brain, the internal surface of the skull. In addition, the expansion spares certain types of volumes such as bony structures or, for the brain, the falx cerebri. It can "wrap" around these structures but cannot cross them directly. We called these elements "barrier." Hence, two types of barriers could be distinguished: 1) Envelope-defined as the external surface of a volume that contains the original structure and is such that no expansion is allowed outside; and 2) obstacles defined as regions that cannot be directly crossed by the expansion, but that can be avoided. They may be partially or completely included within the envelope.

These different regions are shown on Fig. 5. Elements outside of the envelope as well as the interior of the obstacles constituted barriers, which means zones not authorized for expansion. Recognition of the barriers, either for the envelope or the obstacle. was performed slice by slice, by projection of their contours onto the cartesian grid used for the expansion. all points of the grid inside the envelope corresponded to elements of a POINT_AUTHORIZED (i, j) matrix. with a value equal to 1 (= "true"). while the other points were equal to 0 (= "false"). Elements of the POINT_AUTHORIZED (i, j) matrix located inside the envelope and belonging to obstacles were set to a value of 0.

Expansion of the initial structure could therefore not extend beyond the envelope, and it was allowed to skirt around obstacles by means of iterative expansion.

Direct application of the margin m in slice S_C, in presence of obstacles, did not fully mimic the cell spreading process. The reason is that all points belonging to the m expansion [POINT_INTERNAL (i. j) = 1] and outside of the obstacle [POINT_AUTHORIZED (i, j) = 1] were considered as valid, although the "distance" traveled by those "cells" which had to skirt around obstacles could be larger than m (Fig. 6a).

The resulting overestimation of the expanded volume could be avoided by performing iterative expansion in which the number of repetitions n depends on the elementary step p where:

m
n = ----     
p

We, therefore, applied a p margin to the structure. p being 1 mm. The cloud obtained was represented by elements of the POINT_INTERNAL (i, j) matrix equal to 1. These elements were only valid when POINT_AUTHORIZED (i, j) was also equal to 1: otherwise they were attributed a value of 0. The cloud obtained following this test was successively expanded with the same margin of p until the total margin m was obtained after n iterations. This procedure gave more realistic results as shown in Fig. 6b.

The presence of obstacles could result in the formation of separate clouds (influence of underlying and overlying slices), and the volume of expansion could consequently be presented by several separate contours in a slice. In order to account for these situations, the calculation of contours ("vectorization") was not stopped until the whole grid POINT_INTERNAL (i, j) had been explored.

The overall calculation time for typical expansion from GTV to CTV, using an Alpha station DEC 3000/300X, was 1 s for each slice S_C within parallelepiped P.

This methodology has been applied to spherical structures of various radius (r = 0.5, 1, 2, 3, 4, and 5 cm) considering margins between 0 and 3 cm. These spherical structures were defined from circles acquired manually (using a digitizer) in parallel cross-sections of thickness varying from 1 mm in the central region to less than 0.2 mm in the "tangential" portion. The accuracy of the digitalization process was +0.2 min. In addition. to allow for expansion, supplementary cross-sections had to be defined within a margin M on each side of the original structure. They had a thickness of 2 mm.

The "theoretical" volumes of the original structures were calculated from

They were compared to the "computed" volumes. For these volumes, the thickness of the acquired slices was taken into account: the geometrical volume of the sphere was calculated by the sum of the volumes of the prisms of each slice. The volume of the prism was the product of the area (inside the contour) of the polygon (circle) and the thickness of the slice t. Analytical calculation of the area of the polygon as the sum of elementary triangular areas was performed using an arbitrary system of coordinates.

The agreement between theoretical and computed volumes was generally found to be better than 0.1% which is Consistent with the accuracy of the digitalization.

After expansion with margins M of 0.5, 1. 1.5, 2. and 3 cm, a similar comparison was performed between theoretical volume calculated from R = r + M and computed volumes, using the method described above. The volume obtained by 3D expansion was found to be 2 to 5% greater than the spherical theoretical volume,

The precision of the result of the volume expansion depended on several factors. The most important were the following:

If we consider in the current slice S_C the vectorial contour C obtained by joining those points that belonged to the expanded structure and were located at its edge. we obtained necessarily a contour covering an area smaller than the actual contour of the expanded volume. The maximum distance between the actual (unknown) contour and C was equal to the step s. The resulting underestimation of the volume was somewhat overcome by choosing a smoothing algorithm, which tended to use always the more external points.

On the other hand. as discussed previously, the value of D z used to calculate the margin m in S_C. corresponded to the distance between the edge of each slice and not the distance between the centers of the slices. This led to an overestimation of the volume in the (usual) situation where the margin M was larger than the slice thickness t.

The thickness of slices defining the spherical structure and those that were invaded by the expansion, therefore. influenced the volume of the expanded structure. The importance of slice thickness in the result of the expansion could be quantified by expanding a point structure situated in the middle of the serial slices. For example, we performed the expansion of a point with a 3 cm margin (the result should be a sphere) in a series of adjacent slices 1, 2, 3, and 5 mm thick. The volume of the sphere obtained in each series of slices was compared with the theoretical volume

The difference was -0.6, +1.2, +6.1, and +9.9%, respectively. Similar quantitative results were obtained on the additional volume created from cylindrical volumes of various radii using different margins M (Fig. 4). In general, it was found that the percentage error on the additional volume was mostly related to the ratio t/M. It increased linearly with t/M, and it was always smaller than 10% as long as the slice thickness t was smaller than 0.2M. It was around 5% for t = 0.1M. To avoid significant overestimation of the volume of expansion, it is therefore essential for margins M around 10 mm to obtain millimetric slices. These thin slices must be maintained not only through the initial structure but also on both sides longitudinally, within a distance M from top and bottom. respectively. If the acquisition is performed using thicker slices (i.e., 5 mm). some improvement in accuracy can be obtained by transforming each of them into five (or two) slices 1 (or 2.5) mm thick.

In the absence of such 3D algorithm, the margin M is generally applied slice by slice, ignoring the longitudinal expansion. Such a simplification can lead to significant errors, not only because of the lack of expansion in sections adjacent to the initial structure but also because the 3D nature of expansion in underlying and overlying slices is not taken into account. This is illustrated in Fig. 7.

Figure 8 illustrates the difference of expanded volumes obtained by 2D and 3D expansion of the spherical volumes presented above. The difference between 2D and 3D expansion depends on the size of the structure and the value of the margin. Note that reduction of the size of the initial structure accentuates the percentage difference between the volumes obtained by 2D and 3D expansion. For example, for a sphere with a radius of I cm. the application of a 2D margin of I cm resulted in a volume of 21 cm³, while 3D expansion gave a spherical volume of 36 cm³ (underestimation of 42%). The same situation in the case of a sphere with a radius of 3 cm induced an underestimation of the expanded volume of approximately 20%.

If we restrict our comparison to the slices in which the initial structure was defined, the difference is smaller but still significant (around 15%) as illustrated in Fig. 9. In this region, for spherical structures with a radius greater than the margin, the underestimation (%) is accentuated as the margin is increased. For spheres with a small radius, the 2-3D percentage difference of this expanded volume decreases with the size of the margin.

Automatic 3D calculation of the CTV from GTV has been proven to be an invaluable aid to the final estimation of this volume by the physician. It bas been integrated into the treatment planning system ISIS developed and used at Institut Curie. It has been in clinical use since August 1994 for the preparation of treatment plans of all intracranial and prostatitis tumors for which patients enter a protocol of conformal radiotherapy, including, for some of them, a complementary proton treatment performed at the Centre de Protonthérapie d'Qrsay (CPO). After delineation of the GTV and of the various barriers and decision on the margin value, the automatic calculation of the CTV is performed within some seconds. The radiotherapist must then review the contours of the CTV slice by slice and/or using the 3D visualization aids, to confirm its validity or do some manual adjustments. Additional margins can eventually be added to define automatically the Planning Target Volume (PTV) from the CTV, using the same module.

This module can also be used to help in the volumetric interpretation of the 3D dose distribution: For example. can create the internal contours of an hollow organ such as the rectum from its external contours, assuming a Constant wall thickness and using a negative value for margin, It can also be used to create several overlapping "expanded volumes" with different margins. Dose-volume histograms are then calculated and analyzed not only in the different organs, but also in the surrounding "rings" in order to make sure that there is no risk of overdosage or undersosage in the critical region next to the CTV.

In conclusion, this module, which integrates the presence of "barriers" and uses an iterative 3D Contiguous expansion algorithm, is very close to the reality of the tumor cell spreading process. Provided some precautions on the choice of slice thickness, it has become an essential tool in our approach of conformal radiotherapy, allowing a sound calculation of dose-volume histograms accounting for realistic limits of the CTV, and providing additional possibilities for evaluation and optimization of treatment plans.

REFERENCES

1. Brewster, L.; Mageras, G. S.; Mohan, R. Automatic generation of beam apertures. Med. Phys. 20:1337-1342; 1993.
2. Dische, S.; Saunders, M. I.; Williams, C.; Hopkins, A.; Aird, E. Precision in reporting the dose given in a course of radiotherapy. Radiother. Oncol. 29:287-293; 1993.
3. Dowsett, R. J.; Galvin, J. M.; Cheng, E.; Smith, R.; Epperson, R.; Harris, R.; Henze, G.; Needham, M.; Payne, R.; Peterson, M.; Skinner, A. L.; Reynolds, A. Contouring structures for 3-dimensional treatment planning. Int. J. Radiat. Oncol. Biol. Phys. 22:1083-1088; 1992.
4. Goitein, M.; Abrams, M. Multi-dimensional treatment planning: I. Delineation of anatomy. Int. J. Radiat. Oncol. Biol. Phys. 9:777-787; 1983.
5. Goitein, M. Calculation of the uncertainty in the dose delivered during radiation therapy. Med. Phys. 12:608~12; 1985.
6. Herbert, D. E. An extreme value paradigm for the effect of size of target volume on and results in radiation oncology. Med. Phys. 10:589~04; 1983.
7. ICRU REPORT 50 Prescribing, Recording, and Reporting Photon Beam Therapy. Washington; September 1993.
8. Leunens, G.; Menten, J.; Weltens, C.; Verstraete, J.; Van der Schuercn, E. Quality assessment of medical decision making in radiation oncology: Variability in target volume delineation for brain tumors. Radiother. Oncol. 29:169-175; 1993.
9. Mayles, W. P. M.; Chow, M.; Dyer, J.; Fernandez, E. M.; Heisig, s.; Knight, R. T.; Moore, I.; Nahum, A. E.; Shentall, G. S.; Tait, D. M. The Royal Marsden Hospital pelvic radiotherapy trial: Technical aspects and quality assurance. Radiother. Oncol. 29:l8~l9l; 1993.
10. Mohan, R.; Barest, G.; Brewster, L. J.; Chui, C. S.; Kutcher. G. J.; Laughlin, J. S.; Fuks, Z. A comprehensive three-dimensional radiation treatment planning system. Int. J. Radiat. Oncol. Biol. Phys. l5:481-95; 1988.
11. Rosenwald, J. C.; Gaboriaud, G.; Mazal, A.; Pontvert, D. Patient data acquisition for conformal radiotherapy. In: proceeding of EUTECH 94. Genova, 23-24 Feb. 1994, p. 2-4
12. Siddon, R. L. Prism representation: A 3D ray-tracing algorithm for radiotherapy application. Phys. Med. Biol. 30:817-824; 1985.
13. Ten Haken, R. K.; Thornton, A. F.; Sandler, H. M.; LaVigne, M. L.; Quint. D. J.; Fraas, B. A.; Kessler, M. L.; McShan, D. L. A quantitative assessment of the addition of MRI to CT-based, 3-D treatment planning of brain tumors. Radiother. Oncol. 25:121-133; 1992.
14. Urie, M.; Goitein, M.; Doppke, K.; Kutcher, J. G.; LoSasso, T.; Mohan, R.; Munzenrider, J. E.; Sontag, M.; Wong, W. The role of uncertainty analysis in treatment planning. Int. J. Radiat. Oncol. Biol. Phys. 21:91-107; 1991.