Elise Arnaud
email: arnaud@disi.unige.it

Etienne Mémin
email: memin@irisa.fr

Roberto Sosa
email: rsosa@fi.uba.ar

Guillermo Artana
email: gartana@fi.uba.ar

A Fluid Motion Estimator for Schlieren Image Velocimetry

In this paper, we address the problem of estimating the motion of fluid flows that are visualized through a Schlieren system. Such a system is well known in fluid mechanics as it enables the visualization of unseeded flows. As the resulting images exhibit very low photometric contrasts, classical motion estimation methods based on the brightness consistency assumption (correlation-based approaches, optical flow methods) are completely inefficient. This work aims at proposing a sound energy based estimator dedicated to these particular images. The energy function to be minimized is composed of (a) a novel data term describing the fact that the observed luminance is linked to the gradient of the fluid density and (b) a specific div curl regularization term. The relevance of our estimator is demonstrated on real-world sequences.

Introduction

The ability to understand the complexities of fluid flow behavior has large implications in our daily lives and safety as their control and understanding is of the greatest importance in different applications ranging from aero or hydrodynamic studies (air conditioning, aircraft design, etc.) to environmental sciences (weather forecasting, climate predictions, flood disasters monitoring, etc.).

Flow visualization has been a powerful tool to depict flow features. Efforts to develop high-quality flow visualization techniques date back over a century. The analysis of the recorded images consisted firstly to a qualitative interpretation of the streak lines, leading overall global insight into the flow properties but lacking quantitative details on important parameters such as velocity fields or turbulence intensities. Point measurement tools such as hot wire probes or Laser Doppler Velocimetry have typically provided these details. As these probes give information only at the point where they are placed, simultaneous evaluations at different points require to dispose a very large number of probes and the evaluation of unsteady field (most of the flows are unsteady) is almost unachievable.

In an effort to avoid the limitations of these probes, the Particle Image Velocimetry (piv), a non-intrusive diagnostic technique, has been developed in the last two decades. piv enables obtaining velocity fields by seeding the flow with particles (e.g. dye, smoke, particles) and observing the motion of these tracers. An underlying assumption of piv technique is that the motion of these particles follows the motion of the neighboring fluid. This condition is not always satisfied and requires to seed the flow with small sized tracers leading to an increase of the measurement difficulties. Moreover, some phenomena such as natural convection may be influenced by the large amount of seeding particles and the seeding may in return alter results. The setting up of the experiment, adjustment of the seeding concentration and other experimental procedures are in general tedious tasks in many large scale facilities. As a consequence this technique is mainly adapted for test in small closed loops wind tunnels.

Given the various complexities associated to the use of piv, it is important to examine techniques that can be used to generate quantitative measurements of unseeded flows. The techniques that provide useful visualization images and, at the same time, yields high-quality quantitative data about the flow are of particular interest. In general, Shadowgraph, Schlieren and Interferometry fall into this category. These three techniques do not require flow seeding since they are based on index-of-refraction effects. One of the attractive capabilities of the Schlieren technique is that it can be implemented to undertake full scale measurements and outdoor experiments [START_REF] Settles | Full-scale schlieren flow visualization[END_REF][START_REF] Settles | The penn state full-scale schlieren system[END_REF].

The objective of this work is to analyze the ability of a dense motion estimator to extract velocity fields from Schlieren images of fluid flows. To date no satisfying technique exists to perform accurately such velocity measurements. The dense motion estimator we propose here relies on a data model specifically designed for such images. The devised data model has been elaborated on physical grounds. In addition to this constraint, we have also considered a div-curl smoothing function allowing the preservation of curl blobs.

Description of the Schlieren technique

The Schlieren technique [START_REF] Goldstein | Optical systems for flow measurement: Shadowgraph, schlieren, and interferometric techniques[END_REF][START_REF] Settles | Schlieren and shadowgraph techniques: visualizing phenomena in transparent media[END_REF] is an optical method used for fluid flow visualization. Contrary to standard visualization approaches, where a tracer (e.g. solid particle) is followed along the fluid motion or laser-Doppler systems, in which the frequency shift of scattered illumination from such a marker is measured, the Schlieren technique does not require any intrusion in the fluid and prevents any modifications of the considered flow. Such a technique is used to study density fields in transparent media, usually gases or liquids. A typical Schlieren system is described in figure 1. It is based on the fact that a light beam traveling initially in the z direction passing through a medium whose index of refraction varies in x and y direction undergoes a small deviation. For sake of simplicity, the figure 1 presents this phenomenon only in the yz plane. In that case, the light beam has been deviated by an angle α. The Schlieren system is basically a device to observe the angle α as a function of position in the xy plane (respectively the [START_REF] Goldstein | Optical systems for flow measurement: Shadowgraph, schlieren, and interferometric techniques[END_REF] angle in the xz plane). As the light beam deviation depends on the flow density variations, it can be demonstrated that the light pattern obtained with a Schlieren system is determined by the first derivative of the index of refraction such as [START_REF] Goldstein | Optical systems for flow measurement: Shadowgraph, schlieren, and interferometric techniques[END_REF]:

I(s) = K ∂ρ(x, y, z) ∂x + ∂ρ(x, y, z) ∂y dz, (1)
where I(s) is the luminance value of pixel s = (x, y) and ρ(x, y, z) denotes the density of the observed fluid at the physical point of 3d coordinates (x, y, z). The constant K depends on the focal f of the second lens, on the Gladstone-Dale constant C and on a k , the size of the beam cut off by the knife-edge:

K = C f a k . (2)
As described by the equation ( 1), the Schlieren visualization integrates the quantity measured over the length of the light beam. As a consequence, this technique is well suited to the study of almost-2d fields, where no density variation is present in the test section. In that case, the light pattern can be expressed as:

I(s) = K ∆z ∂ρ(x, y, z) ∂x + ∂ρ(x, y, z) ∂y , (3)
where ∆z is the width of the region where the light beam is deflected (supposed small). Since the Schlieren technique is non intrusive and does not require any seeding of particles, this visualization procedure enables studies either for laboratory tests or for full scale models in industrial applications. To illustrate this visualization technique, a sample of images are displayed in figure 2. In particular, figure 2(c) represents a typical image provided by Schlieren systems. Such systems are widely used in experimental fluid mechanics laboratories but to date no satisfying solution exists to analyze image sequences of this nature. Indeed, due to the absence of contrast, no image technique allowing a reliable quantitative evaluation of the visualized fluid flow motion is available. To the best of our knowledge, very few works [START_REF] Fu | Flow visualization using thye negativepositive grid schlieren and its image analysis[END_REF][START_REF] Kegerise | Schlieren image-correlation velocimetry and its application to free-convection flows[END_REF] have been carried out to estimate velocity fields from Schlieren images. All these works rely on correlation methods [START_REF] Settles | Schlieren and shadowgraph techniques: visualizing phenomena in transparent media[END_REF]. These methods suffer from several limitations. Among them, one can cite the fact that the results are sensitive to the size of the correlation window support and the possible lack of spatial coherence of the resulting displacement field. We believe that the use of dense motion estimation using optical flow is an interesting alternative that has not been investigated for the Schlieren imagery. These methods, formalized as the minimization of an energy function, have been already successfully implemented for general fluid flow imagery. We have adapted one of these methods to the Schlieren technique. Nevertheless, before describing the proposed dedicated Schlieren energy function, a brief overview of classical dense motion estimator is given in the next section.

3 Related works on dense motion estimation

Standard optical flow estimation

Dense estimation of the apparent motion aims at recovering a 2d displacement field w defined over the continuous plane domain S. The estimation is based on the knowledge of the luminance function at two consecutive instants denoted I(s, t), s ∈ S.

