Qiuxi Zhu

Francoise Sailhan

Md Yusuf Sarwar Uddin

Valérie Issarny

Nalini Venkatasubramanian

Multi-Sensor Calibration Planning in IoT-Enabled Smart Spaces

Emerging applications in smart cities and communities require massive IoT deployments using sensors/actuators (things) that can enhance citizens' quality of life and public safety. However, budget constraints often lead to limited instrumentation and/or the use of low-cost sensors that are subject to drift and bias. This raises concerns of robustness and accuracy of the decisions made on uncertain data. To enable effective decisionmaking while fully exploiting the potential of low-cost sensors, we propose to send mobile units (e.g., trained personnel) equipped with high-quality (more expensive) and freshly-calibrated reference sensors so as to carry out calibration in the field. We design and implement an efficient cooperative approach to solve the calibration planning problem, which aims at minimizing the cost of the recurring calibration of multiple sensor types in the long-term operation. We propose a two-phase solution that consists of a sensor selection phase that minimizes the average cost of recurring calibration, and a path planning phase that minimizes the travel cost of multiple calibrators which have load constraints. We provide fast and effective heuristics for both phases. We further build a prototype that facilitates the mapping of the deployment field and provides navigation guidance to mobile calibrators. Extensive use-case-driven simulations show our proposed approach significantly reduces the average cost compared to naïve approaches: up to 30% in a moderate-sized indoor case, and higher in outdoor cases depending on the scale.

I. INTRODUCTION

The advent of IoT ecosystems with low-cost sensors and actuators is enabling the widespread deployment of smart spaces in our homes, communities, and cities. Such smart spaces are a valuable source of knowledge and can drive applications to enhance public safety, personal health and community wellbeing [START_REF] Kortuem | Smart objects as building blocks for the internet of things[END_REF], [START_REF] Uckelmann | An Architectural Approach Towards the Future Internet of Things[END_REF], [START_REF] Vermesan | Internet of Things: Converging Technologies for Smart Environments and Integrated Ecosystems[END_REF]. The associated challenges in the realization of this vision are enormous -deployment/operation cost, big data, interoperability, data analytics, and resource-efficient networking, to name a few [START_REF] Atzori | The internet of things: A survey[END_REF], [START_REF] Stankovic | Research directions for the internet of things[END_REF] -and are on the research agenda of next-generation digital sciences. In this paper, we address cost-accuracy issues that arise in the deployment of affordable IoT systems. In particular, the aggregation of relevant knowledge at scale from low-cost smart spaces is problematic since low-cost sensing solutions imply low accuracy and faster degradation/drift. Recent experiences, including by the authors, show that the relative inaccuracy of the connected devices can be alleviated through the automated and in-situ calibration of sensors -made possible due to the linearity of the bias in most cases [START_REF] Stone | Calibration and linear regression analysis: A self-guided tutorial[END_REF]. For example, mobile platforms (including humans that carry calibration devices) can visit the smart spaces, either opportunistically or in a planned manner, to assess the biases of the deployed sensors and compensate for errors [START_REF] Nesamony | On optimal route of a calibrating mobile sink in a wireless sensor network[END_REF]. One can thus generate "sufficiently accurate" knowledge over time through the frequent calibration of the low-cost sensors in the field; however, the low deployment cost might result in increased maintenance costs! Careful planning of the calibration process is therefore essential for the costeffective monitoring of the smart spaces -this is increasingly important as the number and size of smart spaces grow.

In this paper, we address the calibration planning problem. The aim is to develop a plan for the calibration of a large number of inexpensive (and often inaccurate) sensors in a smart space using high-integrity reference sensors that are mobile, such that (a) the deployment and operational costs for calibration are minimized while (b) maintaining a sufficient observation accuracy from the sensor measurements. More realistically, given the knowledge of a sensor's degradation characteristics, we program its calibration with respect to an observed phenomenon so as to maintain an adequate sensing accuracy while minimizing the required effort from the mobile calibrators. We formalize the above as a multi-path planning problem, which we solve via intelligent heuristics, and validate using real-world settings. Our key contributions include:

• A measurement study to help understand the calibration issue of low-cost sensors for environmental monitoring, which motivates the proposed approach for the scalable calibration of IoT-enabled smart spaces (Section II). • Long-term multi-sensor calibration planning as a service that exploits device locality, sensor characteristics, and application needs. • The characterization and formulation of the multi-sensor calibration planning problem (Section III) and discussion of NP-hardness. • A two-phase iterative solution and a family of heuristic methods to enable the cost-effective planning of multisensor calibration in large smart spaces and over longer time periods (Section IV). • The validation of our approach and algorithms leveraging real-world indoor and outdoor smart spaces settings from our ongoing testbeds in Irvine, CA and Paris, France; results show significant cost improvement while maintaining adequate accuracy. • Initial steps towards a prototype of a calibration planning service for IoT smart spaces with: (i) a dashboard to map the deployment of sensor nodes over large spaces, and (ii) a mobile app that the mobile calibrators use along their calibration journeys (Section V). Our first set of experiences in calibration comes from the launch of an urban scale experiment with the city of Paris in 2015 1 . As part of our efforts, we developed a crowdsensing application for mobile phones for monitoring the exposure of the urban population to environmental noise pollution [START_REF] Hachem | Monitoring Noise Pollution Using The Urban Civics Middleware[END_REF]. While crowdsensing may be a cost-effective approach for monitoring urban environmental conditions (e.g., noise) that exhibit high spatio-temporal variability, the relatively low quality of the embedded sensors as well as the uncertain sensing context required dedicated actions. We quickly learned (as anticipated) that only a low percentage of the crowdsensed measurements actually contribute to the analysis of the urban noise pollution. We thus needed to enhance the quality of the observations [START_REF] Issarny | Dos and Don'ts in Mobile Phone Sensing Middleware: Learning from a Large-Scale Experiment[END_REF]. As a first step, we thoroughly studied the bias in smartphone noise sensing against a reference sound level meter, which helped us develop a protocol for calibrating individual handsets [START_REF] Ventura | Evaluation and calibration of mobile phones for noise monitoring application[END_REF]. Through experimental studies, we demonstrated that a calibration protocol can help gather observations that can more accurately map noise pollution at the district level [START_REF]Assimilation of mobile phone measurements for noise mapping of a neighborhood[END_REF]. Unfortunately, the proposed calibration protocol places a high demand on the end-users who must actively participate in the calibration. We subsequently investigated a distributed protocol for the opportunistic multi-party calibration of devices located in the same sensing and communication range [START_REF] Sailhan | Opportunistic multiparty calibration for robust participatory sensing[END_REF].

The second experience in calibration comes from lessons learned while deploying IoT-enabled environmental sensing in the SCALE project [START_REF] Benson | SCALE: Safe community awareness and alerting leveraging the internet of things[END_REF], [START_REF] Uddin | The scale2 multi-network architecture for iot-based resilient communities[END_REF]; here, we deployed inexpensive multi-sensor platforms (called SCALE boxes) at multiple locations worldwide including Irvine, CA; Montgomery County, MD [START_REF] Benson | SCALE: Safe community awareness and alerting leveraging the internet of things[END_REF], and Dhaka, Bangladesh [START_REF] Rahman | Adaptive sensing using internet-of-things with constrained communications[END_REF], [START_REF] Rahman | An adaptive IoT platform on budgeted 3g data plans[END_REF]. A SCALE node is a Raspberry-Pi-based multi-sensor box with middleware (SCALE client) that provides flexible interfaces for data collection for a wide range of sensor types (gas, light, air quality, 1 www.inria.fr/en/centre/paris/news/launch-of-soundcity-mobile-application temperature, seismic, camera, etc.). The SCALE client utilizes a publish-subscribe approach to organize raw data and detected events into topics that are published to a cloud (or edge) data exchange service. IoT applications subscribe to relevant topics and receive updates. The deployments allow us to obtain rich information about the sensor behaviors and especially their respective deterioration across time in various environments. In all cases, we observed a deterioration in the response of the low-cost gas sensors over time.

At UC Irvine, SCALE was used in a smart instrumented building for everyday monitoring [START_REF] Mehrotra | Tippers: A privacy cognizant iot environment[END_REF] -explosive gas sensors, useful in emergency response, were occasionally triggered externally. As part of our testbed, various sensors and monitoring devices (e.g., microphones, surveillance cameras, power meters, BLE beacons -for the location tracking of participants) were deployed in multiple floors of the CS department building. To validate our observations and generate a custom model of sensing deterioration, we developed a measurement platform to reproduce the deterioration process. Here, the custom sensor box was instrumented with 4 pairs of low-cost sensors as a "test node" (Figure 1a). Figure 1b shows the data reported by two pairs of MQ gas sensors 3-4 months after the installment: long-running MQ-2 sensors that were connected to a running node for more than two months generated different values at the same location and show significant difference in sensitivity when exposed to similar stimuli. Such difference has also been observed during/after deployment campaigns [START_REF] Barrenetxea | The hitchhiker's guide to successful wireless sensor network deployment[END_REF], [START_REF] Ni | Sensor network data fault types[END_REF]. State of the Art Approaches: The calibration of a sensor lies in (i) generating a controlled stimuli x and (ii) analyzing the sensor response y. Then, the true (generated) signal x, which serves as ground truth, is compared with the sensor response y so as to obtain the mapping between the sensor response and the true value; the calibration parameters are accordingly adjusted. The calibration is typically performed once in a laboratory so as to avoid the addition of any environmental interference/noise to the high-fidelity/true signal. However, the calibration does not take into account the context wherein the sensor evolves, even though environmental factors, e.g., solar radiation, strongly influence the readings provided by sensors (e.g., air temperature). Recent work focuses on automatically calibrating sensors in the field, without relying on a controlled stimuli and without a well-defined range of conditions. This is commonly referred to as blind calibration and it consists in calibrating a (uncalibrated) sensor using the measurements provided by a neighboring (calibrated) sensor. There are two approaches to the blind calibration of a set of sensors depending on whether the sensors are static or mobile. In the case of (static) WSNs, the common assumption is that the deployment is sufficiently dense so that neighboring sensors provide nearly identical measurements [START_REF] Stankovic | Asynchronous distributed blind calibration of sensor networks under noisy measurements[END_REF], [START_REF] Wang | A deep learning approach for blind drift calibration of sensor networks[END_REF]. Then, the inaccurate readings of an uncalibrated sensor can be compensated using the measurements provided by a nearby (calibrated) sensor. In the case of mobile sensors, multi-hop calibration [START_REF] Hasenfratz | On-the-fly calibration of low-cost gas sensors[END_REF], [START_REF] Saukh | Reducing multi-hop calibration errors in large-scale mobile sensor networks[END_REF], [START_REF] Fu | Multihop calibration for mobile sensing: K-hop calibratability and reference sensor deployment[END_REF] allows mobile sensors to get calibrated "indirectly" using other mobile sensors that have been calibrated by the reference sensor. Multi-party calibration [START_REF] Sailhan | Opportunistic multiparty calibration for robust participatory sensing[END_REF] derives regression models to solve the calibration problem when there are multiple (more than two) participants. Recent work also investigates planning problems, e.g., [START_REF] Fu | Multihop calibration for mobile sensing: K-hop calibratability and reference sensor deployment[END_REF] explores the in-situ placement of reference sensors in the field; [START_REF] Nesamony | On optimal route of a calibrating mobile sink in a wireless sensor network[END_REF] proposes a TSP based path planning algorithm for a mobile calibrator.

