Keywords: Prostate brachytherapy, radioactive seed, CT image, biomedical image segmentation, 3D object location and orientation estimation, mixture model

In this paper, we address the development of an automatic approach for the computation of pose information (position + orientation) of prostate brachytherapy loose seeds from 3D CT images. From an initial detection of a set of seed candidates in CT images using a threshold and connected component method, the orientation of each individual seed is estimated by using the principal components analysis (PCA) method. The main originality of this approach is the ability to classify the detected objects based on a priori intensity and volume information and to separate groups of closely spaced seeds using three competing clustering methods: the standard and a modified k-means method and a Gaussian mixture model with an Expectation-Maximization algorithm. Experiments were carried out on a series of CT images of two phantoms and patients. The fourteen patients correspond to a total of 1063 implanted seeds. Detections are compared to manual segmentation and to related work in terms of detection performance and calculation time. As demonstrated by the results, this automatic method has proved to be accurate and fast including the ability to separate groups of seeds in a reliable way and to determine the orientation of each seed. Such a method is mandatory to be able to compute precisely the real dose delivered to the patient post-operatively instead of assuming the alignment of seeds along the theoretical insertion direction of the brachytherapy needles.

I. INTRODUCTION

ROSTATE cancer is one of the leading cancers in men worldwide with 152 new cases and 23 deaths per 100,000 men reported worldwide per year from 2006-2010 (based on SEER Cancer Statistics Review) [START_REF] Howlader | SEER Cancer Statistics Review, 1975-2010[END_REF]. Low-risk prostate brachytherapy treatment that uses low dose rate radioactive seeds, has emerged as a common and highly effective method to manage localized prostate cancer. The typical implantation procedure is summarized as follows (see Figure 1): based on dose planning, lines of seeds (stranded or loose) are implanted through parallel needles. These needles are inserted into the prostate through the skin of the perineum using continuous transrectal ultrasound (US) guidance and following a preimplantation planning. Once accurate needle placement has been confirmed, the seeds are released through the needles. This process is continued until all seeds have been implanted. In practice, the number of seeds implanted in the prostate commonly ranges from 40 to 100. The goal of a successful operation is to position the seeds in order to get the proper 1 The authors are with TIMC-IMAG laboratory, UJF-Grenoble 1 / CNRS / TIMC-IMAG UMR 5525, Grenoble, F-38041, France. Email: first_name.last_name@imag.fr dose coverage throughout the prostate while limiting the risk for the neighboring organs.

In theory, the seeds are aligned in the needle insertion direction. Figure 2.a illustrates such a planning scheme. However, in practice the seed implantation depends on many biomechanical factors as well as human experience. The seeds may lose their intended position in spite of any special care or effort used when placing the needles and delivering the seeds (Figure 2.b). The examination of images (CT, X-ray, US or MRI) often shows that the seeds are not aligned in the implantation direction especially when using loose seeds. In this later case, the implantation may also result in groups of closely spaced seeds. In this paper, we name such a group a union-seed. Figure 3 shows an example of seed organization in a single CT slice (including one union-seed). In addition, some seeds can migrate out of the prostate as reported in Gao et al. [START_REF] Gao | Effects of seed migration on post-implant dosimetry of prostate brachytherapy[END_REF]. For treatment quality assessment, fluoroscopic images can be acquired immediately after seed implantation. But generally CT data are also acquired one month after the intervention and a CT-based post-implant dosimetry is performed. The delay is such that any inflammatory modification of the prostate has disappeared and the dose computed from the seed positions can be considered as the real delivered dose. Most existing commercial treatment planning software (e.g. VariSeed, Interplant or PSID Brachytherapy software [4]) work under the assumption that all seeds are aligned with the CT axis (as shown in Figure 2.c). However, the American Association of Physicists in Medicine recommends to determine the 3D dose distribution of brachytherapy seeds based on real seed positions and orientations [START_REF] Rivard | Update of aapm task group no. 43 report: A revised aapm protocol for brachytherapy dose calculations[END_REF]. A related clinical question is: "Does taking into account real seed orientation induce significant differences in dose distribution of brachytherapy treatment?". The Dorgipro project we participate in aims to answer this question by comparing the distributed dose calculated using Segmentation, separation and pose estimation of prostate brachytherapy seeds in CT images Huu-Giao Nguyen, Celine Fouard, and Jocelyne Troccaz 1 P standard brachytherapy software (which assumes that the orientation is the planned orientation) to dose distribution obtained when considering the 3D real seed orientation. This paper focuses on the imaging aspects of this research project.

The medical physics aspects are described in a companion paper and will be briefly summarized in section V of this paper. This paper presents the image processing method used to extract the seed positions and orientations, i.e. five parameters (due to the cylindrical shape of the seeds). Our objective is thus to develop a fully automatic software that is able to detect seeds, separate groups of seedsthis stage is also referred to as declustering in the literatureand accurately determine the 5D pose of seeds. In this study, we exploit the high intensity appearance of radioactive seeds in CT images for a solution based on threshold and connected component segmentation [START_REF] Shapiro | Computer Vision[END_REF]; we also consider volume information for outlier removal and computation of the number of seeds in a union-seed. Three declustering methods are considered for union-seed separation: the k-means based method [START_REF] Han | Data Mining: Concepts and Techniques[END_REF] and a modified version of it and the Gaussian mixture model (GMM) with an Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. Finally, the PCA method [START_REF] Jolliffe | Principal Component Analysis[END_REF] for orientation estimation is applied. Reference data coming from the manual segmentation of seeds are used to validate the proposed method.

This paper is organized as follows. We begin by examining the state-of-the-art in prostate brachytherapy seed detection in medical images in Section II. In Section III, we present the proposed solution for seed segmentation (in III.A), union-seed separation (in III.B) and orientation estimation (in III.C) of the prostate brachytherapy seeds in CT images. Experiments and evaluation are reported for datasets generated from phantoms and 14 anonymous patients in Section IV. We then discuss the main contributions and potential extensions of the proposed approach in Section V.

II. RELATED WORK

The accurate localization and orientation estimation of brachytherapy seeds, including the ability to separate unionseeds, is a major challenge and active research field. Numerous studies [START_REF] Corbett | The effect of seed orientation deviations on the quality of 125i prostate implants[END_REF][START_REF] D"amico | An image-guidance system for dynamic dosimetry in prostate brachytherapy using ultrasound and fluoroscopy[END_REF][START_REF] Su | Prostate brachytherapy seed localization by analysis of multiple projections: Identifying and addressing the seed overlap problem[END_REF][START_REF] Wei | Automated localization of implanted seeds in 3d trus images used for prostate brachytherapy[END_REF][START_REF] Singh | Brachytherapy seed localization using geometric and linear programming techniques[END_REF][START_REF] Fallavollita | Registration between ultrasound and fluoroscopy or CT in prostate brachytherapy[END_REF][START_REF] Kuo | Mri-based prostate brachytherapy seed localization[END_REF][START_REF] Lee | Intraoperative 3D reconstruction of prostate brachytherapy implants with automatic pose correction[END_REF][START_REF] Moult | Segmentation of iodine brachytherapy implants in fluoroscopy[END_REF][START_REF] Moult | Implicit active contours for automatic brachytherapy seed segmentation in fluoroscopy[END_REF][START_REF] Dehghan | Ultrasound-fluoroscopy registration for prostate brachytherapy dosimetry[END_REF][START_REF] Kuo | Automatic segmentation of radiographic fiducial and seeds from X-ray images in prostate brachytherapy[END_REF][START_REF] Chng | Prostate brachytherapy postimplant dosimetry: Seed orientation and the impact of dosimetric anisotropy in stranded implants[END_REF][START_REF] Di San Filippo | Declustering n-connected components for segmentation of iodine implants in C-arm fluoroscopy images[END_REF][START_REF] Hu | Fast radioactive seed localization in intraoperative cone beam CT for low-dose-rate prostate brachytherapy[END_REF][START_REF] Di San Filippo | Intraoperative segmentation of iodine and palladium radioactive sources in C-arm images[END_REF][START_REF] Kuo | An image-guidance system for dynamic dosimetry in prostate brachytherapy using ultrasound and fluoroscopy[END_REF] have been published based on different prostate image modalities, including magnetic resonance (MR), ultrasound (US), X-ray and computed tomography (CT) images. As previously mentioned, determining seed positions is useful at different stages of the clinical protocol: during the implantation US and/or X-ray images are necessary for simply monitoring seed deposition or for intra-operative dynamic dosimetry [START_REF] Wei | Automated localization of implanted seeds in 3d trus images used for prostate brachytherapy[END_REF][START_REF] Di San Filippo | Declustering n-connected components for segmentation of iodine implants in C-arm fluoroscopy images[END_REF][START_REF] Hu | Fast radioactive seed localization in intraoperative cone beam CT for low-dose-rate prostate brachytherapy[END_REF][START_REF] Di San Filippo | Intraoperative segmentation of iodine and palladium radioactive sources in C-arm images[END_REF] and iterative correction in case of inaccurate delivery; X-ray images are also used immediately after the implantation for recording purpose; finally, CT data is most often used after one month for dose evaluation. MRI can also be used postoperatively alone or in combination with CT. These modalities have different advantages and drawbacks: US is a non radiating modality as compared to X-ray but seeds are more difficult to detect in US due to resolution, noise and reflection artifacts. Whilst MRI is non radiating, CT is more often used for dosimetric planning and evaluation since it provides useful information about tissue radiological density. X-ray based modalities enhance seed visibility while US or MRI improve the visibility of the prostate and other soft tissue. Clusters of seeds visible in the images may arise either from an inaccurate delivery of loose seeds or from an occlusion of stranded seeds in X-ray projections. These specificities have given birth to a very large collection of methods. Whilst our approach is for the detection of loose seeds in CT post-operative images we give a brief overview of some of the developed methods, with a summary in Table 1.

