Vincent Vidal

Christian Wolf

Florent Dupont

Mechanical Mesh Segmentation and Global 3D Shape Extraction

Keywords: 2-manifold, Triangular Mesh Segmentation, Mesh Decomposition, Line of Feature Edges, Graphical Model

This paper presents a method for segmenting noisy 2-manifold meshes based on a decomposition into local shape primitives maximizing global coherence. This technique works by partitioning the input mesh into regions which can be approximated by a simple geometrical primitive such as a plane, a sphere or a cylinder. The proposed approach is entirely error-driven, convergence-proven, and does not need to specify a number of segments.

The partitioning is guided by robust shape extractions based on RANSAC sampling and by a global graphical model which regularizes the segmented regions. The final decomposition is based on the minimum of the energy associated with this graphical model.

Obtained segmentations on noisy mechanical meshes outperform other approaches in terms of region contour correctness and consistency with mechanical object decomposition. Applications of this work are reverse engineering, mesh structure analysis, mesh feature enhancement, noise removal, mesh compression, piecewise approximation of mesh geometry, and remeshing.

Introduction

In this paper, the NP-hard problem of geometric segmentation [START_REF] Suri | Surface approximation and geometric partitions[END_REF] is addressed, which aims at dividing an input mesh into polygon clusters approximated by a simple shape primitive [START_REF] Digne | Mesh Segmentation and Model Extraction[END_REF]. Existing approaches typically use planar primitives and higher-degree primitives such as spheres, cylinders, cones and tori to segment meshes, by capturing these local shape primitives or pieces of primitives on the mesh surface. Geometric segmentation is particularly suited for mechanical objects which are more likely to exhibit consistent geometric area than natural objects.

This paper presents a method for the robust geometric segmentation in the presence of noise on the vertex positions, which is appropriate for the reverse engineering, analysis, processing and approximation of noisy mechanical meshes. Noisy mechanical meshes come usually from surface reconstruction algorithms, either using 3D points obtained by a scanner device or using multi-view imagery. The noise on the vertex position has therefore several origins, among which are sensors internal noise of the capturing device, environment conditions and surface reconstruction methods.

Related Work

Segmentation of mechanical objects has been widely studied [START_REF] Agathos | 3D mesh segmentation methodologies for CAD applications[END_REF]. There are several approaches to address the geometric segmentation, which consist in several tactics to cluster mesh polygons into semantically consistent regions, among which are hierarchical face clustering [START_REF] Garland | Hierarchical face clustering on polygonal surfaces[END_REF][START_REF] Attene | Hierarchical mesh segmentation based on fitting primitives[END_REF], variational face clustering [START_REF] Cohen-Steiner | Variational shape approximation[END_REF][START_REF] Wu | Structure recovery via hybrid variational surface approximation[END_REF][START_REF] Yan | Quadric surface extraction by variational shape approximation[END_REF], variational face clustering and model selection [START_REF] Li | Variational Surface Approximation and Model Selection[END_REF], RANSAC extraction [START_REF] Schnabel | Efficient RANSAC for Point-Cloud Shape Detection[END_REF][START_REF] Li | GlobFit: Consistently Fitting Primitives by Discovering Global Relations[END_REF], curvature guided model extraction [START_REF] Digne | Mesh Segmentation and Model Extraction[END_REF], and spectral clustering with boundary smoothing [START_REF] Zhang | Variational Mesh Decomposition[END_REF].

In most of approaches, the two key ingredients are the partitioning (labeling) of polygons and the re-estimation (fitting) of shape primitives onto their region. The partitioning step tries to find clusters of polygons with the desired properties such as the total number of regions, the size or the perimeter length of each region (boundary smoothing). The (local) fitting step estimates the parameters of the best shape primitives mostly using least-squares (regression). This fitting can be done over triangles [START_REF] Digne | Mesh Segmentation and Model Extraction[END_REF][START_REF] Cohen-Steiner | Variational shape approximation[END_REF][START_REF] Wu | Structure recovery via hybrid variational surface approximation[END_REF][START_REF] Yan | Quadric surface extraction by variational shape approximation[END_REF][START_REF] Li | Variational Surface Approximation and Model Selection[END_REF] and point clouds [START_REF] Schnabel | Efficient RANSAC for Point-Cloud Shape Detection[END_REF][START_REF] Li | GlobFit: Consistently Fitting Primitives by Discovering Global Relations[END_REF]. For meshes, the square of the Euclidean distance is often integrated over polygons leading to the L 2 metric [START_REF] Cohen-Steiner | Variational shape approximation[END_REF], which compensates for irregular sampling.

In case of the presence of outliers and noise on the vertex positions, a robust fitting is done over inlier points [START_REF] Schnabel | Efficient RANSAC for Point-Cloud Shape Detection[END_REF] or inlier triangles [START_REF] Li | Variational Surface Approximation and Model Selection[END_REF]. Li et al.'s [START_REF] Li | GlobFit: Consistently Fitting Primitives by Discovering Global Relations[END_REF] work suggests that parameter estimation of geometrical primitives should be subject to global constraints to improve the robustness of the fitting. For instance, some plane primitives almost parallel or orthogonal should be forced to be parallel or orthogonal. However, global parameter estimation is much more computationally expensive than local estimation, and such global priors should be handled with care, since only an expert user may fix such constraints.

Model selection, which consists in finding the best description of the data, is a key step in recent segmentation frameworks. Li et al. [START_REF] Li | Variational Surface Approximation and Model Selection[END_REF] have improved the global consistency of the final decomposition by optimizing the number and type of geometrical primitives according to the MDL (Minimum Description Length) criterion. However, their optimization is not global, which may result in poor local minima. The global energy approach of Isack and Boykov [START_REF] Isack | Energy-Based Geometric Multi-model Fitting[END_REF], for robust extraction of lines in 2D point clouds, optimizes the number of consistent lines extracted in a noisy point clouds, and can deal with higher level of noise than approach based on greedy selection of primitives with the most of inliers such as RANSAC.

Contributions

This paper presents a new method for the partitioning of 2-manifold triangular meshes, which is robust to the presence of outliers and noise on the vertex positions. To our knowledge, our framework introduces two novelties :

• the joint segmentation of mesh triangles and mesh edges in order to favor lines of feature edges pass along region boundaries ;

• and the underlying global (graphical) model which sets the label dependencies and the regularization power of each soft constraint.

Our framework optimizes a globally coherent model including the following soft constraints:

• label coherence of neighboring triangles: if an edge is not a feature, adjacent triangles should be associated with the same label (spatial regularization prior) ; on the contrary, in case of a feature edge, adjacent triangle labels should be different ;

• label cost: each type of shape primitives has its own cost (label cost prior), thus allowing the global optimization of the type and number of shape primitives.

Segmentation algorithm overview

The problem of geometric segmentation is clearly continuous. Indeed, each type of shape primitive is described by a set of real-valued parameters. We therefore resort to a discretization of the problem, i.e. manage a set of candidate primitives, which involves the use of three strategies in an iterative process:

• initial discretization is based on RANSAC : according to mesh triangle barycenters and normals, a set of best primitive candidates is extracted ;

• globally coherent partitioning using current set of shape primitives ;

• discretization is refined according to the current partition: if the refitting of a shape primitive onto its inlier triangles leads to a noise variance greater than a predefined maximum variance, new shape primitives are extracted using RANSAC, otherwise the refitted primitive is added to the primitive candidates.

Using our discretization to the geometric segmentation problem, our global model (see section 2) permits to automatically select an excellent subset of shape primitives. Our framework is general, as it encompasses MDL-like criteria: the type and number of shape primitives are globally optimized according to user requirements. In addition, shape primitives are not restricted to approximate connected set of triangles, which permits to decompose a mesh with fewer shape primitives, in particular for mechanical objects where some shape primitives are an extension of others. For instance, the gear model in figure 2, the noisy fandisk model in figure 3, and the mechanical model in figure 8 have a cylindrical proxy associated with a 2-connected region.

Our mesh segmentation model is presented in the next section and the following section deals with the energy function related to our graphical model. Section 4 addresses the minimization of the energy function and the management of shape primitive candidates. Afterwards, experimental results on noisy mechanical meshes and scanned objects are discussed.

Model for Mesh Segmentation and Shape Extraction

In our work, mesh segmentation is coupled with extraction of lines of feature edges and is formulated as an energy minimization problem, where triangles and edges are labeled. Each triangle label corresponds to a shape primitive approximating its region and each binary edge label indicates (or not) the presence of a feature edge along a feature line. In the following, we introduce relevant definitions and notations used throughout the paper, as well as the model variables and its goals.

Definitions and Notations

A feature edge is either a boundary edge or an edge along which a sharp discontinuity of the normal field occurs. A region is a set of triangles associated with the same label. In contrast to frequent definitions, non-connected triangles are allowed to be clustered together. As in previous work [START_REF] Cohen-Steiner | Variational shape approximation[END_REF][START_REF] Wu | Structure recovery via hybrid variational surface approximation[END_REF][START_REF] Li | Variational Surface Approximation and Model Selection[END_REF], a shape primitive approximating a region is called shape proxy or proxy for short.

The following notations are used:

• a capital letter X represents a family of random variables, while a lower case letter x represents a realization of X ;

