Martin Rodriguez-Vega

Carlos Canudas-De-Wit

Hassen Fourati

A Route-based Method for Turning Ratio Estimation: Application to the Grenoble Downtown Traffic Flow and Density Reconstruction

In this paper, we propose a route-based approach to improve the flow and density estimation methods for the case of urban traffic networks. For this proposal, a traffic assignment problem is solved at first, whose outputs are used to estimate the turning ratios for all intersections. Such information is not usually available at each node of the network. Second, the estimated turning ratios are used to better reconstruct the dynamic state of the network: the flow and density. To validate the proposed methods, we use real traffic data collected in the city of Grenoble in France, which include the measurement of Origin-Destination matrices, some turning ratios for validation, mean speeds of road sections and traffic flows at the boundaries of the network.

I. INTRODUCTION

Actually, the development of vehicle mobility in smart cities in one of the main problems in transportation research. Methods to estimate the flow and density of roads in urban networks is of great importance to promote such progress in this area. The widespread use of novel heterogeneous data sources such as Floating Car Data (FCD), Bluetooth (BT) and Wifi vehicle identifiers, among others, has helped in the development of data-based estimation approaches that provide more accurate estimates of traffic states, i.e. speed, flow, and density, in each road of a traffic network [START_REF] Seo | Traffic state estimation on highway: A comprehensive survey[END_REF]. Then, the estimates are used by city operators to implement strategies for congestion reduction or by users to decide about the fastest route to take to reach the destination.

Currently, many of the state-of-the-art methods for traffic states estimation require the knowledge of the Turning Ratios (TRs), which model how flow split at intersection from the inbound to outbound roads. For instance, in [START_REF] Shahrbabaki | State estimation in urban traffic networks: A twolayer approach[END_REF] the authors use the TRs with data from connected vehicles to estimate the flow propagation in an urban network. Similarly, in [START_REF] Rodriguez-Vega | Dynamic density and flow reconstruction in large-scale urban networks using heterogeneous data sources[END_REF], the authors propose a data-based density and flow estimation method, using heterogeneous data sources (flow and TR sensors with FCD) which are suitable for applications in large-scale networks.

In practice, as the number of intersections in a network can become very large, the TRs are quite hard to know everywhere, then to circumvent this problem, the interest is focused on their estimation based on some observed variables. For instance, node-based approaches use the local characteristics of the road adjacent to an intersection to calculate the TRs. The authors in [START_REF] Furth | Model of Turning Movement Propensity[END_REF] use geometrical properties of roads to identify the most likely turns, [START_REF] Krumm | Where will they turn: Predicting turn proportions at intersections[END_REF] uses * Univ. Grenoble Alpes, CNRS, INRIA, Grenoble INP, GIPSA-Lab, 38000 Grenoble, France { martin.rodriguez-vega, carlos.canudas-de-wit, hassen.fourati }@gipsa-lab.fr the connectivity of each outgoing road (dead-ends are lesslikely to be taken than an artery), and [START_REF] Rodriguez-Vega | Dynamic density and flow reconstruction in large-scale urban networks using heterogeneous data sources[END_REF] uses the Functional Road Classification (FRC) [START_REF] D'andrea | A functional road classification with data mining techniques[END_REF] which assigns importance to each road according to its function [START_REF]Highway function classification concepts, criteria and procedures[END_REF]. Although node-based approaches provide an inexpensive way to estimate the TRs, they can introduce biases as they do not take into account the behavior of drivers that travel from Origins to Destinations.

To solve this problem, the Traffic Assignment (TA) can be used as an alternative approach for the estimation of TRs. Such method considers the flow of vehicles in an Origin-Destination (OD) pair, and distributes it into the possible routes that make this connection according to some criteria (e.g. minimize the total traveling time of each driver) [START_REF] Saw | Literature review of traffic assignment: Static and dynamic[END_REF]. In the literature, TA receives a lot of attention, in the case of traffic and infrastructure planning [START_REF] Chiu | Dynamic Traffic Assignment: A Primer[END_REF], and the prediction of congestion patterns [START_REF] Wada | Network throughput under dynamic user equilibrium: Queue spillback, paradox and traffic control[END_REF]. Although TA is a powerful tool to understand path selection and congestion formation, it is computationally intensive, specially for large-scale traffic networks as the number of paths to consider grows exponentially. Some works in the literature propose more efficient algorithms to solve TA problem [START_REF] Rojo | Evaluation of traffic assignment models through simulation[END_REF], [START_REF] Hoang | A linear framework for dynamic user equilibrium traffic assignment in a single origindestination capacitated network[END_REF], but they are restricted to specific networks and may not applicable for general networks. Thus, to the best knowledge of the authors, TA is typically unsuited for real-time flow and density estimation applications in large-scale traffic networks.

In this paper, we propose a new method to estimate the TRs using a route-based TA approach, and use this information later to estimate, in real-time, the flow and density of traffic in a large-scale network of city. The proposed TA solution takes as inputs measured OD matrices and boundary flows of the network. As an optional input, few measured TRs can be included to provide more accurate results, if available. At this step, the calculation is done offline which better reduce the cost computation. Once all TRs are obtained, the real-time flow and density estimator proposed in [START_REF] Rodriguez-Vega | Dynamic density and flow reconstruction in large-scale urban networks using heterogeneous data sources[END_REF] is used, which consists in using speed measurements from FCD in conjunction with the current density estimates to calculate road outflows. We validate the proposed methods using real traffic data collected in the downtown of Grenoble city, in France. The results show that the proposed method produces more accurate estimates when compared to node-based approaches.

The paper is organized as follows: Section II presents the main notions and definitions used throughout the paper. Section III describes the proposed methods and the subprocesses needed for its implementation. Section IV presents the study case and real traffic data used for validation of the method, and Section V discusses the obtained results. Finally, Section VI presents some conclusions.

II. DEFINITIONS, PRELIMINARIES AND MODELS

We define a traffic network as a set {N, E, T } where N = {1, 2, . . . , N } are the nodes corresponding to intersections (or nodes), E = {1, 2, . . . , E} are edges corresponding to the roads (or links), and T ⊂ E×E is a set of road pairs (i, j) that correspond to turns from i to j that can legally be made at an intersection. Each road i ∈ E is characterized by physical properties such as length ℓ i , free-flow speed v max i , number of lanes, etc. The state of road i refers to the values of the inflow φ in i (t), outflow φ out i (t), and density ρ i (t). The entire traffic state is the collection of these quantities in vector form φ in (t), φ out (t), ρ(t) ∈ R E , respectively, which are linked by the conservation equation

ρ(t) = L -1 (φ in (t) -φ out (t)) (1)
where L is the matrix of road lengths, L i,i = ℓ i , and the i-th element of ρ(t) is ρ i (t) (analogously for φ in (t) and φ out (t)).

Denote by E in ⊂ E the sets of roads at the boundaries from which cars enter the network (flow origins), and by E out ⊂ E the set of roads at the boundaries from which cars exit the network (flow destinations), respectively. To describe how flows split at intersections, we assume the existence of the TR parameters, which apply locally at each intersection. A TR r i,j ∈ [0, 1] is the proportion of vehicles that turn from road i to road j, such that (i, j) ∈ T . A simple intersection is shown in Fig. 1. Thus, the incoming flow for each road can be expressed in terms of the outflows of the upstream roads