In contrast to the above efforts, we take a more holistic approach to low-cost calibration in smart spaces at scale. First, we consider the presence of heterogeneous sensor types with varying calibration characteristics. Our proposed approach is application-aware and is able to take into account diverse sensing needs (sensor type, sensing accuracy) presented by the context at hand. We assume the ability to utilize multiple mobile calibrators. We exploit the locality of in-situ sensors (that can calibrate each other) to reduce the cost of mobile calibrators. The multi-step approach, presented next, highlights the overall scheme for a more comprehensive and efficient calibration planning mechanism for IoT-enabled smart spaces.

III. MULTI-SENSOR CALIBRATION PLANNING

Our objective is to enable the cost-efficient calibration planning of a given smart space with multiple types of sensors and diverse applications. The target environment assumes the in-situ deployment of a large number of inexpensive sensors with an adequate number of high-precision mobile calibrators (e.g., trained workers that calibrate the in-situ devices). We propose an effective approach (Figure 2) to the planning of the sensor calibration over a given -possibly long -maintenance period T : (i) First, we partition the space to form calibration spots by exploiting the locality of in-situ sensors; a sensor can calibrate other sensors of the same type at the same spot. (ii) Next, we carry out multiple iterations of calibration during the given maintenance period T . Our challenge lies in determining the number of iterations and the time at which each iteration is executed. Furthermore, in each iteration we determine: the sensors to calibrate, the number of mobile calibrators needed, and the paths taken by each mobile calibrator. The sensor calibration is then carried out accordingly over a number of Ω iterations so that the overall cost over T is kept to a minimum and ensures that all the sensors always comply with the data accuracy requirements. The planning problem then amounts to identifying the least number of iterations and planning the calibration at each iteration ω=1, 2, . . . , Ω so that the number of workers involved is also kept to a minimum. Next, we introduce the notations and assumptions underlying our formalization of the multi-sensor calibration planning problem.

A. Notations and Assumptions

A node n j , j=1, 2, . . . , N , is an IoT device that embeds one or more types of low-cost sensors. A sensor type s k , k=1, 2, . . . , K, refers to the capability of detecting a certain type of phenomenon (e.g., temperature, gas concentration), which usually requires a specific kind (or combination) of lowcost sensor(s). Then, we denote an individual sensor by n j,k , and we introduce the binary sensor presence matrix Q N ×K to characterize the set of available sensors so that Q j,k =1 (resp. 0) if sensor type s k is present (resp. absent) on node n j .

1) Calibration-Related Terms: A reference sensor is a high-quality sensor, whose readings serve as the ground truth and therefore can be used to calibrate other sensors that sense the same phenomenon. We assume that the reference sensors are calibrated offline (e.g., in a lab before each iteration).

A mobile calibrator (or simply a calibrator) m i is a person who carries reference sensors and visits the field to calibrate the deployed low-cost sensors. For simplicity, we assume that a calibrator can be equipped with any of the reference sensors needed across his/her journey. Sensor calibration takes place when a calibrator m i visits a node n j and stays at the spot for long enough to calibrate a sensor n j,k .

Each sensor type s k is associated with a calibration time τ k that ranges from a few seconds to several minutes, depending on the phenomenon detected by the type of sensor and the calibration complexity. We associate each sensor type s k with a calibration period T k , which characterizes the maximum duration, during which the sensors remain valid (i.e. the measurements have sufficient accuracy) once calibrated. The period depends on the usage scenario and may be learned from empirical study.

During operation, each individual sensor n j,k is associated with a time to next calibration (TTNC) F j,k that indicates how soon the sensor needs to be (re-)calibrated. The TTNC matrix F N ×K then represents the TTNCs of all the sensors. F is a function of time, where each F j,k decreases between iterations and is reset to T k when n j,k is calibrated. We further denote F[ω -] (resp. F[ω + ]), the matrix TTNC immediately before (resp. after) the iteration ω. Note that F is non-negative, i.e. F j,k [ω ± ] 0, ∀(j, k), ω∈N + . If node n j does not hold a sensor of type s k , the corresponding TTNC is infinite, i.e.

F j,k [ω ± ]=+∞, ∀ω∈N + , if Q j,k =0.
A sensor selection is a collection of sensors (selected for calibration) represented by a binary matrix Γ N ×K [ω], where Γ j,k Q j,k . If sensor n j,k is selected for calibration at iteration ω (i.e. Γ j,k [ω]=1), then, at the end of this iteration, its TTNC is reset to its calibration period T k , i.e. F j,k [ω + ]=T k ; otherwise, its TTNC stays unchanged (during iteration ω):

F[ω + ] = Γ[ω] • T N + (1 -Γ[ω]) • F[ω -] (1)
Where T N is the nodal calibration period matrix that consists of

N identical rows of [T 1 , T 2 , • • • , T K ]; "•"
F[(ω+1) -]=F[ω + ]-min F[ω + ].
2) Smart-Space-Related Terms: We abstract a smart space as a set of spots that are such that sensor nodes deployed at a spot ν l , l=1, 2, . . . , L, are sufficiently co-located to enable their concurrent calibration by a single calibrator. Note that this may possibly involve leveraging multi-party multi-hop calibration (Section II). We then map a smart space as a directed graph G=(V, E), where each vertex ν l ∈V corresponds to a spot and each edge weight (ν l1 , ν l2 ) denotes the average time taken by a calibrator to move from ν l1 to ν l2 . The graph can be represented as an adjacency matrix G L×L , where G l1,l2 is the weight on the directed edge (ν l1 , ν l2 ) (i.e. the movement cost). As long as each spot is physically accessible, G is a complete digraph and all its edges have finite weights.

The node locations are represented by a binary location matrix D N ×L , where each row represents a node and each column a spot. We set D j,l =1 if node n j is deployed at spot ν l ; or 0 otherwise. A node is deployed at a single spot, so L l=1 D j,l =1. For simplicity, we assume that all the calibrators depart from the same spot, which is referred to as the "depot" in related work on path planning, and is conventionally indexed as the first spot, i.e. ν 1 . We assume that no node is deployed at the depot. Given the node locations D, we can derive the spot selection vector h L from the sensor selection Γ. A spot is selected for iteration ω if any sensor on any node deployed at that spot is selected in Γ[ω] i.e.:

h l [ω] = N j=1 K k=1 Γ j,k [ω]•D j,l (2)
In each iteration, the path of a calibrator is an ordered sequence that starts from the depot and visits a set of nonrepeating selected spots. It can be represented as a binary matrix W L×L , where W l1,l2 =1 if the calibrator visits spot ν l2 immediately after visiting ν l1 ; or 0 otherwise. We require that all mobile calibrators return to depot in the end, i.e.

L l=1 W 1,l = L l=1 W l,1 . The path of calibrator m i in iteration ω is denoted W i [ω]
. Each selected spot is visited exactly once by one calibrator:

M i=1 L l=1 W l,l0,i [ω] = h l0 , l 0 = 2, 3, . . . , L (3)
The other constraints that the matrices {W} should satisfy are discussed in Section IV-B.

B. Definition of the Calibration Cost

We consider three major types of cost: iteration overhead, movement cost, and calibration cost. The iteration overhead C it is the cost of all the activities related to an iteration that are not tied to any specific calibrator, such as preparation, equipment, and the transport to the deployment, etc.

The movement cost C w corresponds to the cost associated with the travel time of the calibrators while moving between spots. The movement time C w,i [ω] of a single calibrator m i in iteration ω can be computed from the map G and the calibrator's path W i [ω]. It equals the sum of the weights on the edges between all the consecutive pairs of spots visited by the calibrator:

C w,i [ω] = L l1=1 L l2=1 W l1,l2,i [ω]•G l1,l2 (4)
The calibration cost C c reflects the time and effort it takes to conduct sensor calibration while staying at the spots. The cost C c,i [ω] of a specific calibrator m i in iteration ω can be computed based on a given selection of sensors Γ[ω]. We indeed know the calibration protocol (and thus duration) that needs to performed per sensor type. We further assume that the calibration that happens at the same spot is done in parallel. Thus, the calibration time that m i spends at spot ν l equals the maximum τ k of all selected sensors at that spot (i.e. Γ j,k =1 and D j,l =1). Then C c,i [ω] equals the sum of the calibration times at all the spots assigned to m i :

C c,i [ω] = L l=1 max j,k D j,l •Γ j,k [ω]•τ k • L l =1 W l,l ,i [ω] (5)
The work load of any calibrator m i at iteration ω is the sum of his/her movement and calibration time:

C i [ω] = C w,i [ω] + C c,i [ω] (6)
For any iteration, we assume that the maximum work load of any calibrator is ĉ, i.e. C i [ω] ĉ, ∀(i, ω).

The total cost of any iteration ω, denoted C[ω], is the weighted sum of: (a) C it , the constant iteration overhead, (b) C w [ω], the total movement time of all mobile calibrators, and (c) C c [ω], the total calibration time at selected sensors:

C[ω] = µ 0 •C it + µ w •C w [ω] + µ c •C c [ω] (7)
where

