Michel Kieffer
email: kieffer@lss.supelec.fr

Distributed bounded-error state estimation

This paper presents a distributed bounded-error state estimation algorithm suited, e.g., to measurement processing by a network of sensors. Contrary to centralized estimation, where all data are collected to a central processing unit, here, each data is processed locally by the sensor, the results are broadcasted to the network and taken into account by the other sensors. A first analysis of the conditions under which distributed and centralized estimators provide the same results is presented. An application to the tracking of a moving source using a network of sensors measuring the strength of the signal emitted by the source is considered * .

INTRODUCTION

State estimation usually involves a state equation describing the evolution of the considered system and measurements provided by sensors [START_REF] Gelb | Applied Optimal Estimation[END_REF]. In a classical centralized estimation scheme, all measurements are gathered at a central processing unit (CPU), responsible for all the signal processing. Alternatively, computations may be performed in a decentralized, or distributed, manner. In such approach, each sensor processes its own measurements together with the data received from other sensors [START_REF] Speyer | Computation and transmission requirements for a decentralized linear-quadratic-gaussian control problem[END_REF], [START_REF] Durrant-White | Toward a fully decentralized architecture for multi-sensor data fusion[END_REF], [START_REF] Rao | A fully decentralized multi-sensor system for tracking and surveillance[END_REF]. The main advantage of the latter approach is its robustness to failures of sensors or of the CPU: in the distributed context, even if one or more sensors fail, other sensors continue to provide state estimates. Distributed estimation techniques are currently gaining a growing interest, thanks to the availability of nodes in wireless sensor networks with increased battery and processing power [START_REF] Kay | The design space of wireless sensor networks[END_REF], [START_REF] Haenselmann | GFDL Wireless Sensor Network textbook[END_REF]. Distributed state estimation is more complex than centralized estimation: the limited computational complexity and transmission capabilities have to be taken into account. This raises several problem, among which one may cite the choice of the type of data to transmit between sensors, of the transmission protocol, of the way these data have to be processed by the sensors, etc.

For what concerns the type of data, measurements or estimates may be transmitted. When the number of node in the network is larger than the dimension of the state vector, transmitting measurements to all sensors may result in a high communication load, except if data are very coarsely quantized, as done in [START_REF] Ribeiro | SOI-KF: Distributed Kalman filtering with low-cost communications using the sign of innovations[END_REF]. Thus, estimates are usually transmitted from one sensor to the others, as in [START_REF] Rao | A fully decentralized multi-sensor system for tracking and surveillance[END_REF], [START_REF] López-Orozco | An asynchronous, robust, and distributed multisensor fusion system for mobile robots[END_REF], Olfati-Saber [2005], [START_REF] Alriksson | Experimental evaluation of a distributed Kalman filter algorithm[END_REF], where decentralized version of the Kalman filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] were proposed. As for centralized Kalman filtering, state perturbations and measurement noise are assumed to be described by Gaussian probability density functions. Moreover, the dy-namical part of the state equation is usually assumed to be linear or is linearized. As a consequence, the best than can be achieved is a an estimate and its associated confidence region: no guarantee may be provided as regards quality of the obtained results.

In the approach proposed in this paper, state perturbation and measurement noise will be assumed to be bounded with known bounds. The aim of distributed state estimation in such bounded-error context is to evaluate, at any node of the network, a set containing all state values which are consistent with the information (model, noise bounds, measurements, data of neighboring sensors) available at that node. For that purpose, set estimates will be exchanged between nodes of the network. The resulting set estimator is guaranteed, since each set evaluated by the nodes contains the actual state value, as in centralized bounded-error state estimation [START_REF] Kieffer | Guaranteed recursive nonlinear state estimation using interval analysis[END_REF][START_REF] Kieffer | Guaranteed recursive nonlinear state bounding using interval analysis[END_REF], provided that the hypotheses on the model and noise bounds are satisfied.

The paper is organized as follows. Section 2 presents the model, before recalling the principle of idealized centralized bounded-error state estimation in Section 3. Its distributed extension is provided in Section 4. The conditions under which the distributed and the centralized algorithms perform similarly is provided. Using interval analysis [START_REF] Moore | Interval Analysis[END_REF], [START_REF] Jaulin | Applied Interval Analysis[END_REF], it will be possible to build an implementable counterpart of the idealized distributed state estimation algorithm. Finally, Section 5 presents the application of the preceding algorithm to the tracking of a source emitting waves which amplitude is measured by the nodes of a wireless sensor network.

BOUNDED-ERROR DISCRETE-TIME STATE ESTIMATION

Consider a system described by a model consisting of a discrete-time state equation

x k = f k (x k-1 , w k , u k ) , (1
) where x k is the state vector of the model at time instant k (the sampling period is T ). The state perturbation vector w k accounts for unmodeled parts of the system and is assumed to remain in some known box [w]. The input vector u k is also assumed known. At k = 0, x 0 is only assumed to belong to some (possibly large) known set X 0 .

Assume that at time k, each node = 1 . . . L of a sensor network has access to a noisy measurement vector y k . The measurement process is described by the observation equations

y k = g k x k , v k , = 1 . . . L
(2) where v k is the measurement noise, assumed bounded in some known box [v]. Usual observation equations are

g k x k , v k = h k (x k ) + v k (3) or g k x k , v k = h k (x k ) • v k , (4
) depending on whether the measurement noise is additive or multiplicative.

Centralized and distributed state estimation aims at evaluating at each time instant k the set of all state values that are consistent with the models (1) and ( 2), the measurements, and the bounds on the state perturbation and measurement noise. The main difference between centralized and distributed estimation is that in the former, all measurements are available at a given point, whereas in the latter, estimation is performed at each node, and data has to be exchanged between sensors.

CENTRALIZED DISCRETE-TIME STATE ESTIMATION IN A BOUNDED-ERROR CONTEXT

This section briefly summarizes centralized state estimation, since it constitutes the reference which distributed algorithms should reach.

When all measurements at time k are available at a CPU, one gets from (1) and (2)

x k = f k (x k-1 , w k , u k ) , y k = g k (x k , v k ) , (5)
with

y T k = ((y 1 k ) T . . . (y L k ) T ) and v T k = ((v 1 k ) T . . . (v L k ) T
). Summarizing the information available at time k, one obtains

I k = X 0 , {[w j ]} k j=1 , {[v j ]} k j=1 , {[y j ]} k j=1 . (6)
Performing a centralized bounded-error state estimation at time k aims at characterizing the set X k|k of all values of x k that are consistent with (1), (2), and I k . One may propose an idealized algorithm [START_REF] Kieffer | Guaranteed recursive nonlinear state bounding using interval analysis[END_REF], alternating, as the Kalman filter a prediction step involving (1)

X k|k-1 = f k (x, w, u k ) | x ∈ X k-1|k-1 , w ∈ [w]
(7) and a correction step accounting for the new measurement using (2)

X k|k = x ∈ X k|k-1 | y k = g k (x, v) , v ∈ [v] L . (8)
The two steps of the idealized algorithm are depicted in Figure 1.

This idealized algorithm requires the evaluation of the direct image of a set by a function in the prediction step ( 7) and the evaluation of the inverse image in the correction step (8). To get an implementable counterpart of this algorithm, first one has to wrap the sets X k|k-1 and

x 1 x 2 X k k -1 -1 | x k-1 y 1 y 2 X k k | X k k | -1 y k , [v] g k f k x k
Fig. 1. Idealized recursive bounded-error state estimator X k|k , which can take any shape, by sets X k|k-1 and X k|k allowing a parametric description. Candidate wrappers are for example ellipsoids [START_REF] Schweppe | Uncertain Dynamic Systems[END_REF], zonotopes [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF], interval vectors or boxes [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF] or union of boxes [START_REF] Kieffer | Guaranteed recursive nonlinear state bounding using interval analysis[END_REF][START_REF] Kieffer | Guaranteed set computation with subpavings[END_REF]. Due to the wrapping effect, some x ∈ X k|k may not correspond to values of the state vector consistent with (6). Depending on the wrapper, prediction and correction steps may require further hypotheses on the state equation. When ellipsoids and zonotopes are considered, (1) and ( 2) have to be linear [START_REF] Schweppe | Uncertain Dynamic Systems[END_REF], [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF]. This condition is not required when considering boxes or unions of boxes, thanks to interval analysis, which provides tools to evaluate outer approximations of ( 7) and ( 8), see [START_REF] Kieffer | Guaranteed recursive nonlinear state bounding using interval analysis[END_REF], [START_REF] Jaulin | Applied Interval Analysis[END_REF].

DISTRIBUTED STATE ESTIMATION

When considering distributed state estimation, all measurements are not directly available at each node of the network. The set estimate X k|k provided by the -th node is thus such that X k|k ⊂ X k|k . Ideally, the data exchange between sensors should allow to get

X k|k = X k|k , = 1 . . . L. (9)
To establish conditions under which ( 9) is satisfied, some notions of graph theory have to be recalled. For more details, the reader is referred to [START_REF] Harary | Graph Theory[END_REF], [START_REF] Bollobás | Modern Graph Theory[END_REF].

Description of the network