φ in j (t) = i|(i,j)∈T r i,j φ out i (t) , ∀j ∈ E \ E in . (2)
Note: The super-indexes in and out for the flows correspond to roads and not intersections. Thus, φ out i is the outflow of road i (that then enters an intersection), and φ in i is the inflow to road i (that comes from an intersection).

To estimate the TRs, we use a route-based approach which is based on a traffic assignment problem using the OD matrix of the network. A route (or path) for a given OD pair (o, d) ∈ E in × E out , is an ordered set P of adjacent roads that begins with o and ends with d, i.e., P = (o, i 1 , i 2 , . . . , d), where i k ∈ E. Multiple paths may exist for each OD pair. Consider the example shown in Fig. 2, and the OD pair (4, 8), for which two possible paths are P 1 = (4, 3, 2, 8) and P 2 = (4, 6, 7, 8). We denote by P the set of all considered paths for all OD pairs. Assume that there is an indexing p = 1, 2, . . . , |P| of the paths. Each path P p ∈ P has an associated steady-state flow f p , and let f ∈ R |P| be the concatenated vector of path flows. These flows are constrained by the given OD flows: define ψ (o,d) as the steady-state OD flow between the pair (o, d), then, the sum of all path-flows connecting o to d must be equal to ψ (o,d) ,

ψ (o,d) = |P| p=1 Ω (o,d) p f p (3)
where

Ω (o,d) p = 1 if path p connects the OD pair (o, d)
, and is 0 else. Written in matrix form, this becomes

ψ = Ωf ( 4
)
where ψ is the concatenated vector of all OD flows, and Ω is called the path-OD incidence matrix.

Similarly, each road i ∈ E has a unique steady-state flow φ in i = φ out i = φ i , which are given by the path flows: the flow in road i is the sum of all the path-flows that contain road i,

φ i = |P| p=1 L i,p f p (5)
where

L i,p = 1 if i ∈ P p 0 else . (6)
In matrix form, this is written as

φ = Lf ( 7
)
where L is the link-path incidence matrix and maps the pathflows f to link flows φ.

III. ROUTE-BASED ESTIMATION METHOD

We propose to use a route-based method to estimate the TRs, which can be done offline using average OD flows during a long aggregation period (e.g. daily). As an 1. Path creation 2. Traffic Assignment 3. TR calculation 4. Flow and density estimation

P, L, Ω f r i,j E, T ψ P, L TR measures ri,j , T * Input demands φ ext (t) Speed measures v i (t) ρ(t), φ in (t), φ out (t)
Fig. 3: Overview of the proposed methodology.

application, the TR estimates are used as an input to a flow and density reconstruction method. The overall diagram of the proposed methodology is shown in Fig. 3 where the main steps are 1) Create a list of paths that connect each origin to each destination. 2) Taking as input the OD flows, solve the TA problem to determine the flow in each path. If available, TR data from measured intersections can be included. 3) Use the computed path flows to estimate the TRs for each intersection. 4) Use the TR values with boundary inflow and speed measurements to estimate the road flows and densities. The details of each step are described below.

A. Path creation

Route-based methods require first the creation of a set of considered paths P between each OD pair. To construct this set, we select the paths with the smallest traveling time in free-flow,

c(P) = i∈P ℓ i /v max i . ( 8
)
This problem is solved by the well-known Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF], which gives a set of K shortest (or in this case fastest) paths for a given (o, d) pair. In practice, this calculations can be time consuming, so we introduce a slightly modified version of the algorithm to speed-up computations.

A new condition is added such that if the new found paths are too slow in comparison to the fastest path, the algorithm stops. Our version will find at least K min paths between o and d, and after K min paths have been found, the algorithm will continue to look for paths until the new path is ϵ times slower than the fastest path. If K max paths are found, the algorithm stops regardless of the length-tolerance condition. This modification allows to have a variable number of paths for each pair, so the total number of iterations is reduced while keeping relevant (not too slow) paths.

Algorithm 1. K-fastest paths

Inputs:

• Road graph {E, T } • Origin link o ∈ E in . • Destination link d ∈ E out . • Minimum number of paths K min • Maximum number of paths K max • Length tolerance ϵ > 1
Initialize:

• Set of optimal paths P = ∅.

• Number of paths found to each link n i = 0, ∀i ∈ E.

• Stack of paths B = {{o}} Steps:

• While B ̸ = ∅ and n t < K max -Find P ∈ B such that c(P) ≤ c(P ′ ) ∀P ′ ∈ B -Let i be the final link in P.

-

n i ← n i + 1 -B ← B \ {P}. -if i = d * P ← P ∪ {P} * if c d ≥ K min and c(P)/min(c(P)) > ϵ • Return P. -if c i ≤ K max * For each j such that (i, j) ∈ T • B ← B ∪ {P ∪ {j}} • Return P.

B. Traffic assignment

The TA problem consists in determining the flows in each path of a considered set P, given a vector of OD flows ψ. We assume that drivers choose their route such that it minimizes total travel time. This is referred to as User Equilibrium (UE), whose solutions are constrained by the Wardrop criteria [START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF]. In the literature, TA problems under UE are commonly written as the following non-linear convex optimization problem [START_REF] Sheffi | Urban Transportation Networks: Equilibrium Analysis with Mathematical Programming Methods[END_REF],

min f i∈E φi 0 t i (x)dx s.t. φ = Lf , Ωf = ψ φ ≤ φ max , f ≥ 0 (9)
where φ max is a vector of road capacities, and t i (φ) is the travel time in road i for a given flow. In this work, we use the commonly used BPR function [START_REF] Zhang | Analytical Model for Travel Time-Based BPR Function with Demand Fluctuation and Capacity Degradation[END_REF] defined as

t i (φ) = ℓ i v max i 1 + β φ φ max i α (10)
where α and β are tuning parameters. This function was first proposed by the Bureau of Public Roads, now the Federal Highway Administration (FHA) in the United States [START_REF] Zhang | Analytical Model for Travel Time-Based BPR Function with Demand Fluctuation and Capacity Degradation[END_REF]. Common choices for the parameters are α = 4, β = 0.15.

This TA problem formulation can be modified to include additional TR information. Suppose some intersections are equipped with sensors such that measured TRs ri,j are available for turns (i, j) in a set T * ⊂ T (see Section IV-B). We can modify [START_REF] Chiu | Dynamic Traffic Assignment: A Primer[END_REF] by including this data to obtain a new optimization problem. From the TR definition, for each measured ratio there is an equation of the form φ i,j -ri,j φ i = 0

where φ i,j is the flow on turn (i, j) ∈ T * . These flows can be written in terms of the path flows,

|P| p=1 L i,p L j,p f p -ri,j |P| p=1 L i,p f p = 0, (12)
where the first term is the sum of flows in all paths that contain both roads i and j, and the second term is the sum of flow in paths that contain road i. This set of equations can be written in matrix form as