They are classified according to different properties:

-the three required abilities: segmentation, 3D orientation estimation and declustering, -the imaging modalities and dimensionality, -when they are used (intra-operative IO, immediate post-operative IPO or post-operative PO), -whether or not they use the planning information, -the type of seeds (stranded or loose -Palladium or Iodine) -and how they were evaluated (phantom study, clinical study). Note that Table 1 does not contain data about commercial software as there is very little information available, they effectively work as "black boxes". To the best of our knowledge no commerical software provides capabilities similar to the one described in this paper.

Recent advances in X-ray images have been reported for position detection and declustering of prostate seeds: for

Author&Ref

D"Amico [START_REF] D"amico | An image-guidance system for dynamic dosimetry in prostate brachytherapy using ultrasound and fluoroscopy[END_REF] 2000 MRI 3D IO (1) N/G (2) X Id C Su [START_REF] Su | Prostate brachytherapy seed localization by analysis of multiple projections: Identifying and addressing the seed overlap problem[END_REF] 2004 X-ray n2D (5) IO or IPO N/G X X N/G C+P Wei [START_REF] Wei | Automated localization of implanted seeds in 3d trus images used for prostate brachytherapy[END_REF] 2006 US 3D IO X L X X D P Singh [START_REF] Singh | Brachytherapy seed localization using geometric and linear programming techniques[END_REF] 2007 X-ray n2D N/G N/G X N/G C+P Fallavollita [START_REF] Fallavollita | Registration between ultrasound and fluoroscopy or CT in prostate brachytherapy[END_REF] 2010 CT 3D IO N/G X N/G C+P Kuo [START_REF] Kuo | Mri-based prostate brachytherapy seed localization[END_REF] 2010 MRI 3D N/G N/G X D P Lee [START_REF] Lee | Intraoperative 3D reconstruction of prostate brachytherapy implants with automatic pose correction[END_REF] 2011 [START_REF] Chng | Prostate brachytherapy postimplant dosimetry: Seed orientation and the impact of dosimetric anisotropy in stranded implants[END_REF] 2012 CT 3D IPO or PO S X X Id P San Silippo [START_REF] Di San Filippo | Declustering n-connected components for segmentation of iodine implants in C-arm fluoroscopy images[END_REF] 2013 X-Ray n2D IO or IPO N/G X X Id C Hu [START_REF] Hu | Fast radioactive seed localization in intraoperative cone beam CT for low-dose-rate prostate brachytherapy[END_REF] 2013 Cone beam CT 3D IO X N/G X X Id/Pd C San Filippo [START_REF] Di San Filippo | Intraoperative segmentation of iodine and palladium radioactive sources in C-arm images[END_REF] 2014

X-Ray 2D IO N/G X Pd C Moult [17] 2012 X-Ray 2D IO N/G X X Id C+P Moult[18] 2012 X-Ray 2D IO N/G X X N/G C Defghan[19] 2012 X-Ray n2D IO or IPO X N/G X Pd C+P Kuo[20] 2012 X-Ray 2D IO N/G X X Pd C+P Chng
X-Ray n2D IO or IPO X L X X Id+Pd C Kuo[25] 2014 X-Ray n2D IO or IPO X N/G X X Pd C+P Proposed method 2014 CT 3D PO L X X X Id C+P Table 1.
Existing brachytherapy seed detection systems and their innovations.