The network of L sensors is represented by a graph G = (V, E). V is the set of L vertices of the graph, each vertex representing a node of the network and E is the set of edges of the graph. An edge {k, } ∈ E connecting two vertices k ∈ V and ∈ V indicates that the two corresponding nodes are able to directly exchange information; the graph is thus undirected. In what follows, it is assumed that G is entirely connected, i.e., that there is always a path from any vertex to any other vertex in G and that each vertex is connected to itself.

The distance between two vertices in G is the number of edges in a shortest path connecting them. Consider a vertex ∈ V, then C ({ }) = {k ∈ V | (k, ) ∈ E} (10) denotes the set of all vertices that are directly connected to , i.e., that are at a distance not larger than one of . More generally, for any W ⊂ V, C (W) ⊂ V is the set of all vertices which are at a distance not larger than one from a given vertex of W. The set

C (C ({ })) = C 2 ({ }) (11
) contains thus all vertices that are at a distance not larger than two of . More generally, C n ({ }) contains all vertices that are at a distance not larger than n of . The eccentricity ε of a vertex ∈ V is the largest distance between and any other vertex in G. Finally, the diameter d of G is the maximum eccentricity of any vertex in G.

Hypotheses and idealized algorithm.

The following measurement processing and communication protocol will be considered.

At time k, each sensor processes its own measurement y k , to get X ,0 k|k . Between time k and k + 1, a first round trip is considered (r = 1) during which each sensor broadcasts its own estimate X ,r-1 k|k to all the sensors of the network (only those which are directly connected to receive the information). Each sensor may communicate in turn using simple time division multiple access, or simultaneously using, e.g., code division multiple access [START_REF] Sklar | Digital Communications: Fundamentals and Applications[END_REF]. Then each sensor receives and processes X s,r-1 k|k , s ∈ C ({ }) to get X s,r k|k . Depending on the sampling time T , more round trips (r > 1) may be considered. Just before time k + 1, each sensor builds a final estimate X k|k . This way of processing and transmitting information leads to the following idealized distributed algorithm.

For each sensor = 1 . . . L,

(1) At time k

X k|k-1 = f k (x, w, u k ) | x ∈ X k-1|k-1 , w ∈ [w] . (12)
X ,0 k|k = x ∈ X k|k-1 | y k = g k (x, v) , v ∈ [v] . ( 13
) (2) Between k and k + 1, for r = 1 to R max (number of round trips) X ,r k|k = s∈C({ }) X s,,r-1 k|k (14)
(3) Just before k + 1

X k|k = X ,Rmax k|k . ( 15
)
As for the centralized algorithm, the idealized distributed algorithm requires the evaluation of the direct and inverse images of a set by a function. Proposition 1 gives some conditions under which the distributed approach gives results similar to the centralized one.

Proposition 1. Consider a WSN of L nodes represented by an entirely connected graph G = (V, E) of diameter d.

Assume that at time

k -1, X k-1|k-1 = X k-1|k-1 for all ∈ V. If the number of roundtrips R max satisfies R max d, then one has at time k X k|k = X k|k (16) for all ∈ V. ♦ Proof. Consider a vertex ∈ V. Since X k-1|k-1 = X k-1|k-1
, after the prediction step ( 12), X k|k-1 = X k|k-1 , where X k|k-1 is provided by (7). The first correction step done at involves only the measurement vector y k to get X ,0 k|k . After the first roundtrip, the estimate at becomes

X ,1 k|k = s∈C({ }) X s,0 k|k = s∈C({ }) x ∈ X s k|k-1 | y s k = g s k (x, v) , v ∈ [v] = s∈C({ }) x ∈ X k|k-1 | y s k = g s k (x, v) , v ∈ [v] = x ∈ X k|k-1 | y C({ }) k = g C({ }) k (x, v) , v ∈ [v] ,
where y

C({ }) k and g C({ }) k
(x, v) are the vector and function consisting of the concatenation of all y s k and g s k (x, v), with s ∈ C ({ }).

After a second roundtrip, the estimate at becomes

X ,2 k|k = s∈C({ }) X s,1 k|k = s∈C({ }) x ∈ X k|k-1 | y C({s}) k = g C({s}) k (x, v) , v ∈ [v] = x ∈ X k|k-1 | y C(C({ })) k = g C(C({ })) k (x, v) , v ∈ [v] = x ∈ X k|k-1 | y C 2 ({ }) k = g C 2 ({ }) k (x, v) , v ∈ [v] .
Similarly, after R max roundtrips, one gets at

X ,Rmax k|k = x ∈ X k|k-1 | y C Rmax ({ }) k = g C Rmax ({ }) k (x, v) , v ∈ [v]
It is now enough to show that C Rmax ({ }) = V in order to prove that X ,Rmax k|k = X k|k . First, one has C Rmax ({ }) ⊂ V. Assume now that there exists some k ∈ V such that k / ∈ C Rmax ({ }). This means that k lies at a distance strictly larger than R max from . Since the diameter d of G is lower than R max the distance between two vertices is necessarily lower than d, which contradicts the initial assumption. Thus any k ∈

V satisfies k ∈ C Rmax ({ }) and C Rmax ({ }) = V.
The result of Proposition 1 mainly shows that when there are enough information exchanges between sensors, the distributed estimate converges at any sensor to the centralized estimate. What is more interesting is that the number of roundtrips needed for convergence depends only on the diameter of the graph associated with the WSN.

When R max < d, the situation is much more complex, since not all sensors will have access to all measurements (or to their contribution to the estimation of x k ). For the first roundtrips at time k + 1, sensor will have to broadcast information about X ,r k+1|k+1 , but also about X ,r k|k , as long as X ,r k|k has not converged to X k|k . Again, the diameter of the graph plays a crucial role. Further analysis is still required.

Practical algorithm

The implementation of the distributed algorithm requires the same wrappers as for the centralized algorithm. The complexity of the wrapper impacts the amount of information that has to be transmitted to the other nodes.

In this paper, a basic version of the algorithm is considered, where sets are wrapped by boxes. This simple description allows to perform the prediction step using basic tools from interval analysis, such as interval evaluations for the prediction step. The correction step is implemented via interval constraint propagation [START_REF] Jaulin | Applied Interval Analysis[END_REF]. The advantage of dealing with boxes is that interval analysis may readily be implemented on chips with reduced computational capabilities, see [START_REF] Piskorski | Efficient 16-bit floating-point interval processor for embedded systems and applications[END_REF]. A more sophisticated version could involve description of sets using subpavings, a prediction step implemented using ImageSp [START_REF] Kieffer | Guaranteed recursive nonlinear state bounding using interval analysis[END_REF] and Sivia Jaulin and Walter [1993] combined with interval constraint propagation for the correction step.

APPLICATIONS

For the application part, a problem of source tracking from readings of signal strength (RSS) data. Centralized approaches have been proposed to solve this problem for acoustic sources [START_REF] Sheng | Maximum likelihood multiplesource localization using acoustic energy measurements with wireless sensor networks[END_REF] and for sources emitting electromagnetic waves, see, e.g., [START_REF] Kontkanen | Probabilistic methods for location estimation in wireless networks[END_REF], [START_REF] Gustafsson | Mobile positioning using wireless networks[END_REF], [START_REF] Gezici | Localization via ultra-wideband radios[END_REF]. In the first case, some knowledge of the decay rate of the RSS (path loss exponent) is needed for efficient nonlinear least squares estimation. In the second case, an off-line training phase is required to allow maximum a posteriori localization. In both cases, a good initial guess of the location of the source facilitates convergence to the global minimum of the cost function. Distributed approaches have also been employed, e.g., in [START_REF] Rabbat | Decentralized source localization and tracking[END_REF], where a distributed version of a nonlinear least squares solver has been presented. When badly initialized, it suffers from the same convergence problems as the centralized approach, as illustrated in [START_REF] Hero | Sensor network source localization via projection onto convex sets (POCS)[END_REF], which advocates projection on convex sets. Nevertheless, the latter requires an accurate knowledge of the source signal strength and of the path loss exponent.

Models

The sensors are assumed to be static, and their location is known and denoted by r ∈ R 2 , = 1 . . . L. The state of the source at time instant k is

x k = (θ 1 k , θ 2 k , ϕ 1 k , ϕ 2 k ) T ∈ R 4 , where θ k = θ 1 k , θ 2 k
T represents the location of the source and

ϕ k = ϕ 1 k , ϕ 2 k T its speed.
The evolution of the source is assumed to be described by

x k =    1 0 T 0 0 1 0 T 0 0 1 0 0 0 0 1    x k-1 + T    0 0 w 1 w 2    . (17)
Since the inputs are unknown, they are considered as bounded state perturbations, thus, w 1 ∈ [w] and w 2 ∈ [w].

The mean power P dB d k (in dBm) received by the -th sensor is described by

P dB d k = P 0 -10n p log d k d 0 , (18)
where P 0 is the power received at a reference distance d 0 , n p is the path-loss exponent, and d k = rθ k , see [START_REF] Okumura | Field strength ans its variability in VHF and UHF landmobile radio service[END_REF]. The received power is assumed to lie within some bounds

P dB (d) ∈ P 0 -10n p log d d 0 -e, P 0 -10n p log d d 0 + e ,
where e is assumed known.

Assuming that P 0 , d 0 , and n p are known, the RSS by sensor = 1 . . . L may be rewritten as