Q(T * )f = 0, with Q k,p = L i k ,p (L j k ,p -ri k ,j k ) (13)
where (i k , j k ) is the k-th element of T * . We introduce this condition as a soft constraint in the objective function of (9), min

f i∈E φi 0 t i (x)dx + γ||Q(T * )f || 2 (14)
under the same constraints, where γ is a tuning parameter to control how relative important the satisfaction of the TR measures are in comparison to the regular TA cost.

C. Turning ratio calculation

The resulting vector of path flows f can be used to estimate the TR values for each intersection. Consider an arbitrary turn (i, j) ∈ T , whose corresponding TR is calculated by

r i,j =   |P| p=1 L i,p L j,p f p   |P| p=1 L i,p f p . ( 15
)
This is the ratio of the total flow of paths containing both roads i, j, and those containing road i. As the former set of paths is always a subset of the latter, this ratio is in [0, 1], unless all the incoming flows are zero.

D. Flow and density estimation

In a previous work [START_REF] Rodriguez-Vega | Dynamic density and flow reconstruction in large-scale urban networks using heterogeneous data sources[END_REF], we proposed a data-based method to estimate flows and densities for all roads, using speed, turning ratio and external input demand data from sensors. In the previous section, (2) allows to write road inflows in terms of the turning ratios, except for the roads at the boundaries of the network E in . For these roads,we have

φ in i (t) = φ ext i (t) ∀i ∈ E in (16)
where φ ext i (t) is the external input demand for road i, and are assumed to be measured. Define the TR matrix R ∈ R E×E whose elements R i,j are the TRs if (i, j) ∈ T or are 0 else. Then (2) and ( 16) can be written in vector form as

φ in (t) = R ⊤ φ out (t) + φ ext (t) (17)
where φ ext (t) is a vector whose i-th entry is φ ext i (t) if i ∈ E in , or is 0 else.

We assume that the speed of each road v i (t) is known, and let V (t) be a diagonal matrix with V i,i (t) = v i (t). Consider the hydrodynamic relation

φ out (t) = V (t)ρ(t). (18)
Substitution of ( 17) and ( 18) into (1) yields the following open-loop density estimator,

ρ(t) = L -1 (R ⊤ -I)V (t) ρ(t) + L -1 φ ext (t). (19)
In [START_REF] Rodriguez-Vega | Dynamic density and flow reconstruction in large-scale urban networks using heterogeneous data sources[END_REF], we proved that under error-free measurements and if the rate of change of the mean speed V (t) is slow enough, this estimator converges asymptotically.

IV. EXPERIMENTAL SETUP

In this paper, we use data provided by the GTL-Ville 1 , which is an experimental platform developed by the ERC Scale-Freeback project [START_REF]Scale-Free Control for Complex Physical Network systems[END_REF] for the real-time collection of traffic data coming from a network of installed sensors in the city of Grenoble in France. It covers an area of 1.4 Km by 1 Km in the city center, which is composed of 463 intersections, 804 road sections, and 1234 turns. The information collected using heterogeneous data sources are:

A. Flow sensors

Two sensor technologies are currently being used to measure the vehicular flow (number of vehicle passages during a period of time) in a selection of roads: magnetic induction loops, and microwave radars. There are 40 flow sensors inside the region of interest as shown in blue, red, and green in Fig. 4. These colors correspond to the classification of the data provided by the sensors:

• In blue, sensors that provide boundary inflow data, located in roads in E in . • In red, sensors that provide boundary outflow data, located in roads in E out . • In green, sensors that provide data for cross-validation, located in roads on the network's interior. Flow data is available in real-time with sampling rate of 1 min.

B. Turning Ratio sensors

To measure the TR parameters for a selection of 12 intersections, Bluetooth (BT) reader devices were located at the adjacent roads. These sensors detect the unique identifiers (IDs) of vehicles equipped with other BT devices and match these IDs with the other sensors to obtain the inbound and outbound road for each detected vehicle. To protect user privacy, only the total number of detected vehicles making each turn during a period of 1 hour is available. The selected intersections where TR sensors were installed are shown as gray circles in Fig. 4 with corresponding label numbers from 1 to 12.

TR data was collected during measuring campaigns, where data for each intersection was collected during a week, between September 21 and October 28, 2020. The TR values are computed as ri,j = Counts(i, j)

k Counts(i, k) (20)
where Counts(i, j) is the number of detected vehicles that

C. Origin-Destination matrix sensors

BT sensors placed in specific locations were used to measure the OD matrices. In total, 11 sensors were located around the boundaries of the network, at the most important roads. OD data was collected during a two-weeks measuring campaign from February 24 to March 11, 2021. As BT devices only detect a sample of the population, only the OD rates can be estimated, using an analogous equation to (20). The location of sensors are shown as gray circles labeled from 1 to 11 in Fig. 5.

In the figure, the colored lines show the computed OD rates between each OD pair. The line color and thickness represents the percentage of vehicles going from each origin to all destinations. The ratios are normalized with respect to the origins, i.e., the sum of rates from one origin to all destinations is 1. As OD pairs are bidirectional (roads can have both directions), the right-hand driving rule is used to separate direction, i.e., for a vertical pair, the north-tosouth direction is to the left, and for an horizontal pair, the east-to-west direction is on top. The data agrees with the expected behavior, as roads with high importance receive a much higher rate than least important ones. For instance, the OD pairs (1,5) and (5,1) are connected by the most important artery in this zone which sees very high flows, and have the highest OD rates in the figure. A similar behavior is seen for pairs (2,9) and [START_REF] Krumm | Where will they turn: Predicting turn proportions at intersections[END_REF][START_REF] Saw | Literature review of traffic assignment: Static and dynamic[END_REF].

D. Floating car data

Floating Car Data (FCD) refers to vehicle traces that are collected for a subset of drivers, usually obtained from GPS devices. For the region of interest, data is obtained through a third party, Tomtom2 , a well-known company for assisted navigation. Due to privacy protection regulations, we do not have access to direct GPS traces. Instead, this information Fig. 6: Speed data in the GTL-Ville. Top: Colored lines correspond to the roads with speed data in real time. Bottom: Time plot of a selected road for one day. is aggregated and only the average speed of all reporting vehicles in a road section during a time period is available, for a subset of roads in the network. This data is obtained in real-time with a frequency of 1 min.

Figures 6 shows the available speed in the GTL-Ville platform. In the map representation, the colored lines show the speed data for each road in real time (green are high speed values and red are low speed values). The bottom of the figure shows a time plot of the speed of one road for one day.

V. RESULTS

For validation purposes, we consider the flow and speed data collected during February 8th, 2021. First, we use Algorithm 1 to calculate the paths between each OD pair, using K min = 2, K max = 10, and ϵ = 1.2. This resulted in a total of 3437 paths.

The OD flow vector ψ was computed using the OD rates as described in Section IV-C, in conjunction with the boundary flows from Section IV-A, averaged during the entire day. As there are more roads at the boundaries than OD rate measurement sites, each road was mapped to the closest OD measure.