The most accurate techniques to address this problem are related to the Horn and Schunck (H&S) optical flow estimator [START_REF] Horn | Determining optical flow[END_REF][START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Mémin | Hierarchical estimation and segmentation of dense motion fields[END_REF][START_REF] Ruhnau | Variational optical flow estimation for particle image velocimetry[END_REF]. Such estimators are formalized as the minimizer of an energy function H composed of a data term H 1 and a regularization term H 2 . The first one describes a consistency assumption of the luminance function along a point trajectory. The standard brightness consistency assumption dI dt = 0 leads to consider the well-know optical flow constraint (ofc):

H 1 (w) = S φ 1 ∇I(s, t) . w(s) + ∂I(s, t) ∂t ds, (4)
where ∇I accounts for the spatial gradient of the luminance function and w(s) = (u(s), v(s)) T is the velocity at point s. The penalty function φ 1 is often chosen as the L 2 norm but better results may be obtained using a robust function that attenuates the effects of areas that do not respect the brightness assumption [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Mémin | Hierarchical estimation and segmentation of dense motion fields[END_REF].

The regularization term captures an a priori on the displacement field. A standard first-order spatial smoothness is usually considered:

H 2 (w) = α S φ 2 [ ∇w(s) ] ds, (5)
where ∇w(s) = ∇u(s) + ∇v(s) with an abuse of notation. Like φ 1 , the penalty function φ 2 authorizes handling local deviations from the smoothness model. The parameter α balances the relative influence of both terms in the functional.

Facing large frame-to-frame displacements, the data term H 1 is not anymore relevant due to its differential nature. To tackle this problem, the brightness consistency assumption has to be expressed in an integrated way, according to the displacement d(s) from time t to t + ∆t instead of the velocity w(s). As we have:

dI dt = lim ∆t→0 I(s + d(s), t + ∆t) -I(s, t) ∆t with d(s) = w(s) ∆t, (6)
by relaxing the constraint on the limit, the integrated version may be readily written as:

I(s + d(s), t + ∆t) -I(s, t) = 0. ( 7
)
To circumvent the high nonlinearity of this form with respect to the displacement field, the solution consists in proceeding to successive linearizations around an increment field dw. This is usually performed within a multiresolution scheme.

A first-order linearization of the first term of ( 7) yields to the following new energy function (where the time increment ∆t has been set to 1 for simplicity):

H(dw) = S φ 1 [∇I(s + w(s), t + 1) • dw(s) + I(s + w(s), t + 1) -I(s, t)] + α φ 2 [ ∇(w(s) + dw(s)) ] ds. (8)
For an interested reader, a state of the art of such techniques as well as their comparison can be found in [START_REF] Barron | Performance of optical flow techniques[END_REF][START_REF] Beauchemin | The computation of optical flow[END_REF].

Dense motion analysis in fluid imagery

As detailed in [START_REF] Corpetti | Dense estimation of fluid flows[END_REF], although estimators based on the energy function [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF] have been used for the velocity estimation of fluid structures, the two main assumptions involved in this function are not well suited to that specific case. First, the brightness consistency assumption involved in the data term is rarely valid for sequences of fluid flows. As a matter of fact, the observed luminance of a fluid structure may exhibits high spatio-temporal variations caused by temperature and pressure variations or due to its inherent deformable nature. The use of the fluid law of mass conservation (also called the continuity equation) as an alternative assumption applied to the evolution of the luminance function has originally been proposed in [START_REF] Schunk | The motion constraint equation for optical flow[END_REF]. Denoting v = (u, v, w) the 3d velocity, the continuity equation reads:

∂ρ ∂t + div(ρ v) = 0, (9)
where div(v) = ∂u ∂x + ∂v ∂y + ∂w ∂z denotes the divergence of the 3d velocity. Making a direct analogy between the density of a fluid particle and its luminance, this law has been integrated in some optical flow schemes [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Béréziat | A generalized optical flow constraint and its physical interpretation[END_REF][START_REF] Kohlberger | Variational dense motion estimation using the Helmholtz decomposition[END_REF]. Nevertheless, let us remark that apart from transmittance images [START_REF] Fitzpatrick | The existence of geometrical density-image transformations corresponding to object motions[END_REF], the use of the continuity equation remains an approximate constraint when applied to the image brightness. For Schlieren imagery we will show that an exact brightness variation model can be devised. This model will be detailed in the next section.

Secondly, concerning the regularization term, it can be demonstrated that a first order regularization is not adapted to fluid phenomena as it favors the estimation of velocity fields with low divergence and vorticity. A second order regularization can advantageously be consider as proposed in [START_REF] Suter | Motion estimation and vector splines[END_REF]:

H 2 (w) = α S ∇div(w(s)) 2 + ∇curl(w(s)) 2 ds, (10)
where div(w) = ∂u ∂x + ∂v ∂y and curl(w) = -∂u ∂y + ∂v ∂x are respectively the divergence and the vorticity of the 2d field w = (u, v). To circumvent the difficulty of implementing second order smoothness constraint, this regularization term can be simplified -in a computational point of view -in two interleaved first-order div-curl regularizations based on two auxiliary variables ξ 1 and ξ 2 approximating the divergence and the vorticity of the flow [START_REF] Corpetti | Dense estimation of fluid flows[END_REF]. Introducing the use of a robust penalty function instead of the quadratic function, we have:

H 2 (w, ξ 1 , ξ 2 ) = α S (div(w(s)) -ξ 1 ) 2 + β φ 2 [ ∇ξ 1 ] + (curl(w(s)) -ξ 2 ) 2 + β φ 2 [ ∇ξ 2 ] ds, ( 11
)
where β is a positive regularization parameter.

Dense estimator dedicated to Schlieren images

Data term

To construct a relevant dense motion estimator for Schlieren image sequences, it is essential to take into account the physical properties of this fluid visualization method. In particular, as previously described, the light pattern at time t is deduced from the density of the fluid (eq. ( 3)). In case of an almost 2d flow, introducing the time variable, we have:

I(s, t) = K∆z ∂ρ(x, y, z, t) ∂x + ∂ρ(x, y, z, t) ∂y . (12)
From that expression, we can deduce:

dI dt = K ∆z ∂ ∂x dρ dt + ∂ ∂y dρ dt A -∇ρ. ∂v ∂x -∇ρ. ∂v ∂y B ( 13
)
where v = (u, v, w) is the 3D velocity. This expression can be modified relying on the continuity equation ( 9) which can be alternatively rewritten after simple manipulations as:

dρ dt + ρ div(v) = 0. ( 14
)
This expression can be advantageously used in the first term A of equation [START_REF] Beauchemin | The computation of optical flow[END_REF].

Using expression [START_REF] Barron | Performance of optical flow techniques[END_REF] we have:

A = -ρ ∂ div(v) ∂x + ∂ div(v) ∂y -div(v) ∂ρ ∂x + ∂ρ ∂y (15) = -ρ ∂ div(v) ∂x + ∂ div(v) ∂y - I K ∆z div(v). ( 16
)
In order to simplify the second term B, let us assume that the two first components of the spatial gradient of the density are of the same order, i.e. ∂ρ ∂x ≈ ∂ρ ∂y with no local favored direction. This assumption does not necessary cancel the possibility that a global preferential direction for the pressure gradients may exist. It may be erroneous to associate directly the flow direction, or the favored pressure gradient direction, as the direction of the local fluid density gradients. Many flows of interest behave as incompressible flows and in these kinds of flows it can be admitted that the pressure gradients that drive the fluid flow may produce only negligible changes in the fluid density. The density field results in general from a complex interaction of the different coupled fields: temperature, pressure, buoyancy forces and velocity. As it is difficult to determine a priori a principal direction for the density gradients, it seemed to us reasonable to admit as a first approach that no direction for density gradients is preferential. Using this assumption, expression [START_REF] Barron | Performance of optical flow techniques[END_REF], and the fact that we are interested in this work on the dense motion estimation of mainly bidimensional fluid flows (i.e. inducing ∂ρ ∂z = 0), we have:

B = - I 2 K ∆z ∂u ∂x + ∂v ∂x + ∂u ∂y + ∂v ∂y . ( 17
)
As a 2D fluid flow is considered, we can also suppose that the apparent 2D velocity is defined by the two first components of the 3D velocity i.e. w = (u, v).

From that hypothesis, we can deduce that div(v) = div(w). Then:

A = -ρ ∂ div(w) ∂x + ∂ div(w) ∂y - I K ∆z div(w) (18)
and

B = - I 2 K ∆z div(w) + ∂v ∂x + ∂u ∂y . ( 19
)
The evolution in time of the luminance is then governed by the expression:

dI dt = - 3 2 I div(w) - 1 2 I ∂v ∂x + ∂u ∂y -ρ K ∆z ∂ div(w) ∂x + ∂ div(w) ∂y .
(20) Finally, as for most flows studied through a Schlieren system, it can be demonstrated that div(v) = 0, i.e. div(w) = 0, the resulting equation reads:

dI dt + 1 2 I ∂v ∂x + ∂u ∂y = 0. ( 21
)
In a similar manner as the standard optical flow estimation ( § 3.1), the expression ( 21) is not relevant for the estimation of large frame-to-frame displacements. An integrated version of this constraint has to be considered. Assuming that the velocity is constant between two instants t and t + ∆t, equation ( 21) is a first order differential equation at constant coefficient (equation of type y (t)m y(t) = p). Choosing I(s, t) as the initial condition, and setting the time interval ∆t to 1, the integrated for of the data model reads:

I(s + w(s), t + 1) exp 1 2 ∂ v(s) ∂x + 1 2 ∂ u(s) ∂y -I(s, t) = 0. ( 22
)
To cope with the non linearity of this constraint regarding to the displacement field, a coarse to fine strategy has to be settled. A first-order linearization of the left term in eq. ( 22) is considered with respect to an increment field dw. Removing the time index for sake of clarity and introducing the following notations I(.) = I(., t) ; I(.) = I(., t + 1), the Schlieren data term can be finally written as:

H 1 (dw) = S φ 1 -I(s) + exp(g(w(s)) I(s + w(s)) + (∇ I(s + w(s)) + ∇g(w(s)) I(s + w(s))) . dw(s)) ds, (23)
where g(w(s

)) = 1 2 ∂ v(s) ∂x + ∂ u(s) ∂y .

Regularization term

As for the regularization term, a second-order div-curl regularizer is considered as it enables the preservation of the fluid structures. To deal with the computational difficulties of second order smoothness functional implementation, the approach proposed in [START_REF] Corpetti | Dense estimation of fluid flows[END_REF] is followed. This leads to a regularization term already described by equation [START_REF] Ruhnau | Variational optical flow estimation for particle image velocimetry[END_REF]. Writing this expression in terms of a function of a velocity field increment to be minimized, we have:

H 2 (dw, ξ 1 , ξ 2 ) = α S (div(w(s) + dw(s)) -ξ 1 ) 2 + β φ 2 [ ∇ξ 1 ] + (curl(w(s) + dw(s)) -ξ 2 ) 2 + λ φ 2 [ ∇ξ 2 ] ds. (24)
This formulation has the very interesting property of allowing the introduction of an a priori information on the divergence and/or vorticity map. In particular, we have seen in the previous paragraph that in the studied experimental images, the divergence of the flow can be considered as null. Such a constraint has to be taken into account in the regularization term also. To that purpose, the term H 2 is modified to consider a constrained minimization implemented through a Lagrangian optimization technique. The new regularization reads:

H 2 (dw, ξ, λ) = α S (curl(w(s) + dw(s)) -ξ) 2 + β φ 2 [ ∇ξ ] + λ (div(w(s) + dw(s))) 2 ds, (25)
where λ denotes the Lagrangian multiplier associated to the constraint div(w(s)+ dw(s)) = 0.

Minimization issues

The incremental estimation of the dense displacement field is conducted through a multiresolution structure that consists in implementing an incremental estimation scheme on a pyramidal hierarchical representation of the image data. At a given resolution level, an incremental displacement field is computed considering that the main component of the displacement is known (supposed null at the coarsest level) and refined by solving:

min dw,ξ,λ H 1 (dw) + H 2 (dw, ξ, λ) (26)
where H 1 and H 2 are defined by equations (23,25). The minimization of the functional is considered through a direct discretization of H 1 and H 2 . The different functions involved in the functional are discretized on the image lattice. A particular attention has been paid for the discretization of divergence and curl operator for which an uncentered discretization scheme has been used. The overall system is constituted by two main sets of variables that have to be estimated. The first one is the motion field w, and the second set comprises the scalar field ξ. The estimation is conducted alternatively by minimizing H 1 + H 2 with respect to dw, λ and ξ respectively. For the motion field, considering the curl estimate ξ as being fixed, the robust minimization with respect to dw is solved with an iteratively re-weighted least squares technique. This optimization is embedded in an efficient multi-parametric adaptive multigrid framework [START_REF] Mémin | Hierarchical estimation and segmentation of dense motion fields[END_REF]. In turn, the motion field dw being fixed, the minimization of the cost function with respect to ξ is in fact equivalent to the minimization of H 2 and is again conducted using an iteratively re-weighted least squares technique.

Experimental results

In this section, experimental results are presented to highlight the relevance of our estimator. Two image sequences are studied. They both have been obtained in a laboratory of fluid mechanics 4 . The first experiment corresponds to a natural convection of a cylinder test in air at rest and the second one corresponds to a forced convection test of a heated cylinder immersed in a free airstream at room temperature. As it can be noticed on figures 3, 5, the obtained images are very difficult to analyze due to low brightness contrasts. It clearly appears that generic methods based on the brightness consistency assumption (correlation approaches, H&S methods) are hardly suited to these images.

The images obtained from the first experiment are shown on fig. 3, as well as the sequence of motion fields and vorticity maps obtained by a dense opticalflow estimator [START_REF] Corpetti | Dense estimation of fluid flows[END_REF]. The difficulties of this sequence lie in the lack of luminance variations and in the large frame-to-frame displacements of the fluid structures. As it can be noticed on the vorticity maps, the emergence of a vortex has been well captured by our estimator, as well as the smaller structures. This result proves the validity of the Schlieren dedicated data term. The impact of the new regularization term (that forces the estimation of a flow with null divergence) is demonstrated on fig. 4. This figure presents a comparison between our method and an optical-flow estimator proposed in [START_REF] Corpetti | Dense estimation of fluid flows[END_REF]. As it can be noticed on the presented displacement fields, this latter generates a motion field with areas of high divergence that are not physically plausible.

The results obtained on the second experiment are shown on fig. 5. These results are displayed in terms of vorticity maps. These pictures show that the moving vertical structures of the fluid flows have been well recovered. We can see in particular the coherent displacement of the lower vortex and the vanishing due to dissipation of the upper vortices. The curl maps also highlight the temporal consistency of the recovered motion fields.

Conclusion

In this paper, we have presented a new method for the estimation of dense fluid motion fields dedicated to images obtained with a Schlieren system. The analysis of the Schlieren images is of great importance in the field of fluid mechanics since this system enables the visualization of unseeded flows. The proposed method is a minimization-based approach where the two terms involved in the cost function have been designed for these images. In particular, the data term has been deduced from the physical relation between the luminance function and the fluid density gradient. The very promising results have demonstrated the interest of using such an approach for the Schlieren image analysis. The following planned step is to validate our approach considering synthetic images produced by a 4 The images have been obtained with a Schlieren system disposed in a Z configuration.

It comprised two spherical mirrors of 35 cm in diameter and the light was cut off with two razor blades disposed in vertical and horizontal positions, thus density gradients in both directions could be detected. The parallel light rays traversed the test section of a low speed wind tunnel with windows in the test section of optical quality to avoid improper light deflections. The images were recorded on a monochromatic digital image camera of 12 bits that enabled fast frame acquisition.

Fig. 1 .

 1 Fig.1. Typical Schlieren system using lenses -figure from[START_REF] Goldstein | Optical systems for flow measurement: Shadowgraph, schlieren, and interferometric techniques[END_REF]

Fig. 2 .

 2 Fig.2. Example of images obtained from a Schlieren system: (a) human thermal plume (image from[START_REF] Settles | Schlieren and shadowgraph techniques: visualizing phenomena in transparent media[END_REF]) ; (b) instantaneous image of bullet and blast (image from[START_REF] Settles | Schlieren and shadowgraph techniques: visualizing phenomena in transparent media[END_REF]) ; air flow ; visualized flames with a color Schlieren system (here displayed in black and white).

Fig. 3 .Fig. 4 .

 34 Fig. 3. natural convection. sequence of (a) images, (b) estimated vorticity maps and (c) estimated displacement fields

Fig. 5 .

 5 Fig. 5. forced convection. images 1, 5, 10, 15, 20 and associated estimated vorticity maps

Direct Numerical Simulation code. From this work, several perspectives can be investigated such as the study of 3D flows (using for example the Schlieren tomography [20]), and the design of dedicated algorithms to track the fluid structures.


Fabio Cuzzolin
email: cuzzolin@cs.ucla.edu

The geometry of relative plausibilities

Keywords: Belief space, Bayesian relatives, relative plausibility, equidistribution, orthogonality

The study of the interplay between belief and probability has recently been posed in a geometric setup, in which belief functions are represented as points of a simplex called belief space. Each belief function b is associated with three different geometric entities, which in turn determine three different Bayesian function related to b. In this paper, in particular, we study the geometry of relative plausibilities, and show how equidistribution can be possibly a trait d'union among all the Bayesian relatives of a belief function.

Introduction

In the last decades a number of different uncertainty measures have been proposed, as alternatives or extensions of the classical probability theory. The theory of evidence is one the most popular formalism, extending quite naturally probabilities on finite spaces through the notion of belief function. Naturally enough, the connection between belief and probability plays a major role in the theory of evidence, and is the basement of a popular approach to evidential reasoning, Smets' pignistic model [START_REF] Smets | Belief functions versus probability functions[END_REF]. In fact, the problem of finding a correct probabilistic approximation of a belief function has been widely studied [START_REF] Kramosil | Approximations of believeability functions under incomplete identification of sets of compatible states[END_REF][START_REF] Bauer | Approximation algorithms and decision making in the Dempster-Shafer theory of evidence-an empirical study[END_REF], mainly in order to find efficient implementations of the rule of combination, as in Tessem [START_REF] Tessem | Approximations for efficient computation in the theory of evidence[END_REF], or Lowrance et al. [START_REF] Lowrance | A framework for evidential-reasoning systems[END_REF]. Voorbraak [START_REF] Voorbraak | A computationally efficient approximation of Dempster-Shafer theory[END_REF], on his side, proposed to adopt the so-called relative plausibility function pl b , i.e. the unique probability that, given a belief function b with plausibility pl b , assigns to each singleton its normalized plausibility. He proved that pl b is a perfect representative of b when combined with other probabilities, pl b ⊕ p = b ⊕ p ∀p ∈ P. Cobb and Shenoy [START_REF] Cobb | On transforming belief function models to probability models[END_REF] described the properties of the relative plausibility of singletons and discussed its nature of probability function which is equivalent to the original belief function. The study of the interplay between belief and probability has recently been posed in a geometric setup [START_REF] Ha | Theoretical foundations for abstraction-based probabilistic planning[END_REF][START_REF] Black | Geometric structure of lower probabilities, Random Sets: Theory and Applications[END_REF], in which belief functions are represented as points of a simplex called belief space [START_REF] Cuzzolin | Geometric analysis of belief space and conditional subspaces[END_REF]. Each belief function is associated with three different geometric entities, namely the line (b, pl b ), the orthogonal complement P ⊥ of the probabilistic subspace P, and the simplex of consistent probabilities P [b] = {p ∈ P : p(A) ≥ b(A) ∀A ⊂ Θ}. These in turn determine three different probabilities associated with b, i.e. the orthogonal projection π[b] of b onto P, the barycenter of P [b], and the pignistic function BetP [b]. We recently showed that if the belief function b is 2-additive then pignistic function, orthogonal projection and intersection of (b, pl b ) with P coincide [START_REF] Cuzzolin | Geometric interplays of belief and probability[END_REF]. However, the analysis of the simplest binary case suggests that the relative plausibility function does not fit in this picture. In this paper we will study the geometry of relative plausibilities in the general case, pointing out its richness and complexity. We will see that their geometry can be described by three planes in the belief space, in turn related to an additional set of three Bayesian functions. As 2-additivity is not able to encompass both the above Bayesian functions and Voorbraak's relative plausibility function, we show that equidistribution can be a trait d'union for all the "Bayesian relatives" of a belief function.

Geometry of belief functions

In the theory of evidence [START_REF] Shafer | A mathematical theory of evidence[END_REF] degrees of (subjective) belief are represented as belief functions. A basic probability assignment (b.p.a.) over a finite set Θ (frame of discernment) is a function

m : 2 Θ → [0, 1] such that m(∅) = 0, A⊂Θ m(A) = 1, m(A) ≥ 0 ∀A ⊂ Θ.
(A) . = 1 -b(A c ) = B∩A =∅ m(B) ≥ b(A).
From pl b we can derive a Bayesian belief function called relative plausibility of singletons [START_REF] Voorbraak | A computationally efficient approximation of Dempster-Shafer theory[END_REF] as follows:

pl b : Θ → [0, 1] x → pl b (x) = pl b (x) y∈Θ pl b (y) . ( 1
)
Consider now a frame of discernment Θ and introduce in the Euclidean space R N , N = |2 Θ | -1, an orthonormal reference frame {x A } A⊂Θ,A =∅ in which each coordinate x A measures the belief value b(A). The belief space associated with Θ is then the set of points B Θ of R N which correspond to a belief function. We assume the domain Θ fixed, and denote the belief space with B. It is not difficult to prove [START_REF] Cuzzolin | Geometric analysis of belief space and conditional subspaces[END_REF] that B is convex. More precisely, if we call

b A . = b ∈ B s.t. m(A) = 1, m(B) = 0 ∀B = A
the unique belief function assigning all the mass to a single subset A of Θ, the belief space B is the convex closure Cl of all the basis belief functions b A :

B = Cl(b A , A ⊂ Θ, A = ∅). ( 2
)
Furthermore, each belief function b ∈ B can be written as a convex sum as follows:

b = A⊂Θ, A =∅ m(A)b A . ( 3
)
Since a probability is a belief function assigning non zero masses to singletons only, the set P of all the Bayesian belief functions is the simplex determined by all the basis functions associated with singletons P = Cl(b {x} , x ∈ Θ). Analogously, we call plausibility space the region PL of R N whose points correspond to admissible plausibility functions. It can be proved [START_REF] Cuzzolin | Geometry of upper probabilities[END_REF] that this is also a simplex PL = Cl(pl A , A ⊂ Θ), whose vertices are

pl A = - B⊂A (-1) |B| b B . ( 4
)
The A-th vertex pl A of the plausibility space is the plausibility vector associated with the basis belief function b A , pl A = pl b A . Each plausibility function pl b can be uniquely expressed as a convex combination of basis plausibility functions (4) as follows

pl b = A⊂Θ m(A)pl A . ( 5
)
Belief functions are then nothing but special points in an Euclidean space. On the other side, given an ordering A 1 , .., A N for the subsets of Θ, each point x = [x 1 , ..., x N ] ∈ R N (where x denotes the transpose of a vector x), N . = 2 |Θ| -1, can be thought of as a function ς : 2 Θ \ ∅ → R s.t. ς(A) = x ord(A) , where ord(A) is the position of A in the ordering. As the Moebius transformation is invertible, for each of these functions ς there always exists another function m ς : 2 Θ \ ∅ → R such that ς(A) = B⊂A m ς (B). Each vector ς of R N can be thought of as a sum function [START_REF] Aigner | Combinatorial theory[END_REF], even though m ς does not in general satisfy the positivity constraint m ς (A) ≥ 0 ∀ A ⊂ Θ. The section x ord(Θ) = 1 of R N corresponds to the constraint ς(Θ) = 1, so that points of this section are sum functions meeting the normalization axiom, A⊂Θ m ς (A) = 1. We call them normalized sum functions (n.s.f.). Analogously to the case of belief functions, we call Bayesian normalized sum function a n.s.f.

ς such that x∈Θ m ς ({x}) = 1. Note that this implies |A|>1 m ς (A) = 0, but not necessarily m ς (A) = 0 ∀|A| > 1.

A geometric interplay of belief and probability

As b(Θ) = pl b (Θ) = 1 for all b.f. b, we can neglect the coordinate x Θ and think of B as a region of R N -1 , in which b Θ = 0 = [0, ..., 0] . Figure 1 shows the geometry of belief and plausibility spaces for a binary frame Θ 2 = {x, y}, where belief and plausibility vectors are points of a plane with coordinates

b = [b(x) = m(x), b(y) = m(y)] pl b = [pl b (x) = 1 -m(y), pl b (y) = 1 -m(x)]
respectively. These two simplices are symmetric with respect to P. Furthermore, each pair of functions (b, pl b ) determines a line which is orthogonal to P, on which they lay on symmetric positions on the two sides of the Bayesian region. The orthogonal projection of b onto P, the pignistic function BetP plement P ⊥ of P. In [START_REF] Cuzzolin | Geometric interplays of belief and probability[END_REF] we proved that the line (b, pl b ) is always orthogonal to P, even in the general case. However, (b, pl b ) does not intersect the probabilistic subspace in general, but there exists a unique Bayesian normalized sum function ς

[b] such that ς[b] = β[b] • pl b + (1 -β[b]) • b (6)
where

β[b] = |A|>1 m(A) |A|>1 m(A)|A| (7)
is a scalar value between 0 and 1 which depends on b. ς[b]'s values on singletons are

m ς[b] (x) = m(x) + β[b] • A⊃x,A =x m(A). (8)
Hence, ς[b] is naturally associated with a Bayesian belief function, assigning the same mass [START_REF] Ha | Theoretical foundations for abstraction-based probabilistic planning[END_REF] to each singleton, namely

σ[b] . = x∈Θ m ς[b] (x) • b {x} . ( 9
)
In [START_REF] Cuzzolin | Geometric interplays of belief and probability[END_REF] we also derived the convex shape of However, it is clear from Figure 1 that Voorbraak's relative plausibility function does not fit in this picture, as 2-additivity is then not general enough to encompass both pl b and the above Bayesian functions. In the following sections we will study its geometry: for sake of simplicity, we will drop the dependence on b from the notations for σ, ς, β.

Geometry of relative plausibility

= x m(x)b {x} +β x (pl b (x)-m(x))b {x} = (1-β) x m(x)b {x} +β x pl b (x)b {x} so that, after defining b = x∈Θ m(x)b {x} = x∈Θ m(x)b {x} +(1-k 2 )b Θ , (11) with k 2 . = |A|>1 m(A), Eq. ( 10
) implies σ = (1 -β) • b + β • pl b ( 12
)
where b ∈ Cl(b Θ , P) . = D and pl b ∈ span(D), the affine (vector) space generated by the simplex D. D is the simplex of discounted probabilities [START_REF] Shafer | A mathematical theory of evidence[END_REF], i.e. belief functions assigning mass to singletons or Θ only. b is in fact the unique belief function assigning to Θ all the mass (1 -k 2 ) which b gives to non-singletons. In the general case, then, the functions on the line (b, pl b ) are associated with other quantities laying on the line ( b, pl b ) in the affine subspace span(D) of discounted probabilities b ↔ b, pl b ↔ pl b , ς ↔ σ which are in the same relative positions on the segment (just compare ( 6) and ( 12)). Again, 2-additivity plays a role, since [START_REF] Cuzzolin | Geometric interplays of belief and probability[END_REF]

Theorem 1. ς[b] = σ[b] iff b is 2-additive, in which case σ[b] is the orthogonal projection π[b] of b onto P.
= x∈Θ m(x)pl {x} + (1 -k 2 )pl Θ = = b + (1 -k 2 )pl Θ pl b = x∈Θ pl b (x)pl {x} + (1 -k 1 )pl Θ = = pl b + (1 -k 1 )pl Θ . ( 13
)
Theorem 2. The line passing through the duals of plausibility of singletons [START_REF] Lowrance | A framework for evidential-reasoning systems[END_REF] and belief of singletons [START_REF] Shafer | A mathematical theory of evidence[END_REF]

)⊥P (φ 1 [b] = π/2) iff ∀x A⊃x,A ={x} m(A) = pl b (x) -m(x) = const.
In this case, recalling the definition of σ (9),

σ(x) = m(x) + β(pl b (x) -m(x)) = m(x) + 1-k 2 y (pl b (y)-m(y)) (pl b (x) -m(x)) = m(x) + 1-k 2 n
and σ assigns the mass originally given by b to non-singletons to all singletons on equal basis. This can be expressed in terms of another probability function

R[b] R[b] . = pl b - b k 1 -k 2 = x pl b (x) -m(x) k 1 -k 2 b {x} . ( 18
(x) -m(x) = A⊃x,A =x m(A) = k 1 -k 2 n for all x.
Replacing this in (18) yields

R[b] = x∈Θ 1 n b {x} .
The value of φ 2 also depends on this mysterious probability (18).

Lemma 1. φ 2 is nil if and only if

1, R[b] 2 = R[b] 2 1, 1 . ( 19
)
where

pl Θ = 1 is the N -dim vector [1, .., 1] . Theorem 4. The angle φ 2 [b] is zero if and only if R[b] is parallel to pl Θ = 1.

Condition (19) has the form

A, B 2 = A 2 B 2 cos 2 ( AB) = A 2 B 2 i.e. cos 2 ( AB) = 1, with A = pl Θ , B = R[b]. This yields cos( R[b]pl Θ ) = 1 or cos( R[b]pl Θ ) = -1, but as both
) never coincide ∀b ∈ B Θ when |Θ| > 2; instead φ 2 [b] = 0 ∀ b ∈ B Θ when |Θ| ≤ 2. Σ b R[b] Θ pl b P /c 1-Σ/ c - - - b b - b pl b pl - b Figure 3: Location of R[b] in P.
The binary case is then really anomalous, a situation signaled by R[b] being parallel to the vacuous plausibility pl Θ . R[b] has in fact a simple geometric interpretation, as Equation (18 Again, a necessary and sufficient condition for The quantity M k (x) is then somehow connected to the geometric condition of orthogonality in the belief space framework. In [START_REF] Cuzzolin | Geometric interplays of belief and probability[END_REF] we have shown that a belief function b is orthogonal to the probabilistic subspace P iff In fact, condition ( 22) is equivalent to

) implies R[b](k 1 -k 2 ) = pl b -b = pl b • k 1 -b • k 2 = pl b • k 1 -b • k 2 + k 1 • b -k 1 • b = k 1 • ( pl b -b) + b(k 1 -k 2 ) from which R[b] = b + k 1 k 1 -k 2 ( pl b -b). (20
φ 3 [b] = 0 can expressed in terms of R[b], as R[b](x) = (pl b (x) -m(x))/(k 1 -k 2 ) = 1 k 1 -k 2 ( k 1 k 2 m(x) -m(x)) = m(x)/k 2 ∀x, i.e.
B⊃y m(B)2 1-|B| = B⊃x m(B)2 1-|B| ≡ n-1 k=1 1 2 k |B|=k,B⊃y m(B) = n-1 k=1 1 2 k |B|=k,B⊃x m(B) ≡ n-1 k=1 1 2 k M k (y) = n-1 k=1 1 2 k M k (x)
for all y = x, and if M k (x) = M k (y) ∀y = x the condition is met. In this case, confirming the intuition given by the binary case (Figure 1), all the Bayesian relatives of b converge to the same probability: Equidistribution is then inherently linked to the geometric notion of orthogonality, and provides a unifying viewpoint for all the geometric entities we encountered in the study of the Bayesian relatives of b: the orthogonal complement of P, the probabilities P [b] consistent with b, and the dual lines (b, pl b ) and ( b, pl b ). The next natural step will be the study of the behavior of all these function under the rule of combination.

(x) = A⊃x m(A) = n k=1 A⊃x,|A|=k m(A) = k M k (x) we get pl b (x) = pl b (x) k 1 = n k=1 M k (x) x n k=1 M k (x) which is equal to 1/n if M k (x) = M k ∀k, x.

Discussion

Appendix

Proof. (Theorem 2) From Equation ( 13) we have that b

-b = (1 -k 2 )pl Θ and pl b -pl b = (1 -k 1 )pl Θ . Hence, β( pl b -b) + b = β[ pl b + (1-k 1 )pl Θ -b-(1-k 2 )pl Θ ]+ b+(1-k 2 )pl Θ = β[ pl b -b + (k 2 -k 1 )pl Θ ] + b + (1 -k 2 )pl Θ = b + β( pl b -b) + pl Θ [β(k 2 -k 1 ) + 1 -k 2 ] but β(k 2 -k 1 )+1-k 2 = 1-k 2 k 2 -k 1 (k 2 -k 1 )+1-k 2 =
-1 = (k 1 -1) 2 ( R 2 -1, R ) (k 1 -1) 2 R] R, R + 1, 1 -2 R, 1
that is equivalent to (after elevating to the square both numerator and denominator) 2 and by erasing the common terms we have as desired.

R[b] 2 ( R[b] 2 + 1, 1 -2 R[b], 1 ) = R[b] 4 + 1, R[b] 2 -2 1, R[b] R[b]

  Subsets of Θ associated with non-zero values of m are called focal elements (f.e.). The belief function (b.f.) b : 2 Θ → [0, 1] associated with a basic probability assignment m is simply b(A) = B⊂A m(B). The basic probability assignment m of a belief function b can be uniquely recovered by means of the Moebius inversion formula m(A) = B⊂A (-1) |A-B| b(B). In particular, a probability function is a peculiar belief function called Bayesian b.f. assigning mass to singletons only: m(A) = 0, |A| > 1. A dual representation of the evidence encoded by a belief function b is the plausibility function, whose value pl b (A) expresses the amount of evidence not against a proposition A: pl b

Figure 1 :

 1 Figure 1: Geometry of the relative plausibility function in the bi-dimensional case Θ = {x, y}.

Proposition 1 .

 1 P [b], and proved how the center of mass of this simplex P [b] is in fact the pignistic function. CITARE REFERENZA The three geometric loci which arise from the analysis of the binary case are in conclusion each associated with a different Bayesian b.f. (b, pl b ) ↔ σ[b], P ⊥ ↔ π[b], P [b] ↔ BetP [b] = P [b]. The trait d'union between those Bayesian functions is again suggested by the binary case [5]. If the belief function b is 2additive (m(A) = 0, |A| > 2) then pignistic function, orthogonal projection onto P and intersection of the line (b, pl b ) with P coincide, BetP [b] = P [b] = π[b] = ς[b] = σ[b].

From

  the definition of relative plausibilities (1) it follows a straightforward geometric interpretation of pl b . Let us call plausibility of singletons the quantity pl b . = x∈Θ pl b (x) • b {x} . (10) pl b can be written as pl b = x pl b (x) • b {x} + (1 -k 1 ) • b Θ with k 1 . = x pl b (x) since b Θ = 0, in order for it to meet the normalization constraint. As 1 -k 1 ≤ 0 pl b is a normalized sum function. Now, by definition (1) pl b = pl b /k 1 . Since b Θ is the origin of the reference frame, pl b lies on the segment ( pl b , b Θ ). This in turn implies that pl b = ( pl b , b Θ ) ∩ P i.e. pl b is the intersection of the line joining the vacuous belief function b Θ and the plausibility of singletons pl b with the probabilistic subspace. From the definition of σ (9) it follows that σ

4. 1 A

 1 three plane geometry In the binary case D = B, b = b and pl b = pl b , so that the plausibility of singletons is a plausibility function (Figure 1). In the general case, though, D = Cl(b Θ , P) is in correspondence with a "dual" simplex D ⊂ PL, namely the set of plausibilities associated with a discounted probability. As pl Θ is the dual of b Θ , and each probability b {x} = pl {x} is the dual of itself, D = Cl(pl Θ , P). The functions b, pl b ∈ D then correspond to two dual quantities in D (recalling Equation (5)) b

  crosses σ[b] too, and β( pl b -b) + b = σ[b] = β( pl b -b) + b. (14) The geometry of the relative plausibility function can then be described in terms of the three planes ( pl b , σ, pl b ), (b Θ , pl b , pl Θ ) and (b Θ , b, pl Θ ) (see Figure 2), where b . = b/k 2 (15) is another probability which is natural to call relative belief of singletons. More precisely, we have just seen how σ is the intersection of both ( b, pl b ) and ( b, pl b ) and lays in the same relative position on the two lines. They then determine a plane ( b, σ, b) = ( pl b , σ, pl b ). Furthermore, by definition pl b -b Θ = ( pl bb Θ )/k 1 while (13) implies pl b = pl b /k 1 = ( pl b -(1 -k 1 )pl Θ )/k 1 so that pl b -pl Θ = ( pl b -pl Θ )/k 1 . The relative plausibility function then lies in the same relative position on the two lines (b Θ , pl b ) and (pl Θ , pl b ), which intersect exactly in pl b . b Θ , pl Θ , pl b , pl b and pl b then determine another unique plane which we can denote with (b Θ , pl b , pl Θ ). Analogously, by definition b -b Θ = ( b -b Θ )/k 2 while (13) yields b -pl Θ = ( b -pl Θ )/k 2 , so that the relative belief of singletons lies in the same relative position on the two lines (b Θ , b) and (pl Θ , b), which intersect exactly in b. b Θ , pl Θ , b, b and b then determine a single plane denoted by (b Θ , b, pl Θ ).

Figure 2 :

 2 Figure 2: Planes and angles describing the geometry of the relative plausibility function.

) Corollary 1 .

 1 The dual line ( b, pl b ) is orthogonal to P iff the Bayesian b.f. R[b] is the uninformative probability P . = x b {x} /n. In fact x∈Θ A⊃x,A =x m(A) = x∈Θ (pl b (x) -m(x)) = k 1 -k 2 so that the condition of Theorem 3 can be written as pl b

  R[b] and pl Θ have all positive components it must be cos( R[b]pl Θ ) = 1 so that R[b]pl Θ = 0. In this case the two lines ( b, pl b ) and ( b, pl b ) coincide. However, R[b] pl Θ means R[b] = α • pl Θ for some scalar α, i.e. R[b] = -α • A⊂Θ (-1) |A| b A (since pl Θ = -A⊂Θ (-1) |A| b A by Equation (4)). On the other side, R[b] = x∈Θ R[b](x)b {x} is a probability and the two conditions are never met together, unless |Θ| = 2. In conclusion Corollary 2. φ 2 [b] = 0 and the lines ( b, pl b ), ( b, pl b

  ) R[b] lies on the line joining b and pl b , outside the segment. Figure 3 illustrates how its location depends on the ratio between the total belief of singletons k 2 and the total plausibility of singletons k 1 . Finally, the angle φ 3 [b] between b and pl b nullifies iff b = pl b , which is equivalent to m(x)/k 2 = pl b (x)/k 1 ∀x ∈ Θ. In other words, relative plausibility and belief coincide iff mass and plausibility are distributed in the same way among the singletons.

R

  [b] = b, with R[b] "squashing" pl b onto b from the outside. In this case the quantities pl b , pl b , pl b , σ[b], b, b, b all lie in the same plane. 5 Equidistribution, orthogonality and Bayesianity This geometric analysis brings to full light the existence of several Bayesian functions linked to the original belief function by means of the balance between the overall plausibility k 1 of the elements of the frame, and the total mass k 2 assigned to them. This flag of the "non-Bayesianity" of b is symbolized by the Bayesian function R[b]. Those functions must be added to the three other probabilities discussed in Section 3. However, 2-additivity falls short of unifying this extended set of Bayesian relatives of b. The binary case again provides us with an intuition of how equidistribution can provide such a common ground. From Figure 1, we can appreciate how belief functions with m(x) = m(y) lay on the bisector of the first quadrant which is orthogonal to P. Their relative plausibility is also equal to their orthogonal projection π[b]. In fact, Theorem 3 can be interpreted in terms of equal distribution of mass among focal elements. Consider the quantityM k (x) of size k subsets to the plausibility of x ∈ Θ. Suppose that M k (x) = const = M k ∀x ∈ Θ for all k = 2, ..., n -1 (trivially true for k = n). Then A⊃x,A =x m(A) = n k=2 |A|=k,A⊃x m(A) = m(Θ) + n-1 k=2 M k which is constant ∀x ∈ Θ.In other words, Corollary 3. If M k (x) = const for all x ∈ Θ, ∀k = 2, ..., n -1 then the dual line ( b, pl b ) is orthogonal to P. It is worth to notice that [5] Proposition 2. If the belief function b is such that for all x ∈ Θ, k = 3, .., n -1 M k (x) = const = M k then pignistic function P [b] and orthogonal projection coincide, P [b] = π[b].

  B⊃y,B ⊃x m(B)2 1-|B| = B⊃x,B ⊃y m(B)2 1-|B| . (22) We can use (22) to prove that Theorem 5. If a belief function b is such that M k (x) = const = M k for all k = 1, ..., n -1 then b is orthogonal to the probabilistic subspace, b⊥P.

Theorem 6 .

 6 If M k (x) = const = M k for all k = 1, ..., n -1 then pl b = π[b] = P [b] = P. To see this let us write the values of the pignistic function P [b](x) in terms of M k (x) as under the hypothesis, yielding P [b] = P. But now, as pl b

From

  the geometrical point of view, each belief function is associated with two different families of Bayesian functions (the Bayesian relatives of b). A first group of probabilities σ[b], π[b], P [b] = BetP [b] coincide under the assumption of 2-additivity. From the geometric study of the relative plausibility of singletons a second family of probabilities b, pl b , R[b] arises which does not meet this condition. The geometry of pl b can be described in terms of three planes and angles, which depend on the balance between total belief and total plausibility of singletons. Equidistribution can provide a common ground for the two families. If subsets of the same size contributes equally to the plausibility of each singleton (M k (x) = const), k = 3, ..., n → P [b] = π[b] k = 2, ..., n → ( b, pl b )⊥P k = 1, ..., n → b⊥P, pl b = P [b] = π[b].

1-k 2 k 1

 21 0 by definition of β[b][START_REF] Cuzzolin | Geometric analysis of belief space and conditional subspaces[END_REF], and[START_REF] Voorbraak | A computationally efficient approximation of Dempster-Shafer theory[END_REF] is met. b -b, b {y} -b {x} = z∈Θ (pl b (z) -m(z))b {z} , b {y} -b {x} = z∈Θ (pl b (z) -m(z)) • [ b {z} , b {y}b {z} , b {x} ] = z∈Θ (pl b (z) -m(z)) • [ b {z}∪{y} , b {z}∪{y} -b {z}∪{x} , b {z}∪{x} ] since it is easy to see by the definition of b A that b A , b B = b A∪B , b A∪B . We distinguish three cases: if z = x, y then |z ∪ x| = |z ∪ y| = 2 and the difference b z∪x 2 -b z∪y 2 goes to zero. If z = x then b z∪x 2 -b z∪y 2 = b x 2 -b x∪y 2 = 2 n-2 -1 -(2 n-1 -1) = -2 n-2 , while if z = y then b z∪x 2 -b z∪y 2 = b x∪y 2 -b y 2 = 2 n-2 . Hence (17) is equal to 2 n-2 • (pl b (y) -m(y)) -2 n-2 (pl b (x) -m(x))and the thesis follows asA⊃x,A =x m(A) = pl b (x) -m(x). Proof. (Lemma 1) After noticing that β[b] = -k 2 we can write pl b -σ = ( pl b -pl b )+( pl bσ) = (1 -k 1 )pl Θ + (1 -β)( pl b -b) = k 1 -1 k 1 -k 2 • ( pl b -b)-(k 1 -1)pl Θ , so that pl b -σ, pl b -σ = pl b -σ, pl b -σ +(1-k 1 ) pl Θ , pl b -σ . But now pl b -σ, pl b -σ = (1-β) 2 pl b -b, pl b -b , being pl b -σ = (1-β)( pl b -b), while pl Θ , pl b -σ = (1 -β) pl Θ , pl b -b where 1 -β = k 1 -1 k 1 -k 2 , so that by definition of R[b] (18) we have pl bσ, pl b -σ = (k 1 -1) 2 R[b], R[b] and (1k 1 ) pl Θ , pl b -σ = -(k 1 -1) 2 1, R[b] . The desired scalar product is then pl b -σ, pl b -σ = (k 1 -1) 2 ( R[b], R[b] -1, R[b] ). Now, cos(π -φ 2 ) = pl b -σ, pl b -σ pl b -σ pl b -σ where pl b -σ = [ pl b -σ, pl b -σ ] 1/2 = (k 1 -1)[ R[b] -pl Θ , R[b] -pl Θ ] 1/2 which is equal to (k 1 -1)[ R[b], R[b] + pl Θ , pl Θ -2 R[b], pl Θ ] 1/2 and pl b -σ = (1 -β)( pl bb) = (k 1 -1) R[b] . So, as φ 2 [b] = 0 iff cos(π -φ 2 [b]) = -1 we can write the condition as


Francesco Amigoni

Simone Gasparini

Maria Gini

Good Experimental Methodologies for Robotic Mapping: A Proposal

A way to significantly advance robotic science is to perform experiments that can be replicated by other researchers and be used to compare different methods. This happens rarely in current robotics research. In this paper we present a methodology for performing experimental activities in the area of robotic mapping. The proposed methodology prescribes a number of issues that should be addressed when experimentally validating a mapping method. We present the application of the proposed methodology to a mapping system we have developed.

I. INTRODUCTION

Efforts are being made to establish standard benchmarks for mobile robots [START_REF] Euron | Survey and inventory of current efforts in comparative robotics research[END_REF], [START_REF]Intelligent systems division -performance metrics and test arenas for autonomous mobile robots[END_REF]. However, benchmarks only allow to compare different robotic systems by considering them as black boxes, namely by comparing their "external" performance. In our opinion, this is not enough to advance significantly the science of robotics. Beyond, and before, benchmarks, what is needed is to give researchers the ability to repeat and to cross-check the experiments done by others. The repeatability and controllability of experiments are distinctive features of any scientific discipline. However, these issues are not sufficiently addressed by current experimental practices in robotics.

In this paper we propose a methodology for performing experimental activities in the area of robotic mapping [START_REF] Thrun | Robotic mapping: A survey[END_REF]. More specifically, we concentrate on mapping methods that operate on segment-based maps. The proposed methodology prescribes a number of issues that, when addressed in the experimental validation of a mapping method, will enable the replication and cross-checking of experiments and the comparison with other methods. We also present the application of the proposed methodology to a specific mapping system we have developed in previous work [START_REF] Amigoni | Building segment-based maps without pose information[END_REF].

The basic idea from which we start is that fully replicable and controllable experiments can enable a widespread comparison of different mapping systems. The methodology presented in this paper is still preliminary. However, we believe it is a first step toward the final goal to foster the good practice of replicating experiments in order to compare different methods and assess their strengths and weaknesses, as usually happens in other scientific disciplines.

The original contributions of this paper, whose nature is more methodological than technical, are (a) a critical survey of the current practices in experimental activities for segment-based mapping methods and (b) the proposal of a F. Amigoni and S. Gasparini are with the Dipartimento di Elettronica e Informazione of the Politecnico di Milano, Milano, Italy, {amigoni,gasparini}@elet.polimi.it.

M. Gini is with the Department of Computer Science and Engineering of the University of Minnesota, Minneapolis, USA, gini@cs.umn.edu.

preliminary methodology for performing these experimental activities. To the best of our knowledge, we are not aware of any other similar effort in the area of robotic mapping.

This paper is structured as follows. In the next Section, we overview the current practice in performing experiments for segment-based mapping methods. In Section III we describe the proposed methodology, which we then apply to our mapping system in Section IV.

II. STATE OF THE ART OF EXPERIMENTS FOR ROBOTIC MAPPING SYSTEMS

We concentrate on mapping systems that build maps composed of 2D line segments. Usually, these systems work by incrementally integrating scans acquired in the environment. For this reason, they are usually based on scan matching methods. Scans are the result of perception actions performed by robots exploiting their sensors (e.g., laser range scanners). Scans are usually represented by line segments obtained by post-processing the sensor data.

In this Section we overview how the experimental activities aimed at validating scan matching methods have been conducted and reported in the literature. Without attempting to be exhaustive, we selected the major methods for building segment-based maps that have been proposed in the last few years. While algorithms are usually well detailed and explained in published papers (even if making code publicly available should be encouraged), the experimental sections often lack information useful to completely replicate the experiments and evaluate the methods. Table I summarizes the basic information about the experimental activities and whether they are provided by the authors in their papers. We have considered only the features that have an impact on the repeatability of the experimental activities and on the comparison between different mapping methods. In the following we discuss the data reported in the table.

Among the methods we considered, almost all experimental activities are conducted in indoor environments with "private" data, not freely available to other researchers. Only in few cases methods are evaluated using public data, available on web repositories such as [START_REF] Howard | The robotics data set repository (radish)[END_REF].

The use of publicly available data can help researchers not only to replicate an experiment but also to evaluate the performance of different methods when applied to the same data. In fact, the publications we examined do not perform any comparison, even qualitative, with other methods, with the exception of [START_REF] Gutmann | Amos: Comparison of scan matching approaches for self-localization in indoor environments[END_REF], [START_REF] Ip | Segment-based map building using enhanced adaptive fuzzy clustering algorithm for mobile robot applications[END_REF] that compare their methods with previous methods they have extended or improved.

Information about the size of the mapped environments is a qualitative gauge of a method, since the ability to cover large environments suggests that the method is reliable and does not diverge even after many iterations. The size is usually explicitly reported or it is inferable from the figures. Only in few cases this information is missing. The number of line segments that constitute the final map is another important information since it is related to the cost of storing the maps produced by the method. Only few authors provide this kind of information. The displacement between scans is also useful to compare the effectiveness of the methods. Some methods are (sometimes implicitly) based on the assumption of (very) small displacements between scans (e.g., [START_REF] Martignoni | Localizing while mapping: A segment approach[END_REF]), which implies a simpler approach to the scan matching problem but the use of a larger number of scans to cover the environment. Conversely, larger displacements between scans require fewer scans but present more challenging problems in matching the scans [START_REF] Amigoni | Building segment-based maps without pose information[END_REF]. Despite its importance, information about displacement is absent in many publications; in some cases it is inferable according to the length of the path covered by the robot and the number of scans [START_REF] Einsele | Real-time self-localization in unknown indoor environments using a panorama laser range finder[END_REF], [START_REF] Garulli | Simultaneous localization and map building using linear features[END_REF], [START_REF] Röfer | Building consistent laser scan map[END_REF], or from the robot positions marked in the maps [START_REF] Gutmann | Amos: Comparison of scan matching approaches for self-localization in indoor environments[END_REF], or from the information contained in the repository if public data are used [START_REF] Amigoni | Building segment-based maps without pose information[END_REF], [START_REF] Howard | Multirobot simultaneous localization and mapping using manifold representations[END_REF].

The obtained final maps are usually evaluated visually with respect to the real environment or to a ground-truth map, if available [START_REF] Cox | Blanche: An experiment in guidance and navigation of an autonomous robot vehicle[END_REF], [START_REF] Ip | Segment-based map building using enhanced adaptive fuzzy clustering algorithm for mobile robot applications[END_REF], [START_REF] Kim | Range and contour fused environment recognition for mobile robot[END_REF], [START_REF] Mazl | Building a 2D environment map from laser range-finder data[END_REF]. Only [START_REF] Ip | Segment-based map building using enhanced adaptive fuzzy clustering algorithm for mobile robot applications[END_REF] evaluates numerically the likeness of the obtained map by considering the geometrical distance between the map and a ground-truth map. Some authors [START_REF] Garulli | Simultaneous localization and map building using linear features[END_REF], [START_REF] Gonzáles-Baños | Navigation strategies for exploring indoor environments[END_REF], [START_REF] Gutmann | Amos: Comparison of scan matching approaches for self-localization in indoor environments[END_REF] use simulated data to better compare the results with respect to an artificial known environment that represents the ground-truth. In some cases the final map is evaluated by calculating the error on the robot pose; for example, by calculating the error between the starting and ending poses of the robot after following a closed path or by comparing the pose of the robot with some ground-truth landmarks [START_REF] Ip | Segment-based map building using enhanced adaptive fuzzy clustering algorithm for mobile robot applications[END_REF]. Note that multirobot mapping systems (e.g., [START_REF] Howard | Multirobot simultaneous localization and mapping using manifold representations[END_REF]) usually need to evaluate the errors on the mutual estimated pose of the robots. However, we do not take into consideration this issue in this paper.

The parameters and their values are vital information in order to replicate the experimental activity. Nevertheless only few authors report the values at least of the main parameters involved. Besides allowing repeatability, the values of the parameters and hints about how they influence the behavior and the performance of a method are helpful to evaluate the method itself. Only [START_REF] Amigoni | Building segment-based maps without pose information[END_REF], [START_REF] Ip | Segment-based map building using enhanced adaptive fuzzy clustering algorithm for mobile robot applications[END_REF] show the behavior of the system as the values of the parameters change. This enables a more thorough understanding of the proposed method, in addition to a better assessment of its reliability and robustness.

Closing a loop path is a challenging problem that a scan matching method must solve [START_REF] Thrun | Robotic mapping: A survey[END_REF]. Accordingly, only few papers do not show any experiment with a loop path.

The processing time for building the map is useful in order to evaluate whether the method can run in real-time or not. This information is not provided by all the authors; few authors report the average time taken between two integrations [START_REF] Amigoni | Building segment-based maps without pose information[END_REF], [START_REF] Einsele | Real-time self-localization in unknown indoor environments using a panorama laser range finder[END_REF], [START_REF] Gonzalez | Map building for a mobile robot equipped with a 2D laser rangefinder[END_REF], [START_REF] Zezhong | Scan matching based on CLS relationships[END_REF], [START_REF] Zhang | Line segment based map building and localization using 2D laser rangefinder[END_REF], whereas some other authors say only that their methods work in real time without any other information [START_REF] Howard | Multirobot simultaneous localization and mapping using manifold representations[END_REF], [START_REF] Ip | Segment-based map building using enhanced adaptive fuzzy clustering algorithm for mobile robot applications[END_REF], [START_REF] Mazl | Building a 2D environment map from laser range-finder data[END_REF], [START_REF] Röfer | Building consistent laser scan map[END_REF].

III. THE PROPOSED METHODOLOGY: DESCRIPTION

From the discussion of the previous Section, we can see clearly the heterogeneity and incompleteness of the exper-imental results provided to validate the proposed mapping systems. In order to contribute to overcome this situation, in this Section we propose a methodology for the experimental activities.

In order for experiments to be repeatable, the following issues have to be addressed when presenting experimental activities:

• The mapping system has to be applied to publicly available data.

• The values of the parameters should be indicated.

• Some experiments in which the mapping system does not perform well should be shown. This can help other researchers to better understand the system. In order to evaluate and to compare the different methods, the following issues have to be addressed when presenting experimental activities:

• All the data about the produced maps should be clearly indicated. These include: the dimensions of the mapped environment, the number of line segments, the time required to build the map, and the displacements between the scans.

• The behavior of the mapping system for different values of the parameters should be shown. This allows one to evaluate the robustness of the system, for example to find out the ranges of values of parameters for which the system works well.

• The map produced following a closed loop path in the environment should be shown, in order to evaluate the ability of the method not to "diverge".

• When a ground-truth map is available (this is not always the case), it should be used to assess the quality of the produced map, by evaluating its distance from the ground-truth map (e.g., according to the Hausdorff metric). With the availability of the above information, other researchers can evaluate the "soundness" of a proposed mapping method. For example, if the map of a 5 m-long empty corridor contains 1, 000 line segments, then it is evident that the mapping system that produced the map is not "sound". However, apart from some extreme cases, evaluating the "soundness" of a mapping method may not be easy: this issue deserves more investigation in future research.

With the availability of the above information, other researchers can compare different mapping methods. When data are publicly available, the mapping systems can be applied to the same data (possibly, after some pre-processing) and the resulting maps can be compared (this could be facilitated by making the resulting maps publicly available as well). The comparison can be done according to several dimensions, including: time, space (e.g., number of line segments), and quality (e.g., distance from a ground-truth map). These dimensions are contrasting with each other; for example, a method that produces very accurate maps could require more time and produce more line segments than a method that produces less accurate maps. A criterion to compare different mapping systems should be based on an application-dependent combination of the above dimensions. This issue deserves more investigation in future research.

Many of the above issues are quite obvious but, as we have seen, only rarely all of them are addressed in published papers. These issues are by no means exhaustive and others can be easily added (e.g., the average length of line segments). However, we think that they are a good starting point for the definition of a complete and rigorous experimental methodology for robotic mapping.

Although we have presented it for segment-based maps, the proposed methodology is more general and could be applied to mapping methods that work with different types of maps. For instance, it could be applied to grid-based mapping systems without any change, except for defining scans appropriately and for substituting the number of line segments with the number of cells. Other issues could also be added, like the dimensions of the cells. Extensions of the methodology should be deeper investigated, since it may not be easy to compare a segment-based mapping system with a grid-based one.

IV. THE PROPOSED METHODOLOGY: APPLICATION

In this Section, we present an example of how to apply the proposed methodology to a mapping system we developed. We first describe our mapping system and then we illustrate some experimental results, according to the proposed methodology.

A. Our Mapping System

The mapping system we use has been extensively described in [START_REF] Amigoni | Building segment-based maps without pose information[END_REF]. Here, we briefly summarize its main features (please refer to the original paper for the details).

The scans and the maps we consider are collections of line segments obtained from 2D laser range data. Our method is exclusively based on the geometrical information contained in the scans to be integrated. In particular, we consider angles between pairs of line segments in the scans as a sort of "geometrical landmarks" on which the matching process is based. We define a INTEGRATE function for integrating two scans into a map. Let's call S 1 and S 2 the two scans and S 1,2 the resulting map. In our method, two points are considered to coincide when they are closer than POINTDISTANCETOL-ERANCE and two angles are considered equal when their values differ of less than ANGLEDIFFERENCETOLERANCE. These are the two main parameters of the method. Function INTEGRATE operates in three major steps, as follows:

a) Determine the possible transformations: This step first finds the angles between the line segments in S 1 and between the line segments in S 2 and then finds the possible transformations (namely, the rotations and translations) that superimpose at least one angle α 2 of S 2 to an equal angle α 1 of S 1 . Recall that angles between pairs of line segments in a scan are the geometrical landmarks we adopt. b) Evaluate the transformations: Every transformation found in the previous step needs to be evaluated in order to identify the best one. To determine the goodness of a transformation t we transform S 2 on S 1 (in the reference frame of S 1 ) according to t (obtaining S t

2 ), then we calculate the approximate length of the line segments of S 1 that correspond to (namely, match with) line segments of S t 2 . The value of a transformation is the sum of the lengths of the corresponding line segments that the transformation produces.

