Image Anal Stereol 2011;30:179-186 Original Research Paper A SEGMENTATION PROBLEM IN QUANTITATIVE ASSESSMENT OF ORGAN DISPOSITION IN RADIOTHERAPY Giovanni Naldi1, Barbara Avuzzi2, Simona Fantini3, Mauro Carrara4, Ester Orlandi3, Elisa Massafra1 and Stefano Tomatis4 1Dipartimento di Matematica, Universita degli studi di Milano, Via Saldini 50, Milan, Italy; 2Prostate program, Scientific directorate, Fondazione IRCCS Istituto Nazionale Tumori, via Venezian 1, 20133 Milano, Italy; 3Department of Radiotherapy, SC RT2, Fondazione IRCCS Istituto Nazionale Tumori, via Venezian 1, 20133 Milano, Italy; 4Medical Physics Unit, Fondazione IRCCS Istituto Nazionale Tumori, Via Venezian, 1, Milan, Italy e-mail: giovanni.naldi@unimi.it, barbara.avuzzi@istitutotumori.mi.it, simona.fantini@istitutotumori.mi.it, mauro.carrara@istitutotumori.mi.it, ester.orlandi@istitutotumori.mi.it, stefano.tomatis@istitutotumori.mi.it (Accepted October 18, 20II) ABSTRACT Radiotherapeutic treatment of cancer is best conducted if the prescription dose is given to the tumor while surrounding normal tissues are maximally spared. With the aim to meet these requirements the complexity of radiotherapy techniques have steadily increased under a strong technological impulse, especially in the last decades. One problem involves the rate of the particular disposition of the structures of interest in a patient. Recently the authors (Tomatis et al., 2010; 2011) have proposed a computational approach in order to represent quantitatively the geometrical features of organs at risk, summarized in characteristics of distance, shape and orientation of such organs in respect to the target. A basic problem to solve before to compute the risk index, is the segmentation of the organs involved in the radiotherapy planning. Here we described a 3D segmentation method by using the clinical computed tomography (CT) data of the patients. Our algorithm is based on different steps, a preprocessing phase where a nonlinear diffusion filter is applied; a level set based method for extract 2D countours; a postprocessing reconstruction of 3D volume from 2D segmented slices. Some comparisons with manually traced segmentation by clinical experts are provided. Keywords: geometrical features, level set, radiotherapy planning, segmentation, 3D reconstruction. INTRODUCTION Radiotherapy is a treatment of cancer and other diseases with ionizing radiation. Ionizing radiation deposits energy that injures or destroys cells in the area being treated, the "target", by damaging their genetic material, making it impossible for these cells to continue to grow. Although radiation damages both cancer cells and normal cells, the latter are, in general, able to repair themselves and function properly. However, Radiotherapeutic treatment of cancer is best conducted if the prescription dose is given to the tumor while surrounding normal tissues are maximally spared. With the aim to meet these requirements the complexity of radiotherapy techniques have steadily increased under a strong technological impulse, especially in the last decades (Connell and Hellman, 2009). The development of computer and of medical imaging provides the possibility of new approaches to improve the balance of coverage of target versus normal tissues, with computer aided optimization and control of planning, delivery and verification. As a fundamental part of the complex process of realizing a radiotherapy treatment, planning is often a demanding task. The planner work load depends on the employed technique and dose constraints but, at the same time, on the spatial distribution of the designed anatomical structures, with a major care needed when critical organs are closer to the target volume. But, how can the particular disposition of the structures of interest in a patient be rated? A complete answer to this basic question has not been provided yet, but it would allow to start a quantitative assessment of the influence of structures geometries on the treatment plan results. A desirable consequence of this evaluation would be a better control for plan optimization and a standardization of the planning procedures and results. In recent papers (Tomatis etal., 2010; 2011) a new computational approach is presented based on what the authors call Expansion Intersection Histogram (EIH), a function defined as the intersection between an organ at risk and the target volume, while the target is expanded in 3D through a scanning procedure. This index is able to represent quantitatively the geometrical features of organs at risk (OARs), summarized in characteristics of distance, shape and orientation of such organs in respect to the target. Fig. 1. Three dimensional view of anatomical volumes manually contoured for a radiotherapy treatment in the head and neck region. Besides the tumour (target) organs shown are spinal cord, brainstem and omolateral parotid gland. A 3D representation of anatomical structures in a typical clinical situation is shown in Fig.1 where a real case was selected in the head and neck region. Although potentially all the contoured volumes could be taken into account, we focused on the high dose target for which the planned dose was 70 Gy and on spinal cord, brainstem and right parotid selected as OARs. in order to quantify the geometrical features of patient volumes, a necessary step is the segmentation of the organ of interest. This task is currently performed by the radiation oncologist who traces the organ contours on each slice of the CT data manually (in some cases interactively). These procedures are time consuming. The problem to develop an automatic or semi-automatic method for several organs (for example for abdominal organs) from CT data is still open due to many factors. in particular, these factors include the high inter- and intra-patient variability of shapes and gray-levels, some factors that limit the image quality (as an example, beam hardening, reconstruction artifacts, and patient movements). Moreover, different organs are often represented with similar gray-levels due to their similar tissue densities. Many contributions have been made to the field of automatic and semi-automatic segmentation (Pham etal., 2000; Bankman etal., 2008) for biomedical applications. However, the complicated structures found in medical imaging offer several unsolved challenges to automated algorithms, including the lack of global morphological characteristics, scanner noise and artifacts, and an incomplete or weak separation between points representing neighboring tissue. Then the same algorithm which gives excellent results for one application, might not even work for another (for example, the segmentation of lungs has different issues than the segmentation of liver). Any of the existent automated algorithms can be shown to fail on certain dataset for reasons specific to each method. Here we focus on a semi-automatic gray-level based segmentation framework inspired by recent works on liver segmentation (Pratissoli, 2009; Campadelli etal., 2010). Our aim is to develop a simple segmentation algorithm based on established, and not critical, anatomical knowledge. Then it can be easily generalized, and adapted to segment different organs, overcoming problems due to the high inter-and intra-patient gray-level and shape variabilities. Moreover we start form the clinical framework, for this we require, - to trace the organ contours on each slice of the CT data (as was done by medical experts), the extracted volumes are then combined to produce the final results. in this way we can make a comparison between the automatic method and the experts results with the ability to easily test the method (and eventually correct). - To have some "analytical control" on each slice contours in order to perform geometrical and morphological operations. We point out that with our real CT data was not possible to proceed with a direct 3D (evolution of a surface) segmentation algorithm due to the difference in resolution within each 2D slice and the distance between slices. Then, in order to consider direct 3D segmentation, we would have had to introduce an interpolation process which may typically introduce artifacts. We illustrate the method by using the bladder as an internal organ of interest (at risk). in any case, our method appears to have good potential to become a general segmentation method (while many existing algorithms are specific to a particular problem). A SEMI-AUTOMATIC SEGMENTATION We address the problem of semi-automatic organs segmentation from computed tomography (CT) images in order to compute a risk index for radiotherapy planning. At present this task is usually done by experts who trace the liver contour using some graphical interface: this procedure is very time-consuming and affected by operator's errors and biases. Our system, inspired by work on liver segmentation (Campadelli er a/., 2010), is a semiautomatic framework which only requires manual intervention only in the initialization phase. The method is based on the following steps: preprocessing, noise reduction, fast marching (FM) contour evolution, postprocessing. PREPROCESSING The CT image is much larger than the region containing the organ of interest (in our case the bladder). Then, each CT slice is resized in order to reduce the computational time of the successive steps. First, we perform the binarization of the gray level image by using a suitable thresholding (the threshold level is determined by observing the typical average gray level of tissues). After binarization, we eliminate the small artifact by using morphological filter. Now, we obtain a suitable mask for the selection of the region of interest. In Fig. 2 we show an example of the procedure, the final region is obtained by a simple superposition between the binary mask and the original image. Next, after the resizing of the image, in order to reduce noise and make regions as piecewise constant as possible, without blurring edges, we use a non-linear diffusion filter (Weickert, 1998). This Partial differential equation (PDF) based method is an adaptive smoothing method which is based on the idea of applying a process which itself depends on local properties of the image. The resized image Ires is considered as the initial data of a non-linear anisotropic diffusion equation. |^=div(c(/,V/,x,r)V/) I{x,0)=Ires, xeQ.,t>0, where I = I{x, t) is the image at the time r, Q c M^ is the spatial domain of the image (div is the divergence operator while V/ the gradient of /). Here the diffusion coefficient c depends on the gradient of the image and it is defined as c = c(V/) = 1 i+m/Kf where A' is a suitable positive parameter which may depend on different organs. We point out that it is possible to show that the non-linear diffusion equation is equivalent to a robust procedure that estimates a piecewise-constant image from a noisy input image (Sapiro, 2001). We chose experimentally the parameter ^ = 8, and the final time r = 100. Fig. 2. Resize operation, (A) original CT image; (B) binary image after thresholding; (C) Binary image after morphological operation (opening), the blue line indicates the boundary of the region of interest; (D) resized image. Fig. 3. Non-linear diffusion filter, (A) resized CT image; (B) filtered image, we can observed that now the image is "piecewise smooth ". In Fig. 3 we report a numerical experiment concerning the non-linear diffusion step (here we impose homogeneous Neumann boundary condition). We observe that initial anisotropic smoothing allows a rough initial segmentation. In fact, as explained by Kawohl (2004), there is a connection between the variational approach for segmentation due to Mumford and Shah, and an anisotropic diffusion approach leading to an evolution type equation by Perona and Malik. Then, the anisotropic smoothing reduces noise and provides a first segmentation. Without this step the level set based algorithm produces "jagged curves". We implemented the anisotropic smoothing filter by suitable convolution operations and by using central finite difference schemes in order to approximate the spatial gradient coupling with a classical forward Euler method. The computational cost is proportional to N"^ if N is the number of involved pixel. After the preprocessing step, we obtained for each CT slice a resized and filtered image Ip. Now, we start a segmentation step based on level set methods (Sethian, 2003). CONTOUR EVOLUTION In order to identify the boundary of the bladder (or the region of the organ of interest) in each 2D slice, and to perform image segmentation we consider a curve evolution approach. In particular, we use a level set based, or implicit active contour, approach which is a PDE-based techniques (Osher and Sethian, 1988; Sethian, 2003; Osher and Paragios, 2003). In level set methods, a contour (or more generally a hypersurface) of interest is embedded as the zero level set of a level set function (LSF). Then the contour is moved by suitable image driven forces to the boundaries of the desired objects. In Fig. 4 we show an example of curve evolution, the black line is the starting contour, while the blue line is the final contour (the boundary of the object). Other lines in figure other represent intermediate stages of the contour evolution. The level set method was introduced (Osher and Sethian, 1988) to approximate the propagating interfaces which occur in a wide variety of settings (as physical entities, they include, for example, burning flames, ocean waves, and material boundaries). In the level set approach a "speed function" for the movement of the zero level set must be defined. Typically, this user defined velocity combines a data term with a smoothing term, which prevents the solution from fitting too closely to noise-corrupted data. In our work, we proposed a new speed function of front propagation suitable for our gray-level based segmentation framework. A closed curve F in R^ is an initial position for a front, the level set method takes the perspective of viewing F as the zero level set of a function (p{x,t = 0), from to R. For example, let = 0) = ±d, where d is the distance from x to F, and the plus (minus) sign is chosen if the point x is outside (inside) the initial curve F. Suppose now that the front moves with a (normal) speed F, we can link the evolution of this function (j) to the propagation of the front F(r) itself by a time-dependent initial value problem, namely the following Hamilton-Jacobi equation. Fig. 4. Example of contour evolution, the black line represents the initial front; red and green lines are some of the fronts at intermediate times; the blue line is the final contour where the initial data at time r = 0 is given. There are several advantages to this level set framework: although (^{x,t) remains a function, the zero level propagating surface ^ = 0 may change topology; a discrete grid can be used together with finite differences to devise a numerical scheme to approximate the solution; the intrinsic geometric properties of the front are easily determined from the level set function; finally, the formulation is unchanged for propagating interfaces in multidimensional case. Suppose we now restrict ourselves to the particular case of a front propagating with a speed F that is either always positive or always negative. We then have a monotonically advancing front. In this case, we can convert our level set formulation from a time dependent PDF to a stationary one. In our two-dimensional case, in which the interface is a propagating curve, let T{x,y) be the time at which the curve crosses the point {x,y). The surface T{x,y) then satisfies the equation F\VT\ = 1 , (1) which simply says that the gradient of arrival time surface is inversely proportional to the speed F of the front. In our approach we introduce, as news, a velocity F related to the image properties in order to have as stationary zero level set the boundary of the organ of interest. The following expression for F is considered 1 FHx,y) = a Ux,y)-ß +ßVUx,y) where Ip is the CT slice after the preprocessing step, ß is the average gray level inside the region bounded by the initial curve F, while a and ß are two positive parameters. Heuristically, the speed F gets smaller and smaller near the edges of the organ. Moreover we have a monotone advancing front. Fq. 1 is numerically approximated by using a suitable entropy method. On a 2D grid with an uniform grid spacing of size h the velocity Fij is computed as = of + miniD^^'T, Of + max(D,./r, 0)2 + min(Dj^r, 0)^ , (2) where D^T = - Tij)/h, Dp'T = (Tij^i - Tij)/h. Then, the scheme for the front propagation follows the Fast Marching Level Set Method proposed by Sethian (2003). The initial curve F is manually traced, this is the only manual intervention in the segmentation process. We point out that It would be possible to have a fully automatic method by computing the initial guess curve. We point out that the fast-marching approach can be applied if the front moves always either outward, F > 0, or inward, F <0, so that the evolving front never crosses the same point more than once. The scheme for the curve evolution is as follows (see Sethian, 2003). 1. Inizialization - Set Tij = oo for all grid points; - set Tij = 0 if the point {i, j) is on the initial curve and mark it as done, - set Tij = h/Fij for all points j) in the band, where band is the set of all 4-neighbors of the points on the front. 2. Marching forward. While any point is left in band, - choose the point z — (imin) imin) in the band with the lowest value of Tij; - mark z as done and move its neighbors to band, if they are not already in it or is done; - recompute T for all the neighbors in band using upwind finite difference method. In Fig. 5 we report an example for a single patient of bladder segmentation, for each plane (slice) we perform the preprocessing step and the level set method in order to extract the edge of the bladder. From the numerical experiments good values for the parameters are a = 1, j3 = 2. Moreover, to stop the evolution of the front we control the value of the speed and stop the algorithm when it falls below the minimum tolerance of 0.3, this tolerance value was found experimentally. Fig. 5. Example of bladder segmentation by using different slice from CT data for a patient. In each picture we show the 2D slice and the final contour (green line) at the end of the fast marching algorithm (the final time was determined experimentally). The computational labor involved in the fast marching step is proportional to N^logN, with N points in each grid direction (Sethian, 2003). 3D POSTPROCESSING The last step is the 3D reconstruction of the organ by stacking the collected planar segmented image. First we compute a binary image from each slice where we assign the value one if the pixel belongs to the organ, zero otherwise. We consider binary image because if we do this then apply morphological operators may have algorithms faster and easier implementation. Then we collect all the planar data in order to obtain the finale 3D organ. The distance between the various slices generally varies from 3 to 5 mm, for this reason to have a spatial resolution of 1 mm in the z direction we use local cubic interpolation scheme. In Fig. 6 we have outiined the process of reconstruction from 2D sections to 3D segmentation. Because of not very high resolution of the initial data, the stair-case effect occurred in the binary volume. Since human organs mostiy have smooth structure, it was reasonable to use some smoothing filter. Then, to smooth the 3D contour we apply a convolution with a Gaussian kernel. At this point we have also completed the post-processing and we have the final segmented organ. Fig. 6. The reconstruction of three-dimensional structure of the organ (here the bladder), at the each slice we consider a binary image obtained after the 2D segmentation, then the whole organ is computed by stacking the planar data. NUMERICAL RESULTS AND CONCLUSIONS To evaluate our semi-automatic segmentation algorithms we compared the computed volumes, Vcom, to the ground truth, Vman, manually traced by experts for bladder in 20 patients. This patient population (66.7% of males) was randomly selected from clinical records at Fondazione IRCCS Istituto Nazionale dei Tumori among non-bladder pathologies, (mainly prostate and rectal cancer 73.3%). Median age was 67 years with a 48-87 years range. The CT data were in DICOM (Digital Imaging and Communications in Medicine) format. Each CT scan involves from 10 to 20 slices, the distance between the various slices varies from 3 to 5 mm, while the size of each slices was of 512 pixels with which corresponds a spatial resolution of 0.9375 x 0.9375 mm^. We consider the following measures (here |y| means the counting voxels measure of the volume V): Overlap Sensitivity 0= 100 5= 100 \vr„„nv, com^ K Krom^Vman Symmetric Volume SVO = 100, —^it—n >-/^{\rcom\ + \rman\) Overlap Overlap Error 0^=1- com I I ^ man ^com^^man Volume Difference VD = _ I ^com I I ^man \ \v„ Each of the previous measures is related to the quality of the segmentation but suffers from some drawbacks, for example overlap can not detect over-segmentation errors while sensitivity can not consider the under-segmentation errors. However, the set of the volume measures can give an overall estimate of the goodness of the proposed method. In Table 1 we report the estimates of the volume measures for all patients. From the Table appears the good performance of our segmentation algorithm. Table 1. Segmentation results by using overlap (O), sensitivity (S), symmetric volume overlap (SVO), overlap error (OE), and volume difference (VD)for 20 patients (pat). In the last line we have computed the mean values (Av). Pat. 0(%) S(%) SVO (%) OE (%) VD (%) 1 87.7 86.7 87.2 12.6 1.1 2 96.5 85.5 90.5 17.3 12.5 3 97.3 80.6 88.2 11.1 20.7 4 92.7 87.2 89.8 18.4 6.3 5 95.4 89.7 94.8 10.9 2.3 6 96.3 88.9 94.0 10.5 1.9 7 81.3 75.6 78.3 25.5 7.5 8 91.4 92.4 91.9 14.9 -1.1 9 85.2 78.3 81.6 21.0 8.7 10 92.9 89.9 91.5 15.7 3.3 11 93.1 90.2 92.0 10.4 4.7 12 95.5 87.0 93.8 9.9 2.1 13 95.5 84.7 89.8 18.4 12.7 14 95.1 79.2 86.5 13.8 20.1 15 93.4 86.2 90.1 12.3 3.3 16 90.3 80.8 85.3 15.6 11.7 17 94.6 79.6 86.2 14.2 19.5 18 89.9 87.8 88.0 11.2 7.9 19 93.6 85.7 89.5 19.0 9.2 20 92.8 84.3 90.1 15.6 8.9 Av. 92.5 85.0 89.0 14.9 8.2 Fig. 7. The comparison between the traced manual volume (blue volume) and the volume computed by the semi-automatic method (green volume). In the four pictures it shows the same numerical experiment from different points of view. For a graphical comparison between manual segmentation and our semi-automatic segmentation we show in Fig. 7 a numerical example where the two volume are overlapped. In order to understand where the biggest errors may occur we look to the 2D segmentation. In fact, there is some difference if the slice is at the top, in the middle or at the bottom of the organs. For a quantitative comparison we consider some metric similar to the one used in clustering algorithms. Let Ccom the close contour computed by the level set approach, and Cman the close contour manually traced by medial experts. We define the distance dc\x,T) = minygr x') between a point x and the closed curve r, here d{x,x') is the usual Euclidean distance in R^. Then we consider the average distance