We consider three methods for calculating the TR values: a) Method 1: Using regular TA: With the set of paths P and the OD flow vector ψ as input, ( 9) is solved to calculate the path flow vector f , and then [START_REF] Sheffi | Urban Transportation Networks: Equilibrium Analysis with Mathematical Programming Methods[END_REF] to estimate the TRs. For this, we used the well-known open-source convex optimization solver CVX in Matlab [START_REF] Grant | Graph implementations for nonsmooth convex programs[END_REF], [START_REF]CVX: Matlab Software for Disciplined Convex Programming[END_REF].

b) Method 2: Using TA with TR measures: This follows an analogous procedure as the Method 1 but uses the modified cost function [START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF] with the TR data described in Section IV-B. The heuristic value γ = 1/|T * | = 1/142 was used.

c) Method 3: Using FRC: For comparison purposes, we also used the TR values produced by the node-based approach we proposed in [START_REF] Rodriguez-Vega | Dynamic density and flow reconstruction in large-scale urban networks using heterogeneous data sources[END_REF], which makes use of the FRC. In this case, we consider 7 different classes, going from 1 (most important, e.g. Highways) to 7 (least important, e.g. local residential roads). To each class c, a weight θ c is given, and each TR r FRC i,j is computed as

r FRC i,j = θ cj k|(i,k)∈T θ c k , (21)
where c i is the FRC of road i, and the superindex is to specify the TR estimation method. The weights θ are computed through a flow-matching optimization problem using the measured boundary flows and imposing the measured TRs as constraints (see [START_REF] Rodriguez-Vega | Dynamic density and flow reconstruction in large-scale urban networks using heterogeneous data sources[END_REF] for details). This method assumes that vehicles tend to turn to the most important road in each intersection, independent from destination.

A. Turning ratio estimation The measured TR values were compared with the estimates provided by the TA based methods. The resulting errors followed a Gaussian distribution with parameters shown in Table I. The route-based TR estimates are seen to be close to the measured values. As expected, the inclusion of TR data in the TA problem decreases the standard deviation of the errors by 4.4%. The errors can be further decreased by increasing the value of γ, although too high values could produce unrealistic results by artificially increasing the total travel time of routes chosen by vehicles.

B. Flow and density estimation

The TR estimates using the three methods previously described were used to estimate the dynamic flow and density in each road of the network using ( 18) and [START_REF]CVX: Matlab Software for Disciplined Convex Programming[END_REF]. For validation, we used the flow measurements for the roads equipped with validation sensors. Figure 7 shows the measured flow values as ground truth, and the dynamic flow estimates using the TR values estimated using the three methods.

In general, the three approaches provide close estimates of the ground truth flow values. Note that for sensors R15 EW and R16 NS, the use of the TA approach (both with and without TR measures) significantly improves the estimation accuracy, especially when TR measures are included. As error metrics, we use the Relative Mean Error (RME)

RME i = T 0 (φ out i (t) -φout i (t))dt T 0 φ out i (t)dt, (22)
and the Relative Absolute Error (RAE) where φ out i (t) and φout i (t) are the real and estimated flows in road i, respectively. To take into account the data from all cross-validation sensors, we use the Cumulative Distribution Function (CDF) for each metric: y = CDF(x) is the fraction of roads y which have an error value less than or equal to x. Figure 8 shows the resulting CDF for both error metrics for all the validation sensor locations.

RAE i = T 0 φ out i (t) -φout i (t) dt T 0 φ out i (t)dt (23)
In general, the TR values obtained using route-based TA approaches outperform the node-based approach, and the inclusion of TR measures in the TA approach improves the accuracy even more: 80% of the roads have an RME below 38%, 29%, and 21%, for Methods 3, 1 and 2, respectively. Comparing the TR estimates from Methods 1 and 3, both have very similar behavior up to an RME of 28% and RAE of 35%. From this point on, the estimates from Method 3 start producing higher errors, except for one location, where Method 1 produces a high jump in RME from 33% to 72%. However, this outlier is not present in the results of Method 2. This is because the high-error location in Method 1 corresponds to a sensor in the interior of the network R16 NS, whose real behavior is not correctly captured by only minimizing traveling time. By inclusion of the local TR information in the optimization problem, more detailed

VI. CONCLUSIONS

In this paper, we have explored the use of route-based traffic assignment methods under user equilibrium, to estimate turning ratios at intersections. These estimates are then used in a computationally efficient and scalable dynamic flow and density estimator. We also consider a slight modification of the Assignment problem such that direct measurements of a set of Turning Ratios can be used when available. To validate these methods, we use real traffic data collected in the city of Grenoble and then compared it against a local node-based approach to estimate Turning Ratios, which we proposed in a previous work.

The results show that in general, the route-based approaches outperform the accuracy of the flow estimates compared to the considered node-based method. Furthermore, the case where Turning Ratio measurements are added to the classical methods provided the best reconstruction, greatly surpassing the other methods. The price to pay for the increased accuracy of route-based methods is the increase in computational cost, as the listing of paths between all Origin-Destination pairs and the TA optimization problems can incur a lot of computations. However, these computations can be performed offline with low-frequency average inputs.

4 Fig. 1 :

 41 Fig. 1: Flow exchange at an intersection.

Fig. 2 :

 2 Fig. 2: Example network. Boundary roads are E in = {1, 4} and E out = {8, 9, 10}.

Fig. 4 :

 4 Fig. 4: Location of flow and turning ratio sensors. Flow sensors in color with text corresponding to the identifier in the database. TR sensor in gray with labeling numbers from 1 to 12.

Fig. 5 :

 5 Fig. 5: Location of sensors for measurement of oD matrices, and the total measured OD ratios, normalized with respect to each origin. The two lines connecting each OD pair follow the right-hand rule to specify direction. enter the intersection through road i and left via road j, during the total duration of the campaign.

Fig. 7 :

 7 Fig. 7: Time series of the real and estimated flows using different TR values.

Fig. 8 :

 8 Fig. 8: Relative error CDF for flow estimation using different methods for TR computation. user route selection can be recovered increasing the accuracy, as seen in the Method 2 results.

TABLE I :

 I Error distribution of TRs using route-based approaches

	Method	Mean Std. Deviation
	Regular TA	0.7%	24.9%
	TA with TR measures	0.7%	23.8%

https://www.tomtom.com/

ACKNOWLEDGEMENTS

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement 694209), see [17].


Richard Moussa

Marie Beurton-Aimar

Guilhelm Savin

Pascal Desbarats

Image segmentation using social agents

Keywords: Image segmentation, Social spiders, Multi-agent system, Articial life. 1 Magnetic Resonance Imaging 1

In the literature, there are many methods for image segmentation. Many methods are often limited to a few number of images or by their sensibility to image artefacts like noise. Thus, there is no generic method for solving problems of image segmentation. In this report, we will introduce a new way to improve segmentation methods which use a model coming from biology: the social spiders and implemented by a multi-agent system. This method developed will be presented with its shortcomings and the way to overcome them. Finally, some perspectives will be shown to improve the method and its results on MRI 1 images.

Introduction