c) Apply the best transformation and fuse the line segments: Once the best transformation t has been found (i.e., the one with the largest value), this step transforms the second scan S 2 in the reference frame of S 1 according to t obtaining S t 2 . The map that constitutes the output of INTEGRATE is then obtained by fusing the line segments of S 1 with the line segments of S t 2 . The main idea behind the fusion of line segments is that a set of matching line segments is substituted by a single polyline, in order to reduce the dimensions of the resulting map.

We have developed methods for integrating a sequence S 1 , S 2 , . . . , S n of n scans by repeatedly calling INTEGRATE. We used one of these methods, called pivot method, in the experimental activities presented in this paper (see [START_REF] Amigoni | Building segment-based maps without pose information[END_REF] for details).

B. Experimental Results

Following the proposed methodology, we applied our mapping system to the "stanford-gates1" data set publicly available in the Robotics Data Set Repository (Radish) [START_REF] Howard | The robotics data set repository (radish)[END_REF]. (Note that the experiments presented in this Section are complementary to those presented in [START_REF] Amigoni | Building segment-based maps without pose information[END_REF].) The data set is a 30-minute tour through the first floor of the Stanford's Gates Computer Science Building. The robot used to collect the data is a Pioneer 2DX with a forward-pointing SICK LMS 200. The laser was running at high speed (75 Hz scans) in the 10 mm, 1 • mode. The data set includes both laser data and odometry data. We considered only laser data (about 115, 000 laser scans). For each scan of the data set, we approximate the points acquired by the laser range scanner by line segments following the approach described in [START_REF] Gonzáles-Baños | Navigation strategies for exploring indoor environments[END_REF] (the distance between a perceived point and its approximating line segment is 25 mm at most). We call scans S x , where x is the time (in seconds) at which a scan has been acquired, according to the timestamps reported in the data set. The mapping system has been coded in C++ and run on a 2 GHz Pentiunm IV computer with Debian Linux 3.1.

