**Segmentation Techniques of Anatomical Structures with Application in Radiotherapy Treatment Planning**

S. Zimeras

*University of the Aegean, Department of Statistics and Actuarial, Financial Mathematics, Karlovassi Samos, Greece* 

#### **1. Introduction**

40 Modern Practices in Radiation Therapy

Ryu, S., R. Jin, J. Y. Jin, Q. Chen, J. Rock, J. Anderson, & B. Movsas. (2008). Pain control by

Sahgal, A., L. Ma, I. Gibbs, P. C. Gerszten, S. Ryu, S. Soltys, V. Weinberg, et al. (2010). Spinal

Scarantino, C. W., D. M. Ruslander, C. J. Rini, G. G. Mann, H. T. Nagle, & R. D. Black. (2004).

Schulman, K. L., & J. Kohles. (2007). Economic burden of metastatic bone disease in the U.S.

Sheng, K., R. Jones, W. Yang, B. Schneider, Q. Chen, G. Sobering, G. H. Olivera, & P. W.

Siewerdsen, J. H., & D. A. Jaffray. (2001). Cone-beam computed tomography with a flat-

Timmerman, R., R. Paulus, J. Galvin, J. Michalski, W. Straube, J. Bradley, A. Fakiris, et al.

Ulin, K., M. M. Urie, & J. M. Cherlow. (2010). Results of a multi-institutional benchmark test

Wagner, T. H., S. L. Meeks, F. J. Bova, W. A. Friedman, T. R. Willoughby, P. A. Kupelian, &

Wiersma, R. D., Z. Wen, M. Sadinski, K. Farrey, & K. M. Yenice. (2010). Development of a

Wu, J. S., R. Wong, M. Johnston, A. Bezjak, T. Whelan, & Cancer Care Ontario Practice

*Symptom Management,* Vol.35, No.3, pp. 292-8, ISSN 0885-3924

*Medical Physics,* Vol.31, No.9, pp. 2658-71, ISSN 0094-2405

*Biology, Physics,* Vol.77, No.5, pp. 1584-9, ISSN 0360-3016

*Medical Dosimetry,* Vol.32, No.2, pp. 111-20, ISSN 0958-3947

*Biology,* Vol.55, No.2, pp. 389-401, ISSN 1361-6560

*Cancer,* Vol.109, No.11, pp. 2334-42, ISSN 1097-0142

pp. 220-31, ISSN 0094-2405

ISSN 0098-7484

ISSN 0360-3016

*Oncology, Biology, Physics,* Vol.77, No.2, pp. 548-53, ISSN 0360-3016

image-guided radiosurgery for solitary spinal metastasis. *Journal of Pain and* 

cord tolerance for stereotactic body radiotherapy. *International Journal of Radiation* 

An implantable radiation dosimeter for use in external beam radiation therapy.

Read. (2011). 3D dose verification using tomotherapy CT detector array. *International Journal of Radiation Oncology, Biology, Physics,* (In press), ISSN 0360-3016

panel imager: Magnitude and effects of x-ray scatter. *Medical Physics,* Vol.28, No.2,

(2010). Stereotactic body radiation therapy for inoperable early stage lung cancer. *JAMA : The Journal of the American Medical Association,* Vol.303, No.11, pp. 1070-6,

for cranial CT/MR image registration. *International Journal of Radiation Oncology,* 

W. Tome. (2007). Optical tracking technology in stereotactic radiation therapy.

frameless stereotactic radiosurgery system based on real-time 6D position monitoring and adaptive head motion compensation. *Physics in Medicine and* 

Guidelines Initiative Supportive Care Group. (2003). Meta-analysis of dosefractionation radiotherapy trials for the palliation of painful bone metastases. *International Journal of Radiation Oncology, Biology, Physics* Vol.55, No.3, pp. 594-605, The definition of structures and the extraction of organ's shape are essential parts of medical imaging applications. These might be applications like diagnostic imaging, image guided surgery or radiation therapy. The aim of the volume definition process is to delineate a specific shape of an organ on a digital image as accurate as possible especially for 3D rendering, radiation therapy, and surgery planning. This can be done, either by manual user interaction or applying imaging processing techniques for the automatic detection of specific structures in the image using segmentation techniques. Segmentation is the process that separates an image into its important features (primitives) so that each of them can be addressed separately. This converts the planar pixel of the image into a distinguishable number of individual organs or tumour that can be clearly identified and manipulated. The segmentation process might involve complicate structures and in this case usually only an expert can perform the task of the identification manually on a slice-by-slice base. Humans can perform this task using complex analysis of shape, intensity, position, texture, and proximity to surrounding structures. In this work we present a set of tools that are implemented on several computer based medical application. Central focus of this work, are techniques used to improve time and interaction needed for a user when defining one or more structures based on segmentation techniques. These techniques involve interpolation methods for the manual volume definition and methods for the semi-automatic organ shape extraction.

#### **2. Radiotherapy Treatment Planning (RTP)**

The goal of radiotherapy treatment planning is to justify an effective treatment that will deliver a precise irradiation dose to the target volume without causing damage to the surrounding normal tissues. Therefore patient positioning, target volume definition and irradiation field placement are very critical steps while planning the irradiation process. Clinically a radiotherapy treatment plan is verified by Virtual Simulators (VS). (Sherouse et. al., 1987) first proposed the concept, often defined as CT-Sim to distinguish it from Sim-CT

Segmentation Techniques of Anatomical

BEV.

**3. Volume rendering** 

Karangelis, 2004).

Structures with Applicationin Radiotherapy Treatment Planning 43

Fig. 2. The six-window layout of the system. On the left side the slices windows, the middle lower window is the OEV, the lower right the Room View and the upper right hosts the