y k (θ k , v ) = h (θ k ) v , ( 19
) with h (θ k ) = A |r -θ k | np , A = 10 P0/10 d np 0 , (20)
and v ∈ [v] = 10 -e/10 , 10 e/10 . The noise, additive in the log domain is multiplicative in the normal domain.

Prediction step: direct image evaluation

Assume that the estimate at time k -1 for node is x k-1|k-1 . The estimate x k|k-1 may be obtained from (17) as follows

x k|k-1 =    1 0 T 0 0 1 0 T 0 0 1 0 0 0 0 1    x k-1|k-1 + T    0 0 [w] [w]    . (21)

Correction step: interval constraint propagation

At sensor , y ∈ y is measured. Some boxes x k is assumed to be available from the prediction step, or as results transmitted by the other sensors to sensor . The parameter vector has to satisfy the constraint provided by the RSS model

y - A r -θ k np = 0. ( 22
)
Using interval constraint propagation, it is possible to reduce the domain for the variables using ( 22). The contracted domains may be written as

                   [y ] = y ∩ A r -[θ k ] np , [θ k,1 ] = [θ k,1 ] ∩ r 1 ± (A/ [y ] ) 2/np -(r 2 -[θ k,2 ]) 2 , [θ k,2 ] = [θ k,2 ] ∩ r 2 ± (A/ [y ] ) 2/np -(r 1 -[θ k,1 ]) 2 . (23)
       [ϕ k,1 ] = [ϕ k,1 ] ∩ ( [θ k,1 ] -[θ k-1,1 ] T + T [w]), [ϕ k,2 ] = [ϕ k,2 ] ∩ ( [θ k,2 ] -[θ k-1,2 ] T + T [w]). (24)
Each sensor will perform this constraint propagation before transmitting its updated estimate to its neighbours.

RESULTS

A field of 50 m×50 m is considered, with its origin at the center. A WSN of L = 25 sensors with communication range of 15 m is spread over this field. The source is placed at θ * = (5 m, 5 m), with characteristics P 0 = 20 dBm, d 0 = 1 m. The measurement noise is such that e = 4 dBm. The path-loss exponent is n p = 2, assumed constant over the field. The sampling time is T = 0.5 s and [w] = [-0.5, 0.5] 2 m•s -2 . Figure 2 illustrates the connectivity of the considered regular WSN and a typical trajectory followed by the source. Each sensor is connected in average with 6.9 sensors of the network.

The simplest algorithm implementation presented in Section 4.3 has been considered: sets are represented by boxes, simple image evaluations using inclusion functions are performed and correction is done by interval constraint propagation. This limits the amount of information to be exchanged between sensors and the computational effort. The localization performance using this algorithm is depicted in Figure 3 for 100 realizations of the source trajectory. The average width of the solution box (top of Figure 3) provided at each time instant decreases very quickly before reaching a floor slightly higher than the minimum width and increases again after about 18 s. At the beginning, the source is close to the middle of the field and many sensors participate to the localization. When the sensor moves near the limits of the field, the number of involved sensors decreases and as a result the localization accuracy worsens. This effect is even more important when the source moves outside the field. A similar behavior is seen for the average norm of the localization error taking the center of the solution boxes at each time instant as estimate. The convergence is quite fast. The number of round trips impacts significantly the quality of the obtained results algorithm. In this paper, we have considered distributed boundederror state estimation applied to the problem of source tracking with a network of wireless sensors. Estimation is performed in a distributed context, i.e., each sensor has only a limited amount of measurements available. A guaranteed set estimator is put at work.

There is still large space for improvements in the considered problem. First, convergence properties have to be more carefully studied. In particular, more general conditions under which the distributed solution coincides with the centralized one have to be determined, in particular when wrappers are used instead of general sets. This type of problem is partly addressed in [START_REF] Yokoo | Distributed Constraint Satisfaction: Foundations of Cooperation in Multi-Agent Systems[END_REF], [START_REF] Bejar | Sensor networks and distributed CSP: Communication, computation and complexity[END_REF]. Robustness to outliers has also to be considered. Extensions of [START_REF] Jaulin | Guaranteed robust nonlinear parameter bounding[END_REF] have to be proposed in a distributed context.

Fig. 2 .

 2 Fig. 2. Trajectory of the source (o); each sensor is represented by a cross (x), distances are in meters Once the estimate for θ k has been improved, it may be useful to use (21) to improve the estimation of the speed of the source as follows    

Fig. 3 .

 3 Fig. 3. Width of the box [θ 1,k ] × [θ 2,k ], and norm of the localization error when the estimate is taken as the center of the solution box, when one roundtrip (top) and three roundtrips (bottom) are allowed; results are averaged over 100 random paths followed by the source

* This work has been partly supported by NEWCOM++


Introduction

The gaze guiding framework described in this technical report was presented at the 15th European Conference on Eye Movements, 24th August 2009, Southampton, UK, as an oral presentation with the title: "Gaze guidance through peripheral stimuli" [START_REF] Lint | Gaze guidance through peripheral stimuli[END_REF].

Gaze guiding can be used in a wide range of applications. Fields of use range from: education, transfer of knowledge from expert to novice, better understanding of art, advertising or numerous applications for driving a motor vehicle, just to name a few. The before mentioned applications can be implemented in the fields of augmented or virtual reality, as well as for attentive displays. We present here a gaze guiding application which can be used in the field of gaze-contingent displays, which are a subcategory of attentive displays. For a general overview paper on attentive displays the reader should refer to [START_REF] Baudisch | Focusing on the essential: considering attention in display design[END_REF].

A general problem of gaze-contingent displays is their limited response latency. We propose a gaze guiding method using near-peripheral stimuli which can decrease this latency by correctly anticipating the users gaze patterns. We also propose a novel and simple way to compute the saliency of the stimuli used in our framework.

Gaze contingent displays can be used in conjunction with the visualization of large data sets on large display devices in cases where interactivity is important. Even with today's powerful graphics hardware, very large data sets viewed at high resolution can cause problems at full screen display which makes a gaze contingent solution feasible in the advent of today's high quality and precise, non-intrusive eye trackers.

We will give at this point short descriptions of concepts from vision science which are important in the design and understanding of our proposed framework. Bottom-up processing, or data-driven processing refers to the type of information processing in the human visual system which can be described as a flowchart where the retinal image stays at the bottom and its subsequent, more processed interpretations are at higher levels. Conversely, top-down processing or hypothesis driven processing refers to processes that work in the other direction, with higher level representations as input and producing lower and lower level representations as output. It is an accepted fact in the pshycophisics literature that most processes in the early stages of visual processing are bottom-up, while there are many top-down processes present at a later level of visual processing -mainly where people have prior existing models of the viewed scene [3, pp. 84-85].

Another interesting subdivision of the human visual system is the classification in the where and what system. The where system consists of larger cell subdivision situated mostly in the periphery on the retina and is responsible for motion perception, depth perception, spatial organization and figure/ground segregation. The where system is color-blind, fast, has low acuity and high contrast sensitivity. The what system has a small-cell subdivision and is mostly placed in the foveal region of the retina. It is responsible for object and face recognition and color perception. The what system is color selective, slow, has high acuity and low contrast sensitivity. For more details on the where and what system, please refer to [START_REF] Livingstone | Vision and Art: The Biology of Seeing[END_REF]Chapter 4].

The main goal of our framework is to guide gaze with stimuli which are not consciously perceived by the subject. We achieve this goal by showing high Figure 1: The first 30 images in our collection resolution patches on a low resolution image viewed by the subject, placed in the peripheral region of the visual field, i.e. we want to stimulate the where system. The user views a set of images and its task is to deblurr all images, successively. If the users gaze is attracted by the current stimulus, an area associated with the stimulus is deblurred and the next stimulus is turned on.

We propose two different strategies for guiding: with (see Figure 5) and without anticipation (see Figure 6). The strategy with anticipation is based on the idea that most of the subjects saccades are preprogrammed and that we can already deblurr the area associated with the stimulus if the direction of the saccade points towards it. In the case of the strategy without anticipation, we only deblurr the area corresponding to the current stimulus when an eye tracker sample is already close to the stimulus.

The position of most our stimuli is given by the position of fixations from a previous pilot study consisting of a free exploration of selected 60 images. A thumbnail overview of our image collection is displayed in Figures 1 and2. An example of a recorded gaze pattern for one image is given in Figure 3. The white circles indicate fixations, with the circle radius proportional to fixation length. The goal of our experiment is to deblur the images in their totality, thus stimuli have to be present even on the periphery of the images. As it can be seen in Figure 3, gaze is mostly attracted to the central parts of the image and this is why we place additional stimuli in the peripheral region, following a spiral distribution.

As already mentioned earlier, each stimulus has its associated area surrounding it which gets deblurred, once the stimulus attracts the subject's gaze. For simplicity, we divide our input images in a grid of cells, and each cell has a stimulus associated with it, as shown in Figure 4. Each cell is labeled with the same color, as well as number of its corresponding stimulus.

INRIA