POINTDISTANCETOLERANCE has been set to 10 mm and ANGLEDIFFERENCETOLERANCE to 0.2 rad. We selected these values by performing some integrations of scans randomly taken from the data set. Then, we tested the robustness of our method with respect to variations of these values. To this end, we considered three randomly selected pairs of scans of the stanford-gates1 data set and we applied our scan integration method varying the values of POINTDIS-TANCETOLERANCE and ANGLEDIFFERENCETOLERANCE. The method has been able to correctly integrate the pairs when the above parameters had values within 3 mm and 13 mm for POINTDISTANCETOLERANCE and 0.16 rad and 1.53 rad for ANGLEDIFFERENCETOLERANCE.

Sometimes, some values for the parameters do not allow to find the correct integration between two scans. Here we report some examples, in order to better clarify the role of the parameters. (In these examples, we do not apply the fusion step of INTEGRATE to the maps.) Fig. 1 shows the integration of scans S 1238 and S 1242 for different values of POINTDISTANCETOLERANCE. In this case, a smaller value of POINTDISTANCETOLERANCE allows our method to integrate the two scans with higher accuracy. In general, smaller values of POINTDISTANCETOLERANCE make the method stricter when evaluating if two points coincide, thus achieving an higher accuracy in integrating scans. Fig. 2 shows the integration of scans S 1536 and S 1540 for different values of ANGLEDIFFERENCETOLERANCE. In this case, a larger value of ANGLEDIFFERENCETOLERANCE allows our method to integrate the two scans with higher accuracy. In general, larger values of ANGLEDIFFERENCETOLERANCE make the method more tolerant when comparing two angles in order to find a possible transformation between two scans, thus increasing the number of transformations to evaluate. We applied our mapping method to build the maps of four large portions of the environment of the stanford-gates1 data set. Fig. 3 shows these maps, while Table II reports some information about them. Finally, we applied our mapping method to integrate these maps to obtain the global map of the environment. Results are shown in Fig. 4 and in the last row of Table II.