Image segmentation is one of the most difficult task in digital image processing. This is an important step for quantitative analysis of brain images and to study many brain disorders. Indeed, structural changes in the brain can be the result of brain disorders. The quantification of these changes, by measuring area of regions of interest, can be used to characterize disease severity or evolution. Manual marking out of cerebral structures in MRI 2 images by an expert is obviously a time consuming process. Moreover, these manual segmentations are prone to large intra/inter-observer variability. Thus, the segmentation remains a challenging task. Image segmentation algorithms subdivide images into their constituent regions, with the level of subdivision depending on the problem to be solved. Robust, automatic image segmentation requires the incorporation and efficient utilization of global contextual knowledge. However, the variability of the background, the versatile properties of the target partitions that characterize themselves and the presence of noise make it difficult to accomplish this task. Considering the complexity, we often apply different methods in the segmentation process according to the nature of the images.

Social spiders belong to spider species whose individuals form relatively longlasting aggregations. Whereas most spiders are solitary and even aggressive toward conspecifics, hundreds of species show a tendency to live in groups and to develop collaborations between each other, often referred to as colonies. For example, these spiders of 5mm in body length are capable to fix silks up to a volume of 100m 3 [6]. This technique is used to trap preys having big forms.

This work shows how it is possible to use the social spiders method to implement an efficiency image segmentation method to treat brain MRI images.

This report is organized as follows: the first section gives a view on image segmentation methods. The second section describes different types of multiagent systems, their utilities for the image segmentation. The third section presents some related works with this subject. The following section presents the method used for segmenting an image before we provide an experimental assessment of the method. Finally, we will conclude and propose some ideas for further works.

Image Segmentation

Image segmentation consists on partitionning an image into a set of regions that covers it. After this process, each pixel/voxel is affected to a region and each region corresponds to a part of the image. The discontinuity between the regions constructs the contour of the object. The segmentation approaches can be divided into three major classes [START_REF] Jähne | Digital Image Processing[END_REF]. The first corresponds to pixelbased methods which only use the gray values of the individual pixels. The second is region-based methods which analyze the gray values in larger areas for detecting regions having homogenous caracteristics, criteria or similitude.

Finally, the edge-based methods detect edges and then try to follow them, this can be done by computing a luminacy function. The common limitation of all these approaches is that they are based only on local information. Even they can only use a part of the information. Pixel-based techniques do not consider the local neighborhood. Edge-based techniques look only for discontinuities, while region-based techniques only analyze homogeneous regions.

Multi-agent system

As we have previously explained, a multi-agent system is a distributed system composed of a group of agents, which interact between them through an environnement. Agents are classified in two categories: cognitive and reactive. Cognitive agents have a global view of the environment, they know the task for which they work. Conversely, reactive agents only know a restricted environment part. They react to environment stimulus and can modify this environment by adding or removing informations. Reactive agents do not know the complex task for which they work: they have a restricted set of simple features and they only apply them. Insects colonies like spiders and ants are an example of reactive agents: each one knows locally what it has to do, but no one knows the more complex task for which they work. Such insects are called social insects.

Domain: Artificial Life Vs Artificial Intelligence

The most part of multi-agent systems are considererd as relevant to Artificial Life Domain. To clarify the field of Artificial Life (AL), we will distinguish it from the Artificial Intelligence (AI): AI takes human cognition as a reference and is interested in designing machines that can simulate human cognition. AL tends to be closer to biological functions and processes over behavior "reflexes" than logical reasoning or cognitive acts. Its field of study is broader, it explores the characteristics of living in general. AL operates the concepts from the current cognitivity and connectivity that the cognitive science share. Therefore, it proposes solutions based sometimes on the symbolic approach, sometimes on the self-organizational approach, or built its models from a combination of both methodologies. One can distinguish two types of work from the AL [START_REF] Bedau | Artificial life: Organization, adaptation and complexity from the bottom up[END_REF]. On the one hand, simulations using only the computer. It is recognized, in this case, that any system can be modeled by computer. That is to say a formal system is adequate to represent a physical system in a likely way. On the other hand, achievements, where the concrete and hardware system is paramount. It is considered here that the physical dimension of a system is irreducible to a symbolic representation.

Related Works

Ramos and Almeida [START_REF] Ramos | Artificial ant colonies in digital image habitats -a mass behaviour effect study on pattern recognition[END_REF] have explored the idea of using a digital image as an environnement for artificial ant colonies. They observed that artificial ant colonies could react and adapt appropriately their behavior to any type of digital habitat. Ramos et al. [START_REF] Ramos | Self-organized data and image retrieval as a consequence of inter-dynamic synergistic relationships in artificial ant colonies[END_REF] also investigated ant colonies based data clustering and developped an ant colony clustering algorithm which he applied to a digital image retrieval problem. By doing so, they were able to perform retrieval and classification succesfully on images of marble samples. Liu and Tang [START_REF] Liu | Adapive image segmentation with distributed behaviorbased agents[END_REF] have conducted similar works and have presented an algorithm for greyscale image segmentation using behavior-based agents that self reproduce in areas of interest. He and Chen [START_REF] He | Artificial life for image segmentation[END_REF] have provided an artificial cell model. In their model, each life is one individual unit, called a cell, which adheres to one pixel in the image. All the cells bear similar structures but mutations may occur during the process of reproduction due to the influence of the environment. Hamarneh et al. [START_REF] Hamarneh | Deformable organisms for automatic medical image analysis[END_REF] have shown how an intelligent corpus callosum agent, which takes the form of a worm, can deal with noise, incomplete edges, enormus anatomical variation, and occlusion in order to segment and label the corpus callosum in 2D mid-sagittal images slices in the brain. Bourjot et al. [START_REF] Bourjot | A new swarm mechanism based on social spiders colonies: from web weaving to region detection[END_REF] have explored the idea of using social spiders as a behavior to detect the regions of the image. The principle is to weave a web over the image by fixing silks between pixels. This method has been selected, implemented (see section 4) and experimented (see section 5) in comparaison with other type of methods.

Segmentation Method

The multi-agent system is composed of an environment and a set of agents. For segmentation purpose, environment is created from a given grayscale picture: it is a matrix of gray pixels. System and agents have a life cycle. Figure 1 shows that a cycle of the system consists in executing the life cycle of each agent. This life cycle is transposed to a step. The number of steps to be executed is given by the user. Algorithm 1 presents a description of figure 1 where for a given number of iterations (steps), each spider computes the three functions. Its complexity is about O(Ite.N ), where Ite is the number of iterations and N the number of spiders. Movement(s, Pixels, scuts)

6:
Skillfixing(s, scuts)

7:
Comeback(s)

8:
end for 9: end while

Environnement

Environment is a matrix of gray pixels. Each pixel is a position for spiders and allows them to access to other pixels. Neighborhood of a pixel p is defined in two ways (figure 2):

• The 8 pixels around p, which are named Neigh.

• All pixel linked to p by a dragline, which are named SCuts.

Set of all reachable pixels are named Access.

Figure 2: Neighborhood of a pixel

Agent

As previously mentionned, spiders are reactive agents. They are defined by an internal state composed of a set of parameters values, a current position and the last pixel where a spider has silked. These spiders have also an ability to move in the environnement, to fix silk3 and to come back 4 . Spiders which detect the same region can be grouped in a set called a colony. Spiders of a same colony share the same set of parameters values.