C = I min -I max I min + I max (1)
It is mainly used in experiments measuring contrast sensitivity. A definition better adapted for contrast in complex images was proposed by Peli [START_REF] Peli | Contrast in complex images[END_REF]. An improved version of the contrast in complex images was proposed by Mantiuk et al. [START_REF] Mantiuk | A perceptual framework for contrast processing of high dynamic range images[END_REF]. We use in our stimulus saliency computations a simplified, yet concise version of contrast measurement by Rizzi et al. [START_REF] Rizzi | A proposal for contrast measure in digital images[END_REF].

For the reader interested in the human visual system in general, good introductory books on the general topic of human visual perception are [START_REF] Livingstone | Vision and Art: The Biology of Seeing[END_REF][START_REF] Marr | [END_REF]. For a neuroscience view on the human visual system, the reader should refer to [START_REF] Hubel | Eye, brain, and vision[END_REF]. Another solid source on visual perception is [START_REF] Palmer | Vision science: photons to phenomenology[END_REF].

Related Work

Dorr et al. [START_REF] Dorr | Dynamic Perception Workshop of the GI Section" Computer Vision[END_REF] examine a system to guide eye movements using two types of stimuli in the periphery. The subjects have to view short videos of natural scenes; the success and "non-detection" of stimuli for guidance is evaluated.

Parkhurst et al. [START_REF] Parkhurst | Modeling the role of salience in the allocation of overt visual attention[END_REF] use a bottom-up model of visual processing based on the two stage processing framework proposed in the literature: early vision and center surround. The authors goal is to provide a framework which can predict fixation locations using a saliency map. They use a two stage multi-resolution approach two compute saliency maps and use them afterwards to evaluate the correlation between areas of high saliency and fixation location while the subjects viewed natural and artificial images. The authors conclude that bottom-up In [START_REF] Parkhurst | Scene content selected by active vision[END_REF] the authors investigate the possible correlation between the local contrast and fixation location for several types of natural images while free viewed by subjects. Analysis is made for image patches of different sizes, but the biggest difference between the contrast of fixated patches compared to random patches is attained at a patch size of 1 • visual angle.

Contrast sensitivity in peripheral vision is scrutinized in [START_REF] Thibos | Characterization of spatial aliasing and contrast sensitivity in peripheral vision[END_REF]. The authors concentrate on two tasks: resolution and detection of stimuli. Several experiments are described and as an important general conclusion the authors state that the detection task is much less affected by the mechanisms that limit resolution in the peripheral vision.

Itti et al. [START_REF] Itti | A Model of Saliency-Based Visual Attention for Rapid Scene Analysis[END_REF] present a saliency-based model of visual attention. Their model includes center-surround operations, a multi resolution Gaussian pyramid as well as a model for the orientation of features in the image. The presented model reproduces well human performance for pop-out tasks as well as detects most regions of high saliency in the input images.

In the fields of computer graphics and image processing, an algorithm for the assessment of image fidelity was proposed by Daly [START_REF] Daly | Visible differences predictor: an algorithm for the assessment of image fidelity[END_REF], called the Visible Differences Predictor. The algorithm tries to asses at which degree any differences physically present in images due to image processing or displaying are also visible differences. It is a complex framework which can be used in many fields for

INRIA

Experimental Setup

We will describe here into more detail our experimental setup. The subjects observe the sequence of images which have only a single deblurred area in the center at the beginning. If the given conditions are satisfied, and the current stimulus is visible and also attracts the users gaze the corresponding area is deblurred and merged with the already deblurred regions. The goal for the user is to deblur the hole image; we have in mean around 10 stimuli per image (depending on the complexity of each given image). After each image is deblurred, the subject is presented a questionnaire regarding several aspects of the process (if the process was annoying, or if the user likes the image visually).

The current stimulus is visible only if the visual angle between its position on the screen and the position of the users current fixation is superior to 8 •which puts the stimulus already in the foveal or para-foveal vision and inferior to 14 • which is the limit of the near-peripheral vision. Eye movements are detected through an ASL-H6 eye-tracker working at 120 Hz. The blinking frequency of our stimuli as a function of distance from the current fixation to the stimulus was adapted from [START_REF] Dorr | Dynamic Perception Workshop of the GI Section" Computer Vision[END_REF] and [START_REF] Dorr | Visibility of temporal blur on a gaze-contingent display[END_REF].

To be able to compare our two proposed strategies we use in random way one of them for each image presented to the subject. This way we can avoid the subjects learning of a particular strategy and have more significant results.

We asses the effectiveness of our proposed guiding method by measuring as principal indicators: the time in which gaze is attracted to every stimulus, and the distance from the last fixation to the current stimulus.

All images where low pass filtered using a simple Gaussian Blurring kernel of 20 pixels in GIMP. The blurring was selected by trial and error in such a way as to keep in most of the cases the contrast between the stimuli and its surroundings as low as to be just detected, but not resolved, see [START_REF] Thibos | Characterization of spatial aliasing and contrast sensitivity in peripheral vision[END_REF].

Guidance

We propose two different guiding strategies. For the strategy with anticipation, a gaze sample has to fall further than half the distance between the previous fixation and the stimulus, as illustrated in Figure 6. By using this strategy we try to anticipate the subjects next fixation and additionally avoid the subjects noticing the current stimulus. In the case of the strategy without anticipation, Figure 5, a gaze sample has two fall into the circle of 2 • visual angle radius surrounding the current stimulus. Because one of the goals of our framework is to keep stimuli as unperceived by subjects as possible in this case the stimulus will be invisible because it falls in the subjects foveal vision. Our system deblurs the area corresponding to the stimulus even if it is invisible, in the case the sample meets the strategy's constraints. We rely on the fact that saccades are ballistic movements and once begun, their trajectory cannot be altered [3, pp. 523]. If the peripheral stimulus was perceived by the subject, there is a good chance that the next fixation will land in the circle of 2 • visual angle surrounding the current fixation.

Stimulus Saliency

The stimuli used in our framework are circular patches with 1 • visual angle radius. They are high resolution patches overlayed on the low resolution image, displayed at salient positions for each image, as described in Section 1.

We propose a novel measure for stimulus saliency based on two metrics: contrast ratio and sum of squared differences between the corresponding high resolution and low resolution patches for each stimulus position. This saliency model is much more simplified than previously proposed models such as the one by Parkhurst et al. [START_REF] Parkhurst | Modeling the role of salience in the allocation of overt visual attention[END_REF] or Itti et al. [START_REF] Itti | A Model of Saliency-Based Visual Attention for Rapid Scene Analysis[END_REF], but has the advantage that it can be computed fast and it serves as an additional tool for the fast evaluation of the results of our effectuated experiments.

To measure the contrast ratio between 2 corresponding low and high resolution patches around all stimuli in the 60 images we use an GNU Octave [START_REF] Eaton | GNU Octave[END_REF] implementation of the contrast measurement method proposed by Rizzi et al. [START_REF] Rizzi | A proposal for contrast measure in digital images[END_REF]. This method computes a local and multi-scale contrast for a given natural image by constructing a pyramid of sub-sampled images and computing the contrast at every level. The pyramid is built by dropping every second row and column from the input image at each level, for N = [log 2 min(m, n)] levels, where I(m, n) is the input image. After the input image I is de-constructed, the contrast for every pixel (x, y) is computed as the average difference between its value and that of the eight nearest neighbors. All these values are afterwards averaged to obtain an average contrast for every given resolution level. A global contrast can be computed by finally averaging the computed contrasts for every given level. The equation for computing the global contrast for the input image I following all steps described earlier is as follows:

C(I) = ∀ level       ∀ pixel 8-neigh |P i -P j | 8 number of pixels number of levels       (2)
Using the contrast measurement method described above, we systematically compute the global contrast for corresponding patches from the high resolution and blurred images from the regions surrounding all stimuli. Having the global contrast for every high resolution as well as blurred patch, we then compute the contrast ratio between each corresponding patch. For each image in our collection we plot the contrast ratios, to be able to asses stimulus saliency based on our contrast metric, as presented in Figure 7. The corresponding image for the computed saliencies is shown in Figure 8. For most cases the proposed contrast ratio measure can be used as a good indicator of stimulus saliency, however there are exceptions where an additional measure of stimulus strength has to be computed in order to be sure about the saliency of the stimulus.

We propose to use the sum of squared differences as a measure for additional investigation of stimulus saliency. For every image patch I around the predefined stimuli for each image, we compute the sum of squared differences (SSD) as:

ssd(I) = ∀ pixel (I blurred (m, n) -I high resolution (m, n)) 2 (3)
where I blurred and I high resolution are the two corresponding 150 x 150 pixel regions around each stimulus from the blurred and high resolution image. An example of this metric for Figure 8 is also presented in Figure 7, alongside the contrast ratio metric.

Conclusion and Future Work

We presented a framework which demonstrates efficient gaze guiding using peripheral stimuli. Our system consists of a personal computer system running the custom software for displaying and free-viewing of natural images by subjects and the interface with an ASL-H6 eye-tracker working at 120 Hz. The images are first displayed in a blurred, i.e. low pass filtered version with a single high resolution patch of 3 • visual angle in the middle and the user deblurs the image step by step by following the stimuli in the peripheral vision with its gaze. The deblurring is done using one of the two proposed strategies for gaze guiding: with and without anticipation. The images viewed by the subjects are divided in a grid and each stimulus is assigned a region including the grid cell containing the stimulus plus other surrounding grid cells. Once the subjects gaze is attracted by a stimulus, its corresponding region is deblurred, i.e. the region is displayed in high resolution. The deblurring process ends for every image once the whole image is displayed in high resolution. The two guiding strategies proposed by us differ mainly by the conditions by which deblurring is effectuated.