(

IO: intra-operative; IPO: immediate post-operative; PO: or post-operative.

(

N/G: information is not given. (

) 5
n2D: multiple 2D images.

example, the partition division on multiple projections of Carm fluoroscopic images [START_REF] Su | Prostate brachytherapy seed localization by analysis of multiple projections: Identifying and addressing the seed overlap problem[END_REF][START_REF] Lee | Intraoperative 3D reconstruction of prostate brachytherapy implants with automatic pose correction[END_REF], the region-based segmentation implicit active contour model [START_REF] Moult | Segmentation of iodine brachytherapy implants in fluoroscopy[END_REF][START_REF] Moult | Implicit active contours for automatic brachytherapy seed segmentation in fluoroscopy[END_REF], the geometric analysis from the graph matching problem [START_REF] Singh | Brachytherapy seed localization using geometric and linear programming techniques[END_REF] or the mathematical morphology analysis [START_REF] Kuo | Automatic segmentation of radiographic fiducial and seeds from X-ray images in prostate brachytherapy[END_REF]. However, these methods may require a sophisticated object-matching algorithm and/or calibration to deal with the substantial distortion of seeds in fluoroscopic images; moreover the orientation of seeds was not considered. The other modalities, US and MRI, have also received a lot of attention in the last decade. For example, the appropriate location estimation of seeds in a target volume of real-time MR imaging [START_REF] D"amico | An image-guidance system for dynamic dosimetry in prostate brachytherapy using ultrasound and fluoroscopy[END_REF] or in IRON images (Inversion-Recovery with On-Resonant Water Suppression) using the Laplacian of a Gaussian technique for blob detection [START_REF] Kuo | Mri-based prostate brachytherapy seed localization[END_REF] was considered. Wei et al. [START_REF] Wei | Automated localization of implanted seeds in 3d trus images used for prostate brachytherapy[END_REF] segmented the seeds from the subtraction map between the background and post-implant US images and then applied a PCA method for orientation detection. Again, these methods do not manage union-seeds and the orientation of the seeds and their detection results are limited by the poor visualization of seeds in US and MR images.

Other approaches are based on the coregistration of different image modalities. For instance, a series of methods [START_REF] Fallavollita | Registration between ultrasound and fluoroscopy or CT in prostate brachytherapy[END_REF][START_REF] Dehghan | Ultrasound-fluoroscopy registration for prostate brachytherapy dosimetry[END_REF][START_REF] Di San Filippo | Declustering n-connected components for segmentation of iodine implants in C-arm fluoroscopy images[END_REF][START_REF] Di San Filippo | Intraoperative segmentation of iodine and palladium radioactive sources in C-arm images[END_REF][START_REF] Kuo | An image-guidance system for dynamic dosimetry in prostate brachytherapy using ultrasound and fluoroscopy[END_REF] propose a volume-to-volume and point-tovolume registration scheme of US images with the implants reconstructed from fluoroscopy. Recently, Hu et al. [START_REF] Hu | Fast radioactive seed localization in intraoperative cone beam CT for low-dose-rate prostate brachytherapy[END_REF] considered the prior knowledge of the US-to-CT transformation via registration and the planning seed position to define an atlas of regions of interest for seed detection. Union-seed separation was considered, however there was no explicit mention concerning the management of seed orientations in these approaches.

Many papers consider stranded seeds whose real orientation is generally quite similar to the planned orientation, making it possible to search for lines of seeds close to the planned orientation. For instance, Chng et al. [START_REF] Chng | Prostate brachytherapy postimplant dosimetry: Seed orientation and the impact of dosimetric anisotropy in stranded implants[END_REF] estimated the seed orientations from the tangent vector to the curve of a seed strand identified in post-implant CT images at each seed position.

Because loose seeds enable the clinician to sculpt the dose to the precise treatment constraints, our objective was to develop a method allowing the use of these loose seeds. Thus, this study addresses the development of an automatic image processing solution for the segmentation, localization and orientation estimation of prostate seeds. Figure 2.d shows an illustration of the expected result of the proposed method for the detection of seed poses in the prostate from the analysis of CT scanner images.

III. METHODS

Figure 4 shows a sketch of the proposed approach. The different steps of this method are further detailed in the following sections, including:

 3D object segmentation and classification for the detection of a set of seeds (single and union-seeds) in section A.  Union-seeds separation in section B.  3D orientation estimation in section C.

A. Seed segmentation using the connected object labeling method and seed classification using the k-means method

This sub-section details the detection of a set of seed candidates and their classification into three groups (outliers, single seeds and union-seeds)see Figure 5. Brachytherapy seeds are small metallic cylinders (typically about 1mm diameter for 5mm length) and appear as high intensity objects. They may produce local artifacts obscuring neighboring tissues. Numerous methods have been developed to segment such small objects in a gray level image including some region-based methods such as watershed transformation [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulation[END_REF][START_REF] Grau | Improved Watershed Transform for Medical Image Segmentation Using Prior Information[END_REF] and level sets [START_REF] Osher | Geometric Level Set Methods in Imaging Vision and Graphics[END_REF][START_REF] Tsai | A Shape-Based Approach to the Segmentation of Medical Imagery Using Level Sets[END_REF]. However these methods require user input by positioning initial seeds or shapes. We choose to use a priori knowledge about Hounsfield values in CT images [START_REF] Robertson | and ultrasound image artifacts from prostate brachytherapy seed implants: the impact of seed size[END_REF], as well as volume information about the seeds in order to limit user interaction.

For instance, the voxel intensities of each type of material in the CT images processed in this study are [min, 400] for phantom material and for patient soft tissues, [0, 1350] for bones and [500, max] for seeds. Here, a threshold-based segmentation method, namely connected component labeling [START_REF] Shapiro | Computer Vision[END_REF] with only an intensity threshold parameter t, is considered to exploit this information for the detection of individual objects in the images. We first threshold the original volume with the threshold parameter t, then each connected component (using 26-connectivity) is assigned a label i and ordered by its size. The location of each component is determined as its center of mass c i . It is evident that the choice of the intensity threshold t is a key issue. It can be heuristically set based on the image characteristics and on the physical characteristics of the X-ray absorption of the seeds (see section IV). Moreover, the proposed method minimizes its influence on the results. Indeed, as the choice of the threshold affects the volume of the obtained connected component, no absolute volume is used, but the relative volumes of the components are compared to each other. Figure 6 shows an example of connected components detected in the CT image of a patient, for a given threshold. The detected objects in CT images of patients are divided into 4 types, including: single seeds, union-seeds, bones and noise.

Comparing the volumes of the detected components to the real volume of the seeds allows the first coarse classification of the objects detected in the images. Let us denote by V real =r 2 l the real volume of a radioactive seed, where r is the radius and l the length of seed. The next step aims at suppressing large and small objects as compared to V real . In practice, the pelvic bones are very large components (with volumes more than 100 times larger than the real volume V real ). Conversely, noise is composed of tiny components (with volumes smaller than a third of V real ). The other components are kept as candidate seeds with two types: single seeds and union-seeds. In practice up to 4 or 5 seeds can be included in a union-seed; however generally only two seeds are grouped. Determining the number of seeds in a union-seed requires an estimation of the seed volume observed in the CT exam; this volume clearly depends on the threshold t.

Without loss of generality, we can assume that single seeds outnumber union-seeds in the remaining candidates (connected components). Indeed, even if union-seeds often appear in loose seed insertions, despite care taken by the clinician, they remain exceptions and most seeds are placed with a reasonable distance between them. We thus investigated the use of an unsupervised partitioning method, k-means clustering [START_REF] Han | Data Mining: Concepts and Techniques[END_REF] to separate the candidate seeds into k groups based on volume analysis. As mentioned, most detected components correspond to the cluster of single seeds; we therefore decided to add two other clusters for smaller (if any) and larger objects (in particular including unions seeds); thus, we set k=3. Figure 7 shows an example of k-means clustering on the volume of seed candidates detected in the CT image of a real patient with 85 radioactive seeds implanted.

Figure 7: Example of the result of the k-means method applied to the volume histogram of connected components detected in CT13 image (see section III). 66 components are assigned to 3 clusters (represented by blue, green, red circle points).

The red square points correspond to the mean volume of each cluster. Here, the largest cluster with 54 members is considered as the cluster of single seeds. The detail of the general k-means method is as follows: we first randomly select k points as the initial cluster centroids. Then, each candidate is assigned to the closest centroid of k clusters based on volume information. We iterate this process until stability is reached. The k-means algorithm shows its computational simplicity in many classification applications. However, the resulting clusters strongly depend on the selection of the initial centroids. To improve the classification of seeds, the three initial centroids of clusters are defined as: the minimum, mean and maximum volume of the detected candidate components. At the end of the process, the largest cluster is selected as the set of single seeds with average volume V mean computed from the cluster. This value of V mean is considered as the observed volume of a seed in the image given by the threshold t. This is a very important element to make the method more robust with respect to the choice of t. The final classification of the components is based on the comparison of their volume to V mean . The number of seeds corresponding to each candidate is calculated (rounded to the nearest integer value) as follows:

n s = V component / V mean (1) where V component is the volume of the candidate component. At the end of this process, components are classified as follows:

-n s =0: the component is considered as an outlier and removed; -n s =1: the component is considered as a single seed; -n s >1: the component is considered as a union-seed composed of n s seeds. Declustering the union-seeds will be necessary before pose determination; it will be presented in the next section. Considering the set of voxels of each union-seed detected by the connected component labeling and the number of seeds n s grouped in this union-seed (Equation 1), this second step aims at separating union-seeds using three unsupervised learning procedures : a k-means, a modified version named kmeans-FS and a Gaussian mixture model method. Figure 8 shows the detail of this step. These three methods enable handling cases that none alone can treat properly. The kmeans clustering method is again used for the voxel locations of each union-seed to separate it. Here, k is given for each union-seed by the computed number n s (in Equation 1) and the method groups voxels based on their proximity in terms of position. However, the resulting seed clusters of k-means also depend on the selection of the initial centroids and on the size and shape differences between the regions shared in the clusters. In practice, the separation of a straight union-seed (see Figure 9.a) is easily achieved by this classical k-mean method. In contrast, the inaccuracies of this method occur when trying to separate groups of parallel seeds (see Figure 9.b) or groups of 4 or 5 seeds that are closely spaced (see Figure 9.c). Hence, we introduce two other methods to improve the separation of union-seeds. The first method consists in choosing the initial cluster of the k-means algorithm by exploiting the orientation information given by the PCA method (denoted k-means-For-Seeds); it is intended to more robustly separate groups of parallel seeds. The second method makes use of the Gaussian mixture model with the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] for better processing of complex groups of seeds.

B. Union-seed separation

1) Introducing the k-means-For-Seeds method

We define the k-means-For-Seeds (k-means-FS) method as follows: first, two main directions {v 1 ,v 2 } of the union-seeds are estimated using the PCA method (see in IIIC). Then, (n s -1) parallel planes are defined by the main direction v 1 of the union-seeds and the distance d=λ 2 /n s between them (Figure 9.b), where λ 2 is the second eigenvalue of the covariance matrix C of the union-seed. These parallel planes divide the union-seed space into n s partitions. Finally, we apply the kmeans clustering algorithm with the initial cluster centers that are the centroids of these n s partitions. In practice, this method presents its strength in solving the problem of parallel seeds. Figure 10 shows an example of the application of the basic kmeans and k-means-FS for a parallel union-seed separation. Note that the seeds displayed in Figure 10 and Figure 11 are the results with their pose information estimated using thePCA method.

2) Using the Gaussian mixture model and the EM algorithm

In some cases, the k-means-FS method does not give optimal results; thus we also consider a Gaussian mixture model (GMM) in order to improve the clustering of complicated union-seeds containing 4 or 5 seeds (Figure 9.c). Each cluster of the GMM is generated by initially choosing a cluster and then drawing from that cluster"s Gaussian distribution (with means μ and covariance C ). The probability given in a mixture of k Gaussians is given by: 𝑝 𝑥 = 𝑤 𝑖 𝑁(𝑥|𝜇 𝑖 , 𝐶 𝑖 ) 𝑘 𝑖=1

(2) where w is the prior probability (weight) of the i th Gaussian, 𝑤 𝑖 𝑘 𝑖=1

= 1 and 0 ≤ 𝑤 𝑖 ≤ 1. The parameters θ={w i ,μ i ,C i } of the i th Gaussian component are learned by the maximum likelihood estimation. The log likelihood function takes the form:

ln 𝑝 𝑋 𝑤, 𝜇, 𝐶 = ln 𝑤𝑖 𝑁 𝑥 𝑗 𝜇 𝑖 , 𝐶 𝑖 𝑘 𝑖=1 = 𝐿(𝜃|𝑋) 𝑁 𝑗 =1
(3) In this maximum likelihood problem, we try to find a set of parameters θ that maximizes L(θ|X) using the standard Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. This algorithm works as follows: we initialize an estimation of a set of Gaussian parameters θ. Here, the means μ i are obtained by the k-means algorithm, the covariance matrices C i are calculated from the distance to the nearest cluster of k-means and all Gaussian weights w i are equally likely. At each iteration of the EM algorithm, we compute the expected values of the unknown data given the observed data and the current model parameters in the expectation (E) step. The maximization (M) step involves optimizing and updating the parameters to be those with maximum likelihood. It can be shown that the log-likelihood was improved at each such iteration. This process is stopped if a local maximum has been reached or some stopping criterion is met, e.g. the default number of iterations. Figure 11 shows an example of the seed separation improvement using the GMM method compared to the k-mean based method when applied to a complex unionseed with 4 closely spaced seeds.

3) Selection of the best separation method

In order to achieve robustness, for each union-seed, we run each of the three methods: the classical k-mean method with different random selections of the initial centroid clusters, the k-means-For-Seeds method and the GMM with EM algorithm. For each of the methods, the cylindrical shapes corresponding to the model of the seeds are positioned as computed. For each solution, the sum of the number of voxels common between these cylindrical shapes and the union-seed detected in the image is calculated. The solution with the largest common volume is selected as the best one.

In our clinical experiments reported in Table 4 of section IV.C, 116 union-seeds in 14 patients were achieved. Among them the best solution was found for 30 cases (25,9%) using the k-means-For-Seeds method, 18 cases (15,5%) using the GMM method and 68 other cases (58,6%) using the classical k-means method.