C w [ω]= i C w,i [ω] and C c [ω]= i C c,i [ω].
Since we assume that each spot is only visited once by one calibrator, the total calibration cost is computed directly from Γ[ω], i.e.:

C c [ω] = L l=1 max j,k D j,l •Γ j,k [ω]•τ k (8)

C. Multi-Sensor Calibration Planning Problem

We now introduce the multi-sensor calibration planning problem to minimize the average cost of operation over the maintenance period T . The problem is formulated as follows: Given the time span T , the map G, the location matrix D and the sensor presence matrix Q of all the nodes, the calibration time τ k and the calibration period T k of all the sensor types, and the initial TTNC matrix F[1 -]; find the total number of iterations Ω, and for each iteration ω=1, 2, ..., Ω, find the time t ω it takes place, the sensor selection Γ[ω], and the number and the paths of calibrations {W[ω]}; such that the average cost of all iterations over the time span T is minimized: where C[ω] is the total cost of iteration ω given by Equation [START_REF] Nesamony | On optimal route of a calibrating mobile sink in a wireless sensor network[END_REF], which depends on the sensor selection Γ and the calibrators' paths {W[ω]}, i.e. C[ω]=C(Γ[ω], {W[ω]}); obviously, it also depends on problem inputs (i.e. G, D, Q, etc.) which are hidden for cleaner expressions. The unnumbered constraints above are related to the definition of sensor selection and TTNC in Section III-A. Constraint [START_REF] Ventura | Evaluation and calibration of mobile phones for noise monitoring application[END_REF] says the iterations need to cover the entire time span of T ; [START_REF]Assimilation of mobile phone measurements for noise mapping of a neighborhood[END_REF] and ( 12) make sure all the spots with selected sensors in Γ are visited in {W}; [START_REF] Benson | SCALE: Safe community awareness and alerting leveraging the internet of things[END_REF] says no calibrator should work for longer than ĉ in any iteration. Additional constraints apply to ensure {W} are valid path(s).

min 1 T • Ω ω=1 C(Γ[ω], {W[ω]}) (9)
s.t. t 1 = 0 Γ j,k [ω] ∈ {0, 1}, ∀ω, ∀j, ∀k W l1,l2,i [ω] ∈ {0, 1}, ∀ω, ∀i, ∀(l 1 , l 2 ) F j,k [ω -] 0, ∀ω, ∀j, ∀k F j,k [ω + ] > 0, ∀ω, ∀j, ∀k Γ j,k [ω] Q j,k , ∀ω, ∀j, ∀k F[ω + ] = Γ[ω] • T N + (1 -Γ[ω]) • F[ω -], ∀ω F[(ω+1) -] = F[ω + ]-min F[ω + ], ∀ω t ω+1 = t ω + min F[ω + ], ∀ω t Ω + min F[Ω + ] T (10)
h l [ω] = N j=1 K k=1 Γ j,k [ω] • D j,l , ∀ω, ∀l (11)
M i=1 L l=1 W l,l0,i [ω] = h l0 , l 0 =2, 3, . . . , L, ∀ω (12)
C i (Γ[ω], W i [ω]) ĉ, ∀ω, ∀i (13)
Algorithm 1: TTNI-driven single-iteration local optimization algorithm for the sensor selection planning problem.

1 function solveSSPSingle (G, D, τ , T, Q, F, µ, ĉ) ; Input : τ -[τ k ]; T -[T k ] k=1, . . . , K; µ -{µ 0 , µ w , µ c };
Refer to Section III-A, III-C for other symbols. Output: Γ -Sensor selection matrix.

2 tMin ← min{F j,k | Q j,k =1∧F j,k >0} ; tMax ← min T k ; 3 tCandSet ← {F j,k | tMin F j,k tMax} ; 4 Initialize minCostAvg ← +∞ ; minGamma ← null ; 5 for each T cand in tCandSet do 6 Γ ← 0 N ×K ; 7 for j in 1, . . . , N ; k in 1, . . . , K do 8 Γ j,k ← (Q j,k =1∧F j,k <T cand ; 9 for j' in 1, . . . , N ; k' in 1, . . . , K do 10 if D j,l =D j ,l , ∀l then 11 Γ j ,k ← (Q j ,k =1∧τ k τ k ) ; 12 H ← {l | ∃(j, k) s.t. Γ j,k =1∧D j,l =1} ; β ← 0 L ; 13 for l in 1, . . . , L do β l ← max j,k (D j,l •Γ j,k •τ k ) ; 14 cost ← µ 0 •C it + µ c •C c (D, Γ, τ ) + µ w •C w (G, solveMPPGreedy(G, ĉ, H, β)) ; 15 if (costAvg ← cost/T cand ) < minCostAvg then 16 minCostAvg ← costAvg ; minGamma ← Γ ; 17 return Γ ← minGamma ;
The sensor calibration planning problem is NP-hard. It tries to minimize the total cost of all iterations while the choices of early iterations can affect and limit the choices of later ones. Also, the cost of each iteration C[ω] involves a movement time C w [ω], which also needs to be minimized, and thus requires an optimization on the paths of the calibrators, which is a variant of the Multiple Travelling Salesman Problem (mTSP) that is known to be NP-hard.

IV. SOLUTIONS AND ALGORITHMS

Our formulation in Equations (9-13) suggests we find Γ[ω] (sensor selection) and {W[ω]} (path plan) simultaneously for all iterations. However, we observe the fact that if we know which spots the calibrators need to visit, we can optimize the paths to visit them accordingly. Hence, instead of attempting to minimize C/T , for each iteration ω, we attempt a twophase local optimization on the single-iteration average cost, C[ω]/(t ω+1 -t ω ), where we decouple the optimization of Γ[ω] and {W[ω]}. Accordingly, for each iteration we have a sensor selection planning phase and a multi-path planning phase. In the selection planning phase, given the initial TTNC matrix F[ω -], we optimize the sensor selection Γ, from which we derive the set of selected spots H={h l | h l =1, ∀l}, which is then used in the path planning phase to decide the number of calibrators and the optimal path(s) to visit the selected spots.

A. Sensor Selection Planning Algorithms

Leveraging the discrete nature of TTNC and the definition of TTNI (time to next iteration), we propose the TTNIdriven local optimization algorithm. The intuition behind this algorithm is to exhaust the possible values of TTNI (i.e. t ω+1 -t ω ) and find the "cheapest" one to fulfill.

The procedure of the TTNI-driven local optimization is shown in Algorithm 1. It involves the following steps: (1) Determine all the possible values of TTNI that could result from any possible sensor selection in this iteration. The minimum TTNI candidate is min{F

j,k [ω -] | Q j,k =1 ∧ F j,k [ω -]>0},
selecting only the sensors that need immediate calibration (Ln 2). The maximum TTNI candidate is min T k , selecting all sensors (Ln 2). All values in F[ω -] between them become TTNI candidates (Ln 3). In the worst case, the number of TTNI candidates is O(N •K) (2) For each TTNI candidate T cand , tentatively assume it to be the desired TTNI and create the minimum selection of sensors to meet the TTNI, i.e. let Γ j,k =1 if Q j,k =1 and F j,k [ω -]<T cand (Ln 5-8); then add all the sensors that are co-located with the selected sensors and that do not induce extra time for calibration (Ln 9-11, because their calibration is done in parallel, if it takes a shorter time). Generating Γ from T cand takes O(N •K+N •L) time. Compute the single-iteration average cost from Γ (Ln 12-14). ( 3) Select the TTNI candidate that gives the minimum average cost (Ln 16), and its corresponding Γ is the output of the algorithm. The worst-case running time excluding the time used to compute or estimate the movement time, is

O(N 2 •K 2 +N 2 •K•L).
If during step [START_REF] Uckelmann | An Architectural Approach Towards the Future Internet of Things[END_REF] we are able to compute the optimal paths of calibrators, we will compute the best cost evaluation for each selection and find the local optima. Unfortunately, multi-path planning is also NP-hard. We then propose two heuristics: a fast nearest-neighbor-based greedy algorithm, and an improved genetic algorithm (GA). During sensor selection, we use the faster greedy algorithm to estimate the movement cost; once the selection is done, we use GA to generate the final path(s) for the iteration.

B. Multiple-Path Planning Algorithms

The multi-path planning problem for a specific iteration ω refers to a sub-problem in our two-phase local optimization solution to the sensor calibration planning problem. The objective is to generate a set of paths {W[ω]} of minimum cost (i.e. movement time C w [ω]) for the selected spots yielded by the sensor selection Γ[ω]. It is a variant of the classic mTSP or VRP: we determine the number of calibrators based on the demand instead of having the number m of travellers given, as in mTSP. Also, evaluating the calibrator workload constraint involves the movement time of individual calibrators, which adds to the complexity of solutions.

Hence, we derive the following mixed-integer-programming (MIP) formulation of the multi-path planning problem based on a flow-based three-index MIP formulation of mTSP [START_REF] Bektas | The multiple traveling salesman problem: an overview of formulations and solution procedures[END_REF], adding appropriate modifications to match our assumptions and constraints: Given a map G, the location of the nodes D, the sensor selection Γ, and the calibration time τ k , ∀k; find W L×L×M (represented in a more general form that could contain calibrators with no assignment) and helper variables

U L×M to min L l1=1 L l2=1 G l1,l2 • M i=1 W l1,l2,i (14) s
.t. W l1,l2,i ∈ {0, 1}, ∀i, ∀(l 1 , l 2 ) W l,l,i = 0, ∀l=2, 3, . . . , L, ∀i L l2=1 W 1,l2,i = 1, ∀i L l1=1 W l1,l,i - L l2=1 W l,l2,i = 0, ∀i, ∀l L l1=1 L i=1 W l1,l2,i = h l , ∀l 2 (15)
u l,i 2, ∀i, ∀l

u l1,i -u l1,i + 1 -(L -1)•(1 -W l1,l2,i ) 0, ∀i, ∀(l 1 , l 2 ) L l1=1 L l2=1 W l1,l2,i • G l1,l2 + Υ l2 ĉ, ∀i (16)
Where W l1,l2,i =1 if calibrator m i visits spot ν l2 immediately after spot ν l1 ; or 0 otherwise. M is the maximum number of calibrators; assuming we always have enough calibrators, L would be an effective upper bound of M to be used in solvers.

The unnumbered constraints are related to the construction of multiple paths. Constraint ( 15) makes sure all selected spots are visited by exactly one calibrator; ( 16) enforces the maximum workload of calibrators, where Υ l is the total calibration time spent at spot ν l , i.e.

Υ l [ω] = max j,k (D j,l •Γ j,k [ω]•τ k ).
This formulation allows us to apply MIP solvers directly. However, the problem is NP-hard; the number of independent variables and the number of constraints in this MIP formulation are both in the order of O(L 3 ), resulting in a huge solution space. It is hard for any MIP solver to optimally solve the problem in a reasonable amount of time [START_REF] Bektas | The multiple traveling salesman problem: an overview of formulations and solution procedures[END_REF]. In particular, we tried two widely used solvers: GLPK (GNU Linear Programming Kit, open-source -https://www.gnu.org/software/glpk/) and Gurobi (commercial software -http://www.gurobi.com/). None of the two solved the problem in less than 48 hours for L 15. To solve the problem at larger scale, we propose two heuristics: a greedy algorithm derived from the nearest neighbor heuristic of traditional TSP, and an improved genetic algorithm (GA) based on the one proposed by Sedighpour, et al. [START_REF] Sedighpour | An effective genetic algorithm for solving the multiple traveling salesman problem[END_REF] for mTSP. For a clean design of the algorithms, after the completion of the sensor selection planning phase, the planning framework computes the set of selected spots H and the calibration time β at these spots from the sensor selection matrix Γ.

1) Nearest-Neighbor-Based Greedy Heuristic: The nearest neighbor algorithm for TSP starts with a tour containing only one spot. At each step, it determines that the next spot to visit as the one that is closest to the last visited spot, and loops until all the spots are visited.

Inspired by this straightforward TSP algorithm, we derive our greedy algorithm for the multi-path planning problem shown in Algorithm 2 as follows: (1) Start with a set of empty paths (i.e. all the calibrators stay at the depot) and the set of all selected spots H. (2) At each step, for every unvisited spot, compute the extra travel and calibration time yield by adding it to the end of the path of each mobile calibrator as long as Algorithm 2: Nearest-neighbor-based greedy algorithm for the multi-path planning problem. 3) Loop until all the spots are visited (Ln 3,[START_REF] Ni | Sensor network data fault types[END_REF]. Note: Our actual implementation of this algorithm caches the movement and calibration time associated with each calibrator to reduce redundant computation, so the worst-case running time of this algorithm is O(L 3 ).

2) Improved Genetic Algorithm (GA): We design our genetic algorithm (GA) based on the mTSP GA solution of [START_REF] Sedighpour | An effective genetic algorithm for solving the multiple traveling salesman problem[END_REF]. Features are added to address the peculiarities of our MPP formulation, i.e. the variable number of calibrators, the workload constraint, and the map (i.e. a directed-graph).