• a subscripted capital letter X i represents a random variable associated with the node i of a graphical model, while a subscripted lower case letter x i represents a realization of X i .

Model Variables

For the joint segmentation of triangles and edges, we propose an undirected graphical model, whose nodes are divided into two categories:

• S (t) represents the nodes associated with mesh triangles ; the neighboring relationship is induced by the mesh: two nodes are neighbors iff the two corresponding mesh triangles share a common mesh edge ;

• S (x) are nodes linked to mesh edges ; the neighboring relationship is induced by the mesh: two nodes are neighbors iff the two corresponding mesh edges share a common mesh vertex.

The corresponding edges of this graphical model are:

• E (t) denotes the set of edges between two neighboring nodes in S (t) ;

• E (x) accounts for the set of edges between two neighboring nodes in S (x) .

The set of edges / ∈ E (t) E (x) , which are between triangle and adjacent edge sites, is not associated with a specific notation. Note that edges are undirected and thus (i, j) and (j, i) represent the same edge. Each node of our graphical model is associated with a discrete random variable. Let T and X be, respectively, a family of random variables over sites S (t) and S (x) . Let G and Y be, respectively, a family of observed variables over sites S (t) and S (x) . Observed variables integrate sensed geometrical features with respect to our global model. g includes the geometry of mesh triangles needed for computing the approximation error between the geometric triangle g i and its proxy P t i , while y is made of all geometrical features used in the model of line of feature edges. P t i (resp. P) is the proxy of the triangle i (resp. the set of all proxies).

Model Objectives

The goals of our graphical model are threefold:

• to label triangles with a label corresponding either to a shape proxy index {0,. . . ,|P|-1} or to -1 ; in the latter case, the kind of proxy is ∅ ; here -1 label means either outliers or cannot be well approximated by current proxy candidates ;

• to label mesh edges by x i ∈ {non-feature=0, feature=1} ;

• to make triangle and edge labeling consistent with each other, by modeling the joint distribution of triangle and edge labels.

Segmentation Pipeline

The whole segmentation process is illustrated in figure 1. A set of shape proxies is initially obtained by a RANSAC discretization of the continuous space of shape primitive parameters. Then, a subset of these shape proxies is assigned to triangle variables by a joint labeling of triangles and edges. More precisely, the labeling process will set triangle variables to a shape proxy index and will set edge variables either to non-feature or feature. Afterwards, the current segmentation is used to update the current discretization by adding new shape proxies and removing unused ones in the set of candidate primitives. The algorithm terminates when the set of candidate proxies and the current segmentation are stable, which occurs at a (good) local minimum of the energy minimized. During the segmentation process, all proxy selection steps are based on the global minimization of the objective function [START_REF] Suri | Surface approximation and geometric partitions[END_REF], presented in the next section, and thus each triangle label change is related to an energy decrease. Based on the current mesh segmentation, a new set of candidates is added to the proxy set. Eventually, when the set of primitives and the mesh segmentation are stable, a final primitive fitting is done to decrease further the energy E. The possibility to remove a shape primitive from the set of candidate primitives is not drawn here for clarity purpose.

Objective Function

The objective function (1) assigns an energy to each possible solution of the global labeling problem and evaluates to lower energy values for more favorable solutions. It is a weighted scalar sum of scale-invariant energy potentials, where five types of energy potentials appear:

• E f eat (x i ; y) is the local evidence of edge i ∈ S (x) being a feature edge or not and is defined in equation (3) ; it is based on machine learning and takes into account the prediction from a statistically learnt model to be less sensitive to the presence of noise on the vertex positions ;

• E dir (x i , x j ; θ, φ)
is a spatial regularization model on couples of neighboring edge sites (i, j) ∈ E (x) , which favors almost aligned couples of feature edges which have similar dihedral angle and is defined in equations ( 4) and ( 5) ;

• E approx (t i ; g i ; P t i ) is called residual error and measures the approximation error between the geometric triangle g i from site i ∈ S (t) and its shape proxy P t i and is defined in equation ( 2) ;

• E reg (t i , t j , x ψ(i,j)
) is a spatial regularization prior on triplet (t i , t j , x ψ(i,j) ) of two adjacent triangles i, j ∈ S (t) and their shared edge ψ(i, j) ∈ S (x) , and is defined in equation ( 6) ; The energy function is given as

• β k |L (L k ) |
E(t; x; g; y; P) = κ i∈S (x) E f eat (x i ; y) +λ (i,j)∈E (x) E dir (x i , x j ; y) +ρ i∈S (t) E approx (t i ; g i ; P t i ) + (i,j)∈E (t) E reg (t i , t j , x ψ(i,j) ) + k β k |L (L k ) | (1)
where κ is a constant adjusting the strength of edge attachment terms of data compared with those of the triangles. λ is a positive scalar used to control the regularization strength of the edge spatial prior. ρ is a positive weight, which permits to monitor the strength of the approximation error with respect to regularization priors. The terms E reg (t i , t j , x ψ(i,j) ) have no visible weight in the equation ( 1) because their weights are internally parameterized. β k costs offer a fine control onto the number of shape proxies of each type, but also permit a more compact description of data through the merging of proxies close enough.

Triangle Residual Error

The triangle residual error E approx (t i ; g i ; P t i ) depends on the label t i of a site i ∈ S (t) , on the geometrical triangle g i and on the shape proxy P t i among all proxy candidates P:

E approx (t i ; g i ; P t i ) = L 2 (g i , P t i ) if P t i = ∅ γ × |g i | else (2)
where L 2 (g i , P t i ) = p∈g i pi (p) 2 dp is a metric which measures the error between the geometrical triangle g i and its shape proxy P t i , and i (p) denotes the orthogonal projection of point p onto proxy P t i . As in VSA [START_REF] Cohen-Steiner | Variational shape approximation[END_REF], the normalization by the region area and the square root have been removed, because they are not useful for optimization purpose. An exact formula for L 2 (g i , P t i ) is given in [START_REF] Cohen-Steiner | Variational shape approximation[END_REF] in case of planar proxies. This formula is also used in case of spherical or cylindrical proxies as a good approximation [START_REF] Wu | Structure recovery via hybrid variational surface approximation[END_REF]. |g i | means the area of triangle g i . γ × |g i | is the cost of a non-associated triangle, it depends on the triangle area as does the L 2 metric. γ is a positive constant which controls the distance of a triangle to its approximating proxy: if γ is large enough, all triangles will be associated to a proxy candidate. On the contrary, if γ is close to zero, a large amount of triangles may remain non-associated along with the segmentation process. To achieve results invariant to scale, the input mesh is normalized such that the largest dimension of its axis aligned bounding box is 1. At the end of the segmentation, the mesh is denormalized.

Feature Edge Model

The feature edge model included in our global model is a CRF for robust extraction of lines composed of feature edges, in the lines of the one proposed by [START_REF] Vidal | Robust feature line extraction on CAD triangular meshes[END_REF]. E f eat (x i ; y) is a dataattached term which gives lower energy for a labeling x i which is consistent with the response of a (learnt) prediction model according to observed geometrical features y.

E f eat (x i ; y) = exp((-1) x i d(y)) (3)
E f eat (x i ; y) directly incorporates classification and confidence of a SVM (Support Vector Machine) -which depends on the obtained prediction model -through the signed distance d(y) of the mapped feature vector associated with edge i to the separating hyperplane (in feature space). It is positive for predicted feature edges and negative otherwise. The feature vector is made of 43 geometrical features such as the edge dihedral angle or the edge mean curvature, and the SVM, a machine learning algorithm, learns how to characterize a feature edge. We employed the same features as in [START_REF] Vidal | Robust feature line extraction on CAD triangular meshes[END_REF]. E dir (x i , x j ; y) is a directional regularization term on couples of neighboring edge sites, which gives lower energy to a labeling (x i , x j ) when x i = x j , while favoring the case of two feature edges only when they are rather aligned. The edge coherence term is given by

E dir (x i , x j ; y = [θ i , θ j , φ ψ (i,j) ])) =    0 if x i = x j -τ if x i = x j = 0 -Γ(θ i , θ j , φ ψ (i,j) ) else (4)
where τ is a positive constant which controls the homogeneity of normal edges. θ i ∈ y is the dihedral angle of edge i. φ ψ (i,j) ∈ y is the turning angle between edges i and j. Γ (cf. equation 5) is a positive functional used in the 2 feature edges case, which reaches larger values for geometrically aligned (within a tolerance) neighboring edges if they have similar dihedral angles.

Γ(θ i , θ j , φ ψ (i,j) ) = exp        -∆     |cos(θ i )-cos(θ j )| σ similarity + 1-cos(φ ψ (i,j) ) alignment            (5)
In equation ( 5), ∆ is a positive constant which monitors the neighboring edge alignment tolerance. The greater ∆ is, the more two neighboring edges must be aligned. σ handles the tolerated variance of ∆|cos(θ i )-cos(θ j )| along feature lines. The similarity term avoids detecting adjacent feature edges with a very different dihedral angle. The alignment term predominates over similarity term, for favoring non-jagged feature line extraction.

Consistency between Triangle and Edge Labels