An agent life cycle, as shown is figure 1, consists in firing each behavioral item according to a probability which depends of parameters values and environmental characteristics at the spider position.

Movement

Movement is computed in two ways: the way where only one colony of spiders is set up in the environment (single region detection) and the way where several colonies try to detect a region (multiple-region detection). In both cases, a spider reaches one neighbor of its current pixel like described in figure 2.

In the first case, spider must choose between move to a pixel in SCuts or move to a pixel in Neigh. This choice is done according to the value of simulation parameter pdraglines. Then, each pixel of the choosen set has the following probability to be choosen:

• if set is Neigh, each pixel has the same probability • if set is SCuts, each pixel p has a probability draglinesp draglines all where draglines p is the number of draglines going to p and draglines all is the number of draglines starting from the current position of the agent.

Probability of a pixel which is in both set is the sum of both probabilities.

In the second way, a weight-function is used to compute the probability to move to a pixel. So, for each pixel p, probability to move to p is

P (M ove(p)) = w(p) a∈Access w(a) (1)
w(p) is constant if p is in Neigh, else if p is in SCuts we distinguish draglines woven by agents colony and those woven by other colonies. Two simulation parameters, attractself and attractother, are used to compute this probability. These parameters are respectively the attraction for draglines woven by agents colony and draglines woven by other colonies. A function F is also used to express draglines-counts influence on weight. So for a pixel p in SCuts, we have the following weight :

self = attractself . F (draglines self ) other = attractother . F (draglines other ) (2)
w(p) = self + other
Algorithm 2 describes the movement function by affecting to the weights of the neighborhood pixels whether the weight of the colony of the spider whether its weight function and finally move the spider to a pixel from the neighborhood according to its weight value. Its complexity is O(N ei) where N ei is the maximun of neighborhoods of a pixel.

Algorithm 2 Movement

Require: S: Spider, Pixels: Matrix of pixels ∈ N 2 , Scuts: Matrix of draglines ∈ N 2 1: weights[neightSize(position(S)) + scutsSize(position(S))]: weights of the neighborhood pixels 2: for i ∈ {0, . . ., neightSize(position(s)) -1} do

Silkfixing

Here, two other simulations parameters are used, reflevel and selectivity. The first one corresponds to the graylevel of the region to detect, and the second one defines the tolerance to fix a dragline with a pixel whose graylevel is not exactly reflevel. Probability to fix a dragline with current pixel (figure 3) follows a gaussian distribution whose mean is reflevel and standard deviation is

Comeback

This behavorial item allows spiders to come back to the last silked pixel. This action is fired depending on the probability value defined as backprobability.

The aim of this action is to detect connected region of pixels: spiders can not go far of pixels they have silked. As it is shown by [START_REF] Bourjot | A new swarm mechanism based on social spiders colonies: from web weaving to region detection[END_REF], disable this item leads to create draglines between two unconnected groups of pixels. Algorithm 4 compares a random number with a probability function and tries to return to last fixed dragline if the condition is satisfied in ligne 2. Its probability is also constant.

Algorithm 4 Comeback

Require: S : Spider 1: p ← random(0, 1) 2: if p < backprobability(colony(S)) then

3:
postion(S) ← lastfixed(S) 4: end if

Dynamic system

Dynamic system is based on stigmergy: each spider lets informations on the environment, those are used by other agents or by itself in a next cycle. Image segmentation emerges from the global task achieved by all the spiders: after a certain number of system iterations, draglines are created. Degree of a pixel defines number of draglines coming in or out of this pixel. In multiple-regions detection, degree is given according to a colony. In this case, global degree is the sum of degree for each colony. Region detected by a colony is composed of pixels having higher degree for this colony.

Problems

The social spiders method raises two important problems :

1. There is a high number of parameters

Computing the number of iterations required or defining a stop condition

The parameters computation can be problematic if it is done empirically. Indeed, to compare the results of social's spiders segmentation among several images, we must be sure that the computation of the parameters will be in an equivalent manner in all cases.

Similarly, it is important that the stop condition meets the same criteria between different acts of segmentation. Otherwise, the results could be evaluated: a number of iterations not big enough lead to a poor qualitative analysis, and a number of iterations too high could increase the simulation time of the method without improving the quality of the results.

Solutions

In this section, we propose solutions to the two problems evaluated above.

Autodetection of parameters

The method of socials spiders has several parameters that will define the regions to be detected. These parameters are presented in Table 1. First, it must know the number of colonies (the number of type of regions) and each of the properties, in particularly the grayscale reference and the selectivity parameter. We have seen that the computation of these empirical parameters could be an obstacle for comparing the results producted by the socials spiders method: if the parameters are not computed in the same way in each case, the comparisons became unreliable because they depends on the computation method. We'll see in this part a method for determining the parameters of the simulation. This method will be used thereafter in all segmentations processes.

Each colony has a grayscale (or intensity) of reference which will serve the spiders of the colony to determine if they should fix a silk. There is two possibilities:

1. Determining "manually" the intensities, which is unattractive except where the intensities are known by the user;

2. Using a method to determine automatically the intensities.

The automatic detection will allow us to obtain optimal parameters without the user knowledge about the image. Thereafter, we will appoint an intensity level of pixels/voxels of this intensity. We will use the histogram of the image to determine the intensity of interest to be detected. Indeed, a region in a image usually implies a peak more or less important in the histogram of the image. Therefore, we will detect successively each peak by eliminating the peak intensities component of future options. The Figure 4 shows the detection of maximum (peak heights) of the histogram: A: detecting a maximum intensity level n from the histogram. B: detecting the neighborhood intensities v such as: Indeed, the maximas caused by noise are on the peaks of the histogram whose slopes are small and will be erased by smoothing them in few iterations. Finally, the maxima representing a region is located on a peak having an important slope which requires a large number of smoothing in order to be eliminated. However, the risk to clear up an interesting maxima is not zero. In order no to fall into aberrations, we add the condition that the number of maximas, and therefore the number of regions to detect, must be greater than 2. Indeed, the image contains at least a substance and a foreground.

We get a serie of maximas M 1 , M 2 , . . ., M k . The number of maximas, k, determines the number of spiders colonies that will be used. Each parameter is the maximum reflevel of the colony. Then, we will partition the intensities in as many classes as maximas detected. To do this, simply find the minimum level of intensity between two maximas. This yields a serie of minimas m 1 , m 2 , . . ., M k-1 so that 0 < M 1 < m 1 < . . . < m k-1 < M k ¡ l where l is the maximum intensity. The k classes intensities obtained are [0; m 1 ],]m 1 ; m 2 ], . . ., ]m k-1 ; l].

The variance of each class provides the parameter selectivity of each column. The number of spiders per colony may be adjusted according to the sum of degrees of class.

We can now determine the essentials parmaters of the socials spiders method. It remains to determine the numbers of spiders per region and the attraction of the silks on spiders.

Stop condition

We have seen that the socials spiders method has a life cycle which is repeated a number of times until the image is segmented. The number of iterations will influence two important points on the result of segmentation:

1. The quality of segmentation; 2. The execution time required to achieve this result.