Volume rendering is the technique according to which a scalar field of data with discrete values, volume data, is selectively sampled in order to generate a useful image in relation to the sampled values. A volume data set is typically a set of discrete samples of one or more physical properties in a limited object space, *V(x), xRn*, in which {*x*} is a set of sampling points; *n* is the dimension of the sampling space, usually *n***=***3*, i.e. 3D volume data; *V*  represents the sampling values, it can be a single-valued or multi-valued function. According to the distribution of sampling points (i.e. the structure of *x*), volume data sets are classified into structured and unstructured data sets. In medical imaging, volume data is usually a structured data set, typically organised as a stack of slices; *V* can be single valued (e.g. the CT Hounsfield value) or multi-valued (e.g. density, T1, T2 in MRI). The resulting data structure is an orthogonal 3D array of voxels, each representing a sampling value. A general pipeline of volume visualisation in medicine can include several steps (Sakas, 1993;

According to the distribution of sampling points, volume data sets are classified into structured and unstructured data sets. In medical imaging, volume data is usually a structured data set, typically organized as a stack of slices. The rendering methods which are proposed in this work are based on the work of (Sakas, 1993): (1) Transparent mode: Digitally Reconstructed Radiographies (DRR) (Sakas, 1993); (2) Surface reconstruction mode: iso-value, gradient (Sakas, 1995) and semi-transparent. DRR images (Cai, 1999) generated from CT digital volumes are often called using the term digitally reconstructed radiographs (DRRs) (Figure 3). The term DRR is used when we refer to those X-ray images that are

where a simulator is modified for CT use and by the late 1990s several designs and clinical assessments of CT virtual simulators have been reported (Sherouse et. al., 1987; Nagata et. al., 1990; Nishidai et. al., 1990; Rosenman et. al., 1991; Sherouse and Chaney, 1991; Perez et. al., 1994; Butker et. al., 1996; Chen et. al., 1996; Michalski et. al., 1996; Purdy, 1996; Ragan et. al., 1996; Rosenman, 1996) . Using VS, the clincal routine is modefied accordingly (Figure 1) (Conway and Robinson, 1997; Valicenti et. al., 1997; Cai et. al., 1997; Karangelis, 2004):


Fig. 1. Current clinical routine for external beam treatment delivery.

The important RTP imaging is the dynamic X-ray image on the Simulator monitor, which is generated from the viewpoint of the beam source, called Beam's Eye View (BEV) image, see the BEV view direction in Figure 2. In the system the corresponding window, called the BEV window, displays the interactive DRR images. Through the BEV window, the physicians can investigate the patient anatomy in the DRR image as they observe it on the Simulator monitor. Observer's Eye View (OEV) window visualizes the patient from the room's eye view, which is the same view as the physicians observes the patient in the Simulator room, see the OEV view direction in Figure 2. As we explained before, all of the images in both BEV and OEV windows are generated with the patient CT scan by the direct volume rendering (Figure 2).

Fig. 2. The six-window layout of the system. On the left side the slices windows, the middle lower window is the OEV, the lower right the Room View and the upper right hosts the BEV.

#### **3. Volume rendering**

42 Modern Practices in Radiation Therapy

where a simulator is modified for CT use and by the late 1990s several designs and clinical assessments of CT virtual simulators have been reported (Sherouse et. al., 1987; Nagata et. al., 1990; Nishidai et. al., 1990; Rosenman et. al., 1991; Sherouse and Chaney, 1991; Perez et. al., 1994; Butker et. al., 1996; Chen et. al., 1996; Michalski et. al., 1996; Purdy, 1996; Ragan et. al., 1996; Rosenman, 1996) . Using VS, the clincal routine is modefied accordingly (Figure 1) (Conway and Robinson, 1997; Valicenti et. al., 1997; Cai et. al., 1997; Karangelis, 2004):

2. Transfer CT data to VS. The physician defines the tumour volume and the organs at risk

3. The simulation plan and the CT data are transferred via DICOM (Digital image and Communication in Medicine) server to the TPS for dose calculation and final treatment

and she/he will place the necessary fields relative to the tumour volume.

5. Perform treatment on the treatment machine (Linear Accelerator or LINAC).

1. Collect patient's CT data including attached aluminium markers.

Fig. 1. Current clinical routine for external beam treatment delivery.

The important RTP imaging is the dynamic X-ray image on the Simulator monitor, which is generated from the viewpoint of the beam source, called Beam's Eye View (BEV) image, see the BEV view direction in Figure 2. In the system the corresponding window, called the BEV window, displays the interactive DRR images. Through the BEV window, the physicians can investigate the patient anatomy in the DRR image as they observe it on the Simulator monitor. Observer's Eye View (OEV) window visualizes the patient from the room's eye view, which is the same view as the physicians observes the patient in the Simulator room, see the OEV view direction in Figure 2. As we explained before, all of the images in both BEV and OEV windows are generated with the patient CT scan by the direct volume

4. Verify patient position on LINAC before irradiation.

plan optimization.

rendering (Figure 2).

Volume rendering is the technique according to which a scalar field of data with discrete values, volume data, is selectively sampled in order to generate a useful image in relation to the sampled values. A volume data set is typically a set of discrete samples of one or more physical properties in a limited object space, *V(x), xRn*, in which {*x*} is a set of sampling points; *n* is the dimension of the sampling space, usually *n***=***3*, i.e. 3D volume data; *V*  represents the sampling values, it can be a single-valued or multi-valued function. According to the distribution of sampling points (i.e. the structure of *x*), volume data sets are classified into structured and unstructured data sets. In medical imaging, volume data is usually a structured data set, typically organised as a stack of slices; *V* can be single valued (e.g. the CT Hounsfield value) or multi-valued (e.g. density, T1, T2 in MRI). The resulting data structure is an orthogonal 3D array of voxels, each representing a sampling value. A general pipeline of volume visualisation in medicine can include several steps (Sakas, 1993; Karangelis, 2004).

According to the distribution of sampling points, volume data sets are classified into structured and unstructured data sets. In medical imaging, volume data is usually a structured data set, typically organized as a stack of slices. The rendering methods which are proposed in this work are based on the work of (Sakas, 1993): (1) Transparent mode: Digitally Reconstructed Radiographies (DRR) (Sakas, 1993); (2) Surface reconstruction mode: iso-value, gradient (Sakas, 1995) and semi-transparent. DRR images (Cai, 1999) generated from CT digital volumes are often called using the term digitally reconstructed radiographs (DRRs) (Figure 3). The term DRR is used when we refer to those X-ray images that are

Segmentation Techniques of Anatomical

(0.8sec to 2.3sec) (Figure 4)

tissues and lung tissues (Karangelis, 2004)

eight (8) voxels (p000, p001, p011, p010, p100, p110, p101, and p111).

*HU*

 

to illustrate the complete range of the HU on a volume 12bit voxels are required.

Fig. 4. Volume rendering modes. On the top row from left to right: isovalue mode,

semitransparent mode and maximum intensity projection. On the lower row X-ray images reconstructed using different tissue ranges. From left to right: full tissue range, muscle

where pµ , pwater are the attenuation coefficient of the material for the specific voxel and the water respectively. When on the above equation one replaces the pµ =pwater, then will receive a HUwater = 0. In addition the air as material corresponds to –1000HU, since pair = 0. The Hounsfield units have no upper limit but usually for medical scanners a range between – 1024 to +3071 is provided. Apparently 4096 different HU values are provided and therefore

Probably the most common methods for reconstructing surfaces directly from voxels are the gradient surface and iso-surface model. This rendering model is widely used in almost all medical data sets to render the surface detected by the gradient operator. (Levoy, 1988) presented this concept, which is effectively used in most medical imaging applications for the last decade. Depending on the data set size and on the size of the resulting 3D image, which influence the number of rays, a 3D view image can be calculated almost in real time

The HU can be estimated using the following formula:

Structures with Applicationin Radiotherapy Treatment Planning 45

where *K*l is the attenuation coefficient and *l* is the wavelength. This interpolation method is known as *bi-linear* interpolation and is commonly applied on image interpolation e.g. when magnifying raster images. In volumetric data linear interpolation is known as *tri-linear*  interpolation and estimates the value of the result voxel considering the neighbourhood of

To reconstruct a DRR the digital CT data of the patient are used. Each voxel acquired using CT has an value that is called the Hounsfield unit which was named after the CT inventor.

> \* 1000 *water water*

(5)

Fig. 3. left image. Reconstruction of CT volume using DRR; right image. Reconstruction using direct volume surface using direct gradient volume estimation (Karangelis, 2004; Zimeras and Karangelis, 2008)

generated with an unrealistic way using direct volume rendering techniques or to those images that are generated from volume data using a better approximation of the physical model. In both cases we try to simulate the attenuation of the X-ray through a medium, in our case through the digital patient's body. This can be achieved by using different transfer functions simulating the classical way of X-ray image reconstruction.

Assuming that a ray with initial energy Io enters a volume, which has a thickness of Δs. and it is composed from different materials. Using discrete values the final energy of the ray is a product of the following equation (Max, 1995):

$$I(s) = I\_0 \exp\left\{-\left(p\_1s\_1 + p\_2s\_2 + \dots + p\_xs\_x\right)\right\} \tag{1}$$

where I(s) is the attenuated energy after leaving the volume, Io the original energy of the ray, pi the linear attenuation coefficient for the corresponding voxel material, si the corresponding distance from the ray entrance point to the ray exiting position of the voxel. In medical data the material type corresponds of course to the different tissue type. For detailed description of the model of the X-ray refer to (Karangelis, 2004). The contrast intensity of the final DRR image on the screen level can be calculated using the equation:

$$I\_p = I\_0 \, \text{\*} \left(1 - L\right) + I\_{buckground} \, \text{\*} \, L \tag{2}$$

where *L ps ps ps* exp 11 22 *x x* . In two dimensions the interpolation scheme will involve four data points' p00, p01, p10 and p11. In this case L becomes:

$$L = \exp\left[F(p\_0 p\_1)\right] \tag{3}$$

where F given by the form

$$F\left(p\_0 p\_1\right) = \bigcap\_{p\_0} K\_l\left(t\right)dt\tag{4}$$

Fig. 3. left image. Reconstruction of CT volume using DRR; right image. Reconstruction using direct volume surface using direct gradient volume estimation (Karangelis, 2004;

functions simulating the classical way of X-ray image reconstruction.

involve four data points' p00, p01, p10 and p11. In this case L becomes:

product of the following equation (Max, 1995):

generated with an unrealistic way using direct volume rendering techniques or to those images that are generated from volume data using a better approximation of the physical model. In both cases we try to simulate the attenuation of the X-ray through a medium, in our case through the digital patient's body. This can be achieved by using different transfer

Assuming that a ray with initial energy Io enters a volume, which has a thickness of Δs. and it is composed from different materials. Using discrete values the final energy of the ray is a

where I(s) is the attenuated energy after leaving the volume, Io the original energy of the ray, pi the linear attenuation coefficient for the corresponding voxel material, si the corresponding distance from the ray entrance point to the ray exiting position of the voxel. In medical data the material type corresponds of course to the different tissue type. For detailed description of the model of the X-ray refer to (Karangelis, 2004). The contrast intensity of the final DRR image on the screen level can be calculated using the equation:

where *L ps ps ps* exp 11 22 *x x* . In two dimensions the interpolation scheme will

 <sup>1</sup> 0

*p*

*l p*

0 1

*I s I ps ps ps* <sup>0</sup> exp 11 22 *x x* (1)

<sup>0</sup> \* 1 \* *<sup>p</sup> background I I LI L* (2)

*L F pp* exp 0 1 (3)

*<sup>F</sup> p p K t dt* (4)

Zimeras and Karangelis, 2008)

where F given by the form

where *K*l is the attenuation coefficient and *l* is the wavelength. This interpolation method is known as *bi-linear* interpolation and is commonly applied on image interpolation e.g. when magnifying raster images. In volumetric data linear interpolation is known as *tri-linear*  interpolation and estimates the value of the result voxel considering the neighbourhood of eight (8) voxels (p000, p001, p011, p010, p100, p110, p101, and p111).

To reconstruct a DRR the digital CT data of the patient are used. Each voxel acquired using CT has an value that is called the Hounsfield unit which was named after the CT inventor. The HU can be estimated using the following formula:

$$H\text{UI} = \frac{\left(\rho\_{\mu} - \rho\_{water}\right)\_{\*}}{\rho\_{water}} \* 1000\tag{5}$$

where pµ , pwater are the attenuation coefficient of the material for the specific voxel and the water respectively. When on the above equation one replaces the pµ =pwater, then will receive a HUwater = 0. In addition the air as material corresponds to –1000HU, since pair = 0. The Hounsfield units have no upper limit but usually for medical scanners a range between – 1024 to +3071 is provided. Apparently 4096 different HU values are provided and therefore to illustrate the complete range of the HU on a volume 12bit voxels are required.

Probably the most common methods for reconstructing surfaces directly from voxels are the gradient surface and iso-surface model. This rendering model is widely used in almost all medical data sets to render the surface detected by the gradient operator. (Levoy, 1988) presented this concept, which is effectively used in most medical imaging applications for the last decade. Depending on the data set size and on the size of the resulting 3D image, which influence the number of rays, a 3D view image can be calculated almost in real time (0.8sec to 2.3sec) (Figure 4)

Fig. 4. Volume rendering modes. On the top row from left to right: isovalue mode, semitransparent mode and maximum intensity projection. On the lower row X-ray images reconstructed using different tissue ranges. From left to right: full tissue range, muscle tissues and lung tissues (Karangelis, 2004)

Segmentation Techniques of Anatomical

daily clinical routine are:

the patient's body.

regions.

ellipse (Figure 7).

Structures with Applicationin Radiotherapy Treatment Planning 47

The first aim of the segmentation process in RTP is to define as accurately as possible the target that will be irradiated. This process is a manual process and usually is performed from an expert oncologist. With the terms of manual process, we mean that the physicians will exam the digital volume, slice-by-slice and using the mouse cursor will illustrate the shape of the tumour. The coordinates of the shape of the tumour are stored from the system for further processing. The target definition process is probably the most time consuming process involving an expert during the RT planning (Ketting et. al., 1997). The following describes a number of reasons proving that the automatic target volume definition is still very difficult to be performed from an expert system. The obvious reasons found in the

1. The irregular tumour properties, like shape, texture, volume, relation with surrounding

2. Location within the human body. A tumour theoretically could grow anywhere within

3. Tumour spreading and variation. Depending from the region the tumour is grown, disease cells might spread to the surrounding region. The disease cell distribution

4. Artifacts in the acquired digital data. When treating elderly patients it is often the case to have prosthesis, usually metallic (heart irritating) like hip prosthesis. These patients cannot be examined in a MRI modality and therefore CT imaging is the only alternative for their RTP. Nevertheless, it is well known that metallic components generate severe artifacts in CT imaging that reduce image quality and blur tumour borders. Thus, the

Classically image segmentation denotes the technique of extraction of images structures (regions or objects) so that the outlines of these structures will coincide as accurately as possible with the physical 2D object outlines. Image segmentation approaches may be performed in one of these ways: Manual segmentation methods include pixel selection, geometrical boundary selection and tracing. Given normal image resolution, selection of

*Linear Interpolation*: Linear interpolation between contours is the first approach used to provide an acceleration tool for the manual contouring. The mechanism of the linear interpolation is applied when between the key contours at least one slice exists. To perform the linear interpolation we create triangles between the contour points of the key contours as described in (Strassman, 2000). For this operation both contour's points must be rotated towards the same rotation direction. The interpolated contour points are created after

*Orthogonal Contour Interpolation*: The orthogonal contour interpolation serves to create a volume combining and interpolating orthogonal drawn contours. Principally the algorithm needs at least 2 orthogonal contours to work. The perpendicular plane to these two contours creates intersection points that are the key points to create the new interpolated contour. In this approach we use the cubic Spine interpolation. In case the lines are completely equal in size and their centres match or have very small distance then the result of the interpolation will be a circle. In any other case that the two vertical lines are unequal the result will be an

calculating the intersection of each triangle side with the intermediate slice (Figure 6)

accurate manual target volume definition is even more necessary nowadays.

cannot be predicted and detected on the digital patient's data.

individual pixels is clearly impractical and rarely used.

A significant part for the radiation field is the virtual light field projection on the patient skin [LuH99]. In physical simulators, a light source is located near the irradiation source. The orientation of the light intensity is diverted through the gantry head using a mirror aperture. The outcome of this process is the exact projection of the radiation field on the patient's skin. The two main axis of the field are indicated as line shadows. In case the radiation field is delineated using shielding blocks, then the light field area is also modified accordingly. This process described above should be performed in a similar manner in the virtual simulation process. In order to realise this principle we take advantage of the convexity of the tetrahedral objects and the Z-depth information derived during polygon scan conversion. Using the conventional Simulator the block shape was drawn manually on the patient's X-ray film, acquired from the BEV direction. Then this shape was digitized and its digital points were saved on the block-cutting machine. Both beam and block geometry is defined using combinations of 3D planes. Each beam is represented as a pyramid. The height of the pyramid is calculated according to the current machine specification; the base of the pyramid represents the irradiation field size projected to the image detector level and each side of the pyramid is assigned to a plane (Figure 5)

Fig. 5. 3D beam object reconstructed with patient's CT data. Green object: tumour

#### **4. Segmentation techniques**

Depending upon the application field, individual-processing steps may be neglected, combined, or reversed. Important to notice is that the final 3D rendering image can be obtained in two ways: either through the intermediate surface representation or through the volume representation (i.e. direct volume rendering segmentation is the process that separates an image into its important features (primitives) so that each of them can be addressed separately. This converts the planar pixel of the image into a distinguishable number of individual organs or tumour that can be clearly identified and manipulated. The segmentation process might involve complicate structures and in this case usually only an expert can perform the task of the identification manually on a slice-by-slice base. Humans can perform this task using complex analysis of shape, intensity, position, texture, and proximity to surrounding structures. All these features are differently qualified depending on the experience of the user. To generate a "complete", segmentation application numerous tools and algorithms must be combined (Kuszy et. al., 1995).

A significant part for the radiation field is the virtual light field projection on the patient skin [LuH99]. In physical simulators, a light source is located near the irradiation source. The orientation of the light intensity is diverted through the gantry head using a mirror aperture. The outcome of this process is the exact projection of the radiation field on the patient's skin. The two main axis of the field are indicated as line shadows. In case the radiation field is delineated using shielding blocks, then the light field area is also modified accordingly. This process described above should be performed in a similar manner in the virtual simulation process. In order to realise this principle we take advantage of the convexity of the tetrahedral objects and the Z-depth information derived during polygon scan conversion. Using the conventional Simulator the block shape was drawn manually on the patient's X-ray film, acquired from the BEV direction. Then this shape was digitized and its digital points were saved on the block-cutting machine. Both beam and block geometry is defined using combinations of 3D planes. Each beam is represented as a pyramid. The height of the pyramid is calculated according to the current machine specification; the base of the pyramid represents the irradiation field size projected to the image detector level and

each side of the pyramid is assigned to a plane (Figure 5)

tools and algorithms must be combined (Kuszy et. al., 1995).

**4. Segmentation techniques** 

Fig. 5. 3D beam object reconstructed with patient's CT data. Green object: tumour

Depending upon the application field, individual-processing steps may be neglected, combined, or reversed. Important to notice is that the final 3D rendering image can be obtained in two ways: either through the intermediate surface representation or through the volume representation (i.e. direct volume rendering segmentation is the process that separates an image into its important features (primitives) so that each of them can be addressed separately. This converts the planar pixel of the image into a distinguishable number of individual organs or tumour that can be clearly identified and manipulated. The segmentation process might involve complicate structures and in this case usually only an expert can perform the task of the identification manually on a slice-by-slice base. Humans can perform this task using complex analysis of shape, intensity, position, texture, and proximity to surrounding structures. All these features are differently qualified depending on the experience of the user. To generate a "complete", segmentation application numerous The first aim of the segmentation process in RTP is to define as accurately as possible the target that will be irradiated. This process is a manual process and usually is performed from an expert oncologist. With the terms of manual process, we mean that the physicians will exam the digital volume, slice-by-slice and using the mouse cursor will illustrate the shape of the tumour. The coordinates of the shape of the tumour are stored from the system for further processing. The target definition process is probably the most time consuming process involving an expert during the RT planning (Ketting et. al., 1997). The following describes a number of reasons proving that the automatic target volume definition is still very difficult to be performed from an expert system. The obvious reasons found in the daily clinical routine are:


Classically image segmentation denotes the technique of extraction of images structures (regions or objects) so that the outlines of these structures will coincide as accurately as possible with the physical 2D object outlines. Image segmentation approaches may be performed in one of these ways: Manual segmentation methods include pixel selection, geometrical boundary selection and tracing. Given normal image resolution, selection of individual pixels is clearly impractical and rarely used.

*Linear Interpolation*: Linear interpolation between contours is the first approach used to provide an acceleration tool for the manual contouring. The mechanism of the linear interpolation is applied when between the key contours at least one slice exists. To perform the linear interpolation we create triangles between the contour points of the key contours as described in (Strassman, 2000). For this operation both contour's points must be rotated towards the same rotation direction. The interpolated contour points are created after calculating the intersection of each triangle side with the intermediate slice (Figure 6)

*Orthogonal Contour Interpolation*: The orthogonal contour interpolation serves to create a volume combining and interpolating orthogonal drawn contours. Principally the algorithm needs at least 2 orthogonal contours to work. The perpendicular plane to these two contours creates intersection points that are the key points to create the new interpolated contour. In this approach we use the cubic Spine interpolation. In case the lines are completely equal in size and their centres match or have very small distance then the result of the interpolation will be a circle. In any other case that the two vertical lines are unequal the result will be an ellipse (Figure 7).

Segmentation Techniques of Anatomical

al., 2001; Gross et. al., 2002).

Structures with Applicationin Radiotherapy Treatment Planning 49

Fully automatic segmentation methods are usually impractical due to image complexity and the variety of image types and interpretation. In addition, low contrast between structures

*Active Contour Model*: Active Contour Models (ACMs) are adaptive contour representations, also known as snakes or deformable models (Terzopoulos and Fleischer, 1998). They are able to recover and represent physical contours of an image, and hence can be used as a model to determine object boundaries in static images as well as for tracking in image sequences. (Grosskopf et al., 1998; 2002) uses the Euler Time Integration to solve the optimisation problem. After initialisation by a user sketch, the contour is deformed to fit the actual object by simulating physical properties of an elastic material or fluid. This method is very reliable to overcome local minima and very fast due to its deterministic character.

So far only a very few techniques propose physicians one or more alternatives for volume definition (Ketting et. al., 1997; Belsh et. al., 1997). Recently (Pekar et. al., 2004) reported a method based on an adaptation of 3D deformable surface models to the boundaries of the anatomic structures of interest. The adaptation was based on a tradeoff between deformations of the model induced by its attraction to certain image features and the shape integrity of the model. Nevertheless, to make the concept clinically feasible, interactive tools were also introduced that allow quick correction in problematic areas in which the automated model adaptation may fail. Currently the tool used to perform the volume definition is the mouse cursor. The user defines a number of digital points on the image level, closing the first and the last point of the contour to generate this way a structure. The connectivity between the *key points* can be linear or higher order. The higher order connectivity can be achieved using interpolation models like Hermitte cubic, Spline curves, Bezier curves (Laurent, 1994; Spath, 1995; Cohen 2001), which are the most common and successfully used techniques for smoothing curves in the systems. Aim of the high order interpolation techniques is to reduce the amount of input points required to describe a smooth shape. Performing a simple comparison between linear and higher order interpolation, it can find out that the amount of points used to illustrate the shape of a structure using linear interpolation requires at least twice as many samples as by using the higher order interpolation algorithms. A common methodology used to combine high order interpolation and image edge properties is the use of active contour models. The active contour models or Snakes can be 2D image curves (Kass et. al., 1987; Blake and Isard, 1998) or 3D polygon meshes (Terzopoulos and Fleischer, 1998), which are adjusted from an initial approximation to the image or volume features by a movement caused by simulated forces. Image features provide the so-called external force. An internal tension of the curve resists against highly angled curvatures, which makes the Snakes movement robust against noise. After a starting position is given, the Snake adapts itself to shape by relaxation to the equilibrium of the external force and internal tension. Snakes have been proven efficient and fast for a number of applications in medicine involving different imaging modalities (McIne and Terzopoulos, 1996; Gross et. a., 1998; Behr et. al., 2000; Sakas et.

The interpolation techniques described above, are usually applied only on a single slice level (2D). The use of high resolution CT data, allows the use of multiplanar reconstructions (MPR) for the sagittal and coronal direction, in relation with the patient anatomy. These two

generally causes many times robust automatic algorithms to fail.

Fig. 6. On the left side a simple case where an interpolated contour is created from the plane intersection with the triangles (Zimeras and Karangelis, 2008)

Fig. 7. The orthogonal contour interpolation. In (a), the two drawn contours and the intersection plane in Z direction. The result of the intersection is shown in (b). After applying the cubic Spline interpolation a new contour is created (c) (Zimeras and Karangelis, 2008).

Fig. 6. On the left side a simple case where an interpolated contour is created from the plane

Fig. 7. The orthogonal contour interpolation. In (a), the two drawn contours and the intersection plane in Z direction. The result of the intersection is shown in (b). After applying the cubic Spline interpolation a new contour is created (c) (Zimeras and

Karangelis, 2008).

intersection with the triangles (Zimeras and Karangelis, 2008)

Fully automatic segmentation methods are usually impractical due to image complexity and the variety of image types and interpretation. In addition, low contrast between structures generally causes many times robust automatic algorithms to fail.

*Active Contour Model*: Active Contour Models (ACMs) are adaptive contour representations, also known as snakes or deformable models (Terzopoulos and Fleischer, 1998). They are able to recover and represent physical contours of an image, and hence can be used as a model to determine object boundaries in static images as well as for tracking in image sequences. (Grosskopf et al., 1998; 2002) uses the Euler Time Integration to solve the optimisation problem. After initialisation by a user sketch, the contour is deformed to fit the actual object by simulating physical properties of an elastic material or fluid. This method is very reliable to overcome local minima and very fast due to its deterministic character.

So far only a very few techniques propose physicians one or more alternatives for volume definition (Ketting et. al., 1997; Belsh et. al., 1997). Recently (Pekar et. al., 2004) reported a method based on an adaptation of 3D deformable surface models to the boundaries of the anatomic structures of interest. The adaptation was based on a tradeoff between deformations of the model induced by its attraction to certain image features and the shape integrity of the model. Nevertheless, to make the concept clinically feasible, interactive tools were also introduced that allow quick correction in problematic areas in which the automated model adaptation may fail. Currently the tool used to perform the volume definition is the mouse cursor. The user defines a number of digital points on the image level, closing the first and the last point of the contour to generate this way a structure. The connectivity between the *key points* can be linear or higher order. The higher order connectivity can be achieved using interpolation models like Hermitte cubic, Spline curves, Bezier curves (Laurent, 1994; Spath, 1995; Cohen 2001), which are the most common and successfully used techniques for smoothing curves in the systems. Aim of the high order interpolation techniques is to reduce the amount of input points required to describe a smooth shape. Performing a simple comparison between linear and higher order interpolation, it can find out that the amount of points used to illustrate the shape of a structure using linear interpolation requires at least twice as many samples as by using the higher order interpolation algorithms. A common methodology used to combine high order interpolation and image edge properties is the use of active contour models. The active contour models or Snakes can be 2D image curves (Kass et. al., 1987; Blake and Isard, 1998) or 3D polygon meshes (Terzopoulos and Fleischer, 1998), which are adjusted from an initial approximation to the image or volume features by a movement caused by simulated forces. Image features provide the so-called external force. An internal tension of the curve resists against highly angled curvatures, which makes the Snakes movement robust against noise. After a starting position is given, the Snake adapts itself to shape by relaxation to the equilibrium of the external force and internal tension. Snakes have been proven efficient and fast for a number of applications in medicine involving different imaging modalities (McIne and Terzopoulos, 1996; Gross et. a., 1998; Behr et. al., 2000; Sakas et. al., 2001; Gross et. al., 2002).

The interpolation techniques described above, are usually applied only on a single slice level (2D). The use of high resolution CT data, allows the use of multiplanar reconstructions (MPR) for the sagittal and coronal direction, in relation with the patient anatomy. These two

Segmentation Techniques of Anatomical

Structures with Applicationin Radiotherapy Treatment Planning 51

Fig. 8. Combination between Manual segmentation methods (between two slices and ACM (voxel reconstruction) on CT slices. In (a) from left to right: coronal, axial and sagittal

The algorithm can automatically trace the organ through the complete volume of cross sections. False contours that are not corresponding to the spine shape and position can be rejected automatically from the system and can be replaced with linear interpolated contours considering as key contours those already found by the system. The boundarytracking methods used, belong to the deterministic approaches and therefore there is the tendency to produce misleading results under some circumstances. To reduce that effect data pre-processing and the gradient volume of the original CT data can be used as input to the segmentation routine. Target volume and critical structure definition is a complex and time-consuming process in radiotherapy. The complexity varies for different anatomic sites. In plan evaluation, both the physicists and radiation oncologists interact closely to subjectively identify the plan most appropriate for the individual patient. In order to reduce the investment of time and effort by the radiation oncology staff, several image analysis

contours. In (b) surface reconstruction (Zimeras and Karangelis, 2008).

images are orthogonal to each other and perpendicular with the axial plane. The navigation though these images help in the observation of complex anatomy. The sagittal and coronal views often offer a better overview of organs 3D shape. Defining volumes in these directions could be of benefit since several organs are aligned along the longitudinal body axis. The problem we have to solve in our case is the generation of a surface and contours from structured closed parallel and non-parallel contours. The data points represent the contour points as they are generated from the user on the different levels of the axial slices or/and on the MPRs. The most common approaches used to reconstruct surfaces form parallel contours and are well established in medical imaging applications (Boissonat, 1988; Meyer et. al. 1992; Payne and Toga, 1994; Bajaj et. al. 1995; Weinstein, 2000). The limitation of these algorithms is that they cannot be applied on non-parallel contours. The problem of the nonparallel contours could be also formulated as generation of surfaces from scatter data, which are very common in industrial applications (Hoppe et. al. 1992; Ament et. al. 1998). For our application, we selected the approach presented from Turk (Turk and O'Brien 1999). Their method adapts earlier work on thin plate splines and radial basis interpolants to create a new technique in generating implicit surfaces. Their method allows direct specification of a complex surface from sparse, irregular scatter samples. The main restriction of the method is the relatively small number of sample data that can be handled. This drawback makes the above approach unsuitable for a number of applications that a large number of sampling points are needed. In this work, we demonstrate the use of implicit function interpolation to reconstruct 3D organ shapes from closed planar parallel and non-parallel contours that have been defined selectively by the user. The total number of contour points will be used as the input data to the implicit surface algorithm with arbitrary order. The number of these sampling points will not exceed the level of few hundred, and therefore the calculation times will be in acceptable ranges despite the complexity of the algorithm (Figure 8).

The output result of the reconstruction algorithm is provided in two forms: as a triangulated mesh or as multiple parallel contours extracted in arbitrary plane directions. For the RTP applications we focus mostly on the reconstruction of contours in the axial direction. The algorithm can be separated into different modules from the point editing until the reconstruction of the surface and contours as follows:


Semi-automatic segmentations methods combine the benefit of bath manual and automatic segmentation techniques. By supplying initial information about the region of interest, the user may guide an otherwise automatic segmentation procedure. Any remaining errors introduced by automatic segmentation methods may be corrected by manual editing. In this work, a boundary tracking algorithm (Karangelis *et al.*, 2001, Zimeras and Karangelis, 2001, 2002) was implemented for the segmentation part.

images are orthogonal to each other and perpendicular with the axial plane. The navigation though these images help in the observation of complex anatomy. The sagittal and coronal views often offer a better overview of organs 3D shape. Defining volumes in these directions could be of benefit since several organs are aligned along the longitudinal body axis. The problem we have to solve in our case is the generation of a surface and contours from structured closed parallel and non-parallel contours. The data points represent the contour points as they are generated from the user on the different levels of the axial slices or/and on the MPRs. The most common approaches used to reconstruct surfaces form parallel contours and are well established in medical imaging applications (Boissonat, 1988; Meyer et. al. 1992; Payne and Toga, 1994; Bajaj et. al. 1995; Weinstein, 2000). The limitation of these algorithms is that they cannot be applied on non-parallel contours. The problem of the nonparallel contours could be also formulated as generation of surfaces from scatter data, which are very common in industrial applications (Hoppe et. al. 1992; Ament et. al. 1998). For our application, we selected the approach presented from Turk (Turk and O'Brien 1999). Their method adapts earlier work on thin plate splines and radial basis interpolants to create a new technique in generating implicit surfaces. Their method allows direct specification of a complex surface from sparse, irregular scatter samples. The main restriction of the method is the relatively small number of sample data that can be handled. This drawback makes the above approach unsuitable for a number of applications that a large number of sampling points are needed. In this work, we demonstrate the use of implicit function interpolation to reconstruct 3D organ shapes from closed planar parallel and non-parallel contours that have been defined selectively by the user. The total number of contour points will be used as the input data to the implicit surface algorithm with arbitrary order. The number of these sampling points will not exceed the level of few hundred, and therefore the calculation

times will be in acceptable ranges despite the complexity of the algorithm (Figure 8).

reconstruction of the surface and contours as follows:

representing our interpolation function.

2002) was implemented for the segmentation part.

3D polygon meshes respectively.

of the contour constraints and filtering of unwanted.

The output result of the reconstruction algorithm is provided in two forms: as a triangulated mesh or as multiple parallel contours extracted in arbitrary plane directions. For the RTP applications we focus mostly on the reconstruction of contours in the axial direction. The algorithm can be separated into different modules from the point editing until the

1. Collection and processing of the given contour points. This step involves the generation

2. Calculation of the implicit functions in 3D. In this step the information produced in step one will be used to solve a linear equation system that will provide the coefficients

3. Evaluating the implicit function over a 2D or 3D grid we extract 2D planar contours or

Semi-automatic segmentations methods combine the benefit of bath manual and automatic segmentation techniques. By supplying initial information about the region of interest, the user may guide an otherwise automatic segmentation procedure. Any remaining errors introduced by automatic segmentation methods may be corrected by manual editing. In this work, a boundary tracking algorithm (Karangelis *et al.*, 2001, Zimeras and Karangelis, 2001,

Fig. 8. Combination between Manual segmentation methods (between two slices and ACM (voxel reconstruction) on CT slices. In (a) from left to right: coronal, axial and sagittal contours. In (b) surface reconstruction (Zimeras and Karangelis, 2008).

The algorithm can automatically trace the organ through the complete volume of cross sections. False contours that are not corresponding to the spine shape and position can be rejected automatically from the system and can be replaced with linear interpolated contours considering as key contours those already found by the system. The boundarytracking methods used, belong to the deterministic approaches and therefore there is the tendency to produce misleading results under some circumstances. To reduce that effect data pre-processing and the gradient volume of the original CT data can be used as input to the segmentation routine. Target volume and critical structure definition is a complex and time-consuming process in radiotherapy. The complexity varies for different anatomic sites. In plan evaluation, both the physicists and radiation oncologists interact closely to subjectively identify the plan most appropriate for the individual patient. In order to reduce the investment of time and effort by the radiation oncology staff, several image analysis

Segmentation Techniques of Anatomical

1 , I(i,j) t (, ) 0 , otherwise background *if T objec Iij* 

segmentation technique is presented below (Figure 11):

for i=1 to Whole\_Voxel

find the starting point

step 1 if there is an object then

change starting point

}while region becomes small

define the threshold labelling

{

 do {

on the rule (Figure 10):

Fig. 10. Thresholding cases

Fig. 11. Pseudo-code for segmentation

else

goto step 1

surfaces from scalar volume data in volume rendering.

Structures with Applicationin Radiotherapy Treatment Planning 53

thresholding is a technique for converting a grayscale or color image to a binary image based upon a threshold value. If a pixel in the image has intensity less than the threshold value, the corresponding pixel in the resultant image is set to white. Otherwise, if the pixel is greater than or equal to the threshold intensity, the resulting pixel is set to black. For that purpose, a low pass filter was used based on the rectangular window or box function based

where I(i,j) image matrix and T is the appropriate threshold. The output is the label "object" or "background" which, due to its dichotomous nature, can be represented as a Boolean variable "1" or "0".Summarizing the above procedures, the pseudo-code for the

> find sharp edge using the binary image (0,1) then define initial directions (UP and RIGHT) calculate contour for the region of interest apply contour filtering mask (Sobel filter)

 reduce the number of points calculate contour measures

The main drawback of the method is that is a binary approach and hence is very sensitive to gray value variations. If the threshold value is not selected properly then the system will fail to detect the appropriate organ shape. Most of the inaccuracies of the segmentation method require the user intervention to optimize the result. To overcome the limitation we calculate a secondary opacity volume from the original CT data that is very often used to visualize

tools are integrated. A function that significantly accelerates the contouring process is the linear interpolation between the original key-contours. The same principle can be applied for defining structures in both planar planes, sagittal and coronal. Organs with large differences in there intensities can be segmented semi-automatically. In terms of user effort the only action required from the user is the selection of an initial point from the algorithm on the original axial slices. The complete 3D geometry of the organ will be traced automatically. Some of the common organs with high sensitivity factor and vital importance are the lungs, the spinal cord and the trachea (Zimeras and Karangelis, 2001; Karangelis and Zimeras, 2002a, b; Zimeras *et al.* 2002). In addition to those organs, the external body contour can be extracted in a similar manner. The contours that are generated semi-automatic can be manipulated and modified at the same manner as those defined manually (Figure 9a).

After segmenting the regions of interest, second task to use a 3x3 mask edge detection operator (Figure 9b) for improving the appearance of the contours within these regions. Edge detection operators examine each pixel neighborhood and quantify the slope, and the direction of the grey-level transition. There are several ways to do this, most of which are based upon convolution with a set of directional derivative masks. For that purpose, a Sobel edge kernel was applied for finding sharp region boundaries, especially for these, which are changing greatly in intensity over short image distances. The two convolution kernels are given by:

$$S\_x = \begin{bmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{bmatrix} \text{ and } S\_y = \begin{bmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{bmatrix} \tag{6}$$

Each point in the image is convoluted with both kernels. One kernel ( *Sx* ) responds maximally to a generally vertical edge and the other ( *Sy* ) to a horizontal edge. The maximum value of the two convolutions is taken as the output value for that pixel.

In the analysis of the objects in images it is essential that we can distinguish between the objects of interest and "the rest." This latter group is also referred to as the background. The techniques that are used to find the objects of interest are usually referred to as segmentation techniques - segmenting the foreground from background. Image thresholding is a technique for converting a grayscale or color image to a binary image based upon a threshold value. If a pixel in the image has intensity less than the threshold value, the corresponding pixel in the resultant image is set to white. Otherwise, if the pixel is greater than or equal to the threshold intensity, the resulting pixel is set to black. For that purpose, a low pass filter was used based on the rectangular window or box function based on the rule (Figure 10):

#### Fig. 10. Thresholding cases

52 Modern Practices in Radiation Therapy

tools are integrated. A function that significantly accelerates the contouring process is the linear interpolation between the original key-contours. The same principle can be applied for defining structures in both planar planes, sagittal and coronal. Organs with large differences in there intensities can be segmented semi-automatically. In terms of user effort the only action required from the user is the selection of an initial point from the algorithm on the original axial slices. The complete 3D geometry of the organ will be traced automatically. Some of the common organs with high sensitivity factor and vital importance are the lungs, the spinal cord and the trachea (Zimeras and Karangelis, 2001; Karangelis and Zimeras, 2002a, b; Zimeras *et al.* 2002). In addition to those organs, the external body contour can be extracted in a similar manner. The contours that are generated semi-automatic can be manipulated and modified at the same manner as those defined manually (Figure 9a).

Fig. 9a. Boundary tracking method. Fig. 9b. Contour filtering

121 000 121

*Sx*

After segmenting the regions of interest, second task to use a 3x3 mask edge detection operator (Figure 9b) for improving the appearance of the contours within these regions. Edge detection operators examine each pixel neighborhood and quantify the slope, and the direction of the grey-level transition. There are several ways to do this, most of which are based upon convolution with a set of directional derivative masks. For that purpose, a Sobel edge kernel was applied for finding sharp region boundaries, especially for these, which are changing greatly in intensity over short image distances. The two convolution kernels are given by:

and

Each point in the image is convoluted with both kernels. One kernel ( *Sx* ) responds maximally to a generally vertical edge and the other ( *Sy* ) to a horizontal edge. The

In the analysis of the objects in images it is essential that we can distinguish between the objects of interest and "the rest." This latter group is also referred to as the background. The techniques that are used to find the objects of interest are usually referred to as segmentation techniques - segmenting the foreground from background. Image

maximum value of the two convolutions is taken as the output value for that pixel.

*Sy*

101 202 101

(6)

where I(i,j) image matrix and T is the appropriate threshold. The output is the label "object" or "background" which, due to its dichotomous nature, can be represented as a Boolean variable "1" or "0".Summarizing the above procedures, the pseudo-code for the segmentation technique is presented below (Figure 11):


Fig. 11. Pseudo-code for segmentation

The main drawback of the method is that is a binary approach and hence is very sensitive to gray value variations. If the threshold value is not selected properly then the system will fail to detect the appropriate organ shape. Most of the inaccuracies of the segmentation method require the user intervention to optimize the result. To overcome the limitation we calculate a secondary opacity volume from the original CT data that is very often used to visualize surfaces from scalar volume data in volume rendering.

Segmentation Techniques of Anatomical

bronchus.

**5. Conclusions** 

Structures with Applicationin Radiotherapy Treatment Planning 55

Fig. 13. Segmentation of different anatomical structures. 2D segmentation of the left and right lung; 3D segmentation of the left, right lung and spinal cord, 3D segmentation of the

In this work, an effective semi-automatic method was presented, based on the *boundary tracking technique*, that improves the time when one or more structures are in use. The implemented algorithms can segment within a few seconds the complete volume of specific organs e.g. lungs, skin, spinal canal, bronchus and brain. The only interaction of the user is to select the starting point in the region of interest and the algorithm will track the object boundaries in 3 dimensions. For each particular organs of interests (lungs, skin, spine canal, bronchus and brain), a different segmentation techniques is proposed, but all of them are based on the *boundary tracking technique.* The 3D-shape is reconstructed out of the segmented regions, in order to illustrate the efficiency of the segmentation techniques. The benefit of the 3D illustration is the visual presentation of the organs. In many medical cases,

In order to speed up the discrete contour drawing mode the user may draw a smaller number of key contours on distinct slices while the intermediate contours are interpolated linear. To perform this, a slope factor **A** is calculated between the key slices. The value of this is given from the equation:

$$A = \frac{(\text{Current\\_slice\\_pos - Precision\\_contor\\_pos})}{(\text{Previous\\_contour\\_pos - Next\\_contor\\_pos})}\tag{7}$$

Since each key contour does not have the same number of segment points the algorithm subdivides each contour into the same number of points (currently 100) in order to simplify the contour interpolation step. Then a starting point is selected -- usually this is the point of the 12o'clock position of the contour (Figure 3). To calculate the X and Y pixel (or voxel) position of the interpolated point we use the following formula:

$$\begin{aligned} \mathbf{X\_{int}} &= \mathbf{X\_{previous\\_cut}} + \mathbf{A^\*(X\_{next\\_cut} - X\_{previous\\_cut})} \\ \mathbf{Y\_{int}} &= \mathbf{Y\_{previous\\_cut}} + \mathbf{A^\*(Y\_{next\\_cut} - Y\_{previous\\_cut})} \end{aligned} \tag{8}$$

For visualizing volumetric medical data, sequences of 2D images are piled up to recreate the three-dimensional structure. Usually, in this step linear interpolation of adjacent slices is needed for the generation of new slices (Figure 12), since one of the problems resulting from image acquisition is the space between slices. This problem occurs because the sampling interval between slices is normally greater than the generated image resolution, and then the volume voxels are not cubic. After interpolation this size distortion is corrected, so that the visualization algorithm could generate correct proportion projections.

Fig. 12. Slice interpolation

Figure 13 illustrates different 2D and 3D reconstruction examples of segmented structures.

In order to speed up the discrete contour drawing mode the user may draw a smaller number of key contours on distinct slices while the intermediate contours are interpolated linear. To perform this, a slope factor **A** is calculated between the key slices. The value of

> (Current\_slice\_pos - Previous\_contour\_pos) (Previous\_contour\_pos - Next\_contour\_pos)

Since each key contour does not have the same number of segment points the algorithm subdivides each contour into the same number of points (currently 100) in order to simplify the contour interpolation step. Then a starting point is selected -- usually this is the point of the 12o'clock position of the contour (Figure 3). To calculate the X and Y pixel (or voxel)

> X X A\*(X X ) int previous\_cnt next\_cnt previous\_cnt Y Y A\*(Y Y ) int previous\_cnt next\_cnt previous\_cnt

For visualizing volumetric medical data, sequences of 2D images are piled up to recreate the three-dimensional structure. Usually, in this step linear interpolation of adjacent slices is needed for the generation of new slices (Figure 12), since one of the problems resulting from image acquisition is the space between slices. This problem occurs because the sampling interval between slices is normally greater than the generated image resolution, and then the volume voxels are not cubic. After interpolation this size distortion is corrected, so that the

Figure 13 illustrates different 2D and 3D reconstruction examples of segmented structures.

position of the interpolated point we use the following formula:

visualization algorithm could generate correct proportion projections.

*A* (7)

(8)

this is given from the equation:

Fig. 12. Slice interpolation

Fig. 13. Segmentation of different anatomical structures. 2D segmentation of the left and right lung; 3D segmentation of the left, right lung and spinal cord, 3D segmentation of the bronchus.

### **5. Conclusions**

In this work, an effective semi-automatic method was presented, based on the *boundary tracking technique*, that improves the time when one or more structures are in use. The implemented algorithms can segment within a few seconds the complete volume of specific organs e.g. lungs, skin, spinal canal, bronchus and brain. The only interaction of the user is to select the starting point in the region of interest and the algorithm will track the object boundaries in 3 dimensions. For each particular organs of interests (lungs, skin, spine canal, bronchus and brain), a different segmentation techniques is proposed, but all of them are based on the *boundary tracking technique.* The 3D-shape is reconstructed out of the segmented regions, in order to illustrate the efficiency of the segmentation techniques. The benefit of the 3D illustration is the visual presentation of the organs. In many medical cases,

Segmentation Techniques of Anatomical

*Berlin, Germany, 366-369*

pp:425-438, 1993.

Darmstadt, 2004.

1990;18:505-513.

Vol: 11, pp: 228-258, 1992

Structures with Applicationin Radiotherapy Treatment Planning 57

G. Karangelis and S. Zimeras (2002a). An accurate 3D segmentation method of the spinal

G. Karangelis and S. Zimeras (2002b). 3D segmentation method of the spinal cord applied

G. Karangelis, S. Zimeras, E. Firle, Min Wang and G. Sakas, (2001). Volume definition tools

G. Sakas, G. Karangelis, and A. Pommert; Advanced Applications of Volume Visualization

G. Sakas; Interactive Volume Rendering of Large Fields, The Visual Computer, Vol: 9,

G. Turk, and J.F. O'Brien; Shape Transformation Using Variational Implicit Functions,

H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, and W. Stuetzle; Surface Reconstruction from Unorganised Points, Computer Graphics, SIGGRAPH' 1992, pp: 71-78, 1992

J Strassman G., Kolotas C., Heyd R.: Navigation system for interstitial brachytherapy,

J. Behr, S.M. Choi, S. Großkopf, H. Hong, S.A. Nam, Y. Peng, A. Hildebrand, M.H. Kim, and G.

Karangelis G. 3D Simulation of external beam radiotherapy, Ph.D. Thesis, University of

M. Kass, A. Witkin, and D. Terzopoulos; Snakes: Active Contour Models, IEEE First Int.

Meyers, S. Skinner, and K. Sloan; Surfaces from contours, ACM Transactions On Graphics,

Michalski, J.M.; Purdy, J.A.; Harms, W.; Matthews, J.W. The CT-Simulation 3D Treatment

N. Amenta, M. Bern, and M. Kamvysselis; A New Voronoi-Based Surface Reconstruction Algorithm, Computer Graphics, SIGGRAPH' 1998, pp: 415-421, 1998 N. Max; Optical Models For Direct Volume Rendering, IEEE Transactions On Visualization

Nagata, Y.; Nishidai, T.; Abe, M.; Takahashi, M.; Okajima, K.; Yamaoka, N.; Ishihara, H.;

Nishidai, T.; Nagata, Y.; Takahashi, M.; Abe, M.; Yamaoka, N.; Ishihara, H.; Kubo, Y.; Ohta, H.;

Kubo, Y.; Ohta, H.; Kazusa, C. CT Simulator: A new 3-D planning and simulating system for radiotherapy: Part 2. Clinical application. *Int J Rad Oncol Biol Phys*

Kazusa, C. CT Simulator: A new 3-D planning and simulating system for radiotherapy: Part 1. Description of system. *Int J Rad Oncol Biol Phys* 1990;18:499-504.

Sakas; Modeling, Visualization, and Interaction Techniques for Diagnosis and Treatment Planning in Cardiology, Computer & Graphics, Vol: 24, pp: 741-753, 2000 J. Boissonat; Shape Reconstruction from Planar Cross Sections, Computer Vision, Graphics

by Stergios Stergiopoulos, CRC Press, ISBN 0-8493-3691-0, 2001

H. Spath; Two Dimensional Spline Interpolation Algorithms, A.K. Peters Ltd, 1995

Computer Graphics, SIGGRAPH'1999, pp. 335-342, 1999

Radiotherapy & Oncology, Vol 56, pp:49-57, 2000

and Image Processing, Vol: 44, pp:1-29, 1988

Conf. Comput. Vision, pp: 259-268, 1987

M. Levoy et al.(1988). Display of surface from volume data. *IEEE CGA*, **8**.

Planning Process. *Front Radiat Ther Oncol* 1996;29:43-56.

And Computer Graphics, Vol: 1, pp:99-108, 1995

on CT data, *Computer Graphics Topics, Vol 14, 28-29.*

*in Computer Sciences 2208, 1295-1297.*

canal applied on CT images *BVM 2002, Conference Proceedings, Springer-Verlag,* 

for medical image applications, *4th MICCAI International Conference, Utrecht, Netherlands, M-A. Viergever, T. Dohi, M. Vannier (eds), Springer-Verlag, Lecture Notes* 

Methods in Medicine. In Advanced Signal Processing Handbook, Theory and Implementation for Radar, Sonar, and Medical Imaging Real-Time Systems, Edited

illustration of the organ involves an anomaly, clinical problem or generally artifacts. Visual representation of the particular organ, in addition with the clinical examinations, could be a powerful tool to the doctors for diagnosis, medical treatment or surgery.

#### **6. Acknowledgement**

The author would like to thank Dr. Grigoris Karangelis and Prof Nikos Zamboglou for the useful scientific help and comments about the progress of this work. Also many thanks to Prof G. Sakas MedCom Company and Städtisches Klinikum Offenbach whom gave me equipments and medical data sets for the implementation of the above work. Finally the author would like to thank the EC and the Marie Curie Fellowship Association for this great opportunity to work with different people. This work had been supported by a Marie Curie Industry Host Fellowship Grant no: HPMI-CT-1999-00005.

#### **7. References**


illustration of the organ involves an anomaly, clinical problem or generally artifacts. Visual representation of the particular organ, in addition with the clinical examinations, could be a

The author would like to thank Dr. Grigoris Karangelis and Prof Nikos Zamboglou for the useful scientific help and comments about the progress of this work. Also many thanks to Prof G. Sakas MedCom Company and Städtisches Klinikum Offenbach whom gave me equipments and medical data sets for the implementation of the above work. Finally the author would like to thank the EC and the Marie Curie Fellowship Association for this great opportunity to work with different people. This work had been supported by a Marie Curie

A. Blake, and M. Isard; Active Contours: The Application of Techniques from Graphics,

B.A. Payne, and A.W. Toga; Surface Reconstruction By Multiaxial Triangulation, IEEE

B.S. Kuszyk, D.R. Ney, and E.K. Fishman; The Current State of the Art in Three Dimensional

Butker, E.K.; Helton, D.J.; Keller, J.W.; Hughes, L.L.; Crenshaw, T.; Davis, L.W. A totally

C. Bajaj, F. Bernardini, and G. Xu; Automatic Reconstruction of Surfaces and Scalar Field from 3D Scans, Computer Graphics, SIGGRAPH' 1995,pp: 109-118, 1995 C.H. Ketting, M.A. Seymour, I. Kalet et al.; Automated Planning Target Volume Generation:

Cai, W. Transfer functions in DRR volume rendering. CARS'99, Paris, France, June 1999, 23-26. Cai, W.; Sakas, G.; Karangelis, G. Volume Interaction Techniques in the Virtual Simulation

Chen, G.T.Y.; Pelizzari, C.A.; Vijayakumar, S. Imaging: The Basis for Effective Therapy.

D. Terzopoulos, and K. Fleischer; Deformable Models, The Visual Computer, Vol: 4, pp: 306-

D. Weinstein; Scanline Surfacing: Building Separating Surfaces from Planar Contours, In

E. Cohen, R.F. Riesenfeld, and G. Elber; Geometric Modeling with Splines: An Introduction,

Computer Graphics and Applications, Vol: 14, pp:28-35, 1994

Radiation Oncology Biol. Phys., Vol: 37, pp:697-704, 1997

Conway,. J.; Robinson, M.H. CT virtual simulation. *Brit J Radiol* 1997;70:106-118.

Graphics and Vision (Graphicon) Moscow, 1999.

Proceeding of 11th IEEE Visualization'2000, 2000

*Front Radiat Ther Oncol* 1996;29:31-42.

Vision, Control Theory and Statistics to Visual Tracking of Shapes in Motion.

Oncologic Imaging: An Overview, Int. J. Radiation Oncology Biol. Phys., Vol:33,

integrated simulation technique for three field breast treatment using a CT

An Evaluation Pitting A Computer-Based Tool Against Human Experts, Int. J.

of Radiotherapy Treatment Planning. International Conference on Computer

powerful tool to the doctors for diagnosis, medical treatment or surgery.

Industry Host Fellowship Grant no: HPMI-CT-1999-00005.

Springer-Verlag, London, 1998

simulator. *Med Phys* 1996;23:1809-1814.

pp:1029-1039, 1995

331, 1988

A.K. Peters Ltd, 2001

**6. Acknowledgement** 

**7. References** 


**4** 

*Japan* 

Tomoki Kimura

**Involved-Field Radiation Therapy (IF-RT)** 

**for Non-Small Cell Lung Cancer (NSCLC)** 

**1.1 Patterns of lymphatic spread of Non-Small Cell Lung Cancer (NSCLC)** 

The pattern and incidence of lymphatic spread of non-small cell lung cancer (NSCLC) is differentiated according to location, size and histologic type of primary tumors. As the primary tumor size increases, the incidence of lymphatic metastasis increases. Ogata et al reported the incidence of mediastinal lymphatic metastasis increased from 24% for tumors under 2 cm in size to more than 40% for tumors larger than 5 cm 1). Hata et al analyzed 192 lymphoscintigraphies in 179 patients to determine the lymphatic drainage from each segmental bronchus into the mediastinum 2). For the right lobes, most of the lymph flowed into the right supraclavicular nodes through the subcarina or right paratracheal nodes. There were few drainages to the left supraclavicular nodes through subcarinal nodes. In contrast, the lymphatic drainage from the left lung was more variable, and four routes were

2. The route runs along the left phrenic nerve through the para-aortic nodes to the left

3. The route runs along the left main bronchus to the left hilar or the left prevascular nodes. From the left hilar nodes, this route divides into two branches. One extends to the right supraclavicular nodes through the right upper paratracheal nodes. The other

4. The route runs under the left main bronchus to the subcarinal nodes. After passing the subcarinal nodes, this route extends to the right supraclavicular nodes along the trachea through the pretracheal nodes and the left upper paratracheal nodes. Some branches extend upwards along the left side of the trachea to the upper paratracheal and

Fig. 1 shows the standard patterns of lymphatic drainage 2). These results suggest that most of mediastinal lymph node metastasis were found ipsilaterally in the right primary lung cancer, and the mediastinal or supraclavicular lymph node metastasis were found bilaterally in the left primary lung cancer. Nohl-Oser examined the location of nodal involvement in 749 patients based on data obtained via mediastinoscopy, scalene lymph node biopsy and surgical specimen (Table 1) 3). The right upper lobe tumors spread to the right upper and

**1. Introduction** 

determined, as follows.

supraclavicular nodes.

supraclavicular nodes.

1. The route through the subaortic nodes.

runs upwards along the left side paratracheal nodes

*Department of Radiation Oncology, Hiroshima University,* 

 *Graduated School of Biomedical Sciences,* 