A chromosome is an integer vector that is made of two parts: a permutation of all the selected spots (1st half) and an assignment mapping the spots to mobile calibrators (2nd half). If the number of selected spots is |H|=L , a chromosome will have length 2L . The assignment (2nd half) is represented by the number of spots visited by each calibrator, so these integers should all be in range [0, L ] and sum up to L . For example, chromosome [2, 4, 5, 6, 3, 2, 3, 0, 0, 0] means L =5 and that m 1 will visit spots ν 2 , ν 4 , and m 2 will visit ν 5 , ν 6 , ν 3 . The fitness is the negative of the total movement time of all the mobile calibrators, and the selection is done by a standard scaled-fitness proportional selection.

The initial population is composed of randomly generated individuals. The permutation is performed by a uniformly random permutation generator, and the assignment is done by uniformly and randomly picking an integer and subtracting it from the total number of selected spots until none is left. The crossover is done by applying a standard "order crossover" on the first half of the chromosome.

Because of the variable number of mobile calibrators, we design three helper functions that apply to chromosomes: (a) compress: shift all zeros in the assignment section to the end and non-zeros values to the beginning; (b) split: check if any assignment ( 2 spots) leads to an overloaded calibrator, randomly split it into two calibrators, and loop until none is found; (c) merge: check if there exists a pair of assignments that can be merged into one without overloading the calibrator; then merge the first pair found. Among the three, compress and split are applied to every newly-generated chromosome during population initialization, mutation, and crossover, while merge is applied as one type of mutation.

Apart from "merge", there are three other types of mutation: (a) two-point swap, (b) segment reversal, and (c) 3opt local optimization. Every time a mutation is triggered, we randomly pick one of the four types of mutation functions. (a) and (b) are straightforward. 3-opt [START_REF] Boese | Cost versus distance in the traveling salesman problem[END_REF] is a local optimization for TSP, which tries to break a tour into three segments by removing three edges, and reconnect the three segments into a new but shorter tour. 2-opt is a commonly used local optimization for TSP on undirected graphs, but an odd numbered opt is required for digraphs to avoid reversing any segment, which makes it faster to compute the new movement time.

Finally, the tunable parameters such as the population size, elite-keeping size, and the termination conditions are assigned by the framework according to the problem size (i.e. L ).

V. VALIDATION

We evaluate the performance of our proposed multi-sensor calibration planner using realistic data derived from testbeds and present the steps we are taking towards a usable system compares the performance of multi-path planning algorithms for the number of spots L 12: GA provides close-to-optimal solutions but takes 8-10 sec for L=12 (20-60 sec for L=60); the greedy heuristic is much faster (approx. 0.5 sec for L=60), which makes it suitable as an estimator during the sensor selection planning. GLPK and Gurobi (i.e. the MIP solvers) could not terminate within 48 hours for L 15 so we aborted.

Towards a Usable System: Creating a calibration planning service is not trivial. In addition to the algorithm for planning, we are implementing a service to enable flexible calibration. To support usability, we have implemented a toolkit that serves to generate 3-dimensional maps and provides navigation guidance to the mobile calibrators. We create the 3-D navigable maps using OpenStreetMap (www.openstreetmap.org), which provides a basemap of the outdoor environment, e.g., roads (lines), junctions (points), buildings (polygons). The basemap is further edited (Figure 8a) so as to (i) detail the inner structure of the buildings of interest (if any), and (ii) layer a number of features including the paths accessible in the indoor environment. To design the inner structure of the building, we rely on JOSM (josm.openstreetmap.de). Following, we extract a representation of the connection graph from the resulting map using GraphHoper (www.graphhopper.com). This graph is composed of edges corresponding to streets, roads as well as indoor-pedestrian paths, and nodes representing junctions.

The next step consists in placing on the map, the spots that should be visited along with the nodes and the related sensors.

The wireless connections among sensing nodes are also modeled to inform the propagation of calibration parameters. The map are exploited by the Android app that we implemented for assisting the mobile workers. The app computes the shortest path (avoiding obstacles) between any 2 spots and provides this information to the multi-sensor calibration planner. For instance, with 60 spots and 77 nodes equipped by 188 sensors mapped over an urban area of 10 km 2 , the resulting connection graph is generated in 2.967s. This graph, which contains 8768 edges (e.g., streets) and 5661 vertices (e.g., junctions), is further used to compute the paths between any two spots; such a computation takes in average 0.0195s for an average distance separating two spots of 1387 meters. Relying on the multisensor calibration planner, the mobile app further provides navigation guidance to the calibrators (Figure 8b).

VI. CONCLUSION

Effective deployment of IoT remains a challenging task with a multitude of pitfalls: once deployed, sensors are subject to drifts, bias and thereby fail to provide the expected and meaningful data. Cost-effectively planning the on-site calibration of IoT devices helps in the sustainable long term operation of deployments. Our work focuses on collaborative and distributed maintenance of an IoT-based system. We build on our experiences in deploying sensors over smart spaces and on addressing the calibration tasks to enhance the quality of the gathered observations. We frame multi-sensor calibration planning as an optimization problem where we propose a two-phase iterative local optimization approach to determine (a) how many calibration iterations are necessary, (b) which sensors should be calibrated at each of these iterations, and (c) the number of mobile calibrators that are required (as well as their respective calibration paths), such that the average cost of all iterations is minimized and under the constraint that the calibrators should not be overloaded. Our evaluation shows that the proposed algorithms solve the calibration planning problem effectively compared to naïve/simple solutions. Future Work: Introducing a service to facilitate large-scale calibration in IoT deployments is promising; while our initial studies demonstrate the value and feasibility of this approach, long-term studies with an end-to-end system will help adapt and fine-tune the calibration process for different applications and external dynamicity.

Fig. 1 .

 1 UCI testbed in Donald Bren Hall (Dept. of Computer Science), showing (a) a node deployed in our lab, where we reproduced the sensor deterioration process; (b) readings from two pairs of low-cost sensors on the test node, 3 months after their installment II. MOTIVATION AND BACKGROUND EXPERIENCE The widespread deployment of low-cost IoT technologies and the increasing popularity of corwdsensing enable fine grained monitoring of physical phenomena that are spatiotemporally distributed. Experiences in crowdsensing and IoT deployments indicate that, in general, measurements gathered from low-cost sensors deviate from the ground truth.

Fig. 2 .

 2 Fig. 2. Calibration planning over a long maintenance period.

  {W[ω]} are valid path(s): constraints in IV-B apply.