The goal of the consistency term is to favor equal (resp. different) triangle labels between two neighboring triangles when their shared edge is not (resp. is) a feature. The consistency term is defined as follows:

E reg (t i , t j , x ψ(i,j) ) =        µ if t i = t j and x ψ(i,j) = 0 ν if t i = t j and x ψ(i,j) = 1 ζ if t i = t j and x ψ(i,j) = 0 η if t i = t j and x ψ(i,j) = 1 (6)
where µ, ν, ζ, and η are positive constants which penalize, respectively, the label difference of two triangles sharing a normal edge, the label difference of two triangles sharing a feature edge, the label equality of two triangles sharing a normal edge, and the label equality of two triangles sharing a feature edge. Note that when µ = ν ≥ 0 and ζ = η = 0, the consistency term corresponds to the Potts/Ising model [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF], independent of the type of an edge, which is used by the PEaRL algorithm [START_REF] Isack | Energy-Based Geometric Multi-model Fitting[END_REF]. Therefore, our proposed consistency term offers richer interactions. However, because we want to favor equality (resp. inequality) of the labels t i and t j of two neighboring triangles only when the shared edge is not (resp. is) a feature, we restrict the parameter setting to µ = ν > 0, ζ = 0 and η ≥ µ, ν, ζ. Intuitively, this setting means that the consistency term: gear 26 planes, 2 cylinders Figure 2: Results of our method on a gear (with feature edges): (left) one color per shape primitive index and (right) one color per shape primitive type (blue = plane, green = cylinder). Our method can segment non-noisy mechanical models by increasing the cost related to the approximation error. Two cylinder proxy associated with non-connected regions are successfully found by our framework.

• penalizes the difference of neighboring triangle labels,

• and penalizes even more the equality of neighboring triangle labels when their shared edge is a feature.

Minimization and Proxy Management

In general, in discrete optimization, better guarantees can be obtained in the binary case (two possible values for each variable) than in the multi-label case which is often solved using successive binary decisions between keeping the current label or taking the new one.

In the binary case, where only two shape proxies P 0 and P 1 are available, when E reg (t i , t j , x ψ(i,j) ) is submodular, the objective function (1) can be globally minimized using the graph cuts framework, especially Kolmogorov and Zabih's st-graph construction method [START_REF] Kolmogorov | What Energy functions can be Minimized via Graph Cuts[END_REF]. The submodularity of E reg depends on the parameters of the functional and means that