We also presented a novel method to determine stimulus saliency which computes in parallel the contrast ratios and the sum of squared differences INRIA Figure 7: Stimulus saliency measures for a selected image. We measure the contrast ratio and sum of squared differences between corresponding patches between corresponding stimuli in a high resolution, respectively a blurred image. The proposed measure can be used to asses stimulus saliency also for other similar applications.

As future direction we could include in the computation of stimulus saliency the eccentricity of stimuli or the local lightness adaptation level of the subject for every given stimulus. In the current implementation of our framework, all blurred images are precomputed with a given Gaussian kernel. Thus, the stimuli, which are the high resolution patches appearing on the blurred images have all an equal low pass filter applied. A possible improved system could take into account local adaptation levels of the subjects visual system and apply a locally adaptive blurring such that the stimuli stay at the detection threshold.

Acknowledgments

I would first like to thank Prof. Noëlle Carbonell who gave me the opportunity to work on this project during my post-doc stay here at LORIA. Sadly, Noëlle passed away in April 2009; she was a wonderful person and a dedicated scientist. I will miss her very much; I enjoyed a lot the lively scientific as well as nonscientific discussions with her.

A lot of our thanks go to Daniel Gepner, who worked earlier on this topic and performed the pilot studies to acquire subject scan paths while free-viewing the Figure 8: Image with the position of stimuli highlighted corresponding to the stimulus saliency measured in Figure 7 images included in our framework. I would also like to thank Olivier Christmann and Jerome Simonin for being such nice and helpful colleagues and good friends during my time at Loria. Last but not least I would like to thank my dear friend Cristian Popi, for all the interesting and inspiring discussions we had in Nancy.

Figure 2 :

 2 Figure 2: The next 30 images in our collection

Figure 3 :

 3 Figure 3: Image showing scan-paths and fixations used for stimulus placement

Figure 4 :

 4 Figure 4: Our stimuli are placed at the numbered positions marked by the red circles. All corresponding surrounding cells are labeled with the same number

Figure 5 :

 5 Figure 5: A schematic description of the mechanism used for our without anticipation strategy. The center shaded circle represents the area of 1 • visual angle occupied by the stimulus. The outer circle, of 2 • visual angle is the area in which if an eye-tracker sample falls the deblurring is activated

INRIA

Figure 6 :

 6 Figure6: A schematic description of the mechanism used for our with anticipation strategy. The center shaded circle represents the area of 1 • visual angle occupied by the stimulus. If the subjects gaze is attracted towards the stimulus from the previous fixation, deblurring occurs when an eye-tracker sample is at least at half the distance between the previous fixation and the stimulus

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

RR n°0123456789


Nadeem Salamat
email: nsalamat@univ-lr.fr

El-Hadi Zahzah
email: ezahzah@univ-lr.fr

SPATIAL RELATIONS ANALYSIS BY USING FUZZY OPERATORS

Keywords: Spatial Relations, Force Histogram, Polygonal Approximation, Temporal Complexity, Fuzzy aggregation operators

 and largely developed by Matsakis [12] who used a quantitative representation of relative position between 2D objects. Fuzzy Allen relations are used to define the fuzzy topological relations between different objects and to detect object positions in images.

Concept for combined extraction of topological and directional relations by using histogram was developed by J.Malki and E.Zahzah [3], and further improved by Matsakis [14]. This algorithm has high computational and temporal complexity due to its limitations of object approximations. In this paper we apply fuzzy aggregation operators for information integration along with polygonal approximation of objects. This approach gives us anew, with low temporal and computational complexity of algorithm for the extraction of topological and directional relations.

INTRODUCTION

Space relations has a remarkable importance in computer vision and image analysis as in Content based image retrieval, similarity based image retrieval [START_REF] Wang | Quantitative Representation of Spatial Relations for Similarity-based Image Retrieval[END_REF], to identify forms, manage data bases, support spatial data in artificial intelligence (AI), cognitive science, perceptual psychology, geography particularly geo-information system (GIS), indexation and comparing objects scene and model are major applications of space relations. P.Matsakis et all [START_REF] Wendling | Linguistic description of relative positions in images[END_REF] discussed different applications usable in satellite, robotics and for developing linguistics expressions. Isabelle Bloch [START_REF] Ralescu | directional relative position between objects in image processing: a comparison between fuzzy approaches[END_REF] discussed different applications in fields of medical imaging, video and vision. Different approaches for finding spatial and topological relations have been developed according to the need for applications and object representations. Qualitative methods for directional and topological relations includes Max J. Egenhofer's method of four intersections [START_REF] Franzosa | Point Set Topological Relations[END_REF]. These methods are considered most important in GIS community. Directional relations are defined on relative frame of reference and absolute frame of reference. In relative frame of reference position of a simple object is made with respect to an oriented line or an ordered set forming a vector to some intrinsic properties of reference object. Methods like angle histogram introduced by K.Miyajima and A.Ralescu [START_REF] Miyajima | Spatial Organization in 2D Images,Fuzzy Systems[END_REF] and statistical method developed by MinDeng.Zalimli [START_REF] Deng | Zalim li, A Statistical Model for Directional Relations Between Spatial Objects[END_REF] depends upon relative frame reference frame. Method introduced by MinDeng.Zalimli is fast but we loose the structural information of objects. Computational complexity for handling the information related to objects increases with increase in dimensions. Matsakis [START_REF] Desachy | Représentation de La Position Relative d' Objets 2D au Moyen D Un Histogramme de Forces[END_REF] introduced 1D representation of 2d objects by the union of longitudinal sections which is the extension of angle histogram. It processes vector and raster data, net objects as well as fuzzy objects [START_REF] Desachy | Représentation de La Position Relative d' Objets 2D au Moyen D Un Histogramme de Forces[END_REF][START_REF] Bailey | Qualitative analysis of sketched route maps: Translating a sketch into linguistic descriptions[END_REF][START_REF] Matsakis | F-Histograms and Fuzzy Directional Spatial Relations, Proceedings[END_REF]. Combining the different knowledge based information increases the efficiency of the system. Approaches for combining informations has been developed for both qualitative representations such as [START_REF] Li | Combining Topological and Cardinal Directional Relation Information in Qualitative Spatial Reasoning, Adaptive and Natural Computing Algorithms[END_REF]and quantitative manners like P.Matsakis. The derivation of combined topological and directional relations by using force histogram [START_REF] Nikitenko | Combined Extraction of Directional and Topological Relationship Information from 2D Concave Objects, in fuzzy modeling with spatial informations for geographic problems[END_REF] was first introduced by J.Malki and E.Zahzah [START_REF] Mascarilla | Indexation et recherche d' image fondées sur les relations spatiales entre objets[END_REF], then P.Matsakis raised some problems regarding fuzziness of some relations like meet and meet_by and some others which exist at segmentation level. Extraction of these relations depend on the fuzzy membership functions µ instead of depending on function φ. In case of longitudinal section consist of segment or union of segments then they are treated by histograms f and F while in case of histogram of Allen relations µ comprises three histograms. When objects are convex and have the connected boundary (when µ acts as φ histogram or f histogram), calculations are very simple but in contrary if objects are concave or objects have the disconnected boundaries then fuzzification of segments is complex (when µ acts as F histogram). As a result temporal and computational complexity increases. Approximating the object by its boundary, length of longitudinal sections can be computed as distance between the intersecting points of oriented line and object boundary. The degree of fuzzy membership function depends upon three values x, y and z. In tuple (x, y, z), the pair (x, z) are the lengths of longitudinal sections and y is the difference between maximum value of intersecting points of object B and minimum value of object A, i.e. y ∈ R or y ∈ Z. By this approach of object approximation, temporal complexity decreases from n √ n to N log(N ) where n is number of pixels of objects under consideration and N is the number of vertex of object polygons. Temporal complexity for the said algorithm is not given but in general temporal complexity of force histograms is discussed in [START_REF] Nikitenko | Applying Soft Computing in Defining Spatial Relations, Understanding the Spatial Organization of Image Regions by Means of Force Histograms[END_REF] for different object types. We assume same temporal complexity for objects because segmentation level problems raised by P.Matsakis forced the object as raster data and in addition to this algorithm for fuzzification of longitudinal sections increases temporal and computational complexity. These problems no more exist if we consider objects by their boundary, then need for P.Matsakis's algorithm remain for objects having disconnected boundaries i.e. objects with holes and for concave objects. Each segment is separated by a certain distance for objects with holes or concave objects. This internal distance has a significant impact on directional and topological relations because values of fuzzy membership functions are distance sensitive(In this case y justifies the grade of a fuzzy membership value). This difficulty is coped off with application of fuzzy disjunctive operators. These fuzzy operators have been developed to summarize information distributed in different sets by grades of fuzzy membership values. This paper is structured as follows. First of all we describe Allen relations in space. In section 2 we describe different fuzzy Allen relations defined by P.Matsakis and changes in mathematical formulation of fuzzy histogram of Allen relations due to object approximation.

