Segmentation Techniques of Anatomical Structures with Application in Radiotherapy Treatment Planning

Cancer is the leading cause of death in economically developed countries and the second leading cause of death in developing countries. It is an enormous global health encumbrance, growing at an alarming pace. Global statistics show that in 2030 alone, about 21.4 million new cancer cases and 13.2 million cancer deaths are expected to occur, simply due to the growth, aging of the population, adoption of new lifestyles and behaviors. Amongst the several modes of treatment for cancer available, Radiation treatment has a major impact due to technological advancement in recent times. This book discusses the pros and cons of this treatment modality. This book "Modern Practices in Radiation Therapy" has collaged topics contributed by top notch professionals and researchers all around the world.


Introduction
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. 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).

Volume rendering
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), xR n , 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 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 I o 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 where I(s) is the attenuated energy after leaving the volume, I o the original energy of the ray, p i the linear attenuation coefficient for the corresponding voxel material, s i 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 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 (p 000 , p 001 , p 011 , p 010 , p 100 , p 110 , p 101 , and p 111 ).
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: where p µ , p water are the attenuation coefficient of the material for the specific voxel and the water respectively. When on the above equation one replaces the p µ =p water , then will receive a HU water = 0. In addition the air as material corresponds to -1000HU, since p air = 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) 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

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).

www.intechopen.com
Segmentation Techniques of Anatomical 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 daily clinical routine are: 1. The irregular tumour properties, like shape, texture, volume, relation with surrounding regions. 2. Location within the human body. A tumour theoretically could grow anywhere within the patient's body. 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 cannot be predicted and detected on the digital patient's data. 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 accurate manual target volume definition is even more necessary nowadays.
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). 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; 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;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 www.intechopen.com 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: 1. Collection and processing of the given contour points. This step involves the generation of the contour constraints and filtering of unwanted. 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 representing our interpolation function. 3. Evaluating the implicit function over a 2D or 3D grid we extract 2D planar contours or 3D polygon meshes respectively.
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  was implemented for the segmentation part. 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 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 Karangelis and Zimeras, 2002a, b;. 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: Each point in the image is convoluted with both kernels. One kernel ( x S ) responds maximally to a generally vertical edge and the other ( y S ) 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 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. 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: 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: 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.

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, 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.

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.