The obtained single seeds are then processed in order to compute their orientation as described in section III.C.

C. Orientation estimation using the PCA method

In this third step, we aim to estimate the orientation of each 3D object detected by the connected component labeling or union-seeds declustering. Numerous approaches have been proposed to estimate the 3D object orientation in point clouds. Among the most popular, the 3D Hough transform [START_REF] Qiu | Needle segmentation using 3d hough transform in 3d trus guided prostate transperineal therapy[END_REF][START_REF] Borrmann | The 3d hough transform forplane detection in point clouds: A review and a new accumulator design[END_REF][START_REF] Zana | A Multimodal Registration Algorithm of Eye Fundus Images Using Vessels Detection and Hough Transform[END_REF] focuses on the definition of the 3-dimensional Hough Space of each point (see Figure 12.a). The computational complexity is a major drawback of the Hough transform approximated by O(s p-1 n), where n is the number of points, p is the number of parameters and s is the number of samples along one Hough

Figure 10. Example of the separation result for a parallel union-seed (detected voxels visualized in blue) with two methods: classical kmeans (yellow cylinders) and k-means-FS (cyan cylinders). In this case, the best solution is achieved using k-means-FS. Figure 11. Separation result of a union-seed with 4 closely spaced seeds using 2 methods: k-mean-FS (yellow cylinders) and GMM with the EM algorithm (cyan cylinders). The original detected union-seed voxels are visualized in red. In this case, the best solution is given by the GMM and EM algorithm.

dimension. Another category of 3D orientation estimation is based on finding minimal enclosing boxes [START_REF] Lahanas | Optimized bounding boxes for threedimensional treatment planning in brachytherapy[END_REF][START_REF] O'rourke | Finding minimal enclosing boxes[END_REF]. Approximate minimum-volume bounding box methods (see Figure 12.b) were shown to be particularly efficient with a complexity of O(n+1/ε 4.5 ) compared to the Hough transform approaches, where ε is the approximation parameter. Such minimum bounding box models are however not well-suited to our problem because they require a heuristic parameter ε for grid search of the bounding box. This parameter ε has no physical meaning related to the CT image acquisition.

In this work, we focused on a solution for 3D orientation estimation that would improve both aspects, by investigating the Principal Components Analysis (PCA) method [START_REF] Jolliffe | Principal Component Analysis[END_REF]. The PCA method is the simplest and most robust mathematical procedure for compressing and extracting the description of a set of correlated observations by rejecting low variance features. Considering p-dimensional feature vectors (in our case, 3D), the PCA method is the projection of this data onto q principal components. The first principal component v 1 is the feature space along which projections have the largest eigenvalue λ 1 of the covariance matrix C of the point cloud. This is chosen as the orientation of the object (Figure 12.c). The second principal component v 2 is the direction which maximizes the variance among all directions orthogonal to the first one. The second direction v 2 is exploited for the separation step of union-seeds.

III. EXPERIMENTS AND EVALUATION

A. Experimental setup

As mentioned in Section II, various methods have been developed for different types of seeds. In Europe, the most frequently used isotope for permanent prostate seed implantation is iodine-125 (see the edition 2011 of Radiotherapy in Practice -Brachytherapy [START_REF] Hoskin | Radiotherapy in Practice-Brachytherapy[END_REF]). That is why we focused on such type of seeds. The validation of the proposed method, described in this section, was done using CT images of brachytherapy seeds implanted into: 1) two specially created phantoms and 2) data from 14 real patients provided by the Grenoble University Hospital. 1) Radioactive iodine-125 seed.

The clinical team of the Grenoble University Hospital uses BEBIG IsoSeed®I-125 seeds. A seed is made of a cylindricalshaped ceramic material, saturated with radioactive iodine-125 compound and a gold marker located in the center, all enclosed by a laser-sealed titanium tube. The outer physical dimensions of the seed are l=4.5±0.2mm length and r=0.4±0.02mm external radius. The iodine-125 isotope emits photons at a maximum energy of 35keV and has a half-life of 59.46 days. This information is provided by the manufacturer [START_REF]IsoSeedI-125 (I25.S06) Instructions for Use[END_REF]. Note that, the iodine-125 was taken off for the case of phantom to avoid the risk of radioactive contamination.

2) CT Images

The 3D CT images were obtained using a GE Lightspeed RT16 scanner with the default X-ray tube parameters:120kVp, 380-440mA.s. The slice thickness was 0.625mm with 16 frames/sec for each slice. The image reconstruction matrices were 512x512 archived in DICOM 3.0 format with 16-bit gray-level intensities. These acquisition parameters were experimentally determined by the radiophysicists using phantoms so that the seeds could be seen on 3 to 5 slices [START_REF] Meneu | Consideration of seeds orientation in prostate brachytherapy and dosimetric analysis[END_REF].

3) Evaluation

In this paper, the evaluation of the proposed approach is given in three terms: location and orientation detection performance and calculation time. A fully automatic software based on the proposed solution was built on the open-source framework CamiTK [START_REF] Fouard | CamiTK: a modular framework integrating visualization, image processing and biomechanical modeling[END_REF] (using C++ with VTK and ITK libraries). This software also provides some post-processing tools to verify seed by seed in case of errors or inaccuracy. The user can modify the seed position to better fit the image data; this manually edited seed location is the reference information to which the automatic detection is compared.

For this, we used the Euclidean distance between their centroids c and c ref and the dot product of their orientation vectors v and v ref (Equation 4).

where ||v i || is the magnitude of the vector v i , the unit of distance Δd is expressed in mm, and Δθ is in degrees. Here, the time evaluation is on a computer of 3.4GHz Intel Core i7-2600 CPU.

B. Phantoms

The radiophysicists of the clinical team created two phantoms (see Figure 13) wherein the seeds were precisely positioned with different orientations on the surface of the slab. Both phantoms had the same physical dimensions of 9x9x0.5cm3. The reference position of the seeds in the phantom based on manual detection was also provided for each case.

1) Single seed in a slab phantom

The first phantom (cf. Figure 13.a) that we considered had a single seed located in certain "pre-defined" orientations to evaluate the pose detection of the proposed method. The 3D image size of this phantom is 512x512x41 and the voxel size is 0.199x0.199x0.625. We created 11 orientation cases for this phantom. Figure 14 illustrates the different orientations of the seeds positioned in this phantom. Table 2. The details of the orientation differences Δθ and distances Δd between the detected pose of the seeds and their reference value for t=1500 (in degrees and millimeters). Table 2 shows the values of the reference orientations and the comparison between the results of the proposed method and the reference with respect to orientation Δθ and distance Δd when using an intensity threshold of t=1500. With this value, the greatest Δθ orientation error of 1.69 o is observed for the case (i) (for about 70°) and the other cases are equal to or smaller than 0.83 o . The best detection was obtained with the reference orientation close to 0° (seed perpendicular to the CT acquisition plane) where Δθ=0.2° and Δd =0.02mm.

Additionally, to investigate the dependency of the proposed method to the choice of t, we also ran this experiment with 10 choices of t in an interval of [800,1700]. The mean value and the standard deviation of the orientation error are Δθ=0.96°±0.4. The mean value and standard deviation of the distance error are Δd=0.08±0.04 mm).

2) Multiple seeds in a multi-slab phantom

We also report a second experiment about the detection of 73 seeds implanted in a more complex phantom composed of 9 slabs (Plastic water®-LR, Medi-Test, Saclay, France). The placements of the seeds look more similar to what could be a real implantation for a patient. Nine same size slabs are pressed together by 4 screws. Holes were drilled into the slabs to place the radioactive seeds in different orientations. Figure 13.b shows the design of the central slab of the phantom. The 3D image size of this phantom is 512x512x73 and the voxel size is 0.217x0.217x0.625 mm 3 . In this experiment, we have tested the proposed method with 10 choices of t in the interval of [800,1700]. Figure 15 shows the distribution of the 73 seeds detected in a CT image when using the intensity threshold t=1500.

Compared to the reference data, the mean and standard deviation of the orientation error are Δθ=1.32°±0.9; the mean and standard deviation of the position error is Δd=0.13±0.07mm for 10 choices of the intensity threshold t. Some seeds are perfectly located (Δθ=0°, Δd=0mm), e.g., the seeds positioned vertically in the central slab when using an intensity threshold of t=1500. In contrast, the most inaccurate detections were found for the oblique seeds, where Δθ=1.8° and Δd=0.65mm.

C. Patient experiment

We report 14 cases of radioactive seed detection on real patient images. The CT scanner images were taken one month after the implantation procedure in the Grenoble University Hospital. Data were anonymized before export and processing. Table 4 details the information of the scanner images and the number of seeds implanted in the planning of the brachytherapy treatment. It should be reminded that, some seeds may migrate, therefore the number of detected seeds may be different even if perfectly successful. In this experiment, we ran the proposed method with different values of the intensity threshold t for all 14 patient cases. Again, the choice of the intensity threshold t was based on the HU values of each material in the CT images: in these data, soft tissues were in the range [min,400], bones were in [0,1350] and seeds were in [500,max]. We considered the intensity threshold in the interval of [1350,2000] with a step of 100 (beginning with 1400).

Seed detection errors, which we will call "false detections" (FD), can be classified into three categories as follows:

* FD1: the undetected seeds because of their migration out of the imaged region or because of the limitations of the method (false negative). * FD2: the noisy objects detected in image such as calcifications that can be recognized as potential seeds if their intensity and volume are compatible with the seeds or unionseeds characteristics (false positive).

* FD3: correspond to the errors due to the incorrect separation of union seeds; that may be due to a wrong estimated number of seeds.

The sum of all types of false detections, when using the proposed method with different values of the intensity threshold t in the range [1350,2000], is plotted in Figure 16.

Here, the best results were obtained with the thresholds t=1500 or 1600 where the sum of false detections (FD) are minimum for all 14 patient cases. When a low threshold t (e.g., t=1350) is chosen FD2 increases with extra noise. In contrast, the number of wrongly separated seeds (FD3) is bigger when a high threshold t (e.g. t=1900) is used; this is because of the loss of shape information of the seeds. The total number of false detections, summed from the 14 cases, for each choice of threshold t is given in Table 3.

From this experimental result, we suggest choosing the intensity threshold in Hounsfield Units in [1500,1700], for which the radioactive seeds can be most successfully separated from the other material.

Table 4 details the results obtained for the 14 prostate cases using the proposed method with an intensity threshold of t=1500. Compared to the reference data, the maximum (respectively minimum) orientation error Δθ is 3.180°±0.9 (resp. 0.680°±0.2) and the maximum (resp. minimum) distance error Δd is 0.50±0.16mm (resp. 0.15±0.09mm). The number of false detections is also reported in detail for each case. Compared to the number of implanted seeds, 8 seeds in total could not be found (FD1) in CT3, CT5, CT7, CT8 and CT9 images. Some existing objects of FD2 type were also detected in CT3 (cf. Figure 17), CT5, CT7, CT10 and CT11. Considering FD3 false detections (wrong union-seed separation), two groups of 4 real implanted seeds were detected as union-seed of 5 seeds; this occurs in CT3 and CT5 images.

Figure 18 also shows the detail of each type of false detection for a single patient (CT14) and for different values of the threshold t. The mean calculation time of the proposed method was 9.7s over 126 runs for 14 patient cases with different intensity thresholds t. It is therefore a very fast solution when compared to the half day that was required for the careful manual segmentation of the reference data or compared to the average 30 minutes for conventional postprocessing of images treated by our clinical team after detection of seeds with their commercial system. Figure 19 shows some examples of the seeds detected in patients. 2 In this example a smaller ROI would have allowed to avoid such a false detection ; however, we kept it because the visualization in this figure was easier than when false detections are in the middle of the other seeds.

IV. DISCUSSION