19 H

 19 .del (minSp) ; 20 return {W} ← convertPathVecToMatrix (pathSet) ; the calibrator is not overloaded. Pick the spot-calibrator pair that induces the least additional movement time. Note that the spot could be added to an old calibrator (Ln 5-13) or a new calibrator (Ln 15-18). (

  (c) Average cost vs. number of sensors; emergency condition. Fig.6. Impact of the number of sensors on the sensor selection in an indoor environment containing 100 nodes.

  Movement time vs. number of spots.

  Number of calibrators vs. number of spots.

  Running time (sec) vs. number of spots Fig.7. Scalability of the multi-path planning solvers and proposed heuristic algorithms.

Fig. 8 .

 8 (a) Modeling indoor space in a building and (b) Mapping of a indoor/outdoor space: Spots (green markers) are placed; nodes (red markers) compose a spot and communicate through wireless links (blue lines).

  is the element-wise multiplication of matrices. Typically, there is no need to run an iteration if no sensor needs immediate calibration; also, special needs and unexpected changes could be easily addressed by altering the TTNC matrix. Hence, the time to next iteration (TTNI) after ω is the minimum TTNC of all the sensors, i.e. t ω+1 -t

ω = min F[ω + ]; thus, the TTNC matrix immediately before the next iteration is

  1 function solveMPPGreedy (G, ĉ, H, β) ; Input : G -Map; ĉ -Maximum workload; H -Set of selected spots; β -Vector of calibration time at selected spots. Output: {W} -Set of paths.

	2 Initialize pathSet ← [ ] ;
	3 while H is not empty do
	4	minInc ← +∞ ; minSp ← minPath ← null;
	5	for each path in pathSet ; last ← path[-1] do
	6	oldTime ← getMoveTime (path) + l∈path β l ;
	7	for each sp in H do
	8	dCw ← G[last, sp] + G[sp, 1] -G[last, 1] ;
	9	if oldTime + dCw +β sp ĉ then
	10	if dCw < minInc then
	11	minInc ← dCw ;
	12	minSp ← sp ; minPath ← path ;
	13	if minInc is finite then minPath.append (minSp) ;
	14	else
	15	for each sp in H do
	16	if (dCw ← G[1, sp] + G[sp, 1]) < minInc then
	17	minInc ← dCw ; minSp ← sp ;
	18	newPath ← [minSp] ; pathSet.add (newPath) ;

TABLE I EXPERIMENTAL

 I SETUP FOR THE PERFORMANCE EVALUATION.

		Input Data	Indoor Normal Emergency	Outdoor
	Map	L Num. of Spots G Pairwise Dist. a	60 53 sec	63 73 min
		K Num. of Types	10		8
	Sensor	τ Calib. Time	1-30 min	0.25-1 min
		T Calib. Period	14-91 d	7-91 d	28-123 d
	Nodes	N Num. of Nodes	100 (varies if independent variable)
		Sensor Presence	50 (varies if independent variable)
		ĉ Max. Workload	2 hours		4 hours
	User Req.	µ Coefficients	C it =10000, µ 0 =1 µw=5, µc=1
		T Maintenance P.	360 days

a Pairwise distance correspond to the shortest traveling time.

Acknowledgement:

The authors thank all members of the ISG and DSM groups at UCI for the valuable discussions. We also thank Inria for supporting the MINES associate team. This project is funded under NSF award CNS-1450768 and NIST award 70NANB17H285.

for calibration planning This includes modules to facilitate the modeling of the (indoor/outdoor) environment and to provide navigation guidance to the calibrators (via an Android app). Experimental Setup: We conduct a series of evaluations using three sets of input data. The two first involve the instrumented building at UC Irvine (Figure 4a) we discussed in Section II, which is used for everyday monitoring (normal condition) and for supporting emergency operations when needed. The desired calibration frequency is a parameter that is learned from empirical study and that depends on the context. We consider the actual deployment in our real testbeds and generate additional spots (with sensors and nodes) using a similar pattern. Using the service we present in Section V, we synthesized the third data set that relates to an outdoor urban environment in which sensors are placed to monitor noise and air quality. The parameters are summarized in Table I.

We compare our sensor selection with two naïve sensor selection strategies that aim at "always selecting all sensors" (regardless of their TTNC) and "only selecting the minimal set of sensors" (i.e. those we must calibrate because their TTNC reaches 0). In addition, we also investigate two simple selection strategies. The former consists in selecting all sensors that are co-located with the sensors that form the minimal set ("local"). The latter consists in only selecting the sensors that can automatically calibrate with each other without human intervention and that henceforth do not induce additional calibration time ("local bounded"). For multi-path planning, we evaluate the performance of two MIP solvers (GLPK and Gurobi), our two path planning heuristics (NN-based greedy and GA), and a naïve strategy that sends one calibrator to each spot and that should give the highest cost. Algorithm running time is evaluated on the OpenLab cluster of Dept. Computer Science at UCI, where each computing node has 2x Quad-core Intel Xeon 3.0GHz CPU E5450 CPUs. Indoor vs. Outdoor Results: In an indoor environment (Figure 3), our algorithm (TTNI-driven sensor selection and GAbased multi-path planning) always result in a lower average cost for N ranging from 5 to 200. Compared to the naïve sensor selection strategies, such as "selecting all sensors" and "selecting the minimal set of sensors" (still considering GA-based multi-path planning), our algorithm combination provides up-to 30% improvement in the long-term average cost. Note that even though our algorithm does not always end up with the lowest cost per iteration (Figure 3b), it makes a fair trade-off between the cost and the time (between iterations). Figure 3c shows that the time spent to select sensors is shortless than 1 sec for a reasonably complex building incorporating 200 nodes and 60 spots. The same trend also applies in the outdoor environment (Figure 5), where the distances between spots are significantly longer (Table I). As the spatial span of the setup grows, the difference among the algorithms becomes more dramatic (note the different y-scale in Figures 3a,3b and5a,5b). Certain naïve approaches are very sensitive to this change (Figure 5b, "minimal" and "local bounded"), while our algorithm shows stable performance in both settings. Normal vs. Emergency Condition Results: Having demonstrated the effectiveness of our algorithm in indoor and outdoor settings, we further study the performance of our approach in an emergency scenario where the calibration requirements of certain sensor types are increased (Figures 6a and6c, note the difference in y-scale). When calibration is required more frequently for some sensors, the naïve/simple approaches suffer from a big increase in average cost, especially when the sensors are deployed densely ( Q/(N •K)>0.5), while the performance of our algorithm and "local bounded" are less affected. We also notice that unlike the simpler approaches ("local bounded" or "minimal"), our TTNI-driven sensor selection algorithm avoids the desynchronization of the periodical calibrations, while the "local bounded" strategy does so with a small number of sensors (Figure 6b). Scalability Results for the Multi-Path Planning:


Simon-Frédéric Désage
email: simon-frederic.desage@univ-savoie.fr

Gilles Pitard
email: gilles.pitard@univ-savoie.fr

Maurice Pillet
email: maurice.pillet@univ-savoie.fr

Hugues Favrelière
email: hugues.favreliere@univ-savoie.fr

Jean-Luc Maire
email: jean-luc.maire@univ-savoie.fr

Fabrice Frelin
email: fabrice.frelin@univ-savoie.fr

Serge Samper
email: serge.samper@univ-rennes1.fr

Gaëtan Le Goïc
email: gaetan.le-goic@u-bourgogne.fr

Extended visual appearance texture features

Keywords: Soft metrology, visual inspection, BTF, ABRDF, Texture, Haralick features, BHF, surface characterization

The research purpose is to improve surface characterization based on what is perceived by human eye and on the 2006 CIE report. This report denes four headings under which possible measures might be made: color, gloss, translucency and texture. It is therefore important to dene parameters able to discriminate surfaces, in accordance with the perception of human eye. Our starting point in assessing a surface is the measurement of its reectance (acquisition of ABRDF for visual rendering), i.e. evaluate a set of images from dierent angles of lighting rather than a single image. The research question is how calculate, from this enhanced information, some discriminating parameters. We propose to use an image processing approach of texture that reects spatial variations of pixel for translating changes in color, material and relief. From a set of images from dierent angles of light, we compute associated Haralick features for constructing new (extended) features, called Bidimensional Haralick Functions (BHF), and exploit them for discriminating surfaces. We propose another framework in three parts such as color, material and relief.

This report denes four headings under which possible measures might be made: color, gloss, translucency and texture. However, when a customer, and thus imitating his behavior, a human controller, inspect the visual appearance of a product, they can not dierentiate the four parties. Considering color and translucency as known and controlled by the industry, we have chosen to work on the gloss and texture to quantify the impact on the human perception. This paper is then an echo of a soft metrology approach, 2 because this paper presents some measurements to compare with the visual human behavior. A known diculty is the interdependence of these measurements because translucency can inuence color, which may inuence gloss, and texture is probably a combination of all three.

3 Our starting point is the study and the imitation of human behavior in the surface quality inspection (Figure 1). From surface quality inspection, the reectance problem is a simple problem between three parts : light, the studied surface and the observer. Indeed, the formalization of human quality inspection consists of three main phases : Exploration, Evaluation and Decision. [START_REF] Baudet | The visual inspection of product surfaces[END_REF] There are only two equivalent interactions dening the human perceptual system with the visual surface appearance (Figure 2). The rst interaction corresponds to the one between visual surface appearance and the human eye, and the second one between scanned surface and the human brain. We propose to follow these interactions to construct objective visual texture features. For imitating the rst interaction, a goniophotometer is used for optic datas acquisition and the exploitation of most advanced visual surface representation allows to the link with human eye. The rst part describes the use of digital surface appearance and the measurement instrument. The second part of article recalls classical image processing approach, which allows to human-like texture classication.

The interest of this paper, described in the third part, is the combination of two previous approaches to describe better surfaces and their gonio-appearances in relation to human perception. The problem relates to the characterization of complex surfaces, such as for surfaces having the same hue, the same relief and dierent material.

DIGITAL SURFACE REPRESENTATION

Some recent works have shown dierent representations of visual surface as appearance measurement. Here, we would like to highlight the book "Visual Texture" 5 of Mr Haindl and Mr Filip is remarkably clear about the reectance problem, and which lists these recent works. We remind taking direct their denitions of the General Reectance Function (GRF) (Equation 1 and Figure 3), their reectance functions Taxonomy (Figure 4) and the model corresponding to our measurement system : the Surface Reectance Field (SRF) (Figure 5.)

Y GRF r = GRF (λ i , x i , y i , z i , t i , θ i , ϕ i , λ v , x v , y v , z v , t v , θ v , ϕ v , θ t , ϕ t ) (1)
where r = [r 1 , ..., r 16 ] is the multi-index with corresponding partial indices. All possible values of the index will be denoted by •, e.g., a color input spectrum in the RGB space Y •,r2,...,r16 = [Y R,r2,...,r16 , Y G,r2,...,r16 , Y B,r2,...,r16 ] and the missing index by φ, e.g., a monospectral input Y φ,r2,...,r16 .

GRF describes the incident light with spectral value λ i ; illuminating surface location x i , y i , z i in time t i ; under spherical incidence angles ω i = [θ i , ϕ i ] and observed at time t v from surface location x v , y v , z v under spherical reectance angles ω v = [θ v , ϕ v ] and spectrum λ v ; here ω t = [θ t , ϕ t ] are the corresponding transmittance angles where ω = [θ, ϕ] are the elevation and azimuthal angles, respectively. The model height parameters z i , z v indicate that radiance along light rays is not constant but depends on the height. The GRF function (Equation 1) is too complex to be accurately measured or modeled, hence some simplifying assumptions are inevitable in any practical application. The taxonomy of simplifying assumptions can be divided into two subgroups based on the possibility of neglecting a surface texture. We also share the denition of a visual texture, which is a resolution-based relative notion. The scale of observation is important because any natural surface material is textured and our perception of surfaces as textured or smoothly homogeneous (i.e., non-textured) only depends on the corresponding surface resolution.

The same surface observed from a distance can be categorized as smooth, while its close observation may reveal rough-textured surface. The GRF simplifying taxonomy (Figure 4) is obviously not exhaustive, but it allows the reader to quickly locate themselves in relation to the range of possibilities, and locate the model used. Two particular features may correspond to an industrial system, with limited degrees of freedom because to meet industrial requirements, it is necessary to x extent possible the maximum settings.

In our vision problem, we have three sources of settings: the light source, the study area and the viewpoint.

In the general cases of inspection systems, the study position of the object is xed and at least the position of another between the light source and the viewpoint. If the light source position is xed, a Surface Light Field (SLF) model-like function is measured and verify the equation 3, as the gure 6 shows. On the contrary, if the viewpoint is xed, a Surface Reectance Field (SRF) model-like function is measured and verify the equation 2, as the gure 5 shows. Between to multiply the light sources or viewpoints, the calibration protocol such as nancial stress have occurred in our choice. In other words, a light source costs less than a photonic sensor and is easier to calibrate, so we opted for a system complying with the surface light eld model. The Surface Reectance Field (SRF) 5 has then as 5 dimensions (5D) with a xed viewpoint and a free light position :

Y SRF r = GRF (λ, x, y, θ i , ϕ i ) (2)
The Surface Light Field (SLF) 5 is a similar function of SRF, with a xed incident light and a free observer :

Y SLF r = GRF (λ, x, y, θ v , ϕ v ) (3)
To be rigorous physical part of the system used, we took into account assumptions. We only mention the following assumptions, useful to simplify the SRF model:

A1 light transport takes zero time (t i = t v and t v = φ) A2 reectance behavior of the surface is time invariant (t i = t v = const., t i = t v = φ) A3 interaction does not change wavelength (λ i = λ v , i.e., λ v = φ) A4 constant radiance along light rays (z i = z v = φ) A5 no transmittance (θ v = ϕ t = φ) (no
simultaneous reectance and transmittance) A6 incident light leaves at the same point

x i = x v , y i = y v (x v = y v = φ) A13 xed illumination (θ i = const., ϕ i = const.); For SLF only A14 xed viewing angle (θ v = const., ϕ v = const.); For SRF only 2.

Acquisition

In industrial terms of quality control, the Commission Internationale de l'Eclairage (CIE) recommends lighting conditions such as the type of light and its position on the surface (Figure 7). The standard illuminant is D65 1 (Figure 8). There are other experimental measurement protocols to identify the specularity type 6 of the color of specular surfaces. We propose to use a similar method for gloss measurement with a recent visual texture representation called Bidimensional Texture Function (BTF) 8 (Figure 10), because it is the closest function of controller observation.

One diculty for controller is to repeat the same light path and the same observation function, because he can not realise the whole BTF. In pratice, each human inspection is a part of BTF. Our bias is to achieve a similar and regular inspection for each surface. We set the position of viewpoint and the relative position of light sources to be used to perform the structured lighting sequence (Figure 11) for the chosen SRF model (Figure 5). This behavior is caused by the surface relief, such as self-shadowing, self-occlusions, subsurface scattering and other complex eects occurring in the material structure which are not represented in the true BRDF models. The PTM method models illumination dependence of individual pixels using the following pixel-wise biquadratic formula :

Reconstruction

SRF v (r, i) ≈ a 0 (r)u 2 x + a 1 (r)u 2 y + a 2 (r)u x u y + a 3 (r)u x + a 4 (r)u y + a 5 (r) (4)
where u x , u y are projections of the normalised light vector into the local coordinate system r = (x, y). The set of n i pixels is considered as reectance data, where i = 1, ..., n i is the illumination position index and v is the actual view position index v = 1, ..., n v . The n p = 6 polynomial coecients a 0 -a 5 are tted in each pixel by means of singular value decomposition. From these coecients,the approximated function is known.

Extraction

Once reconstructed, each ABRDF can be used for calculate the surface appearance corresponding for any incident light position from superior hemispheric space. Whereas ABRDF has not the physic sense like a true BRDF, the shape of ABRDF is may be seen as the Phong model (Figure 13 We propose a surface appearance framework from ABRDF slice in three parts: color, material and relief.

The shape of ABRDF slice is a feature of material, respectively the power/amplitude for the color (Dependent of each wavelength), and the direction for the relief related to the optical law of Snell-Descartes (Figure 13 (a)).

The amplitude corresponds to the maximum of function and the function direction is the direction of maximum such as the specular direction. We can extract three types of image for each information type from ABRDF slice. This framework (Figure 14) can help to discuss about visual texture denition. A photography, as image texture, is a visual texture slice and it is a local result of color, material and relief combination. There is surface texture such as only relief texture.

Figure 14: Texture framework One image combines these three types of information whereas we can distinguish three sub-textures from ABRDF slice (Figure 13). For quality inspection, the advantage is to guarantee to nd directly visual anomaly and their causes. The right image is the one with the anomalies if they exist. Then, we must characterize each image photography or numerical representation for each visual texture component.

TEXTURE CLASSIFICATION

Some works have shown the useful of Haralick features for texture classication.

11 The diculty to get a good classication was to have the good image of the product surface, i.e. an image with a uniform or well directed lighting and a uniform or well directed viewing. For quality inspection, the right image is the one with the anomalies if they exist, even if lighting is not uniform. Of course, it is necessary to have a suitable resolution.

So, it is necessary to identify the image contents. Each method is used for one image.

Co-occurrence matrix

Several works are identied as using gray-level co-occurrence matrix (GLCM) method for defect detection.

12 The co-occurrence matrix method is a statistical method for textural defect detection. There are other statistical methods such as histogram properties, local binary pattern (LPB) or autocorrelation. There are also other methods: structural methods, lter-based or model based methods 12 . 13 In this paper, we focus on the co- occurrence matrix method to ease the demonstration.

Figure 15: The 4 directions considered with a step of 1 pixel

As classical processing, Haralick features used to describe a single image. We do a recall of the generic method to compute a GLCM and applied Haralick features, 11 while GLCM method is one of the most well-known and widely used texture features. GLCM measures the spatial dependency of two gray-levels, given a displacement vector. The classical displacement has a size of 1 pixel, but it can be adapted in function of texture frequency.

There is 4 directions in an image (with a square matrix unit). There are many interpretations of co-occurrence matrix. The most common method is to use sum of 4 directions matrix.

m[I, θ](i, j) = Card {p, p + t} I 2 ; f (p) = i, f (p + t) = j and dir(t) = θ (5)
where I denotes the image, f is its associated function, θ the transition direction, i, j are gray-levels, p is the start pixel and t denotes transition step to get the second pixel.

To be rotation invariant in the square matrix of pixels, we propose to use the sum of the 4 directions matrix such as :

M [I] = 4 nt=1 (m[I, θ] θ=45 * (nt-1) (i, j) i [0..2 N ],j [0..2 N ] ) (6)
where N is the gray-levels number. 11 There are also other features 14 set out more recently to describe a co-occurrence matrix. For principle, our demonstration will be limited to typical and the most used in the litterature attributes, such as energy, contrast and homogeneity. One must know that some attributes have been the subject of intuitive analysis in their contribution to characterization of textures.

14

• Feature number 01 : Angular Second Moment or Energy

f 01 = N -1 i=0 N -1 j=0 {M [I, t](i, j)} 2 (7)
• Feature number 02 : Contrast

f 02 = N -1 n=0 n 2 N -1 i=0 N -1 j=0 M [I, t](i, j) (8)
• Feature number 05 : Inverse Dierence Moment or Homogeneity

f 05 = N -1 i=0 N -1 j=0 1 1 + (i -j) 2 M [I, t](i, j) (9)
The approach is to process an image of a surface (Figure 17 the one with the anomalies if they exist. Hence, the idea is to stick to the controller behavior and collect all the useful images of surface. However, it is necessary to process this new and enriched data relative to a single image. We present the (extension) adaptation of Haralick to BTF and some results of Bidimensional Haralick Functions (BHF).

The ideal case will be to compute co-occurrence matrix from exact ABRDFs (or maybe their coecients if they are exacts too) and their Haralick features, but it is [too (for the moment)] complex. So the idea is to do the inverse, i.e., to compute co-occurrence matrix for acquired images sequence and their Haralick features, for reconstructing a bidimensional function of Haralick features.

The data volume does not allow to present each result. We have chosen to give the preview of theoretical cases to present the generic method.

The demonstration will be focus on simple and interesting cases. Before to process complex cases where surface has many dierent subtextures, we process with uniform textures. The interesting cases are ones which have goniovariance, i.e. an anisotropic appearance.

Received incident light function

Some cases, such as a lambertian and plane surface, require no new feature since they are isotropic, and this is true if the lighting incident is isotropic, i.e., compensated for the surface to receive the same amount of light regardless of the impact. Otherwise, for these, it will be sucient to calculate Haralick features once to know the magnitude of the associated hemispherical BHF.

Considering the previous A4 assumption about constant radiance, there are two possible interpretations.

The rst interpretation is each lighting source emits the same quantity of light (Figures 18 (a

Study cases

Considering the previous framework : Material, color and relief, the demonstration show dierent cases with some states for each parameter. As shown previous (Figure 19), we consider two relief. It is normal to consider equally two color levels and two material levels. For color levels, there are a case with the same graylevel and a case with a graylevels gradient, as dierences between Figures 20 (a) and (b). For material, the dierence is made between the behavior : the rst case is Lambertian, i.e., totally diusing, and the second is only specular as a metal, such as Figure 13 show generic BRDF decomposition. We take homogeneity for example, because it is a known parameter as relevant for quality assessment, particularly in the perceptual quality of printed surfaces.

15 It is interesting to observe that the transformation of the Lambertian or specular surfaces BHF are not the same as the surfaces are subject to the same changes (Figure 26). These features of transformation can then be considered as a signature. Otherwise, the results of the cases Clearly, the fourth polynomial coecient seems pointless to classify dierent cases.

Figure 28: Polynomial decomposition coecients corresponding to respective BHF

CONCLUSION AND PERSPECTIVES

We have proposed a combined method for surface characterization through an image processing approach and scanning their complete visual texture, i.e., for the surface appearance seen under an illumination hemisphere.

The method provides enriched features compared to the conventional method. Visual trend of these new features is that they can help distinguish visual texture into three sub-textures related relief, material and color.

In this sense, this trend should be conrmed with the application of the method on a larger variety of surfaces and taking into account a signicant number of features (at least all Haralick features, for example). The application could be a design of experiments with real samples.

There are three perspectives of this work for visual inspection. The rst perspective is when visual inspection is applied to aesthetic eld, there are two challenges:

1 -Have a Repeatable and Reproducible (R&R) method. [START_REF] Maire | Gage r2&e2: an eective tool to improve the visual control of products[END_REF] 2 -Have a human perception-like processing with visual attributes.

We can use some elements of this paper for establishing features of visual data. This data can be used for trying to understand the human visual behavior, because we can verify uniformity and reproducibility of visual exploration and evaluation data. The second perspective is for automatic defect classication. In fact, some works have shown the useful of Haralick features for automatic defect classication from dierent co-occurrence matrix of one image 14 . 17

  quality of visual appearance is [still] an industrial matter. The main diculty is to link objective measurements and subjective human aesthetic indicators. The CIE (International Commission on Illumination) provided a technical report CIE 175:2006 A framework for the measurement of visual appearance .

Figure 1 : 1 Figure 2 :

 112 Figure 1: The vision scheme from CIE 1

Figure 3 :

 3 Figure 3: The General Reectance Model 5

Figure 4 :

 4 Figure 4: The General Reectance Taxonomy 5

Figure 5 :

 5 Figure 5: Surface Reectance Field Model

5

 5

Figure 6 :

 6 Figure 6: Surface Light Field Model

5

 5

7

 7

Figure 7 : 1 Figure 8 :

 718 Figure 7: The acquisition system 1

Figure 9 :

 9 Figure 9: (a) The brilliant eect (b) The black light eect (c) The mirror eect

Figure 10 : 5 Figure 11 :

 10511 Figure 10: The BTF model 5

Figure 12 :

 12 Figure 12: (a) The acquired images sequence (b) The measured points cloud for one pixel (c) The ABRDF function

Figure 13 :

 13 Figure 13: (a) The Phong BRDF model (b) Phong-based Surface appearance framework (c) Adapted Phong model to ABRDF (d) ABRDF-based Surface appearance framework

Figure 16 :

 16 Figure 16: (a) Some transitions for 0 • direction (b) GLCM for 0 • direction (c) Some transitions for 4 directions (d) GLCM sum for 4 directions

Figure 17 :

 17 Figure 17: (a) (resp.(d)) One example image (b) (resp.(e)) The associated GLCM (c) (resp.(f )) The associated features value

  ), (b) and (c)), and the second interpretation is the surface receives as much light for each source (Figures18 (d)). Only the second interpretation is isotropic. Each interpretation denes one received lighting function, which inuences the measured appearance behavior of surface.

Figure 18 :

 18 Figure 18: (a) Illumination received by a unit area for a light source disposed to the surface normal (b) Respectively in (a) to a source having a direction inclined with respect to the normal (c) Slice of received illumination function for constant sources (d) Slice of compensated illumination function for a constant radiance

Figure 19 :

 19 Figure 19: (a) Plane relief (b) Surface plane and Lambertian (c) No plane relief (d) Lambertian surface with no plane relief

Figure 20 :

 20 Figure 20: (a) Plane and uni-color (b) Plane and color gradient (c) No plane relief and uni-color (d) No plane and color gradient

Figure 21 :Figure 22 : 6 (Figure 25 :

 2122625 Figure 21: The study cases scheme

4 and 8 ,Figure 26 :Figure 27 :

 82627 Figure 26: Projection of BHF for Homogeneity in the X-Y reference

  The third perspective is to use the dierent well-known methods for editing BTF-like functions, such as Polynomial extended Lafortune Reectance Model (PLM RF), Reectance eld factorization (PCA RF) or Hemispheric harmonics (HSH) 5 . 10

ACKNOWLEDGMENTS

We thank our SYMME partners in MESURA project as well as Savoie Mont Blanc Industries and Conseil General 74 to enable us to carry out this research by giving us resources.


Safaa Jarifi
email: safaa.jarifi@enst-bretagne.fr

Dominique Pastor
email: dominique.pastor@enst-bretagne.fr

Olivier Rosec
email: olivier.rosec@francetelecom.com

Cooperation between global and local methods for the automatic segmentation of speech synthesis corpora

This paper introduces a new approach for the automatic segmentation of corpora dedicated to speech synthesis. The main idea behind this approach is to merge the outputs of three segmentation algorithms. The first one is the standard HMM-based (Hidden Markov Model) approach. The second algorithm uses a phone boundary model, namely a GMM (Gaussian Mixture Model). The third method is based on Brandt's GLR (Generalized Likelihood Ratio) and aims to detect signal discontinuities in the vicinity of the HMM boundaries. Different fusion strategies are considered for each phonetic class. The experiments presented in this paper show that the proposed approach yields better accuracy than existing methods.

Introduction

This paper deals with the problem of automatic segmentation of speech corpora for concatenative TTS synthesis systems. In the development process of such systems, the segmentation of large databases constitutes a key task. Obviously, the optimal segmentation to use in these systems is the manual one. Nevertheless, an accurate automatic segmentation saves a lot of human effort and time in creating new synthetic voices and thus drastically simplifies the personalization of a TTS speech synthesis.

Up to now, the HMM approach [START_REF] Brugnara | Automatic segmentation and labelling of speech based on Hidden Markov Models[END_REF][START_REF] Nefti | Segmentation automatique de la parole en phones[END_REF] is the most widely used for automatic segmentation and it is considered as the most reliable. This approach is linguistically constrained because it needs the true phonetic sequence associated to the recorded utterances in order to estimate the HMM sequence. Then it applies a forced alignment between this HMM sequence and the speech signal. However, this approach has still some limitations for building voices for TTS systems based on the principles of unit-selection and concatenative synthesis. The main limitation is that HMMs model well steady areas but are not really suited to detect locally the transitions between phonemes in a speech signal. For this reason, in order to guarantee a good quality of synthetic voices, a manual checking is applied to the HMM segmentation before synthesis.

Brandt's GLR algorithm [START_REF] Obrecht | A New Statistical Approach for the Automatic Segmentation of Continuous Speech Signals[END_REF] is another suitable approach for

This study is supported by France Telecom segmenting speech signals. Nevertheless, it produces insertions and omissions because it is linguistically unconstrained.

With respect to the foregoing, the purpose of this paper is to combine global and local automatic segmentation algorithms. To achieve this, we choose in this paper three automatic segmentation algorithms. The first is the HMM segmentation. The second uses a boundary model which is estimated on a small database and which is used to refine the HMM segmentation marks. The third one is Brandt's GLR method which was modified in order to avoid omissions and insertions. These algorithms are described in section 2. In section 3, two combination methods are proposed and evaluated on French and English corpora dedicated to speech synthesis.

Segmentation algorithms

Segmentation by HMM

This approach generally consists of two steps. The first step is the training phase that aims at estimating the acoustic models. In the second step, these models are used to segment the speech signal by the Viterbi algorithm. This one applies a forced alignment between the models associated to the known phonetic sequence and the speech signal.

Note that the training is a decisive step because the accuracy of the obtained segmentation heavily depends on the quality of the estimated models. One solution to perform well the training step is to resort to an iterative training [START_REF] Kim | Automatic segmentation combining an HMM-based approach and spectral boundary correction[END_REF]. The phone labels resulting from the previous iteration are used for initializating and re-estimating the HMMs via the Baum-Welch algorithm. After a few iterations, mismatches between segmentation marks produced by an HMM approach and marks obtained manually are considerably reduced as shown in [START_REF] Jarifi | Coopération entre méthodes locales et globales pour la segmentation automatique de corpus dédiés à la synthèse vocale[END_REF]. Another method to train the models is to use a representative small speech database manually labeled and segmented [START_REF] Matousek | Automatic Segmentation for Czech Concatenative Speech Synthesis Using Statistical Approach with Boundary-Specific Correction[END_REF]. We estimate first the models using this small database. Then we segment the whole corpus with the models. As the estimation of the models on the small corpus is accurate, the processing offers better results than the iterative training. For that reason this strategy is used in this paper.

Refinement by boundary model

The main idea of this method is to train models for boundaries on the basis of a small database segmented and labeled manually. Then, these models are used to refine an initial segmentation [START_REF] Wang | Refining Segmental Boundaries for TTS Database Using Fine Contextual-Dependent Boudary Models[END_REF].

For each boundary of the training database, we create a super vector as mentioned in figure 1 by concatenating the acoustic vectors of size Nc associated to the (2N + 1) frames around the boundary. Because the number of labeled data is limited in practice, the boundaries are clustered into classes using a classification and regression tree (CART). Then a Gaussian model is estimated for each class.

The second step aims at refining each boundary of every segment given a labeled sentence and its initial segmentation. For that purpose, we seek, in a certain vicinity of each boundary, the time instant that maximizes the likelihood of its super vector in comparison with the Gaussain model of this transition.

Brandt's GLR algorithm

The aim of this method is to detect discontinuities in speech signals. Speech signals are assumed to be sequences of homogeneous units. Each unit or window w is a finite sequence w = (yn) of samples that are assumed to obey an AR model : yn = P p i=1 aiyn-i + en. In this equation, p is the model order, which is assumed to be constant for all units and en is a zero mean white Gaussian noise with variance equal to σ 2 . Such a unit is thus characterized by the parameter vector Θ = (a1, . . . , ap, σ). Let w0 be some window of n samples and Θ0 the corresponding parameter vector. The authors of [START_REF] Brandt | Detecting and estimating parameters jumps using ladder algorithms and likelihood ratio test[END_REF][START_REF] Obrecht | A New Statistical Approach for the Automatic Segmentation of Continuous Speech Signals[END_REF] attempt to decide whether w0 should be split into two subsegments w1 and w2 or not. In fact, a possible splitting derives from the detection of some jump between the parameter vectors Θ1 and Θ2 of w1 and w2 respectively. Brandt's GLR method decides that such a jump has occurred by comparing :

Dn(r) = nlogσ0 -rlogσ1 -(n -r)logσ2
to a predefined threshold λ. Note that Dn is merely the GLR. In the equation above, r is the size of the time interval covered by w1, whereas σ1 and σ2 are the noise standard deviation estimates of the models characterized respectively by the parameter vectors Θ1 and Θ2. Thus, the change instant corresponds to arg(maxr(Dn(r)) ≥ λ).

As mentioned before, the basic Brandt's GLR method is an algorithm capable of detecting discontinuities of speech signals without any further knowledge upon the phonetic sequence. As this algorithm is linguistically unconstrained, it makes insertions and omissions. However, for TTS synthesis, the phonetic sequence of every utterance is known. For this reason, we decide to take into account this information by using the boundaries produced by a segmentation algorithm that uses the phonetic sequences as the HMM segmentation. More precisely, let (U0, U1, . . . , UL) be the boundaries obtained by such initial algorithm. For i in {1, . . . , L -1}, we seek a speech discontinuity between Vi =

(U i-1 +U i ) 2 and Vi+1 = (U i +U i+1 ) 2
by applying a modified Brandt's GLR method : to avoid omissions and insertions, the use of the threshold is replaced here by the maximization of the GLR on [Vi, Vi+1].

Evaluation of the three algorithms

In this section, we present the experimental results obtained by the three segmentation algorithms described above on French and English corpora. Both acoustic databases were recorded by a professional native female speaker and sampled at 16 kHz. The French corpus "corpusFR" and the English corpus "corpusEN" respectively contain 7300 and 8900 sentences.

The segmentation by HMM uses the HTK toolkit [START_REF] Young | The HTK Book for HTK V3.2.1[END_REF] for the acoustic analysis, the training and the segmentation steps. It considers mixtures of 2 Gaussian density and the acoustic vector contains 39 coefficients which are the 12 MFCCs (Mel Frequency Cepstral Coefficients), the normalized energy, and their first and second derivatives. The HMM segmentation obtained by the use of a small training database is called hereafter HMMSeg. Twenty iterations of the Baum-Welch algorithm are applied to train the HMM. The refinement by boundary model is applied to the HMM segmentation and the parameters N , M and e of figure 1 are fixed to 2, 20 ms and 30 ms respectively. These parameters were adjusted on the French corpus in [START_REF] Jarifi | Brandt's GLR method & refined HMM segmentation for TTS synthesis application[END_REF]. The segmentation thus obtained is called RefineSeg. The segmentation obtained with Brandt's GLR method initialized by the HMM segmentation is denoted BrSeg. The model order is set to 12 and the minimal length of w1 and w2 are equal to 10 ms.

The training phases of HMMSeg and RefineSeg were carried out with the number SizeAlg of training sets equal to 100, 300 and 700. For each set of learning sentences randomly chosen, the test set was built by considering the remaining sentences in the corpus under study. Moreover, to illustrate the consistency of the results a cross-validation procedure was used where 3 trials were done for each value of SizeAlg. All the accuracies are calculated at a tolerance equal to 20 ms. This value is commonly considered as an acceptable limit to guarantee a good quality of a synthetic voice. The results presented here are obtained by averaging the accuracies using this cross-validation procedure.

Table 1 presents the accuracies of each algorithm with respect to SizeAlg for each corpus. Table 2 shows the limit of performance of each algorithm. This limit corresponds to using the whole database for the training of HMMSeg and RefineSeg. According to these results, we can make the following remarks :

-the refinement by boundary model gives the best results when the corpus size equal 300 or 700 for "corpusFR" and 700 for "corpusEN". This is normal because the boundary models are well trained ; -the modified Brandt's GLR method is inaccurate at 20 ms in comparison with the other algorithms ; Because the important measure in TTS systems is the accuracy at 20 ms, it seems reasonable to say from table 1 that the refinement by boundary model is the most accurate algorithm. Nevertheless, we should not forget that the algorithms are suited to different phonetic classes. In fact, during these tests, it turned out that Brandt's GLR method detects well some boundaries like silence/speech and voiced/unvoiced transitions. Thus, depending on the classes to the right and the left of a transition to detect, it seems relevant to take into account that the algorithms do not perform equally. This is the purpose of combination methods.

Combination of several segmentations

Principles of the merging processes

Segmentation algorithms behave differently according to the transitions they are asked to detect. The main idea here is to take into account the different behaviors of segmentation algorithms so as to favor more some segmentation marks than others given a certain type of transition to detect. We thus propose a processing that merges K boundaries produced by K different algorithms.

Let {c1, . . . , cT } be a set of T phonetic classes. By using a small database segmented manually, we start by estimating the segmentation accuracy α k (ci, cj ) at a tolerance of 20 ms for every algorithm indexed by k, k = 1, 2, . . . , K and every pair (ci, cj ) of classes, (i, j) ∈ {1, . . . , T } 2 .

Then, let s be an unknown transition time instant to estimate given two classes c ℓ (s) and cr(s). Given t k (s), k = 1, 2, . . . , K, the K estimates of s returned by the K available algorithms, the merging phase consists in estimating s on the basis of these K estimates and the accuracies α k (c ℓ (s), cr(s)).

The first solution that we propose consists in choosing the algorithm that offers the best accuracy computed during the training for the type of transition under consideration. Following the terminology used in [START_REF] Pigeon | Authentification multimodale d'identité[END_REF], we adopt a linear hard combination given by :

thard (s) = P k∈A t k (s) Card(A) ( 1
)
where A is the set of algorithms k that maximize α k (cg(s), c d (s)), k ∈ {1, . . . , K}. Note that A is not restricted to one element. For example, if the transitions between two classes i and j are absent from the training database, the segmentation accuracy α k (i, j) is not defined and thus, we impose α k (i, j) to be equal to 1 for each k. In this case, Card(A) = K and the equation (1) becomes a simple average of the K time instants produced by the K algorithms.

Still following [START_REF] Pigeon | Authentification multimodale d'identité[END_REF], another solution is to apply a soft combination for the K boundaries performed by the K algorithms. Thus the estimated boundary for the transition s is computed as the barycenter of the K time instants as :

tsoft (s) = P K k=1 α k (c l (s), cr(s))t k (s) P K k=1 α k (c l (s), cr(s)) . (2)
Similarly to the hard combination, α k (i, j) is set to 1 for each k if the transitions between the classes i and j are absent from the training database. Obviously, this case is equivalent to a hard combination.

Experimental results

In this section, we present the results obtained by the application of the two combination methods (see the equations ( 1) and ( 2)) to the triplet (HMMSeg, RefineSeg, BrSeg).

For the French corpus, the fusion was achieved by using 12 classes : unvoiced plosives, voiced plosives, unvoiced fricatives, voiced fricatives, oral vowels, nasal vowels, diphtongues, nasal consonants, liquid consonants, semi vowels, pauses and silences. For the English corpus, 10 classes were used : vowels, voiced/unvoiced plosives, voiced/unvoiced fricatives, nasal consonants, liquid consonants, semi vowels, pauses and silences. Let SizeComb denote the number of sentences in the training database used for the merging process. The estimation of the accuracies for the soft and hard combinations is achieved by using two different values of SizeComb : 100 and 300. The chosen corpora are different from those used for the training of HMMSeg and RefineSeg. Thus, the accuracies given in this section are computed at a tolerance of 20 ms and evaluated on all the sentences of the database except those employed for training HMMSeg, RefineSeg and the soft and hard combinations. As in section 3, the results presented here are obtained by averaging the accuracies using a cross-validation procedure.

The hard and soft combinations are compared to the so-called isobarycenter method and the optimal soft combination. The isobarycenter method averages the three time instants obtained by HMMSeg, RefineSeg and BrSeg. The optimal soft combination is the soft combination when the accuracies α k (ci, cj ) are estimated on the whole corpus.

The results of the combination methods are given in tables 3 and 4. For every pair (SizeComb,SizeAlg), the accuracies obtained by the four combination methods are larger than the best result of the line associated to SizeAlg in table 1.

Similarly to table 2, table 5 presents the results obtained by using the whole database for training HMMSeg and RefineSeg and estimating the accuracies for the combination methods. These results can be regarded as the limit performance of the merging methods. The difference between this ideal case and the accuracies presented in Tables 3 and4 is not very important, which illustrates that using more learning data hardly improves the results.

In order to get further insight into the behavior of the three combination methods, we analyzed the ability of each combination method to correct errors made by HMMSeg, RefineSeg and BrSeg. By error, we mean a segmentation mark further that 20 ms from the manual boundary. For that purpose, two configurations were studied : the first one where the 3 marks are located on the same side relatively to the manual boundary and the second one where these marks are on both sides of the manual boundary. The results obtained on corpusFR and corpusEN are presented in table 6 and show that in the first configuration, all the combination methods lead to similar performance. However, in the second configuration, the accuracies yielded by the soft and isobarycenter combination methods are respectively 95.07% and 95.71% for corpusFR and 95.44% and 96.30% for corpusEN. These results are significantly higher than the rate obtained by the hard combination method (79.14% for corpusFR and 84.27% for corpusEN).

Conclusion

In this paper we have analyzed the performance of three automatic segmentation algorithms for French and English corpora : a global method based on HMM and two local methods that aim at detecting a transition in the vicinity of a boundary (Refinement by boundary model and Brandt's GLR method).

We have also proposed two methods capable of merging the boundaries produced by the different segmentation algorithms. The experimental results of these two methods applied to two languages show a clear improvement of the accuracy at 20 ms. Furthermore, these methods are simple and not computationally expensive. They seem to be a good prospect regarding the segmentation problem for TTS synthesis applications.

FrameFIG. 1 -

 1 FIG. 1 -Elements of a super vector