In section 3 we discuss different fuzzy operators, section 4 describes experimental results. In section 5 temporal complexity is calculated, section 6 describes different affine transformations and section 7 concludes the paper.

Allen relations in space

Allen in [START_REF] Allen | Maintaining Knowledge about Temporal Intervals[END_REF], introduced the well known 13 mutually exclusive exhaustive interval relations based on temporal interval algebra. These relations are arranged as

A = {< m, o, s, f, d, eq, d i , f i , s i , o i , m i >}. where {< m, o, s, f, d >} resp. ({< d i , f i , s i , o i , m i ,

FUZZY HISTOGRAM OF ALLEN RELATIONS

In real applications, small errors in crisp values can change the entire result when gradual changes of topological relations occur over time. To cope these problems fuzzification was introduced, it comprises the process of transforming crisp values into grades of membership for linguistic terms of fuzzy sets. Fuzzification process of Allen relations does not depend upon particular choice of fuzzy membership function, trapezoidal membership function is used due to flexibility in shape change.

Let r ( i, j) is Allen relation and r is the distance between r ( i, j) and its conceptional neighborhood. We consider a fuzzy membership function µ : r -→ [0, 1]. The fuzzy Allen relations defined by P.Matsakis [START_REF] Nikitenko | Combined Extraction of Directional and Topological Relationship Information from 2D Concave Objects, in fuzzy modeling with spatial informations for geographic problems[END_REF] are

f b (I, J) = µ (-∞,-∞,-b-3a/2,-b-a) (y) f m (I, J) = µ (-b-3a/2,-b-a,-b-a,-b-a/2) (y) f O (I, J) = µ (-b-a,-b-a/2,-b-a/2,b) (y) f f (I, J) = min(µ (-(b+a)/2,-a,-a,+∞) (y), µ (-3a/2,-a,-a,-a/2) (y), µ (-∞,-∞,z/2,z) (x)) f s (I, J) = min(µ -b-a/2,-b,-b,-b+a/2 (y), µ (-∞,-∞,-(b+a)/2) (y), µ (-∞,-∞,z/2,z) (x)) f si (I, J) = min(µ (-(b+a)/2,-a,-a,+∞) (y), µ (-3a/2,-a,-a,-a/2) (y), µ (z,2z,+∞,+∞) (x)) f f i (I, J) = min(µ -b-a/2,-b,-b,-b+a/2 (y), µ (-∞,-∞,-(b+a)/2) (y), µ (z,2z,+∞,+∞) (x)) f d (I, J) = min(µ (-b,-b+a/2,-3a/2,-a) (y), µ (-∞,-∞,z/2,z) (x)) f di (I, J) = min(µ (-b,-b+a/2,-3a/2,-a) (y), µ (z,2z,+∞,+∞) (x)) f Oi (I, J) = µ (-a,-a/2,-a/2,0) (y) f mi (I, J) = µ (-a/2,0,0,a/2) (y) f a (I, J) = µ (0,a/2,∞,∞) (y)
where a = min(x, z), b = max(x, z), x is the length of longitudinal section of argument object A, and z is the length of longitudinal section of reference object B. Most of relations are defined by one membership function and some of them by the minimum value of more than one membership functions like d(during), d i (during_by), f (f inish), f i (f inished_by) and by the definition of fuzzy set in fuzzy set theory, sum of all the relations is one. This gives the definition for fuzzy relation equal, then the definition of Histogram of fuzzy Allen relations stated as in [START_REF] Nikitenko | Combined Extraction of Directional and Topological Relationship Information from 2D Concave Objects, in fuzzy modeling with spatial informations for geographic problems[END_REF]:

( r A F r (q, A q (v), B q (v))dv = x + z w c k=1 m k i=1 n k j=1 [x k i z k j (a k -a k-1 )]r(I k i , j k j ) (1)
where

x = m k i=1 x k i , z = n k j=1 z k j and w = c k=1 m k i=1 n k j=1 [x k i z k j (a k -a k-1 )]
In this case m k and n k represents the total number of segments of argument object A and reference object B and c represents number of loops in algorithem described by P.Matsakis [START_REF] Nikitenko | Combined Extraction of Directional and Topological Relationship Information from 2D Concave Objects, in fuzzy modeling with spatial informations for geographic problems[END_REF]. Here terms multiplied by the Allen relation represents the sum of object areas. In our approach, fuzzification of segments is no more required. Fuzzy Allen relations are computed for each segment. Since each fuzzy Allen relation is the grade of a fuzzy membership value, for simplicity fuzzy Allen relation for each segment is a fuzzy set and fuzzy operators or fuzzy aggregation operators are used to combine different values of fuzzy grades. Hence at the next stage fuzzy aggregation operator is used for integration of different fuzzy Allen relations. Fuzzy operators and fuzzy aggregation operators are discussed in next section. In polygonal approximation of objects, fuzzy Allen relations are calculated for a limited number of segments considering that a complex 2d extended object is the union of small and simple 2d objects like triangles, parallelograms and trapeziums. We find spatial relations only on the boundary of these simple objects and within this region spatial relations do not change so simply generalize the given relations. Definition of fuzzy histogram of Allen relations given by P.Matsakis [START_REF] Nikitenko | Combined Extraction of Directional and Topological Relationship Information from 2D Concave Objects, in fuzzy modeling with spatial informations for geographic problems[END_REF]. "Histogram of fuzzy Allen relations represents the total area of subregions of A and B that are facing each other in given direction θ" Area of subregions of object A and B represented by dark gray color in figure 2 represents a histogram of fuzzy Allen relation. By polygonal approximation, area of objects can be calculated by trapezoidal rules and finally definition of histogram of fuzzy Allen relations can be rewrite as:

( r A F r (q, A q (v), B q (v))dv = (x + z) n k=1 r(I k , J k ) ( 2
)
where z is the area of reference object and x is area of augmented object in direction θ which can be calculated by trapoziodal rules, n is total number of sagments treated and r(I k , J k ) is any Allen relation.

FUZZY OPERATORS AND TREATMENT OF LONGITUDINAL SECTIONS

In real word, objects are complex and they may have different boundary types, they are further divided into convex and concave object classes according to shape. During the decomposition process of an object into segments, there can be multiple segments depending on object shape and boundary which are called longitudinal section. We know that different segments of a longitudinal section are at a certain distance and these distances might effect end results. Low temporal complexity and high precision rates are main objectives of different approaches of finding fuzzy directional relations. Different problem raised by P.Matsakis [START_REF] Nikitenko | Combined Extraction of Directional and Topological Relationship Information from 2D Concave Objects, in fuzzy modeling with spatial informations for geographic problems[END_REF] exists at segmentation level, if we define object by its boundary, then all of these problems can be successfully handled by detecting the object boundary and longitudinal sections is the distance between two intersecting points. But need for the fuzzification process developed by P.Matsakis is still there when object has the disconnected boundary or object is concave. In this case there exist number of 1D segments of concave object or object having disconnected boundary. Each segment and its distance from other segment has its own impact on fuzzy Allen relations of whole object. On the other hand if we follow the process of fuzzification of longitudinal sections, we will be away from our objective. To cope with this, we use fuzzy operators. In this case each fuzzy Allen relation is a member of fuzzy set, therefore fuzzy set theory can be used for information integration. In literature of fuzzy set theory there exist variety of operators such as fuzzy T -norms,T -conorms and so on, which can be used for fuzzy integration of available information. Some mostly used operators are [START_REF] Jo | Fuzzy logic intergration for landslide hazard mapping using spatial data from boeun,korea, symposium on geospatial theory[END_REF] 1.

µ (OR) (u) = max(µ (A) (u), µ (B) (u)); 2. µ (AN D) (u) = min(µ (A) (u), µ (B) (u)); 3. µ (SU M ) (u) = 1 -Π 2 i=1 (µ (i) (u)), 4. µ (P ROD) (u) = Π 2 i=1 (µ (i) (u)) 5. µ (γ) (u) = [(µ (SU M ) (u)] γ * [(µ (P ROD) (u)] 1-γ where γ ∈ [0, 1]
When fuzzy operator OR (respectively AN D) is used, only one fuzzy value contributes for the resultant value which is maximum (respectively minimum). For other operators both values contribute. Fuzzy algebraic sum (product)operator makes the resultant set larger than or equal (less than) the participating value. For the fuzzy γ operator resulting value changes from minimum to maximum values depending on the choice of γ. In our case each Allen relation has a fuzzy grade so we use these operators for integration of fuzzy grade values. Our objective is to accumulate the best available information. We consider that there exist number of segments and each segment has a fuzzy Allen relation with segment of other object. Now suppose longitudinal section of object B has two segments such that z = z 1 + z 2 where z 1 is the length of first segment and z 2 is the length of second segment and z is length of longitudinal section. Let µ (A) (y 1 ) defines the value of fuzzy Allen relations with the first segment and µ (B) (y 2 ) represents value of fuzzy Allen relations with the second segment where y 1 and y 2 are the distances between object A and two segments of B. Now Fuzzy OR operator is used to get consequent information obtained from two sets of fuzzy Allen relations. (Fuzzy γ operator can also be used to make possible the contributions of two fuzzy values but in this case finding compensations values of γ is a problem). For more than two segments composition of operator can be used.

EXPERIMENTS AND INTERPRETATION

For the experiment purpose 360 directions are considered (angle increment is 1 degree) and lines are drawn by the well known Bresenham digital line algorithm. Instead of considering all the v values, we consider only those lines which passes through vertex of polygon. Longitudinal sections are computed and all pairs of segments can be treated simultaneously. Fuzzy Allen relations are computed for each segment, if there exit longitudinal section (More number of segments for an object) then fuzzy aggregation operator is applied to obtain the resultant fuzzy Allen relation of each fuzzy Allen relation. Each relation is associated with the gray scale value like bef ore with black and white represents af ter. We use the same notations as P.Matsakis, except changing the boundary color of each relation for better visualization of relations. Opposite relations have the same boundary color such as m(meet) and (meet_by) relations have the yellow boundary color. Object A has the light gray color while object B is represented by dark gray color. The thirteen histograms that represent directional and topological relations are plotted in the same diagram. For a given angle all the histograms are represented in the layers and each vertical layer represents total area of objects in that direction. Here histograms are not normalized. All relations are symmetric in nature except d(during) and di(during_by).

f AB b (θ) = f AB a (θ + π), f AB mi (θ) = f AB m (θ + π), ; f AB oi (θ) = f AB o (θ + π), f AB si (θ) = f AB f i (θ + π) f AB f (θ) = f AB s (θ + π),

TEMPORAL COMPLEXITY

Finding the exact temporal complexity is a tough job, major aim of this study is to find time length or amount of memory required by the algorithm. We need a language that expresses the computational time as a function of N , grows on the order of N . Five symbols for comparing rates are used such as o ,O ,θ ,Ω and ∼ (Asymptotically equal,irreverently, tiddle). In fact asymptomatically equality is an formalism of idea to find the conditions that two functions have same growth rate.i.e.

limit n→∞ ( an b n ) = 1 and a n = O(b n ) if | an b n |is bounded
We find temporal complexity of a algorithm by asymptotic analysis. For this purpose we find a function which represents upper bound of our function. In our case time constraint depends upon length of line, contour length and number of polygon vertexes. We summarize time for all 360 directions when line length is fixed to 1000,1200,1400 pixels. Following tables represent different computations where L is length of line and N represents number of polygon vertexes.

CONCLUSION

It is shown that histogram of fuzzy Allen relations associated with pair of objects carry a lot of information. To deal with concave objects or objects having disconnected boundaries, fuzzy operators are used. Use of these operators is simple so polygonal approximation of objects and application of fuzzy aggregation operator simplifies the algorithm given by P.Matsakis [START_REF] Nikitenko | Combined Extraction of Directional and Topological Relationship Information from 2D Concave Objects, in fuzzy modeling with spatial informations for geographic problems[END_REF]. This approach decrease its temporal and computational complexity due to avoiding the fuzzification process developed by P.Matsakis. Certainly this approach of using fuzzy operator will open new fields of applications for fuzzy aggregation operators. Here we calculated all the directions for experimental purpose, in practice we performed only limited number of directions according to the requirement of application. Allen relations are used for describing the relative object position in image understanding.

  >}) are the relation meet, overlap, start, finish, during (resp the inverse relations of the cited ones). The relation eq is the equality spatial relation. All the Allen relations in space are conceptually illustrated in figure extracted from 1. These relations have a rich support for the topological and directional relations.

Fig. 1 .

 1 Fig. 1. Black segment represents the reference object and gray segment represents argument object. figure extracted from[START_REF] Desachy | Représentation de La Position Relative d' Objets 2D au Moyen D Un Histogramme de Forces[END_REF]

Fig. 2 .

 2 Fig. 2. Gray area represents the histogram of Allen in direction θ

  and for d(during) and di(during_by) it will be f AB d (θ) = f BA di (θ).

Fig. 3 .Fig. 4 .

 34 Fig. 3. (a)Explanation of gray level value associated with a relation (b)Object pair representation. (c) Corresponding histogram

Fig. 5 .

 5 Fig. 5. Convex and concave pair of objects and their corresponding fuzzy histogram of Allen relations

Fig. 6 .

 6 Fig. 6. (a)Graph of some known functions,(b)Graphs of data given in table No.1(c)Graphs of data given in table No.2

Table 1 .

 1 contour of 1300 pixels

	N L 1000 1200 1400
	24	63 67.14 73
	25 66.5 73	78
	26	68 74.3 80
	27	72	79 84.34
	28 73.20 82 86.20

Table 2 .

 2 contour 3300 of pixels

N L 1000 1200 1400 25 97.5 103 113.5 26 102 107 115.5 27 109 115.5 123 28 113 119 126 29 120 125 129.4 30 124 129 135

If we analyze keenly the data in table 1 and table 2, there will be a certain symmetry between the different values of cost function (Time function) and the number of polygons vertexes.

  Now if we observe the graphical representation of data,(graph fig.6(a) ) of data in table 1 and graph fig.6(b) of data in table2) each time graph for a fixed length of line and given objects sizes (length of contours). It seems that graph is displaced by a constant value of T . It seems that graphically function f (n) = nlog(n)) representation the upper bound of our graphs. (Graphes given infig.5) So histogram of fuzzy Allen relations are of order O(N log(N ))

table No.1(c)Graphs of data given in table No.2


Nadeem Salamat
email: nsalam01@univ-lr.fr

El-Hadi Zahzah
email: ezahzah@univ-lr.fr

Fuzzy histogram for internal and external fuzzy directional relations

Keywords: Spatial Relations, Direction Relations, fuzzy internal cardinal directions, Scale Space, Medical Images

Spatial relations have key point importance in image analysis and computer vision. Numerous technics have been developed to study these relations especially directional relations. Modern digital computers give rise to quantitative methods and among them fuzzy methods have core importance due to handling imprecise knowledge information and vagueness. In most fuzzy methods external directional relations are considered which are useful for small scale space image analysis but in large scale space image analysis internal cardinal directions are also important. In this paper a new function is introduced which can be equally used for both internal and external fuzzy directional relations. Trigonometric function is used to determined the directions. This approach has core importance in medical images. A practical example in medical images is discussed.

Introduction

One of the fundamental tenet of modern sciences is that the phenomenon can not be claimed to be well understand until it can be characterized in quantitative manners. Given our veneration for what is precise and quantitative. Advent of digital computers has resultant in a rapid expansion in the use of quantitative methods through most of human knowledge. Un precise knowledge information and vagueness gives rise to fuzzy methods. In the last decade of 19 th century much attention has been paid to spatial relations between objects and directional relations has a key point importance among them. These relations are used in computer vision and image analysis [START_REF] Ralescu | Directional Relative Position Between Objects in Image Processing: A Comparison Between Fuzzy Approaches[END_REF], [START_REF] Héde | Using Relative Spatial Relationships to Improve Individual Region Recognition, European Workshop on the Integration of Knowledge, Semantics and Digital Media Technologies[END_REF], content based image retrieval (CBIR) [START_REF] Matsakis | The Use of Force Histograms for Affine-Invariant Relative Position Description[END_REF], similarity based object retrieval [START_REF] Wang | Histograms: Efficient Representation of Spatial Relations Between Objects of Arbitrary Topology[END_REF], identify forms, manage data bases, support spatial data in artificial intelligence (AI), cognitive science, perceptual psychology, geography particularly geoinformation system [START_REF] Egenhofer | Categorizing Binary Topological Relations Between Regions,Lines, and Points in Geographic Databases[END_REF], [START_REF] Egenhofer | A Critical Comparison of The 4-Intersection and 9-Intersection Models for Spatial Relations: Formal Analysis[END_REF], [START_REF] Roop | Extended Spatial Objects[END_REF], comparing objects scene and model and indexation space are major applications of space relations. Satellite and robotics [START_REF] Wendling | Linguistic Description of Relative Positions in Images[END_REF], developing linguistics expressions and medical imaging applications [START_REF] Nikitenko | Applying Soft Computing in Defining Spatial Relations, Understanding the Spatial Organization of Image Regions by Means of Force Histograms[END_REF]. Number of approaches has been developed to study the directional relations qualitatively and quantitatively. All of quantative methods are sensitive to distances between objects which is proportional to distribution range of angles between objects. Some methods provide only distribution rang of directions among the objects such as statistical method [START_REF] Deng | Zalim li, A Statistical Model for Directional Relations Between Spatial Objects[END_REF]. This method supports topological relation disjoint,meet and metric relations N ear, F ar and F arof f . Another class of methods which capture shape information along with the distribution range of directions like angle histogram [START_REF] Miyajima | Spatial Organization in 2D Images,Fuzzy Systems[END_REF], force histogram [START_REF] Matsakis | F-Histograms and Fuzzy Directional Spatial Relations[END_REF]. Both angle histogram and angle distance histogram supports meet relation along with disjoint and metric relations N ear, F ar and F ar of f . Force histograms strictly support disjoint relation and other methods support overlapping objects like R-histogram [START_REF] Wang | Histograms: Efficient Representation of Spatial Relations Between Objects of Arbitrary Topology[END_REF] and R * -histogram [START_REF] Wang | R-histogram: quantitative representation of spatial relations for similarity-based image retrieval[END_REF]. These are quantitative methods but not fuzzy moreover labeling process is to complex. Main objective of this paper is to introduce a method which can compensate the overlapped objects like the R -histogram. It is a quantitative representation of the relative objects position. All of the above cited methods are used in small scale space but in large scale space internal object position is also important. It is needed to develop a single method which can be used in small and large scale spaces unlike separate qualitative method internal 5 -intersection, internal 9 -intersection and internal 13 -intersection models introduced by Y. Liu et al. [START_REF] Wang | Qualitative Spatial Reasoning About Internal Cardinal Direction Relations[END_REF]. For external directions separate methods are used. Another important fact that meet is a crisp point relation and only one overlapping pixel can change the entire phenomenon and stop the method of force histograms from working. Paper is structured as follows, in section 2 proposed function and calculation of histogram is described, in section 3 concerns about experiments. Section 4 concludes the work.

Oriented lines, Used Function and calculation of histogram

This section is concerned with object decomposition, use of function and treatment of different parts of object and finely calculation of histograms.

Oriented lines, segments and longitudinal sections

Let A and B are two objects. (v, θ) ∈ R, where v is any real number and θ ∈

[0, 360] or θ ∈ [-π, π]. ∆ θ (v)
is an oriented line at orientation angle θ as explained in figure 1(a). A ∩ ∆ θ (v) is the intersection of object A and oriented line ∆ θ (v). It is denoted by A θ (v), called segment of object A and its length is x. Similarly for object B where B ∩ ∆ θ (v) = B θ (v) is segment and z is its length. y is the difference between the maximum value of B ∩ ∆ θ (v) and minimum of A ∩ ∆ θ (v)(for details see [START_REF] Matsakis | A New Way to Represent the Relative Position between Areal Objects[END_REF]). In case of polygonal object approximations, triplet (x, y, z) can be calculated from intersecting points of line and object boundary, only oriented lines which passes through vertexes of polygon are taken.

Use of function

To define the object position a function is defined.

f : R + × R × R + -→ [0, 1] mathematically f (x, y, z) =      0 if y ≤ -z/2 2y+z z if y ∈ ( -z/2 0 ) 1 ify ≥ 0 where (x, y, z) ∈ R + × R × R + .
This function can easily be modeled as fuzzy membership function.

Treatment of points

If both the objects A, and B are point objects or oriented line ∆ θ (v) passes through the vertex of polygons. Length of both segments will be zero. i.e x = 0,z = 0. If both the points are same then y = 0 in this case f (x, y, z) = 0 and if both points are not same then either A is before B or after. In first case y < 0 and f (x, y, z) = 0 and in second case y > 0 and f (x, y, z) = 1. Third possible case if x = 0 and z = 0 then value off (x, y, z) depends on y. If y ∈ (-z 2 , 0) then f (x, y, z) ∈ (0, 1) and for y > 0 f (x, y, z) = 1 otherwise f (x, y, z) = 0

Treatment of segment

It is considered that the length of segments x and z are two points on real line and they form the support of a fuzzy membership function, y is taken as its evaluation point. Calculation of x,y and z are explained in above cited figure 1(b). Then given function can be modeled as fuzzy membership functions, for this purpose trapezoidal membership function can be used.

f (I, J) = µ (-z/2,0,∞,∞) (y)

Treatment of longitudinal sections

In decomposition process of an object into segments, there can be multiple segments depending on object shape complexity which is called longitudinal section. Each segment of a longitudinal section is at a certain distance and these distances might effect end results. Different fuzzy T -norms,T -conorms are used fuzzy integration of available information. Here for cause of simplicity Standard t -conorms is used. For standard t -conorm only one value contributes.

µ (SU M ) (u) = 1 -Π 2 i=1 (µ (i) (u)) 4) Fuzzy Algebraic Product: µ (P ROD) (u) = Π 2 i=1 (µ (i) (u)) 5) Fuzzy γ Operator: µ (γ) (u) = [(µ (SU M ) (u)] γ * [(µ (P ROD) (u)] 1-γ
where γ [0, 1] When fuzzy operator OR (respectively AN D) is used, only one fuzzy value contributes for the resultant value which is maximum (respectively minimum). For other operators both values contribute. In this case each relation has a fuzzy grade, objective is to accumulate the best available information. In case of longitudinal section, there exist number of segments and each segment has a fuzzy relation with segment of other object. Now suppose longitudinal section of object A has two segments such that x = x 1 + x 2 where x 1 is the length of first segment and x 2 is the length of second segment and x is length of longitudinal section. Let µ 1 (y 1 ) defines the value of fuzzy directional relations with the first segment and µ 2 (y 2 ) represents value of fuzzy directional relation with the second segment where y 1 and y 2 are the distances between object A and two segments of B. µ(y) = µ 1 (y 1 ) µ 1 (y 1 ) where is a fuzzy operator. Existence of longitudinal section is explained in above cited figure 3(f).

Histogram calculations

The above cited function is used to detect the object. If the argument object is within the object and in the lateral half of segment then it will assign a value in (0, 1). If objects are disjoint then it will assign a value 1 and if argument object is in the first half of an object then it will assign the value zero. Let A and B be two objects, and (v, θ) are two reals (I, J) are the one of argued and reference object segments. f (I, J) represents the argument weight that A is in θ direction of B.

Histogram will be sum of all the arguments and represented as:

+∞ -∞ (F (θ, A θ (v), B θ (v))dv = n j=1 f (I j , j j ) ( 1
)
where n is total number of segments and histogram is represented by

F AB (θ) = +∞ -∞ (F (θ, A θ (v), B θ (v))

Experiments and interpretation

For experiment purpose, first of all its sensitivity to shape and distance is verified for disjoint objects. In this example position of object A is changed and proposed histogram is compared. Second set of examples is concerned with comparing different fuzzy methods for directional relations when objects are disjoint but when object A inside the object B, no fuzzy method exist, we show that proposed histogram F can be used for internal fuzzy directional relations.

f rightof = cos 2 (θ) if θ ∈ [-π 2 , π 2 ] 0 otherwise
Results are compared in table 1. Another important thing is that directional relations are concerned with the number of visual points. In quantitative methods final decision depends upon number of visible points. When objects are disjoint, range of directional distribution are limited but in case argument object is inside the reference object then every point of reference object is visible from some points of argument object. This phenomenon will change the process of reasoning. And object is visible in all directions. All the directions are calculated and for conclusion we take only two maximum directions. As in case of figure 4

Conclusion

Fuzzy directional relations have key point importance in image analysis. The most developed fuzzy approaches deals with disjoint and meet topological relations. Only R -histogram deals with overlapping topological relation between object which is not fuzzy. We have developed a function. This function is sensitive to internal cardinal directions and equally sensitive to external cardinal directions. P roposed histogram can be successfully used to detect object position within the object or outside the object. This approach is simple and easy to implement. When the objects are disjoint this approach of object detection is similar to the F 0 (constant histogram forces introduced by Matsakis [8]. This method can in the same manner be employed on both large and small scale image analysis.

Fig. 1 :

 1 Fig. 1: (a)Oriented line,(b) Treatment of segment,(c) Treatment of longitudinal section.( Matsakis [9])

1 )

 1 Fuzzy OR: µ (OR) (u) = max(µ (A) (u), µ (B) (u)) 2) Fuzzy AND: µ (AN D) (u) = min(µ (A) (u), µ (B) (u)) 3) Fuzzy Algebraic sum:

  dv, for normalization divide by the max value of F. • Position invariance: It is invariant and in general F θ (A, B) = F θ (B, A) + π • Shape dependance: Histogram does not depend on object shape different objects may have same histogram. • Distance dependance: Histogram is also sensitive to distance between objects

  In this figure objects are represented at different positions and fig2(c) shows that histogram is sensitive to shape and distance. fig 2(a) represents the objects as in source [6], fig. 2(b) represents objects when object A is translated by (30, 0) and fig 2(c) represents their histograms where black solid line represents histogram between objects given in fig2(a) and dotted line represents histogram between objects given in 2(b).

Fig. 2 :

 2 Fig. 2: (a)(b) objects A,B at different positions(d) Corresponding histograms

Fig. 3 :

 3 Fig. 3: (a) kidney image with a stone inside(source:www.sodahead.com). pairs of objects with object A at different positions

  (a), Stone lies below_lef t side of the kidney and in case of figure 4(b) stone lies at below_right side of the kidney. Results are verified in table 2

  .

Fig. 4 :

 4 Fig. 4: Pairs of objects with object A at different positions inside the object B

Table 1 :

 1 Comparison table when objects are disjoint

	1 y r	otherwise
	because for this case, every pair (v, θ), y is always negative.
	Resulting that force histogram have zero value and this
	restricts the application of force histogram for objects inside
	each other i.e. for topological relation overlap and contain

or contained_by. Following different cases are considered when object A is inside object B. Results are compared in table

2

.

Table 2 :

 2

	Object above	right	below	lef t
	4(a)	.13	.15	.45	.28
	4(b)	.1	.3	.48	.12
	4(c)	.18	.14	.2	.48
	4(d)	.3	.48	.15	.07
	4(e)	.5	.3	.07	.13

Calculation of directional relations when object A is inside B by the proposed fuzzy histogram

Dir. relations object A is inside B