At each iteration, the spiders will weave between the pixels that will be used to determine to which region is the pixel. If the number of iteration is not enough, the number of pixels that do not belong to any region will be important and the result will be of poor quality. On the contary, if the number of iterations is too large, the spiders will only increase after a certain time the silks already existing without providing any new information. This last point will have as effect a longer execution time for an identical quality result.

As for the parameters, the computation of the number of iterations is an important point to obtain comparative results that are credible. Rather than fixing a number of iterations, it is possible to determine a stop condition to be verified before an iteration.

Definition Let β be the number of silks fixed between pixels whose degree is zero during an iteration.

The result of socials spiders segmentation depends on the silk which will be fixed between the pixels. During an iteration, when β tends towards zero, we can consider that the system stabilizes and the spiders only reinforce the existing silks.

Figure 6 shows the evolution of β along the segmentation process. It is possible to detect when β remains at zero. This moment determines to stop the simulation. We can improve this condition by adding two parameters:

• A threshold that determines when β may be considered invalid;

• The number of authorized β zero before stopping the simulation.

The threshold may be determined by the number of spiders. Indeed, at each iteration, each spider has the possibility to fix a silk. The number of silks fixed during an iteration is bounded by the number of spiders.

Results

In this section, we will present the results obtained by our simulation on MRI images. First, we will present the image used to improve our results. Then, we will describe the methods used for comparaison. Finally, we will show the results of the simulations and argue them.

BrainWeb

BrainWeb5 is a simulator of MRI brains images. It helps to validate technical medical analysis as segmentation. Indeed, the image is obtained from a model which we know its different components. A slice of this model is presented figure 7.

Figure 7: Slice of the 3D model We can distinguish ten items in this model are listed in Table 2 The parameters of the model to be used in this study are:

• Modality: T1

• Slice thickness: 1mm

• Noise: 3%

• Intensity non-uniformity: 20%

• Dimensions: 181 x 217 x 181vx

First, we will segment a 2d slice of this image (Figure 8) and then attempt the same process on a 3D image. This slice choosen to our tests is the 90th slice of the image. This image is used to determine if the social spiders method is able to segment images having complex forms. The result obtained will allow us to decide if it is feasible to segment images in three dimensions.

Comparison methods

In this section, we will compare the social spiders method with other segmentation methods. These comparisons will allow us to determine whether the social spiders method brings something positive compared to traditional segmentation methods.

We will use these comparisons on two other methods:

• A classification method by thresholding: the otsu method;

• A region-based method: the region growing method.

To compare these methods, we need to establish criteria to be used on all test images. We will compare the results on several points:

1. The number of regions;

2. Correspondence between the regions of the model and the segmentation result;

3. Execution time.

Definition Let Γ be an image and ∆ its segmentation result. We call Γ i the region i of the image and ∆ j the region j of the result.

The number of regions will allow us to determine whether the method considered detects a number of regions close to reality. In the case of noisy images, it is possible that some methods detect regions with insignificant size. That is why we will add to the total number of regions, the number of regions having insignificant size.

Definition A region will be considered insignificant if its size is less than 10 pixels.

The computation of the number of regions will be done on the segmentation method result on which a labeling is added to the connected components to consider the regions connected.

The correspondence between initial image and its result will enable us to determine if the regions identified by the method correspond to the regions defined in the initial image. This is only possible in the case of images for which we know lots of informations. This match will be appointed accuracy.

To compute the accuracy, it is necessary to determine which region ∆ j matches the most the region Γ i . This region is determined by:

n i =
the total number of pixels the voxels of Γ i n j = the total number of pixels the voxels of ∆ j

n i j = (3)
the total number of pixels the common voxels between ∆ j and Γ i Thus, it is possible to compute δ i,j = n i j n i and γ i,j = n i j nj representing respectivly the proportion of pixels / voxels of Γ i belonging to ∆ j and the proportion of pixels / voxels of ∆ j belonging to Γ i . We have two ways to choose the region that corresponds to ∆ j corresponding the most to Γ i :

1. ∆ k as the value of δ i,k is maximum: in this case, we prefer the fact that Γ i and ∆ j have a maximum of pixels / voxels in common;

2. ∆ k as δ i,k + γi, k is maximum: same as above, but we add the requirement that ∆ k must have a minimum of pixels / voxels in other regions than Γ i .

In our results, we indicate two points, accuracy δ and accuracy δ+γ , which corresponds respectively to the two choices of ∆ k described above. In both cases, the accuracy will be the average values for all regions of the model.

Region growing

The Region growing method as described by Shapiro and Stockman [START_REF] Shapiro | Computer Vision[END_REF] consist on building a region from one chosen pixel and then adding recursively neighbors whose grayscale difference with the original pixel is below to a threshold. This method tries to grow an initial region by adding to this region pixels that do not belong to any region but to the pixels neighboorhood already in the region and whose grayscale is sufficiently close to the area. When it is not possible to add pixels, we create a new region with a pixel that has not been selected yet, then we grow the region.

The method ends when all the pixels were chosen by a region. To automatize the process of this algorithm, the threshold was choosen as the local minima of the image histogram.

Otsu

Otsu [START_REF] Otsu | A threshold selection method from gray level histograms[END_REF] has developped a multi-level thresholding method. Its aim is to determine, for a given number of regions, the optimum values of different thresholds based on the variance of subdivisions created.

The basic method consists on separating the foreground from the backgroud. In this case, we search the optimal threshold to split the pixels in two regions. For a threshold t, it is possible to compute the between-class variance σ 2 (t). This measure is derived from the average intensity µ 1 , µ 2 and µ of classes [0; t], [t + 1; L] and [0; L] where L is the maximum intensity.

The Equation 4shows us the computation of σ 2 , where w 1 and w 2 represent the proportion of pixels in the class [0; t] and [t + 1; L] compared to the total number of pixels.

σ 2 (t) = w 1 (t)(µ 1 (t) -µ) 2 + w 2 (t)(µ 2 (t) -µ) 2 (4)
The Otsu method shows that the optimal threshold t * is obtained for a betweenclass variance. The method consists on computing the variance for all possible thresholds (t ∈ {1; . . .; L -1}) and determining its maximal value.

This method could be extended easily to the computation of M classes with M -1 thresholds {t 1 ; t 2 ; . . .; T M -1 -1} (t 1 < t 2 < . . . < t M -1 ). The betweenclass variance is defined then as follows:

σ 2 (t 1 , . . . , t M -1 ) = M M -1 w k (µ k -µ) 2 (5)
where w k represent the proportion of pixels in the class [t k-1 ; t k ]6 , µ k the intensity average of this same classe and µ the intensity average of the class [0; L].

For each M-1-uplet, we compute thresholds of the between-class variance. The optimal thresholds, (t * , . . ., t * M -1 ), correspond to the maximum value of the between-class variance.

Chen et al. [START_REF] Chen | A fast algorithm for multilevel thresholding[END_REF] propose an algorithm that minimizes the number of necessary computation to obtain a faster algorithm. This method had been implemented for the evaluations tests.

Experimentations

This section will explore if the socials spiders method supports complex images such as MRI images. For the 2D case, the results are presented in the table 3 and the figure 9.