E reg (0, 0, 0) + E reg (1, 1, 0) ≤ E reg (0, 1, 0) + E reg (1, 0, 0) (7)
and E reg (0, 0, 1) + E reg (1, 1, 1) ≤ E reg (0, 1, 1) + E reg (1, 0, 1). ( 8) Submodularity (cf. [START_REF] Kolmogorov | What Energy functions can be Minimized via Graph Cuts[END_REF]) is an important concept in discrete optimization. Its meaning is the equivalent of the concept of convexity, restricted to discrete sets. Equation ( 7) requires ζ ≤ µ which is always true. Equation ( 8) requires η ≤ ν, which is generally false, since we would like to penalize more the case of triangle label equality with a feature edge. Consequently our energy function is not submodular. Fortunately, there exist minimization techniques able to deal with non-submodular terms such as BP (Belief Propagation) [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems[END_REF] or QPBO (Quadratic Pseudo-Boolean Optimization) [START_REF] Rother | Optimizing Binary MRFs via Extended Roof Duality[END_REF]. While both methods give approximated results on nonsubmodular functions, the convergence of BP is not guaranteed in the presence of cycles in the graphical model, which leads us to choose QPBO whose convergence is guaranteed. In addition, QPBO gives optimal results on submodular functions.

In the multi-label case (at least 3 proxies), successive binary decisions between keeping the current label or choosing the new label (one loop over all possible labels for one iteration) are taken using QPBO. If the objective function ( 1) is submodular, a strong local minimum is found by successive QPBO, that is to say a minimum within a factor of the global optimum noisy fandisk 13K triangles 8 planes, 2 spheres, 6 cylinders noisy fandisk 52K triangles 8 planes, 2 spheres, 6 cylinders

Figure 3: Results of our method (with feature edges): (left) one color per shape primitive index and (right) one color per shape primitive type (blue = plane, red = sphere, green = cylinder). Segmented models show that our method is robust to the sampling density once there are enough proxy inliers. Our framework permits to extract a cylinder over 2 nonconnected regions while VSA-like approaches will use two cylinders instead.

(see [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], section 6). QPBO guarantees that after each binary decision the global energy does not increase, which proves the convergence of our segmentation algorithm via the minimization approach presented in section 4.1.

A last advantage of QPBO, in both binary and multi-label cases, is that it experimentally finds a good local minimum when the number of non-submodular terms is low.

Minimization

The objective function ( 1) is hard to minimize, because it is not submodular for selected parameters. In the proposed method, (1) is iteratively minimized by α-expansion move with label cost [START_REF] Delong | Fast approximate energy minimization with label costs[END_REF]: starting from an initial configuration, and for each available triangle label, a binary decision is done for each triangle between keeping the current label or the new one. More precisely, for each binary decision ( 1) is approximately minimized with QPBO [START_REF] Rother | Optimizing Binary MRFs via Extended Roof Duality[END_REF].

Noting that in general the number of normal edges is largely dominant on the number of feature edges, and that for normal edges, consistency terms are submodular, we finally opt for a two steps minimization procedure: first of all, the line of feature edges are extracted using the method in [START_REF] Vidal | Robust feature line extraction on CAD triangular meshes[END_REF] ; then, shape primitive extraction and global selection steps are iterated until convergence: the selection step will use our graphical model, while observing edge labels. Because QPBO works well when the number of non-submodular terms is low, an excellent local minimum is generally found. Note that the method for feature line extraction is now independent from the global model and thus another line extraction method may be used.

Management of Proxy Candidates

The management of proxy candidates begins with the initial proposal of shape proxies, the discretization step: starting from a triangle soup, the largest set of triangles supporting either a planar, a spherical or a cylindrical proxy is found by a RANSAC. Then, a shape primitive fitting is done on these inliers triangles using least-squares and the obtained shape proxy is added to the set of proxy candidates. The discretization continues on the remaining triangles up to no more shape proxy can be extracted. Once the initial proxy set is found, it will be modified by proxy additions (section 4.2.1) and proxy removals (section 4.2.2).

Refined Discretization according to Mesh Segmentation

After a labeling step, each segmented region associated with a shape proxy is analyzed as follows:

1. current proxy parameters are re-estimated in the least-square sense over all region triangles ;

2. parameters are checked to make sure they are meaningful: if a spherical or cylindrical proxy has a radius greater than a specified maximum radius, the best planar proxy is extracted instead of current proxy, using RANSAC ;

3. the last step is an intra-region variance monitoring of the newly estimated proxy: if the current new proxy is associated with a variance below a maximum variance, this new proxy is added to the set of proxy candidates ; otherwise, the best proxy of each type is extracted by a RANSAC and added to the list of candidate proxies ; in the case of a planar proxy with too high variance, the two best planar proxies are added instead of the best planar proxy.

For regions labeled as ∅, i.e. non-associated, the best proxy of each type is extracted by a RANSAC and added to the list of proxy candidates. Note that at this point, no decision to use a newly proposed proxy has been fixed. This will be done in the next joint segmentation step that will make the final decision among old and new candidates. One may expect to propose additional proxy candidates based on the merging of neighboring regions. However, assuming that RANSAC extraction has found all proxies associated with intermediate size region, this possibility is not necessary.

Shape Proxy Removals

A high number of proxy candidates may degrade the performances of our segmentation algorithm, because the chosen algorithm for discrete optimization is linearly dependent on the number of triangle label candidates. As proposed by Li at al. [START_REF] Li | Variational Surface Approximation and Model Selection[END_REF], a shape proxy candidate not supported by any triangle of the initial mesh during two successive steps of proxy discretization and global selection, should be removed. This heuristic works well in practice and experimentally outperforms similar heuristics using different values for the number of successive whole steps. However, we found by experimentation that sorting in decreasing order the set of candidates according to their number of support points allows to run only one α-expansion (with label cost) for the selection step (cf. section 4.1).

Experimental Results

To demonstrate the efficiency of our method, we applied it to several mechanical models with noise on their vertex positions. More specifically, there are 11 normalized 2-manifold meshes among which 5 are obtained by adding 0.5% Gaussian noise intensity to mesh positions and the 6 remaining models are composed of 4 scanned objects and 2 mechanical objects. The normalization process consists in mesh scaling to make the maximal dimension of a mesh bounding box equal to one, and as mentioned earlier is needed to make our segmentation results scale-invariant. 0.5% thus refers to the maximal bounding box dimension.

The presented results have been obtained on an Intel Core I7-2760QM (2.4 GHz) with 8GB RAM, with a basic RANSAC sampling (not the one in [START_REF] Schnabel | Efficient RANSAC for Point-Cloud Shape Detection[END_REF]).

Models

Discussion

To illustrate the ability of our mesh segmentation method to capture the right decomposition of a mechanical object, we applied it onto simple noisy meshes, whose ideal decomposition is obvious: cylinder, screw, joint and sphere. Table 1 and figure 6 confirm that the obtained decompositions are the right ones. To evaluate the ability to segment more complex mechanical objects, we applied our method on a gear model (see figure 2), on the noisy fandisk models (see figure 3) and on two scanned mechanical objects (see figure 4). We also applied it to two scanned natural objects for which no feature edge has been extracted (see figure 5). Visual results show that our method is able to extract a compact surface description made of shape proxies with a large support. Regions approximated by a proxy may be non-connected. If the end-user needs connected regions, a simple post-processing step which extracts connected components can be used. In addition, region complexity is kept low and contour smoothness is encouraged by our method, because of the spatial regularization prior which favors compact regions with as few as possible jagged contour edges. Note that if crest lines are available for natural objects, their use in our algorithm as feature lines input will make the final segmentation respecting the minima-rule [START_REF] Hoffman | Salience of visual parts[END_REF], i.e. a crest line will be at the junction of two different regions.

The joint segmentation of mesh triangles and edges permits to attract region boundaries towards lines of features edges, which improves the global coherence of final results, and also avoids overflowing of one region label across the line of feature edges. In addition, our method is not sensitive to isolated feature edges, that is, region boundaries pass along relevant line of feature edges if there is a line of feature edges close enough.

Our geometric segmentation can cope with noise on the vertex positions, and is robust to irregular (non-uniform) sampling, in particular to remeshing. The extracted mesh decompositions in figure 3 for the fandisk model are indeed very similar for two different sampling rates. The natural model decomposition into simple shape primitives is therefore not sensitive to the density of mesh vertices if a sufficient number of inliers for each mesh region is present (for RANSAC). However, if the region boundaries are affected by a mesh pre-processing, the final segmentation may slightly differ.

Comparisons with previous work

We compare our method to HFC [START_REF] Attene | Hierarchical mesh segmentation based on fitting primitives[END_REF] and VSA [START_REF] Cohen-Steiner | Variational shape approximation[END_REF][START_REF] Wu | Structure recovery via hybrid variational surface approximation[END_REF] approaches using our own implementation based on the fitting over mesh triangles (in original HFC paper, the fitting is done over points). We use our own implementation, because the proxy fitting influences the final segmentation ; thus by using the same fitting algorithm for the 3 methods we are able to compare the 3 frameworks without either the fitting method or fitting implementation bias. HFC and VSA need the user to set a number of segments to work, whereas our method does not need it. As can be seen in figure 7, our method gives better results in terms of region complexity and contour smoothness due to the absence of spatial regularization prior in HFC and VSA methods. The worst region boundaries are obtained for HFC, because of its greedy nature which can get stuck into a poor local minimum. Note in the same figure that VSA may suffer from small island configuration, while it does not experimentally appear for our method. In figure 8 (a), we can see that our framework is able to recover shape proxy associated with a non-connected region (the red and white cylinders), while HFC and VSA are not able to do so. Several previous works applied their method to the bunny model [START_REF] Attene | Hierarchical mesh segmentation based on fitting primitives[END_REF][START_REF] Wu | Structure recovery via hybrid variational surface approximation[END_REF][START_REF] Yan | Quadric surface extraction by variational shape approximation[END_REF][START_REF] Li | Variational Surface Approximation and Model Selection[END_REF]. On average, the segmented bunny with our method without feature line has a more compact description with region boundaries less complex. Our method seems also to preserve the symmetries on the bunny model.

Implementation Details and Settings

Parameters have been experimentally fixed with the values presented in table 2. Penalizing less a proxy type is similar to taking into account a priori knowledge about more probable proxy type. Cylindrical proxies are usually more present in mechanical objects than spherical proxies, that is why the cost of a spherical proxy is greater than the cost of a cylindrical proxy. Penalizing more a proxy will also favor the merging of two close proxies of the same type, which is the case with the chosen label costs. By increasing all label costs, an end-user can obtain another level of detail for the mesh segmentation. However, the best way to control different levels of detail is by the weight ρ used to penalize the residual error between triangles and their proxy.

Another relevant a priori knowledge is the maximum radius size of kept spherical and cylindrical proxies, which is set to 7 × M BB for both. M BB is the maximal bounding box dimension. Indeed, without such an information, a planar area may be approximated by a cylindrical proxy whose radius is very large.

Let mBB be the minimal dimension of the mesh bounding box. Concerning the RANSAC thresholds needed for the initial primitive extraction, we experimentally set the distance threshold to 4.10 -2 × mBB and the deviation threshold to 14.07 degrees. In the case of RANSAC extraction based on the current segmentation, we experimentally set the distance threshold to 8.10 -3 × mBB and the deviation threshold to 11.48 degrees. For the shape proxy (c) VSA. Our method automatically finds the best number of proxies, 8, while for VSA and HFC we experimentally found 10 and 16 as the best solution. Our framework describes the mesh geometry with fewer proxies and permits to extract a cylinder over 2 non-connected regions while VSA-like approaches will use two cylinders instead.

initial estimate, needed for RANSAC, we used those proposed in [START_REF] Schnabel | Efficient RANSAC for Point-Cloud Shape Detection[END_REF]. For the planar proxy fitting, we use the classical method based on Principal Component Analysis [START_REF] Attene | Hierarchical mesh segmentation based on fitting primitives[END_REF]. For the fitting of spherical and cylindrical proxies, starting from the initial estimate, we use conjugate gradient minimization for triangle soups with less than 90 triangles, and otherwise we use the Attene et al.'s fitting method [START_REF] Attene | Hierarchical mesh segmentation based on fitting primitives[END_REF], a computationally efficient approach. With a high number of triangles, Attene et al.'s fitting method is sufficiently accurate, while for a low number of points gradient descent is much more accurate.

For the α-expansion move with label cost [START_REF] Delong | Fast approximate energy minimization with label costs[END_REF], which is an energy minimization algorithm able to minimize an energy function in the multi-label case, we use our implementation based on the QPBO graph cut technique given in Rother et al. [START_REF] Rother | Optimizing Binary MRFs via Extended Roof Duality[END_REF].

Conclusion and future work

In this paper, we have presented a graphical model for the geometric segmentation of 2manifold triangular meshes, coupled with the extraction of lines of feature edges. The obtained results compare favorably to previous work: our method is not sensitive to noise on the vertex positions, not sensitive to irregular surface sampling and is able to capture the right decomposition of simple mechanical objects, without specifying the number of regions.

As future work, we will investigate learning of our model parameters using scanned mechanical objects whose decomposition into shape primitives is known. We would like to include a finer control on the authorized noise variance to cope with meshes with non-uniform noise on their vertex positions. We wish to propose a hierarchical version of our graphical model, to favor consistent decomposition at several scales. We would like to improve the graphical model of line of feature edges, to make it more adapted to scanned objects by taking into account extended dihedral angle integrated along geodesic lines at different scales. Finally, we wish to take into account the area of the extracted regions into the RANSAC score function to favor proxies associated with large area first instead of a large number of inlier triangles.

Figure 1 :

 1 Figure 1: Pipeline of our method. A first RANSAC model extraction initializes the set of geometrical primitives. A global selection is then completed minimizing the objective function (1). Based on the current mesh segmentation, a new set of candidates is added to the proxy set. Eventually, when the set of primitives and the mesh segmentation are stable, a final primitive fitting is done to decrease further the energy E. The possibility to remove a shape primitive from the set of candidate primitives is not drawn here for clarity purpose.

  penalizes the proxy use of type k ∈ {∅, plane, sphere, cylinder} by the positive cost β k ; L (L k ) is the set of indices of shape proxy of type k used in the current configuration L = k L k ; |.| is the number of elements in the set (cardinality) ; β k |L (L k ) | serves as an automatic model selection method.

Figure 4 :

 4 Figure 4: Results of our method: feature lines extracted using dihedral angle have been used for the blade model, while no feature line has been used for the rockerArm model: (left) one color per shape primitive index and (right) one color per shape primitive type (blue = plane, green = cylinder). Feature line information avoid region overflowing.

κ = 1 λ

 1 = 0.2 ρ = 11150 β plane = 0.102 β sphere = 0.110 β cylinder = 0.106 β ∅ = 0.102 γ = 10 -2 τ = 0 ∆ = 15 σ = 10 µ = ν = 6.10 -3 ζ = 0 η = 8.10 -3

Figure 5 :

 5 Figure 5: Results of our method (without feature edge): (left) one color per shape primitive index and (right) one color per shape primitive type (blue = plane, red = sphere, green = cylinder). Our method works well on natural objects: the symmetries of the bunny model are preserved in its decomposition and the fingers of the hand model are associated to either one or two cylinders depending on the finger' bending.

Figure 6 :Figure 7 :

 67 Figure 6: Results of our method on noisy mechanical objects (with feature edges): (left) one color per shape primitive index and (right) one color per shape primitive type (blue = plane, red = sphere, green = cylinder). Our method finds the right decomposition of simple mechanical objects.

Figure 8 :

 8 Figure8: Shape approximation of a mechanical model using (a) our method, (b) HFC and (c) VSA. Our method automatically finds the best number of proxies, 8, while for VSA and HFC we experimentally found 10 and 16 as the best solution. Our framework describes the mesh geometry with fewer proxies and permits to extract a cylinder over 2 non-connected regions while VSA-like approaches will use two cylinders instead.

Table 1 :

 1 Statistics on the segmented models: number of triangles and edges, the running time, and the number of planar, spherical and cylindrical primitives in the final segmentation.

		#t (K) #e (K) Time #p #s #c
				(s)
	cylinder	2.36	3.54	3	2	0	1
	screw	2.48	3.72	2	4	0	3
	joint	6.02	9.02	13	9	0	3
	sphere	9.80	14.70	5	0	1	0
	fandisk	12.99	19.48	31	8	2	6
	blade	25.20	37.80	373	5	0	11
	rockerArm	80.35	120.53	988	3	0	8
	hand	11.69	17.54	66	0	1	12
	bunny	69.67	104.50	377	2	6	2
	mech	1.86	2.78	5	6	0	2
	gear	5.31	7.97	56	26	0	2

Table 2 :

 2 The parameters of our model (cf. equations 1, 2, 4, 5, and 6). ρ = 26650 (resp. 30000) in figure2(resp. figure4).

	hand	0 plane, 1 sphere,
		12 cylinders
	bunny	2 planes, 6 spheres,
		2 cylinders


Romain Guigourès

Dominique Gay

Marc Boullé

Fabrice Clérot

Fabrice Rossi

Country-Scale Exploratory Analysis of Call Detail Records Through the Lens of Data Grid Models

published or not. The documents may come

Introduction

Telco operators' activities generate massive volume of data, mainly from three sources: networks, service platforms and customers data bases. Particularly, the use of mobile phones generates the so called Call Detail Records (CDRs), containing information about end-point antenna stations, date, time and duration of the calls (the content of the calls is excluded). While this data is initially stored for billing purpose, useful information and knowledge (related to human mobility [START_REF] Wang | Human mobility, social ties, and link prediction[END_REF][START_REF] Becker | Human mobility characterization from cellular network data[END_REF], social interactions [START_REF] Vieira | Characterizing dense urban areas from mobile phone-call data: Discovery and social dynamics[END_REF] and economic activities) might be derived from the large sets of CDRs collected by the operators.

domain, e.g., [START_REF] Frías-Martínez | An agent-based model of epidemic spread using human mobility and social network information[END_REF] suggest a model for epidemic spread. Mobile phones may also provide other types of data (e.g., the Nokia Mobile Data Challenge [START_REF] Laurila | From big smartphone data to worldwide research: The mobile data challenge[END_REF]), like applications events, WLAN connection data, etc. For instance, [START_REF] Jiang | A review of urban computing for mobile phone traces: current methods, challenges and opportunities[END_REF] pre-processed phone activities of one million users to obtain information about their approximative temporal location, then mined daily motifs from the spatio-temporal data to infer human activities. Finally, smart phones are or will be equipped with accelerometers and/or gyroscopes providing data about physical activities of users: [START_REF] Lockhart | The benefits of personalized smartphone-based activity recognition models[END_REF] suggest a complete system of activity recognition based on smartphone accelerometers with potential application to health monitoring. Research work related to data grid models: We are not coclustering data (objects × attributes) like pioneering work of Hartigan [START_REF] Hartigan | Direct clustering of a data matrix[END_REF]. Data grid models are related to the work of Dhillon et al. [START_REF] Dhillon | Information-theoretic co-clustering[END_REF] who have proposed an informationtheoretic coclustering approach for two discrete random variables Y 1 and Y 2 : the loss in Mutual

Information M I(Y 1 , Y 2 ) -M I(Y M 1 , Y M 2
) is minimized to obtain a locally-optimal grid with a user-defined number of clusters for each dimension. This is limited to two variables and requires to choose the number of clusters per variable. Going beyond 2D matrices, recent significant progress has been done in multi-way tensor analysis [START_REF] Sun | Beyond streams and graphs: dynamic tensor analysis[END_REF][START_REF] Kolda | Scalable tensor decompositions for multi-aspect data mining[END_REF]. Dealing with k-adic data, (also known as co-occurrence data, like contingency table), [START_REF] Peng | Temporal relation co-clustering on directional social network and author-topic evolution[END_REF] suggest a coclustering method for social network and temporal sequence (with pre-discretization of time). The Information Bottleneck (IB) method [START_REF] Tishby | The information bottleneck method[END_REF] stems from another informationtheoretic paradigm: given the joint probability P (X, Y ), IB aims at grouping X into clusters T in order to both compress X and keep as much information as possible about Y . IB also minimizes a difference in Mutual Information: M I(T, X) -βM I(T, Y ), where β is a positive Lagrange multiplier. Wang et al. [START_REF] Wang | Information bottleneck co-clustering[END_REF] build upon IB and suggest a coclustering method for two categorical variables. Extending IB for more than two categorical variables, Slonim et al. [START_REF] Slonim | Agglomerative multivariate information bottleneck[END_REF] have suggested the agglomerative multivariate IB that allows constructing several interacting systems of clusters simultaneously; the interactions among variables are specified using a Bayesian network structure. To the best of our knowledge, our summarization approach is the only one to combine the following advantages: it is parameter-free, scalable and can be applied to mixed-type attributes (categorical, numerical, thus multiple types of time dimensions without pre-processing). Therefore, the same generic method can be used to analyze network graph, temporal sequence and mobility data.

Impacts on economic strategy

Besides the high-level knowledge extracted from country-scale data and confirmed by local sociologists from the University of Bouaké in Ivory Coast, these studies have also a strong impact on future economic development strategy, mainly in two identified branches:

-Network planning strategy: In 2014, there are around 20M inhabitants in Ivory Coast and the mobile service penetration rate is 84% -with a still growing mobile phone market in a context of demographic growth. The analyses of the first two case studies and the resulting map projections (that can be seen as the network of calls available at various granularities, see Sections 5.1 and 5.2) are considered as an additional input for network planning and investment; for instance to help network designer in answering questions about how many and where the next antennas have to be set while preserving the quality of service at a reasonable cost.

Exploratory analysis through data grid models

Data grid models aim at estimating the joint distribution between K variables of mixed-types (categorical as well as numerical). The main principle is to simultaneously partition the values taken by the variables, into groups/clusters of categories for categorical variables and into intervals for numerical variables. The result is a multidimensional (K-d) data grid whose cells are defined by a part of each partitioned variable value set. Notice that in all rigor, we are working only with partitions of variable value sets. However, to simplify the discussion we will sometime use a slightly incorrect formulation by mentioning a "partition of a variable" and a "partitioned variable". In order to choose the "best" data grid model M * (given the data) from the model space M, we use a Bayesian Maximum A Posteriori (MAP) approach. We explore the model space while minimizing a Bayesian criterion, called cost. The cost criterion implements a trade-off between the accuracy and the robustness of the model and is defined as follows:

cost(M ) = -log(p(M | D) posterior ) ∝ -log(p(M ) prior × p(D | M ) likelihood )
Thus, the optimal grid M * is the most probable one (maximum a posteriori) given the data. Due to space limitation, the details about the cost criterion and the optimization algorithm (called khc) are available in appendix of [START_REF] Guigourès | Country-scale exploratory analysis of call detail records through the lens of data grid models[END_REF]. Hereafter, we focus on the tools for exploiting the grid and their applications on largescale CDRs data. The key features to keep in mind are: (i) khc is parameter-free, i.e., there is no need for setting the number of clusters/intervals per dimension; (ii) khc provides an effective locally-optimal solution to the data grid model construction efficiently, in sub-quadratic time complexity (O(N √ N log N ) where N is the number of data points). Figure 1 illustrates the input data and output results of khc on an examplary mobility data stemming from CDRs. Fig. 1: From 3D mobility data, stemming from CDRs data, to data grid.

Data grid exploitation and visualization

Because of the very large number observations in CDRs data, the optimal grid M * computed by khc can be made of hundreds of parts per dimension, i.e., millions of cells, which is difficult to exploit and interpret. To alleviate this issue, we suggest a grid simplification method together with several criteria that allow us to choose the granularity of the grid for further analysis, to rank values in clusters and to gain insights in the data through meaningful visualizations.

Dissimilarity index and grid structure simplification. We suggest a simplification method of the grid structure that iteratively merge clusters or adjacent intervals -choosing the merge generating the least degradation of the grid quality. To this end, we introduce a dissimilarity index between clusters or intervals which characterize the impact of the merge on the cost criterion.

Definition 1 (Dissimilarity index). Let c .1 and c .2 be two parts of a variable partition of a grid model M . Let M c.1∪c.2 be the grid after merging c .1 and c .2 . The dissimilarity ∆(c .1 , c .2 ) between the two parts c .1 and c .2 is defined as the difference of cost before and after the merge:

∆(c .1 , c .2 ) = cost(M c.1∪c.2 ) -cost(M ) (1)
When merging clusters that minimize ∆, we obtain the sub-optimal grid M (with a coarser grain, i.e. simplified) with minimal cost degradation, thus with minimal information loss w.r.t. the grid M before merging. Performing the best merges w.r.t. ∆ iteratively over the K variables without distinction, starting from M * until the null model M ∅ , K agglomerative hierarchies are built and the end-user can stop at the chosen granularity that is necessary for the analysis while controlling either the number of clusters/cells or the information ratio kept in the model. The information ratio of the grid M is defined as follows:

IR(M ) = cost(M ) -cost(M ∅ ) cost(M * ) -cost(M ∅ ) (2)
where M ∅ is the null model (the grid with a single cell).

Typicality for ranking categorical values in a cluster. When the grid is coarsen during the hierarchical agglomerative process, the number of clusters per categorical dimension decreases and the number of values per cluster increases.

It could be useful to focus on the most representative values among thousands of values of a cluster. In order to rank values in a cluster, we define the typicality of a value as follows.

Definition 2 (Typical values in a cluster). For a value v in a cluster c of the partition Y M of dimension Y given the grid model M , the typicality of v is defined as:

τ (v, c) = 1 1 -P Y M (c) × c j ∈Y M c j =c P Y M (c j )(cost(M |c \ v, c j ∪ v) -cost(M )) ( 3
)
where P Y M (c) is the probability of having a point with a value in cluster c, c\v is the cluster c from which we have removed value v, c j ∪v is the cluster c j to which we add value v and M |c \ v, c j ∪ v the grid model M after the aforementioned modifications.

Intuitively, the typicality evaluates the average impact in terms of cost on the grid model quality of removing a value v from its cluster c and reassigning it to another cluster c j = c. Thus, a value v is representative (say typical) of a cluster c if v is "close" to c and "different in average" from other clusters c j = c. Notice that this measure does not introduce any numerical encoding of the categories of the categorical variable under study.

Insightful visualizations with Mutual Information. It is common to visualize 2D coclustering results using 2D frequency matrix or heat map. For KD coclustering, it is useful to visualize the frequency matrix of two variables while selecting a part of interest for each of K -2 other variables. We also suggest an insightful measure for co-clusters to be visualized, namely, the Contribution to Mutual Information (CMI) -providing additional valuable visual information inaccessible with only frequency representation. Notice that such visualizations are also valid whatever the variable of interest.

M I(Y M 1 ; Y M 2 ) = J1 i1=1 J2 i2=1 M I i1i2 where M I i1i2 = p(c i1i2 ) log p(c i1i2 ) p(c i1. )p(c .i2 ) ( 4
)
where M I i1i2 represent the contribution of cell c i1i2 to the mutual information, p(c i1i2 ) is the observed joint probability of points in cell c i1i2 and p(c i1. )p(c .i2 )

is the expected probability in case of independence, i.e., the product of marginal probabilities.

Thus, if M I i1i2 > 0 then p(c i1i2 ) > p(c i1. )p(c .i2 ) and we observe an excess interaction between c i1. and c .i2 located in cell c i1i2 defined by parts i

1 of Y M 1 and i 2 of Y M 2 . Conversely, if M I i1i2 < 0, then p(c i1i2 ) < p(c i1. )p(c .i2
), and we observe a deficit of interactions in cell c i1i2 . Finally, if M I i1i2 = 0, then either p(c i1i2 ) = 0 in which case the contribution to MI and there is no interaction or p(c i1i2 ) = p(c i1. )p(c .i2 ) and the quantity of interactions in c i1i2 is that expected in case of independence between the partitioned variables.

The visualization of cells' CMI highlight valuable information that is local to the K -2 selected parts and bring complementary insights to exploit the summary provided by the grid.

Exploration results

Each application of khc (available at http://www.khiops.com) for the various case studies data is achieved within a day of computation on a commodity computer -which confirms the efficiency of the method.

Analysis of call network between antennas

The application of data grid models on the CDRs provides a segmentation with 1150 clusters, that corresponds to nearly one antenna per cluster. This is due to the large amount of data -471 millions CDRs. Indeed, the number of calls is so high for each antenna that the distribution of calls originating from (resp. terminating to) each antenna can be distinguished from each other. In order to Fig. 2: Evolution of the information kept in the data grid model w.r.t. the number of clusters using the ascending hierarchical post-processing -from optimal data grid M * (100%) to the null model M ∅ (0%). obtain a more interpretable segmentation, we apply the post-processing introduced in the Section 4.1. Figure 2 shows the information ratio (see definition 1) versus the number of clusters for all intermediate models obtained during the ascending hierarchical post-processing. Interestingly, the resulting Pareto curve shows that very informative models are obtained with few clusters. In our study, we decrease the number of clusters until keeping 60% of the model informativity -corresponding to 20 clusters, an admissible number for the interpretation. Throughout the simplification process, both partitions of source and target antennas stay identical. Thus we consider only the partition of source antennas for the rest of the study. Those clusters are projected on a map of Ivory Coast in Figure 3. Antennas are identified using dots, which color matches with the cluster they belong to.

The first observation is the strong correlation between the clusters and the geography of the country. Indeed, antennas from a same cluster are close to each other. The size of the clusters is almost the same in terms of area and match with the administrative zones of the country, except for Abidjan, the economic capital, which is split into four clusters. This is due to the high concentration of antennas in the city (32% of the ivorian antennas) and the dense phone traffic (34% of the calls).

We use the typicality (see defini-Fig. 3: Twenty clusters displayed on Ivory Coast map. There is one color per cluster.

tion 2) to rank the antennas of each cluster. The place, where the antenna with the highest typicality is located, is used to label the cluster. On the map in Figure 3, the size of the dots are proportional to the antenna typicality. Most typical antennas are located in the main cities of Ivory Coast. This phenomenon has already been observed in [START_REF] Blondel | Regions and borders of mobile telephony in Belgium and in the Brussels metropolitan zone[END_REF] and [START_REF] Guigourès | Segmentation of towns using call detail records[END_REF]: the discovered clusters match with the area of influence of the main cities of a country. We observe few exceptions: the cluster of the city of Sassandra contains the antennas of the city of Divo, while Divo is almost 4 times bigger than Sassandra (population wise) and is the sixth Ivorian city. Antennas in Divo are 40% less typical than the ones in Sassandra, meaning that allocating them to another cluster would be less costly for the criterion. Actually, calls emitted from Divo are significant in direction to other regions of Ivory Coast whereas calls from Sassandra are more internal to its region. In more formal terms, the calls distributions of the antennas in Divo are closer to the marginal distribution than to its cluster's distribution. This observation is not really surprising because Divo has experienced a recent growth of its population, due to migrations within the country [START_REF] Gnabéli | La production d'une identité autochtone en Côte d'Ivoire[END_REF]. Divo is also located in an area specialized in the intensive farming, that attracts seasonal workers from other parts of Ivory Coast. Now, focusing on the segmentation of Abidjan: the city is divided into four parts with a strong socioeconomic correlation. The first cluster -in red in Figure 3 -covers central Abidjan, including the Central Business District (le Plateau), the transport hub (Adjamé) and the embassies and upper class area (Cocody). The second cluster -in light green -is located in the South of the city. The covered neighborhoods are mainly residential areas and ports. Note that this cluster and the previous one are separated by a strip of sea, except for its North part that is included in the previous cluster. This very localized neighborhood matches with the party area of Abidjan. Finally, the last two clusters group antennas located in two areas with a similar profile: these are lower class neighborhoods. These clusters are separated not only because they are located in different parts of the city but especially because their call distribution differs: Abobo in dark blue and Yopougon in grey in the Figure 3.

Traffic between clusters. Now, we Fig. 4: Analysis of Excess of calls between clusters of antennas analyze the distribution of calls between clusters of antennas using the contribution to the mutual information. We suggest to visualize the lacks and excesses of calls between the clusters, compared to the expected traffic in case of independence. Whatever the granularity level of the clustering, we observe a strong excess of calls from the clusters to themselves and weaker excesses and lacks between clusters. Studying the traffic within the clusters has a limited interest. We only focus on the inter-clusters traffic. To visualize the traffic between clusters, we use a finer clustering than previously. Here, we have 355 clusters for 95% informativity (see Figure 2). Figure 4 depicts the excesses of traffic between clusters -highlighted with red segments. The end points of the segments are drawn at the positions of the most representative antennas of the associated clusters (i.e with the highest typicalities). The opacity of a segment is proportional to the value of the contribution to mutual information and its width is proportional to the number of calls between clusters. The biggest cities -like Bouaké, San Pedro and Man -are clearly marked on the map: they are regional capitals, a fact that is confirmed and highlighted by the call traffic visualization. The case of Bouaké is particularly interesting: although it is not the country capital, its national influence seems bigger than the one of Yamoussoukro, the actual capital. Yamoussoukro is twice smaller than Bouaké (population wise) and is a quite recent city where there is no major economical activity, contrary to Bouaké. This fact can explain our observation. We also observe that excess of traffic between major cities is a rare phenomenon. Cities are more like phone hubs, except in the West of the country around Soubré. This area is not a densely populated area but corresponds to a region with important migration flows. Finally, in Abidjan, we observe important excesses of traffic within neighborhoods, but not between neighborhoods.

Temporal Analysis of the Calls Distribution

From previous section, we learnt that the correlation between source and destination antennas is very high. The evolution of the calls distribution over time might be the same for both sets of antennas. Therefore, to track the evolution of traffic over time, we only study the evolution of the originating calls: one call is described by the emitting antenna and a day count (stemming from the date).

Again, the clustering of antennas resulting from the optimal data grid is also too fine for an easy interpretation (1051 clusters of antennas and 140 intervals for the day count). We coarse the grain of the grid with our hierarchical post-processing so that the informativity of the model is 80%, with ten clusters of antenna and twenty time segments. Since, missing values are abundant in this data, i.e., some antennas emitted no call during some time periods, consequently, we obtain time segments that are strongly correlated with missing data. For the same reason, antennas are grouped together because they experienced an absence of calls during one or several similar periods. In the Figure 5, the colored antennas belong to clusters having experienced simultaneous absences of calls. We observe that the green, orange, light blue and purple clusters are located in localized area. The missing data appear during short periods for these clusters. This grouping might be due to localized technical issues on the network. The antennas of the yellow cluster are spread over the country. These antennas are grouped because they have been activated at the same date. This use case provides a better understanding the dysfunctions in the network over the year.

Analysis of output communications w.r.t. week day and hour

Our objective is to build simultaneously a partition of the antennas, a partition of the week days and a discretization of the hour, i.e., a triclustering. For the same reasons as previously, we only keep the emitting antennas. At the finest level, we obtain a triclustering with 806 clusters of emitting antennas, 7 clusters of days and 22 time segments. Again, these results must be simplified to ease the interpretation. However, we fix the numbers of clusters of days and time segments, since they are acceptable for the analysis and we only reduce the number of clusters of antennas. With four clusters of antennas, we keep 51% of the informativity of the model.

Antennas are displayed on the map Fig. 6: Clusters on the map of Ivory Coast. Dots are antennas. There is one color per cluster.

of Figure 6. We also build a calendar (see Figure 7) for each cluster with days in columns and time segments in lines. The color of the cells indicates the excesses (red) or the lacks (blue) of traffic emitted from the corresponding cluster. The lacks and excesses are measured using the contribution to the mutual information (see definition 3) between the cluster and the cross product of the cluster of weekday and the time segment:

M I(X M 1 ; X M 2 × X M 3 )
, with X M 1 the partition of the antennas, X M 2 the partitions of the weekdays and X M 3 the discretization of the time. Now we focus on the analysis of each cluster of antennas that we can easily label manually: Abidjan -Le Plateau (yellow) This cluster covers exactly the Central Business District of Abidjan. In the calendar of Figure 7, we observe an excess of calls from the Monday to the Friday, between 8-9am and 4-5pm. The rest of the time, there is a low lack of traffic emitted from this area. In other words, during the office hours, the phone traffic is higher than expected and lower the rest of the time. This is expected and representative of this type of area: a non-residential business district.

Economic zones (red) The antennas of this cluster are located either in the commercial areas of the cities or in areas with a strong economic activity, like plantations or mines. In Abidjan, these antennas are located in industrial zones (South and North-West), the shopping districts (North of the business district) and the universities and embassies neighborhood (East). The traffic in these areas is mainly in excess from the Monday to the Saturday between 9 am and 5 pm. The correlation is very strong between the working hours and the calls traffic on these areas.

Urban residential areas (blue)

The antennas belonging to this cluster are mainly located in the cities like Abidjan, Bouaké and Yamoussoukro. If we focus on Abidjan, we observe that the cluster covers the residential neighborhood located in the West and in the North-East of the city. At a finer level of partition of the antennas, this cluster would be split according to the socioeconomic class of the neighborhood: the upper class neighborhood in the East of the city is separated from the lower class neighborhoods, located in the North and the West. The calendar shows lacks of calls during the office hours and excesses the weekend, the night and the early morning during the week. This is correlated with the presence of people in residential areas. Note that the excesses of calls start around 8 pm, while it stops around 5 pm in the Central Business district or in economic areas. This time lag is due to the cheaper price of calls after 8 pm.

The countryside (green) The antennas of this cluster are spread over the country, except in Abidjan and other cities in general. The calendar for this cluster is quite similar to the one of the urban residential areas, except that the excess periods are limited to the early evening and the whole Sunday. The selected group of users mainly connects to the antennas located in the East of Abidjan after 6 pm during the working days, while they rarely connect to the same antennas before 6 pm the same days. Then, it can be assumed that the selected cluster of users is composed people living in the same area. This hypothesis is reinforced by the socioeconomic nature of this part of Abidjan: it is a residential area. The contributions to mutual information of the other clusters of antennas are smaller. Three areas experience excesses of traffic before 6 pm and lacks after 6 pm. They correspond to the business district (Le Plateau), the embassies and universities neighborhood and the industrial zone located in the West of the city. The common feature of all these areas is their economic activity during the day. To sum up, we can assume that the users of the selected cluster are similar in that they live in the same area and work during the week in three localized area of Abidjan. Similar observations stand for several other clusters of users -thus we are able to summarize users' mobility behavior.

Conclusion

Motivated by two key points of economic development strategy of a telco in emerging countries, we have instantiated a generic methodology for exploratory analysis of CDRs data. Our method is based on a joint distribution estimation technique providing the user analyst with a summary of the data in a parameterfree way. We have also suggested several tools for exploring and exploiting the summary at various granularities and highlighting its relevant components. We have demonstrated the applicability of the method on graph data, temporal sequence data as well as user mobility data stemming from country-scale CDRs data. The results of the exploratory analysis are currently considered as valuable additional input to improve network planning strategy and pricing strategy.

Definition 3 ( 1 and Y M 2 (

 312 Contribution to mutual information). Given the K -2 selected parts c i3...i K , the mutual information between two partitioned variables Y M from the partition M of Y 1 and Y 2 variables induced by the grid model M ) is defined as:

Fig. 5 :

 5 Fig. 5: Antennas activity clusters projected on Ivory Coast map. Colored clusters show inactivity periods while grey clusters indicate antennas whose traffic is complete over the period.

Fig. 7 :

 7 Fig. 7: Calendars of excesses (red) and lacks of calls emitted from each of the four clusters of antennas, in function of the weekday and the daytime.

-

  Yield management pricing strategy: a part of the pricing policy, called Bonus Zone, established in Ivory Coast offers discount prices (from 10% to 90%) to calling users depending on the location and hour of the emitting call. Maps and calendars resulting from the last two case studies on temporal distribution of output calls (see Section 5.3) and on mobility data (see Section 5.4) that are available at various granularities, provide valuable information to economic analysts in order to design optimized spatio-temporal pricing policy in the context of Bonus Zone.

	i.e. approximatively 55 millions records. The records are described by the four
	following variables: anonymized user ID (50000 categorical values); connexion
	antenna (1214 categorical values); time of call (minute precision); date fo call
	(from 2012/12/12 to 2012/12/24).
	From this data, we consider the user trajectories (identified by user ID) inside
	the network for the following analysis:
	1. Analysis of user mobility w.r.t. week day and hour. We consider the user ID,
	antennas, week day and hour. This data set can be considered as a set of
	spatio-temporal footprints, where each user ID is associated with a sequence
	of antenna usage over the time dimension. Here again, the time dimension
	is represented by two variables and the data of the whole period is folded up
	to week day and hour.
	namely commu-
	nication data and mobility data:
	Case studies on communication data. Communication data consists in 471 mil-
	lions mobile calls and covers a 5-month period (from 2011, December 1st to 2012,
	April 28th). The records are described by the four following variables: emitting
	antenna (1214 categorical values); receiving antenna (1216 categorical values);
	time of call (with hour precision); date of call (from 2011/12/01 to 2012/04/28).
	From this data set, we consider three subsets for:
	1. Analysis of call network between antennas. Considering emitting antennas,
	receiving antennas and the calls made between antennas, the data set can
	be seen as a directed multigraph where nodes are antennas and links are the
	calls between antennas.
	2. Analysis of output traffic w.r.t. date of call. We consider emitting antennas
	and the number of days for each call from referral to first day of recording.
	This data set can be considered as a temporal event sequence spanning over
	the whole observation period, where the time is the number of days passed
	and the events are the emitting antenna IDs.
	3. Analysis of output traffic w.r.t. week day and hour of call. We consider emit-
	ting antennas, the day of the week (stemming from the date and considered
	as a numerical variable) and the hour of the day for each call. Here the time
	dimension is represented by two variables and the data of the whole period
	are folded up to week day and hour.
	Case studies on mobility data. Mobility data consists in mobility traces of 50000
	users over a 2-week period (from 2012 December 12th to 2012 December 24th),

Data description & studies. The CDRs data under study come from the Orange D4D challenge 4 (Data For Development

[START_REF] Blondel | Data for development: the D4D challenge on mobile phone data[END_REF]

). We consider several case studies on two anonymized CDRs data sets from Ivory Coast,

Romain Guigourès was with Orange Labs when this work began.

http://d4d.orange.com/en/home

User mobility analysis w.r.t. week day and hour

Among the 50000 anonymized users, we focus on mobile users characterized by a frequent use of a large set of distinct antennas: after filtering, 6894 users are under study. For these 4-d data (user, antenna, week day and hour), khc operates a tetra-clustering: as a result, users with the same mobility profile are grouped together, i.e., users who have connected to similar groups of antennas, on similar days of the weeks at similar time periods. At the finest grain, we obtain 237 clusters of users, 218 clusters of antennas and three time segments ,while week days remain as singletons. Again, the granularity prevent us from an easy interpretation, and we simplify the model. We keep 50% of informativity, that enables a reduction of the numbers of clusters of users and antennas to 40, and the numbers of groups of week days and hour segments to two. The week is divided in two parts: the working days and the weekend. For the hour dimension, the split occurs around 6 pm. The intervals are 0 am -6 pm and 6 pm -12 am. Note that the bound at midnight is artificial, because the day start as this time. The cut at 6 pm is the last in the hierarchy of the time segmentation. Then it would have been more relevant to consider a day from 6 pm to 6 pm the next day. Nevertheless, it is easier to have an interpretations on a "usual" time period between 0 am and 12 pm. Therefore we keep the following segmentation: 0 am -6 pm, 6 pm -12 pm. To illustrate the characterization of users' behaviors in terms of mobility provided by the grid, we focus on a group of users. The maps of Figure 8 shows the excesses and lacks of traffic in Abidjan during the week, for both periods of the day and for the selected group of users. The colors correspond to the mutual information M I(X M 1 ; X M 2 × X M 3 × X M 4 ) where X M 1 is the partition of antennas; X M 2 , the partition of the weekdays; X M 3 the discretization of the daytime; and X M 4 , the selected partition of the users.


Vishnupriya Raghavan
email: vishnupriyaraghavan@gmail.comdolphy.abraham

Dolphy M Abraham

Madhumita G Majumder
email: madhumita.gm@alliance.edu.in

Quality Improvements for Ensuring e-Retailing Success in India: Constructs and Frameworks

This extended abstract presents a review of various constructs and evaluation frameworks proposed in the literature for e-Retailers. Our study shows that the existing frameworks apply to either the technical or the nontechnical elements of an e-Retailing site and not both of them. Therefore, a comprehensive framework covering all aspects of quality is what is required. Furthermore, for the long run sustenance and growth of e-Retailing, it is necessary to focus on the service provided and not only the technical aspects. In the Indian context, where e-Retailing is beginning to enter a rapid growth phase, evaluation methods and metrics which are appropriate are necessary. Our extended abstract highlights the key issues that will help define these constructs

Introduction

E-Retailing has grown as a service sector globally. The e-Retailing market stood at $198.8 billion in 2011 in the United States, and €200.52 billion in Europe (Center for Retail Research, 2012). E-Retailing in India has started to gain momentum. The industry stood at a valuation of $948.6 million in 2011 and is expected to reach $260 billion by the year 2025 (First Data, 2011). Quality improvement in e-Retailing can ensure a robust growth of the sector and encourage standardization and improvement of overall business practices. For evaluating the measurement of quality improvements a broad approach suggested by Practical Software and Systems Management [START_REF] Statz | Measurement for Process Improvement[END_REF] is employed to understand the comprehensiveness of the constructs and frameworks available in literature.

Practical Software and Systems Management [START_REF] Statz | Measurement for Process Improvement[END_REF] identified four major areas of quality improvement derived from the balance score card approach. These areas are: • Financial (F): That deal with financial goals from the project.

•

Customer Satisfaction (CS): The satisfaction for both internal and external customers is important. These relate to things like mean time to failure, response time, price/performance etc.

•

Internal Business Processes (IBP): These goals relate to practices and methods for product and service development, management of people in the organization etc. 1 and2 (column 2) which lists the areas covered by the study in terms of the four suggested dimensions, that most of the measures focus on three aspects of quality improvement i.e., customer satisfaction, financial and internal business processes and learning and growth has been mostly ignored. Few researchers have addressed the e-Retailing quality dimensions in developing countries which can be applied to the Indian context. For example studies have focused on core service dimensions (Malhotra et al, 2004); system availability and fulfillment [START_REF] Kim | E-service Quality Perceptions: A Cross-Cultural Comparison of American and Korean Consumers[END_REF]); e-Commerce infrastructure [START_REF] Okoli | A Framework for Assessing E-Commerce in Sub-Saharan Africa[END_REF]; strategic alliance and innovative business strategies [START_REF] Li | A Holistic Framework of E-Business Strategy: The Case of Haier in China[END_REF]; cultural adaptation (Zahedi & Bansal, 2011;[START_REF] Singh | Measuring Cultural Adaptation on the Web: A Content Analytic Study of US and Japanese Web sites[END_REF]. Researchers suggest that developing countries have typical issues like infrastructure, technology acceptance etc. (Malhotra et al, 2004) and hence require different quality dimensions.

Conclusion and Future Research Opportunities

This research provides a list of measures for evaluating e-Retailing quality. It highlights the need for a comprehensive framework that encompasses the major aspects of e-Retailing quality i.e., customer service, financial, internal business processes, and learning and growth. However, this paper suffers from the limitation in that it lacks a proposed framework and an empirical analysis.

As e-Retailers expand their markets across the globe future research on e-Retailing quality improvements must also include country specific issues. Extensive studies are required to examine if country-specific websites reflect national cultural values (Zahedi & Bansal, 2011). Future research should focus on development of an evaluation framework for e-Retailing quality improvement that addresses technical, non-technical and country specific issues.

Table 1 .

 1 Important measures for e-Retailing quality improvement.Keeping the above framework in mind, we focus on measures for e-Retailing service quality improvement as suggested by researchers over the years. Table1and 2 provides a list of measures employed in various studies to evaluate an e-Retailer's quality. Table1lists important constructs and table 2 lists important frameworks identified by researchers. The last column of table1i.e. 'suggested measures' lists the measures that the researchers have suggested for future investigation.It can be observed from Tables

	•	Learning and Growth (L&G): This relates to people related capabilities of
		the organization like technical skills of staff, staff growth in terms of
		numbers, domain knowledge, morale and turnover.
	Author(s) & Year		Major Areas Existing Measures		Suggested Measures
	Baty & Lee (1995)		CS, IBP, F	Wide reach, multimedia, better	Semantic representation, advanced
					market functionality, information	navigation	and	product
					availability, shopping agents , better	comparability features
					search, reduced transaction costs
	Hawkins, Mansell,		IBP, F		Provision of wider product basket	Customer perception issues, develop
	& Steinmueller (1999)				(complementary) e.g. linking cash to	and harness virtual communities,
					smart cards			disintermediation	to	re-
								intermediation
	Srivastava & Mock (1999)	CS		Trust, business practices, transaction	Trust should be incorporated as an
					integrity, information protection,	essential construct
					legal environment
	Casati & Shan (2000)		IBP		Different vendors like SAP, ATG,	Integrated platform for a seamless
					Oracle			process like Enterprise Application
								Interface (EAI)
	Rust & Kannan (2003)		CS,IBP, F	Strategic, cost reduction, supply	Revenue expansion, information
					chain efficiency, brand equity	flows,	customer	profitability,
								customization
	Van der Merwe &		CS,IBP		Design,	content,	navigation,	Conversion of buyers decision
	Bekker (2003)				reliability			making process into technical
								aspects
	Croom & Johnston (2003)	CS,	IBP,	Cost	efficiency,	process	Greater	awareness	of	and
			L&G		conformance, internal customer	importance to internal customer
					satisfaction			satisfaction
	Burt & Sparks(2003)		CS,IBP, F	Cost reduction, low inventory,	ROI is not very clear for traditional
					branding			retailers
	Desai, Richards &		CS, IBP		Trust, process efficiency		Employment of EDI, and other
	Desai (2003)							standards that can iron out
								deficiencies causing trust issues
	Klischewski & Wetzel, (2003)	CS, IBP	XML based representation	customer orientation

Table 2 .

 2 Major frameworks for e-Retailing quality improvement

	Author(s) & Year	Major Areas	Title of Framework	Measures Employed
	Bressolles (2006)	CS, IBP	NetQual	Ease of site use, design, reliability, security
	Chiou, Lin & Perng	CS, IBP, F	4PsC	Place, product, price, promotion and customer
	(2010)			relation.
	Yoo & Donthu (2001)	CS, IBP	Sitequal	Ease of use, design, speed of order processing,
				security of financial information
	Barnes & Vidgen (2003) CS, IBP	Webqual	Usability, quality of information, quality of
				interaction
	Wolfinbarger & Gilly	CS, IBP	E TailQ	Site design, customer service, reliability, security,
	(2003)			privacy.
	Parasuraman, Zeithaml	CS, IBP, F	E S Qual/ E RecSQual Efficiency, fulfillment, system availability and
	& Malhotra (2005)			privacy. E RecS Qual is a measure for non-frequent
				customers
	2 e-

Retailing in the Development Context According

  to A. T.Kearney (2012), India's retail sector ranks 5 th in the world based on the global retail development index. This growth in e-Retailing is fuelled by increased broadband connectivity, rising living standards, busy lifestyle and traffic, much wider product range, and convenient processes. As the industry matures it would require quality improvement to keep up with the rising customer expectations and maturing industry standards. Indian e-Retailers can greatly benefit from studies on e-Retailing available from developed countries. Most e-Commerce research in India is based on adoption factors and employs qualitative measures without empirical analysis[START_REF] Vaithianathan | A Review of E-Commerce Literature on India and Research Agenda for the Future[END_REF]). e-Retailers like flipkart.com and snapdeal.com have benefited by following the business models and best practices of established e-Retailers like amazon.com and groupon.com.