The results show that our mapping system is robust and produces maps composed of a limited number of line segments (according to the dimensions of the environment). Note that the majority of the line segments are of little importance (they are due to noise in data acquisition and in scan integration): filtering out those shorter than 300 mm

TABLE II INFORMATION ABOUT MAPS A, B, C, D, AND THE GLOBAL MAP

the number of line segments in the four maps drops to 150, 178, 51, and 273, respectively. The system can also perform in real-time, requiring (on average) less than 1 s to integrate a scan. The error in closing a 50 m loop (upper corridor of the map A in Fig. 3) is about 0.5 m, that is reasonable given the dimensions of the environment. Note that if, differently from our case, a mapping system involves casual choices, mean and variance of closing-loop error should be reported.

We expect that the information presented in this Section and the fact that we have used publicly available data (following the proposed methodology) will allow other researchers

V. CONCLUSIONS

In this paper, we have presented a preliminary version of a methodology for performing good experimental activities in the field of robot mapping systems. The proposed methodology has been applied to a mapping system to show how it can improve the repeatability and controllability of experimental results and the comparison of different mapping systems.

Future work will address a more complete definition of the proposed methodology. For example, some ways to evaluate the "soundness" and the robustness of a mapping method and to compare different mapping methods are needed. An interesting approach has been adopted in [START_REF] Minguez | Metric-based iterative closest point scan matching for sensor displacement estimation[END_REF], where the authors compare their method with other two methods based on the Iterative Closest Point framework. They match each scan against itself with an initial random displacement. In this way, the ground-truth is known and the accuracy of each method can be evaluated and compared. Although ICP methods work with point-based maps, this experimental approach can be easily redeployed for segment-based maps. Finally, the methodology should be "put at work" by applying it to different mapping systems in order to assess its potential.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Integration of scans with POINTDISTANCETOLERANCE = 9 mm (left) and with POINTDISTANCETOLERANCE = 2 mm (right)

Fig. 3 .

 3 Fig. 3. Maps A (top left), B (top right), C (bottom left), and D (bottom right)

Fig. 4 .

 4 Fig. 4. Global map

TABLE I FEATURES

 I REPORTED IN EXPERIMENTAL ACTIVITIES FOR SOME SCAN MATCHING METHODS

VI. ACKNOWLEDGMENTS

The authors would like to thank Igor Giussani for producing some of the experimental results and Brian Gerkey for providing the stanford-gates1 data set.