The execution time to be given comes from the simulation of the methods on a machine equipped with two Intel (R) Xeon (R) E5345 (8 cores having 2.33GHz) and 8GB of RAM. The operating system of this machine is a Linux kernel 2.6.21 x86 64. The Socials Spiders method produced the highest number of significant regions, that closest to the actual number of regions. However, it is also the The accuracy displayed by the spiders method is quite correct compared to the other methods. However, a number of regions rather high can be explained by a number of voxels non-detected more important, leading to disconnection of the regions.

Unfortunately, this accuracy has a very important consuming time over the execution time. Therefore, as well as for the results on two dimensions images, we get a bad time result. The Otsu method, although with a precision of a few percent lower, produces a good result for a shorter time computation of about 4000 times.

In addition, the number of regions produced by the spiders methods is greater by a factor of 10 compared to the other methods.

Discussion

In this report, we have presented a new method of segmentation that showed various problems of implementation. Indeed, this method requires a large number of parameters depending on the image to be processed. In addition, the method did not really have a stop condition which impose the user to fix in advance the number of iterations required.

First, we solved the problem of computation of the parameters by proposing a method to compute automatically these parameters. This method can determine the most important parameters of the method as follows:

• The number of colonies;

• The grayscale reference of each colony;

• The selectivity parameter of each colony.

The other parameters were set as they have less importance in the segmentation result.

Then, we also saw a method based on the number of silks fixed at an iteration to determine a stop condition to the method. This stop condition has allowed us to eliminate the problem of the number of iterations needed to produce a good segmentation.

We have made comparaisons between the results of the social spiders method, the region growing method and the otsu method. These comparisons focused on the accuracy of methods, the number of regions produced and the time processing of the methods. They are not exhaustive comparisons in the sense that all aspects of segmentation are not taken into account.

Through these comparisons, we have put forward some drawbacks on the socials spiders method. Particularly, we have seen that this method produced a significant number of areas and that the execution time was particularly long.

However, the socials spiders method is based on an architecture that is ideal to be parallel. This method is composed of a group of agents that can be spread over several processors. The fact that they all run the same algorithm suggests using a SIMD7 architecture and therefore using this method on graphical processors (GPU8 ).

We have seen that the second drawback of the socials spiders method is producing an oversegmentation of the image. It is possible to solve this problem by initiating the algorithm by merging insignificant regions to ratach them to significant regions or merging them as appropriate. Merging operation is not an expensive treatment time processor, therefore, adding a post-treatment to the social spiders method adds a negligible computation time compared to the method itself.

Some solutions have been considered to improve the method. The spiders seem sensitive to the topology of the image, so it is possible to guide the movement of spiders with a gradient or a laplacian. Indeed, these measures will provide informations on the possible presence of contours. It would be then possible to use spiders in two ways:

1. Gradient would be repellent which would partition a colony of spiders in a region;

2. On the contrary, the gradient could have an attractive effect. In this case, spiders would be used to detect the contours of regions.

Conclusion

Finally, we saw a first approach to image segmentation in three dimensions using the socials spiders method. The method has been used without other modification than the extension of the neighborhood. Therefore, there is no optimization related to adding an extra dimension. However, the BrainWeb segmentation model showed that the method was capable to produce a result with good accuracy.

Optimizations are needed to reduce the number of regions produced and specially to reduce the time process of the method. In the previous section, solutions have been proposed to resolve or at least reduce these defects.

Figure 1 :

 1 Figure 1: System overview

3 :

 3 weights[i] ← w(p)(colony(S)) 4: end for 5: for i ∈ {neightSize(position(S)), . . ., Size(weights) -1} do 6: weights[i] ← WeightFunction(S) 7: end for 8: Choose a pixel in the neighborhood according to its weight and move the spider to the pixel choosen

1 selectivity

 1 . Algorithm 3 consists on choosing a random number and evaluating it with a gaussian function to add or not a dragline. Its complexity is constant.

Figure 3 : 3 :

 33 Figure 3: Fix a dragline with current pixel

  probability to weave a dragline from the current position of the insect attractself attracts toward the dralines of the same colony attractother attracts toward the dralines of the other colony backprobability probability to turn back w(p) weight of a pixel with a direct neighborhood saturationvalue maximun weight of a pixel with an indirect neighborhood

Figure 5 :

 5 Figure 5: Histogram smoothing

Figure 6 :

 6 Figure 6: Number of silks fixed during an iteration

Figure 8 :

 8 Figure 8: 2D slice: 181 X 217 px

Figure 9 :

 9 Figure 9: 2D segmentation of a brain slice: a) Social spiders, b) Region growing, c) Otsu

Table 1 :

 1 Colony parameters.

  .

	Element	Pourcentage	Form	Element	Pourcentage	Form
	background	42.2%		muscle skin	8.7%
	csf	5.2%		skin	10.2%
	grey matter	12.7%		skull	5.1%
	white matter	9.5%		glial matter	0.1%
	fat	2.1%		connective	4.2%

Table 2 :

 2 BrainWeb Model: the various components.

Table 3 :

 3 2D resultsmethod that produced the greatest number of insignificant regions, meaning that the spiders have oversegmented the image. This oversegmentation of the image has not been distented by the usage of a complex BrainWeb model. The execution time of the Socials Spiders method of spiders is too high compared to the two other methods. The number of silks fixed is less important because the image is composed of a greater number of regions, each region thus representing a smaller proportion of the image. Spiders make more time (number of iterations) to find pixels corresponding to their colony. The table4shows the results of different segmentations of the image BrainWeb in three dimensions.

		Region	Region > 10px	Accuracyσ	Accuracy σ+γ	Time
	Social spiders method	4361	434	60.8%	54.9%	50 s
	Region growing method	908	664	46.2%	21.1%	0.5 s
	Otsu method	1875	567	90.7%	89.3%	0.2 s
		Region	Region > 10vx	Accuracyσ	Accuracy σ+γ	Time
	Social spiders method	45575	13391	87.3%	87.3%	8 h
	Region growing method	7866	1634	75.4%	46.7%	102 s
	Otsu method	8497	865	83.9%	83.9%	10 s

Table 4 :

 4 3D results

Magnetic Resonance Imaging

Weave a dragline between two pixels

Return to the last fixed pixel

http://www.bic.mni.mcgill.ca/brainweb/

t 0 = 0 and t M = L

Single Instruction Multiple Data

Graphical Processing Unit

This method allows us to obtain a list of maximas in the histogram. However, this method is particularly sensitive to noise at Step B. Therefore, we will introduce a method to reduce the noise effects on the histogram by smoothing it and allowing to obtain a maximum number of representative of the number of classes intensities of the image.

Histogram smoothing

Histogram smoothing will reduce the peaks caused by noise (see figure 5). The method proposed here is done by iteration. At each iteration, the degree of intensity n becomes the average degrees of intensity n -1, n, n + 1. We propose a method that detect automatically the optimal settings for the segmentation of the image. The method determines the number of settlements and the parameters reflevel, selectivity and the number of spiders in each colony.

First, We will determine the maximas in the histogram as described above. Then, the histogrtam will be smoothed until fewer maximas are significant.