In this study, an efficient and fast approach for the pose estimation of brachytherapy seeds in CT images has been presented. The key methods used in this work may be summarized as follows:

 Classification of detected objects based on a priori intensity and volume information  Estimation of 3D objects based on the extraction of principal components.  Separation of groups of seeds using k-means, a modified k-means clustering method and GMM method with an EM algorithm.  Individual seed orientation estimation using the PCA method. The results herein were quickly obtained with only small differences compared to the reference data for both phantoms and patients. For example, the proposed method achieves the orientation error Δθ=0.96°±0.4, the distance error Δd=0.08±0.04 mm and the calculation time 0.67s±0.3 for a phantom in section IV.B.1. These evaluation values can be compared to Δθ=2°, Δd=0.3mm and calculation time 9s of a previous related work for CT images (Chng's method [START_REF] Chng | Prostate brachytherapy postimplant dosimetry: Seed orientation and the impact of dosimetric anisotropy in stranded implants[END_REF]). This comparison in terms of detection performance and calculation time pointed out the relevance of our contributions. However these comparisons must be interpreted carefully since the phantoms and computer systems were different.

In clinical practice, the role of the human operator can be limited to a verification task and modifications for misdetections if any. In fact, the migrated seeds (false detection type: FD1) always lead to a significant uncertainty in the post-implant dosimetry calculation. The proposed method described in this paper works without the prior knowledge about the number of implanted seeds. Therefore finding less seeds than expected by the planning should warn the operator about possible migration or misdetections. On the other hand, the false positive detections (FD2) are often due to calcifications very frequent in the prostate gland of men with benign prostatic hypertrophy and prostate cancer. The distinction between calcifications and brachytherapy seeds may be still a challenging task even for a very skilled clinician. Thus, the verification by the clinician in this case is always necessary but may not be sufficient.

The advantage of the clustering approaches (k-means or GMM) used in this work for the separation of union-seeds is their small computational complexity compared to other approaches such as RANSAC [START_REF] Moradi | Needle path detection for brachytherapy dosimetry based on lateral power imaging and template matching[END_REF] or Hough transform [START_REF] Zana | A Multimodal Registration Algorithm of Eye Fundus Images Using Vessels Detection and Hough Transform[END_REF]. However, there was still a few false detections due to inappropriate seed separation (here, FD3): less than 0.19% (see Table 4). Therefore, a manual correction step may also be necessary. Some approaches based on the morphology analysis or 3D template matching could be considered in the future to improve the performance of this separation task. Overall of 1063 seeds implanted in 14 patient cases in Table 4, the false positive (FD2+FD3) percentage of our approach is 1.03%. This value is lower compared to the state-of-the-art methods reported in the Table 2 of the publication of San Filippo et al. [START_REF] Di San Filippo | Intraoperative segmentation of iodine and palladium radioactive sources in C-arm images[END_REF], where Moult et al. [START_REF] Moult | Segmentation of iodine brachytherapy implants in fluoroscopy[END_REF] was 2.2%, San Filippo et al. [START_REF] Di San Filippo | Declustering n-connected components for segmentation of iodine implants in C-arm fluoroscopy images[END_REF][START_REF] Di San Filippo | Intraoperative segmentation of iodine and palladium radioactive sources in C-arm images[END_REF] were 1.7%. These results open the door to accurate dose calculation and procedure quality assessment. As already mentioned in the introduction, a dosimetry study has been launched in parallel in our institution on a series of patients including the 14 patients of this paper [START_REF] Meneu | Consideration of seeds orientation in prostate brachytherapy and dosimetric analysis[END_REF]. Its aim is to evaluate the impact of accurate pose evaluation onto dose distribution for prostate brachytherapy treatment. It will be published separately but a few elements may be summarized As concerns the Dose Volume Histograms no significant difference could be demonstrated between dose evaluation using or not orientation information but our number of patients was quite small. In a very recently published paper [START_REF] Fekete | Quantifying the effect of seed orientation in postplanning dosimetry of low-dose-rate prostate brachytherapy[END_REF] concerning a study using 5 fluoroscopic images and a CT of 287 patients, the authors demonstrate small but significant dose difference evaluated on organs at risks. In our much more limited study we also exhibited significant local dose differences that could have a clinical impact. A more extensive clinical study is necessary to draw useful conclusions. In the context of this study which aim was to evaluate a new method, man-machine interaction has been restrained to what was strictly necessary. It is clear that the routine use of the approach would require special care for assisting the operator in the verification task, in particular we envision to orient his/her screening towards suspicious detections. Work remains to be done for scoring the detection.

In future work, the potential of an automatic choice of an optimal intensity threshold will be further explored from the analysis of quantitative and geometric information of false detections. Robust solutions of seed localization for different prostate image modalities such as ultrasound or image registration for improved evaluation with respect to the anatomy of patient are also among the key issues that should be addressed.

V. CONCLUSION

This paper has presented an automatic, accurate, robust and fast approach for the automatic localization of brachytherapy loose Iodine-125 seeds in CT post-operative images. It was evaluated on phantom and patient data. The key originalities of this work lie in the ability to separate groups of seeds and to determine their orientation for improved evaluation of dose distribution to the patient. Based on intensity information and observed volume of the detected objects the method was able to accurate determine the 5D pose of seeds with very few false detections. Further work both concerns technical aspects and a larger clinical evaluation.

Figure 1 .

 1 Figure 1. Prostate brachytherapy implant technique (source http://www.prostatespecialist.co.uk)

Figure 2 .

 2 Figure 2. Seed distribution in prostate brachytherapy: (a) Seeds aligned with the insertion direction as defined in the planning (from http://cancer.uc.edu/cancerinfo/TypesOfCancer/ProstateCancer/InterstitialBrachtherapy.aspx); (b) Real distribution of seeds one month after the implantation; (c) and (c1) Seed distribution as handled by existing commercial software; (d) Seed distribution detected by the proposed method.

Figure 3 .

 3 Figure 3. Example of implanted seeds as visible in a single slice of a CT image.

  seed; Pd : Palladium 103; Id: Iodine125. (4) P: Phantom; C: Clinical data.

Figure 4 .

 4 Figure 4. The main step of our proposed approach.

Figure 5 .

 5 Figure 5. The detail of seed segmentation and classification step.

Figure 6 .

 6 Figure 6. Example of connected component detection in a 3D CT image with 190 objects obtained using an intensity threshold of 700.

Figure 8 .

 8 Figure 8. The detail of union-seed separation step.

Figure 9 .

 9 Figure 9. Seed separation. A case where the standard k-means algorithm works well (a) and two more difficult cases: (b) k-means-For-Seeds method with the choice of the initial centroid to improve the clustering performance of k-means, where the red point are the centroids of the k partitions, x is the centroid of the union-seed, the dotted-lines are used to determine the k-1 parallel planes w.r.t. the distance d= λ 2 /k. (c) Separation result obtained using the GMM with the EM algorithm.

Figure 12 :

 12 Figure 12: Different approaches for 3D orientation estimation of a 3D cloud of points: (a) Hough transform (image is taken from [33]). (b) Minimal bounding box. (c) Principal Components Analysis.

Figure 13 .

 13 Figure 13. Water equivalent phantoms: (a) a single seed in a slab phantom (cf. section III.B.1); (b) seeds located on the surface of the central slab of the other phantom (cf. section III.B.2)

Figure 14 .

 14 Figure 14. Detection of a seed with different orientations. The overlaid curves are the detected contours for different intensity thresholds.

Figure 15 .

 15 Figure 15. 3D CT scanner image of the second phantom and the distribution of the 73 seeds: theoretical positions of seeds are painted in green and detected seeds are in red.

Figure 16 .

 16 Figure 16. Number of false detections obtained for the different choices of the intensity threshold t for the 14 cases.

Figure 17 .

 17 Figure 17. Example of FD2 false detection in CT3 with t=1500 (cyan cylinder at the top of the image)2.

Figure 18 .Figure 19 .

 1819 Figure 18. Number of false detections for different values of t for a patient example (CT14 case).

Year Imaging Modality Image Dimensionality When? Using Planning Stranded \ Loose Position Detection Orientation Estimation Declustering Seed Types

Table 3 .

 3 Cumulated

number of false detections for the 14 cases as a function of threshold value t. Here, the total number of implanted seeds is 1063 for 14 cases.

ACKNOWLEDGMENTS

This work is supported by grants from INSERM (Dorgipro project PC201114) and from ANR (CAMI labex ANR-11-LABX-0004).

The authors would like to thank Francois Meneu, Jean-Yves Giraud, Dr. Camille Verry and Pr. Michel Bolla from the Radiotherapy department, Grenoble University Hospital for providing their expertise and the clinical images used in this study. The authors also thank PhD. Cecilia Hughes and PhD. Nikolai Hungr for critically editing the manuscript. Table 4. Detailed results for 14 CT scanner images of patients. The intensity threshold is t=1500. The minimum intensity value (in HU) over the 14 cases is -3024. The 3 types of false detections are mentioned. The detected seeds are classified as: single, group of two, three, four or five seeds. The false detection percentage is in the last line.


Sébastien Bougleux
email: sebastien.bougleux@greyc.ensicaen.fr

Abderrahim Elmoataz
email: abder.elmoataz@greyc.ensicaen.fr

Mahmoud Melkemi
email: mahmoud.melkemi@uha.fr

Discrete Regularization on Weighted Graphs for Image and Mesh Filtering

We propose a discrete regularization framework on weighted graphs of arbitrary topology, which unifies image and mesh filtering. The approach considers the problem as a variational one, which consists in minimizing a weighted sum of two energy terms: a regularization one that uses the discrete p-Laplace operator, and an approximation one. This formulation leads to a family of simple nonlinear filters, parameterized by the degree p of smoothness and by the graph weight function. Some of these filters provide a graph-based version of well-known filters used in image and mesh processing, such as the bilateral filter, the TV digital filter or the nonlocal mean filter.

Introduction

In many computer vision applications, it is necessary to filter and to simplify images or meshes. In the context of image processing, smoothing and denoising constitute important steps of filtering processes. Among the existing methods, the variational ones, based on regularization, provide a general framework to design efficient filters. Solutions of variational models can be obtained by minimizing appropriate energy functions. The minimization is usually performed by designing continuous partial differential equations (PDEs), whose solutions are discretized in order to fit with the image domain. A complete overview of these methods can be found in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] [START_REF] Chan | Image Processing and Analysis -variational, PDE, wavelets, and stochastic methods[END_REF][3] [START_REF] Alvarez | Axioms and fundamental equations of image processing[END_REF] and references therein.

Another important problem of computer vision is mesh smoothing or denoising. This process is dedicated to noise removal, causing minimal damage to geometric features. Most of mesh smoothing methods are based on the discrete Laplace-Beltrami regularization or on the discrete curvature regularization [START_REF] Taubin | Geometric signal processing on polygonal meshes[END_REF] [START_REF] Ohtake | Mesh regularization and adaptive smoothing[END_REF].

Variational Beltrami flows have also been used to denoise and regularize data defined on manifolds [START_REF] Sochen | Variational Beltrami flows over manifolds[END_REF]. Other mesh smoothing methods, based on feature preserving, were mostly inspired by anisotropic diffusion in image processing [START_REF] Desbrun | Anisotropic feature-preserving denoising of height fields and bivariate data[END_REF][9] [START_REF] Hildebrandt | Anisotropic filtering of non-linear surface features[END_REF].

Inspired by continuous regularization of images, we propose a general framework based on a discrete regularization on weighted graphs of arbitrary topology. This framework unifies the regularization of images and meshes. Let G w = (V, E) be a weighted graph which consists of a set of vertices V , a set of edges E ⊂ V ×V , and a similarity weight function w defined on edges. Let H(V ) be a Hilbert space defined on the vertices of G w . We formalize the discrete regularization of a function f 0 ∈ H(V ) by the following minimization problem:

min f ∈H(V ) E p (f, f 0 , λ) = v∈V ▽ v f p + λ f -f 0 2 H(V ) , (1)
where p ∈ [1, +∞) is the smoothness degree, λ is the fidelity parameter, and ▽f represents the weighted gradient of the function f over the graph. The solution of problem (1) leads to a family of nonlinear filters, parameterized by the weight function, the degree of smoothness and the fidelity parameter. There exist two main advantages of using this framework, which can be considered as a discrete analogue of continuous regularization on weighted graphs. Firstly, the regularization is expressed directly in discrete settings. Secondly, filters are computed by simple and efficient iterative algorithms, without solving any PDEs. Second, the topology of graphs can be arbitrary. Since the proposed approach is general, any discrete data set can be transformed into a weighted graph, by using a similarity measure between data. Thus, we can consider any function defined on these data as a function defined on the vertices of the graph.

The family of filters we propose includes graph-based versions of well-known filters used in image and mesh processing. If w = 1, they correspond exactly to the digitized PDE filters, introduced in the context of image restoration on grid graphs [START_REF] Osher | Digitized PDE method for data restoration[END_REF] (L 2 digital filter for p = 2 and TV digital filter for p = 1). If w = 1, our the filters behave like weighted TV regularization or weighted L 2 regularization. In particular, if p = 2 and λ = 0, the choice of the weight function w allows to find the bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] and the nonlocal mean filters [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. In this particular case, we also show that the discrete regularization is linked to spectral graph theory [START_REF] Chung | Spectral graph theory[END_REF] and to Markov matrix filtering [START_REF] Coifman | Geometries of sensor outputs, inference, and information processing[END_REF]. We can quote other existing methods, developed in the context of image filtering, that can be considered as discrete regularizations on unweighted graphs [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF] [START_REF] Darbon | Exact optimization of discrete constrained total variation minimization problems[END_REF]. These regularizations yield to Markov random fields where only binary variables are involved in the minimization.

The solution of problem ( 1) is obtained by defining local differential geometry operators on weighted graphs. The idea of using differential geometry on graphs, in a regularization process, has also been proposed in other context, such as semisupervised data learning [START_REF] Zhou | Regularization on discrete spaces[END_REF] and image segmentation [START_REF] Bougleux | Image smoothing and segmentation by graph regularization[END_REF].

In this paper, we firstly define differential geometry operators on weighted graphs in Section 2. Section 3 presents the regularization problem (1) and the associated family of filters. Section 4 analyzes the obtained filters and gives relations to existing methods. In Section 5, we show some regularization examples, in the context of image denoising, image segmentation and mesh denoising.

Differential Geometry on Weighted Graphs

In this section, we recall some basic definitions on graphs, and we define local differential geometry operators which can be considered as discrete versions of continuous differential operators. Analogue definitions and properties have also been used in the context of semi-supervised learning [START_REF] Zhou | Regularization on discrete spaces[END_REF], and differential calculus on graphs [START_REF] Bensoussan | Difference equations on weighted graphs[END_REF] [START_REF] Friedman | Wave equations for graphs and the edge-based laplacian[END_REF].

Preliminary Definitions

A graph G w = (V, E) consists of a finite set V of N vertices and a finite set E ⊆ V × V of edges. We assume G w to be undirected, with no self-loops and no multiple edges. Let uv be the edge that connects the vertices u and v. An undirected graph is weighted if it is associated with a weight function w : E → R + satisfying w uv = w vu , for all uv ∈ E, and w uv = 0 if uv ∈ E. The weight function represents a similarity measure between two vertices of the graph.

Let H(V ) denotes the Hilbert space of real-valued functions on vertices. A function

f : V → R in H(V ) assigns a vector f v to each vertex v ∈ V . Clearly, f can be represented by a column vector of R N , f = [f 1 , . . . , f N ] T .
The function space H(V ) is endowed with the usual inner product f, h H(V ) := v∈V f v h v , where f, h ∈ H(V ). Similarly, one can define H(E), the space of real-valued functions on edges, endowed with the inner product f, h H(E) := v∈V u∼v f uv h uv , where f and h are two functions in H(E).

Weighted Gradient and Divergence Operators

Let G w = (V, E) denotes a weighted graph. The difference operator d :

H(V ) → H(E) of a function f ∈ H(V ) on an edge uv ∈ E, is defined by: (df ) uv := √ w uv (f v -f u ), ∀uv ∈ E. (2)
The directional derivative (or edge derivative) of a function f ∈ H(V ) at a vertex v along an edge e = uv, is defined as ∂ v f u := (df ) uv . This definition is consistent with the continuous definition of the derivative of a function, e.g., if

f v = f u then ∂ v f u = 0. Moreover, note that ∂ v f u = -∂ u f v , and ∂ v f v = 0. The weighted gradient operator ▽ of a function f ∈ H(V ) at a vertex v is the vector operator defined by ▽ v f = (∂ v f u : u ∼ v) T .
The local variation of f at v, is defined to be:

▽ v f := u∼v (∂ v f u ) 2 = u∼v w uv (f v -f u ) 2 . ( 3
)
It can be viewed as a measure of the regularity of a function around a vertex. The amplitude of the graph gradient is defined by ▽ f := ▽f, ▽f

1/2
H(E) . The adjoint operator of the difference operator, denoted by

d * : H(E) → H(V ), is defined by df, h H(E) := f, d * h H(V ) , with f ∈ H(V ) and h ∈ H(E).
Using the definitions of the inner products in H(V ) and H(E) (see Section 2.1), and definition (2), we obtain the expression of d * at a vertex v:

(d * h) v = u∼v √ w uv (h uv -h vu ). ( 4
)
The divergence operator, defined by -d * , measures the net outflow of a function in H(E) at each vertex of V .

A Family of Weighted p-Laplace Operators

The weighted p-Laplace operator, ∆ p :

H(V ) → H(V ) with 1 ≤ p < +∞, is defined by ∆ p f := d * ( ▽ f p-2 df
). Substituting ( 2) and ( 4) into the definition of ∆ p f , we obtain the expression of ∆ p at a vertex v:

(∆ p f ) v = u∼v γ uv (f v -f u ) , (5)
where γ is the function defined by

γ uv := w uv ▽ v f p-2 + ▽ u f p-2 .
The operator ∆ p is nonlinear, with the exception of p = 2. Furthermore, ∆ p is positive semi-definite:

f, ∆ p f H(V ) = f, d * ( ▽ f p-2 df ) H(V ) = df, ▽ f p-2 df H(E) = v∈V ▽ v f p-2 u∼v (df ) 2 uv = v∈V ▽ v f p ≥ 0. ( 6
)
The definition of ∆ p can be considered as the discrete analogue of the p-Laplacian in the continuous case. When p = 2, ∆ 2 represents the weighted Laplace operator on G w , and (5) reduces to:

(∆f ) v := (∆ 2 f ) v = 2 u∼v w uv (f v -f u ). (7)
When p = 1, ∆ 1 represents the weighted curvature operator on G w , and expression (5) reduces to:

(κf ) v := (∆ 1 f ) v = u∼v w uv 1 ▽ v f + 1 ▽ u f (f v -f u ). ( 8
)
In practice, to avoid zero denominator in [START_REF] Desbrun | Anisotropic feature-preserving denoising of height fields and bivariate data[END_REF], the local variation (3) is replaced by its regularized version:

▽ v f ǫ := ▽ v f 2 + ǫ 2 , with ǫ → 0 fixed.

p-Laplace Regularization on Weighted Graphs

In this section, we present the discrete regularization problem (1) and associated filters. Let G w = (V, E) be a weighted graph. The regularization of a given function f 0 ∈ H(V ), using the weighted p-Laplace operator, consists in seeking for a function f * ∈ H(V ) which is not only smooth enough on G w , but also close enough to f 0 . This optimization problem can be formalized by the minimization of a weighted sum of two energy terms:

f * = min f ∈H(V ) E p (f, f 0 , λ) := v∈V ▽ v f p + λ f -f 0 2 H(V ) . ( 9
)
The first energy in ( 9) is the smoothness term or regularizer, meanwhile the second is the fitting term. The parameter λ ≥ 0 is a fidelity parameter, called the Lagrange multiplier, which specifies the trade-off between the two competing terms. Both energy functions in E p are strictly convex functions of f . In particular, by standard arguments in convex analysis, the problem ( 9) has a unique solution, for p = 1 and p = 2, which satisfies:

∂E p (f, f 0 , λ) ∂f v = ∂ ∂f ▽ v f p + 2λ(f v -f 0 v ) = 0, ∀v ∈ V . ( 10
)
Using equation ( 6) to compute the derivative of the first term in E p , the system of equations ( 10) is rewritten as:

(∆ p f ) v + 2λ(f v -f 0 v ) = 0, ∀v ∈ V . ( 11
)
The solution of problem ( 9) is also the solution of the system of equations [START_REF] Osher | Digitized PDE method for data restoration[END_REF]. This is a nonlinear system, with the exception of p = 2 (see Section 2.3). Substituting the expression of the p-Laplace operator into [START_REF] Osher | Digitized PDE method for data restoration[END_REF], we obtain:

2λ + u∼v γ uv f v - u∼v γ uv f u = 2λf 0 v , ∀v ∈ V . ( 12
)
Among the exiting methods to solve the system of equations [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF], we use the Gauss-Jacobi iterative algorithm. Let t be an iteration step, and let f (t) be the solution of equation ( 12) at the step t. The corresponding linearized Gauss-Jacobi algorithm is given by:

     f (0) = f 0 γ (t) uv = w uv ▽ v f (t) p-2 + ▽ u f (t) p-2 , ∀uv ∈ E f (t+1) v = 2λ + u∼v γ (t) uv -1 2λf 0 v + u∼v γ (t) uv f (t) u , ∀v ∈ V (13)
where γ (t) is the function γ at the step t. The weights w uv are computed from f 0 , or can be given a priori. We define the function ϕ at an iteration t of algorithm (13) by:

ϕ (t) vu = γ (t) uv 2λ + u∼v γ (t) uv if u = v, and ϕ (t) vv = 2λ 2λ + u∼v γ (t) uv
Then, an iteration of the regularization algorithm ( 13) is rewritten as:

f (t+1) v = ϕ (t) vv f 0 v + u∼v ϕ (t) vu f (t) u , ∀v ∈ V . ( 14
)
At each iteration, the new value f (t+1) , at a vertex v, depends on two quantities, the original value f 0 v , and a weighted average of the existing values in a neighborhood of v. This shows that the proposed filter, obtained by iterating [START_REF] Coifman | Geometries of sensor outputs, inference, and information processing[END_REF], is a low-pass filter which can be adapted to many graph structures and weight functions.

Filter Analysis and Examples

The case of an arbitrary p is not considered in this paper. In the sequel, we discuss the cases of p = 2 and p = 1, and we give some existing results related to our filter. When p = 2, it follows from equation ( 11) that the solution of problem ( 9) is based on the weighted Laplace operator defined by equation [START_REF] Sochen | Variational Beltrami flows over manifolds[END_REF]. Equation [START_REF] Osher | Digitized PDE method for data restoration[END_REF] reduces to ∆f * + 2λ(f * -f 0 ) = 0. In this case, the iterative filter ( 13) is linear on the graph structure, and the coefficients given by the function γ do not have to be updated at each iteration because they depend on the function w. When p = 1, it follows from equation ( 11) that the solution of problem ( 9) is based on the weighted curvature operator defined by equation [START_REF] Desbrun | Anisotropic feature-preserving denoising of height fields and bivariate data[END_REF]. Equation ( 11) reduces to κf * + 2λ(f * -f 0 ) = 0. In this case, the iterative filter ( 13) is nonlinear, and the coefficients given by the function γ are adaptively updated at each iteration in addition of updating the function f .

Regularization of Functions on Discrete Data

The family of filters presented in Section 3 can be used to regularize any function defined on discrete data by constructing a weighted graph, and by considering the function to be regularized as a function defined on graph vertices.

Let V = {v 1 , . . . , v N } be a finite data set such that v i is a vector of R m . There exist several popular methods that transform the set V with a given pairwise similarity measure w into a graph G w = (V, E). Constructing similarity graphs consists in modelizing local neighborhood relationships between data points. Among the existing methods, we can quote the ǫ-neighborhood graph where two points u, v ∈ V are connected by an edge if v-u < ǫ, ǫ > 0. Another important graph is the k-nearest neighbors graph where two points u, v ∈ V are connected by an edge if u is among the k nearest neihgbors of v. This definition leads to a directed graph because the neighborhood relationship is not symmetric. In order to make this graph symmetric, let nn k (u) be the set of k-nearest neighbors of the point u. Then, a point v is connected to

u if u ∈ nn k (v) or v ∈ nn k (u).
Let f ∈ H(V ) be a function defined on each point of the data set V . Similarities between data points are estimated by comparing their features. Features generally depend on the function f and the set V . Every point v ∈ V is assigned with a feature vector denoted by F f (v) ∈ R q . We propose to use one of the two general weight functions:

w uv = exp - F f (v) -F f (u) 2 h 2 , ∀uv ∈ E, (15)
and,

w uv = exp - u -v 2 2σ 2 exp - F f (v) -F f (u) 2 h 2 , ∀uv ∈ E. ( 16
)
where σ and h are two parameters depending on the variations of u -v and

F f (v) -F f (u) over the graph.
The graph structure, associated with one of the above weight functions, describe a general family of filters. This family is linked to several filters defined in the context of image and mesh processing. When p = 2, filter (13), associated with the weight function [START_REF] Darbon | Exact optimization of discrete constrained total variation minimization problems[END_REF], is equivalent to the bilateral filter, introduced in the context of image denoising [12] [21]. It is a nonlinear filter that combines geometric and range filtering. Bilateral filtering is also used to denoise meshes [START_REF] Fleishman | Bilateral mesh denoising[END_REF]. It is obtained by using the scalar feature F f (v) = f v for all v ∈ V . Using the same parameters, filter [START_REF] Chung | Spectral graph theory[END_REF] can also be considered as a discrete nonlocal mean filter, introduced in the context of images [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. Indeed, it is obtained by using the weight function [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF] with the feature vector

F f (v) = [f u : u ∈ B v,s ] T and B v,s a bounding box of size s centered at v.
When λ = 0 and w is constant, filter ( 13) corresponds exactly to the digitized PDE filters proposed in the context of image restoration [START_REF] Osher | Digitized PDE method for data restoration[END_REF]. If p = 1, it is the TV regularization. If p = 2, it is the L 2 regularization. In general, if the weight function is not constant, filter [START_REF] Chung | Spectral graph theory[END_REF] corresponds to the weighted L 2 regularization and the weighted TV regularization on arbitrary graphs.

Relationships with Spectral Graph Filtering

We consider the regularization problem (9) for p = 2 and λ = 0, and we show that it can be expressed in terms of spectral graph theory [START_REF] Chung | Spectral graph theory[END_REF]. From expression [START_REF] Coifman | Geometries of sensor outputs, inference, and information processing[END_REF], the filter reduces to:

f (t+1) v = u∼v ϕ vu f (t) u , ∀v ∈ V , (17)
where ϕ vu = w uv / u∼v w uv , ∀uv ∈ E. As we have ϕ vu ≥ 0 and u∼v ϕ vu = 1, ϕ vu can be interpreted as the probability of a random walker to jump from v to u in a single step [START_REF] Coifman | Geometries of sensor outputs, inference, and information processing[END_REF]. Let P be the N × N Markov matrix defined by: P (v, u) = ϕ vu if the edge uv ∈ E, and P (v, u) = 0 otherwise. Let F be the matrix form of the function f . With these notations, expression ( 17) is rewritten as:

F (t+1) = P F (t) = P t F (0) . ( 18
)
An element P t (v, u), vu ∈ E, describes the probability of transition in t steps. The matrix P t encodes local similarities between vertices of the graph and it diffuses or integrates this local information for t steps to larger and larger neighborhoods of each vertex.

The spectral decomposition of the matrix P is given by P φ i = λ i φ i , with 1 ≥ λ 1 ≥ . . . ≥ λ i ≥ . . . ≥ λ N ≥ 0 the eigenvalues of P , and φ i its eigenvectors. The eigenvectors associated with the k first eigenvalues contain the principal information. Thus, an equivalent way to look at the power of P in filter ( 18) is to decompose each value of F on the first eigenvectors of P . Moreover, the eigenvectors of the matrix P can be seen as an extension of the Fourier transform basis functions with eigenvalues representing frequencies [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data[END_REF]. It defines a basis of any function f in H(V ), and the function f can be decomposed on the k first eigenvectors of P as: f = i=k i=1 f, φ i H(V ) φ i . This can be interpreted as a filtering process in the spectral domain.

Applications

The family of filters proposed in Section 3 can be used to regularize any function defined on the vertices of a graph, or on any discrete data set. Through examples, we show its efficiency in the case of image denoising, image simplification, polygonal curve denoising and surface mesh denoising. We also compare several results obtained for different weight functions and regularization parameters. Image Smoothing and Denoising : Let f 0 be a gray level image of N pixels, 1 illustrates the regularization of a noisy image f 0 (Fig. 1(b)) on a grid graph of 8-adjacency (Fig. 1(c),(d) and (e)), and on a ǫ-neighborhood graph (Fig. 1(f)). The weight function is chosen such that Fig. 1(c) corresponds to the unweighted L 2 regularization, Fig. 1(d) to the unweighted TV regularization, Fig. 1(e) to the bilateral filter, and Fig. 1(f) to the nonlocal mean filter. The use of a non-constant weight function implies an anisotropic diffusion which better preserves image features.

f 0 = [f 0 1 , . . . , f 0 N ] T with f 0 : V ⊂ Z 2 → R. Figure
In the case of color images, the function f 0 : V ⊂ Z 2 → R 3 is a mapping from the vertices of the graph to a vector of color channels. The regularization can be applied on each channel leading to an iteration of filter (13) rewritten as: [START_REF] Bensoussan | Difference equations on weighted graphs[END_REF] where γ (t) uv depends on the norm of the p-Laplace operator defined by

  f c1 f c2 f c3   (t+1) v = pλ + u∼v γ (t) uv -1   pλ   f c1 f c2 f c3   0 v + u∼v γ (t) uv   f c1 f c2 f c3   (t) u    ,
▽ v f = ▽ v f c1 2 + ▽ v f c2 2 + ▽ v f c3 2
. This norm takes into account the inner correlation aspect of color vector data (in the case of p = 2). The above iteration can be extended to any vector-valued function R m → R k .

Image Simplification : One can simplify an image by first considering a fine partition of this image (or over-segmentation), where the pixel values of each region of the partition is replaced by the mean or the median pixel value of this region. The partition can be structured by a region adjacency graph (RAG), where each vertex represents a region and where edges are linking adjacent regions. Figure 2(b) and(c) illustrate a fine partition of an image and its associated RAG. Let G w = (V, E) be a RAG. Let f 0 : V ⊂ Z 2 → R m be a mapping from the vertices of G w to the mean or median value of their regions. Then, the simplification is achieved by regularizing the function f 0 on the subgraph of G w composed of the edges of E for which w uv > µ > 0. Mesh Smoothing and Denoising : By nature, polygonal curves and surface meshes have a graph structure. Let V be the set of mesh vertices, and let E be the set of mesh edges. If the input mesh is noisy, we can regularize vertex coordinates or any other function f 0 : V ⊂ R 3 → R m defined on the graph G w = (V, E). Results of filter [START_REF] Bensoussan | Difference equations on weighted graphs[END_REF] are given in Fig. 3 for polygonal curves, and in Fig. 5 and Fig. 4 for surface meshes.

Conclusion

We propose a general discrete framework for regularizing real-valued or vectorvalued functions on weighted graphs of arbitrary topology. The regularization, based on the p-Laplace operator, leads to a family of nonlinear iterative filters. This family includes the TV digital filter, the nonlocal mean filter and the bilateral filter, both widely used in image processing. Also, the family is linked to spectral graph filtering. The choice of the graph topology and the choice of the weight function allow to regularize any discrete data set or any function on a discrete data set. Indeed, the data can be structured by neighborhood graphs weighted by functions depending on data features. This can be applied in the context of image smoothing, denoising or simplification. We also show that mesh smoothing and denoising can be performed by the same filtering process. The main ongoing work is to use the proposed framework in the context of hierarchical mesh segmentation and point cloud clustering.

Fig. 1 .

 1 Fig. 1. Application to image denoising. Regularizations performed with t = 5. (c) Unweighted L2 regularization (13) on an 8-adjacency graph. (d) Unweighted T V regularization (13) on an 8-adjacency graph (p = 1). (d) Regularization (13) with p = 2 and weight function (16) and scalar feature F f = f on an 8-adjacency graph. (e) Regularization (13) with weight function (15) and feature vector F f (v) = [fu : u ∈ Bv,7] T on an ǫ-neighborhood graph with ǫ = 4.

Fig. 2 .

 2 Fig. 2. Image simplification. (a) Original image I. (b) The region map R obtained by an energy partition of I. (c) The RAG Gw associated with R, weighted with weight function (15), h = 8 and F f (v) = f . The function considered on the vertices corresponds to the mean value of a region. (d) The mean region map associated with R. (e,f) Regularization of Gw with filter (13).

Figure 2 ( 1 Fig. 3 .

 213 Fig.3. Polygonal curve denoising. Edges of polygons are weighted with equation[START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF]. The regularization of vertices position is performed in 10 steps. When λ = 0, the regularization introduces shrinkage effects. They are reduced using a value of λ = 0.

original mesh normal noise p = 2 ,Fig. 4 .

 24 Fig. 4. Mesh denoising by regularizing the position of vertices (t = 10 and w =(15)).

5 Fig. 5 .

 55 Fig. 5. Mesh denoising. (a) Original Stanford Bunny with |V | = 35949. (b) Noisy Bunny with normal noise. (c) Regularization of vertex coordinates on mesh edges in 8 iterations. We use the weight function (15) with the scalar feature F f (v) = fv.