Middle East and North Africa Countries' Agricultural Export Potentials under Trade Reforms

Keywords: Integration &lt, F1 -Trade &lt, F -International TRADE REFORM, MENA COUNTRIES, MEDITERRANEAN COUNTRIES,

Introduction

Middle East and North Africa (MENA) 1 countries, amid a globalized and changed world, could significantly contribute to tackle the world's upheavals of today -illegal immigration, economic stagnation, political unrest, public's uproar -providing that their own course to economic growth and stability is ensured. Since 1995, the EU has acknowledged this important universal role of MENA countries, by the so called Barcelona agreement and attempted to expand its trade and relations to this particular region. Strengthening the vital economic sectors of MENA countries via free trade safeguards a steady development course for the whole region. Agriculture, a major economic component and employment provider for all MENA countries, constitutes, perhaps, an important economic sector 2 . At large, MENA countries concentrate on exporting agricultural products, in particular fruits-vegetables-olive oil or cotton and wheat. Though intra-MENA trade constitutes a significant part of their transactions, "extra" trade is taking place mainly between MENA countries and EU, roughly half of MENA countries exports target the EU and a significant part of their imports originates from EU too (Achy and Sekat, 2003).

Forthcoming trade reforms, either in a global, WTO, or a regional level, EU-Med Partnership (Kniper and dell' Aquila, 2004) can substantially influence MENA countries' export orientation, and subsequently, the pace of their economic growth. EU-Med Association Agreements, of great importance for MENA countries, were Trade liberalization is strictly related to changes in employment and economic growth rates and in relative prices [START_REF] Martin | The Social Impact of Euro-Mediterranean Free Trade Areas: A First Approach with Special Reference to the Case of Morocco[END_REF]. The effects of the liberalization on MENA countries have attracted the interest of several scholars [START_REF] Augier | The welfare implications of trade liberalization between the Southern Mediterranean and the EU[END_REF]Nabli and Veganzones-Varoudakis, 2004;[START_REF] Abu-Qarn | The validity of the ELG hypothesis in the MENA region: cointegration and error correction model analysis[END_REF]Siliverstovs and Herzer, 2007), while the agricultural export performance have not considered in particular. However, in non-MENA countries, trade liberalization and the impact on agricultural sector have been introduced [START_REF] Hertel | Applied General Equilibrium Analysis of Agricultural and Resource Policies[END_REF][START_REF] Hertel | Agriculture and nonagricultural liberalization in the Millenium Round[END_REF][START_REF] Gohin | The phasing out of export subsidies: Impacts on EU agriculture[END_REF][START_REF] Keeney | GTAP-AGR: A Framework for Assessing the Implications of Multilateral Changes in Agricultural Policies[END_REF]Valenzuela et al, 2006). The premise in this work is that agricultural exports form the underpinnings to build up smoothly and gradually the rest of the economic sectors, as it can refrain huge labor force and foster regional development. Thus, a particular focus is given on agricultural exports.

The purpose of this work is to assess the likely impacts due to undergoing trade reforms, WTO or EU-Med agreements, upon the MENA countries with a particular emphasis on To understand the potential effects and to trace out probable directions of the expected changes, three distinct scenarios were defined from a broader spectrum of anticipated trade reforms. Certainly, the examined scenarios do not represent the real pace of WTO or EU-Med partnership reforms in terms of details in application and needed time span to be completed. Nevertheless, the scenarios can serve the goals and the focus of this work to provide a policy supporting base. The three scenarios to be examined are the following: first, global trade reform (scenario -1) elimination of all tariffs on agricultural imports and subsidies on agricultural exports throughout the world is assumed; second, EU and MENA trade agreement (scenario -2) removal of all trade barriers on agricultural products is achieved; third, MENA special provisions (scenario -3) agricultural exports of MENA countries face no trade barriers but MENA countries maintain their barriers to agricultural imported goods (non-reciprocity), an option that can be granted to developing countries in a WTO agreement.

In the next section the methodology and data are briefly presented, while in section three the results of the various trade reform scenarios are given. Finally, in the sections to follow, the effect of bilateral trade agreements between MENA countries and EU are explored. The paper ends up with the main conclusions and implications.

Methodology and data

Computable general equilibrium model (CGE) 4 , a widely followed approach in trade analysis, and GTAP global database (1997) were used [START_REF] Mcdougall | Global Trade, Assistance, and Protection: The GTAP 4 Data Base[END_REF]; [START_REF] Hertel | Structure of GTAP[END_REF] The (CGE) model was developed by Diao et al. (2001), and details are not provided here as can be easily found in the relative literature. The CGE model is global in the sense that all regions of the world are included, and production and consumption decisions in each region follow behavior that is consistent with economic theory. Trade flows among regions are multilateral and world prices are determined by world market clearing conditions or, in other words, excess demand for each commodity in the world is zero. The general equilibrium feature of the model means that resources can move among sectors, securing consistent changes among all sectors (Diao and Somwaru, 2001). Thus, adjustments in the livestock sector, for example, are consistent with adjustments in the feed grains sector.

The assumption that labor and capital are mobile between agriculture and nonagricultural sectors is introduced. Relaxing this assumption would slow the supply response from countries having a comparative advantage in world agricultural markets causing world agricultural prices to rise more than the predicted increase in this analysis.

Moreover, the assumption that labor is fully employed places upward pressure on prices, since abundant labor is likely in MENA countries, supply response can occur at a lower cost.

A "base" scenario was developed initially to represent a stylized view of agricultural production and trade in the world under current trade policies. First, a global trade reform scenario was simulated (scenario -1). In this scenario, all tariffs and export in welfare and gross domestic product. The analysis is based on the 1998 levels of applied agricultural tariffs, domestic support and export subsides, and the use of tariff rate quotas. Caveats need to be noted. First, tariff rates and tariff equivalent rates are based on the 1998 data. Since tariff reductions have been undertaken by many countries after 1998, and since the bound rates are much higher than the applied rates in many cases, the analysis may overestimate the extent of tariff reduction that would take effect after 2000 for some countries. In this situation, the analysis may overestimate the extent of all import barriers. Though this is a very extreme and unrealistic scenario, assumed that all tariffs on agricultural imports and subsidies on agricultural exports worldwide are eliminated, it can provide valuable information on the direction of anticipated changes in a completely free trade for agricultural commodities.

From a theoretical point of view, restricting imports in the import-protecting countries causes domestic consumers to face higher food prices than world prices and to employ more resources in agriculture. Eliminating import tariffs will induce a demand rise for agricultural imported goods and a supply drop of domestic produced products, placing upward pressure on world agricultural prices. This upward pressure in turn induces agricultural exporting countries to increase production. This is in line with the results, as the level of world agricultural prices rises by 11.6 per cent relative to the level of world non-agricultural prices, the worldwide agricultural production increases by 1.15 per cent and the trade flows in value and volume increase by 40 and 23 per cent, respectively.

The well-accepted equivalent variation (often referred to as the willingness to pay) is used to measure the social welfare gains or losses due to trade liberalization. One-time welfare effects are considered. The one-time effects are measured by using the statusquo (pre-reform) prices as the base, and address the question: what income would be equivalent to the change brought about by the trade liberalization (Varian, 1984). The Welfare effects of trade reform suggest that consumers can be worse off if the country's terms of trade deteriorate following liberalization. That is, if the prices of the goods a country exports fall relative to the prices of goods the country imports, then consumers can be adversely affected since their expenditures on imported goods increase while their income from exported goods falls.

As table 1 shows, most countries, but MENA countries, experience an increase in welfare measured by changes in GDP and equivalent variation that accounts for the welfare gains or losses due to agricultural trade liberalization. The negative effect of global trade reform on MENA welfare is mainly caused by erosion in the preferential treatment by the EU. EU has granted special trade preferences to the MENA countries.

After worldwide reform, MENA countries may experience a welfare loss because they suffer a decline in demand for agricultural goods that would have been exported to the EU. In other words, MENA countries will experience deterioration in their terms of trade as a result of global reform in the agricultural sector (namely, net importers of agricultural products). 1).

Trade flows by country/region will be substantially influenced as provided in table 2.

Australia, New Zealand, the EU, the U.S., and Japan will benefit the most from export growth. Value changes in exports will reach almost fifty per cent for the EU and the U.S.. MENA countries' exports will rise by 29.57 per cent. Changes in agricultural import values depict a very different story from that of exports as the highest increase in imports occurs in MENA countries, followed by Japan and the EU. As the paper focuses on MENA countries, table 3 presents changes in commodity trade flow for twelve agricultural commodity/aggregates.

For MENA countries, vegetables, fruits, and olive oil products are of paramount importance and changes in exports / imports in those sectors would play a significant role in the agricultural economy of the MENA region. Under scenario -1 exports and imports of vegetables and fruits, in value terms, increase by 28 and 148 per cent, respectively. Almost the same pattern is followed by olive oil products where the increase in imports counterbalances the increase in exports.

Looking at export / import changes of agricultural products under this scenario, the effects are devastating for the economy. The imports increase outweighs the exports increase for all agricultural products. Thus, an abrupt liberalization of the agricultural trade may bring the whole economy to disarray, since the agricultural sector will suffer the most. In this respect, pursuing a gradual liberalization in agricultural markets may ease the negative effects.

EU and MENA Trade Agreement (Scenario -2)

Under this scenario, the EU and MENA countries sign a bilateral regional agreement, abolishing all current trade barriers on agricultural products. This scenario tries somehow to depict the current EU-Med agreement roadmap in an abstracting way. In accordance with the Barcelona agreement, and the follow-up agreements, this is the ultimate goal that may come into effect in 2010 or afterwards.

From an economic point of view, this scenario can be classified as a regional or bilateral agreement. Two main changes in the trade are expected, the "trade creation" and the "trade diversion". Trade creation occurs if the agreement permits efficient producers in one member country to sell into a previously protected neighbouring market without affecting the exports of more efficient non-members. When trade-creation occurs, capital and other factors of production are reallocated toward more efficient uses, raising the returns to those factors and improving the overall economic welfare of members.

Countries outside a trade-creation agreement could benefit as well, if the efficiency and welfare gains in member countries generate trade and growth opportunities for nonmembers. Trade-diversion, on the other hand, causes importers to switch from more efficient suppliers outside the agreement to less efficient suppliers within the agreement, distorting the allocation of resources and harming non-members of the agreement. Under this scenario world prices of agricultural products increase by 4.2 per cent while worldwide trade in value and volume increases by 1.46 and 0.92, respectively. As expected, welfare in EU and MENA countries, would increase minimally. Trade flows between EU and MENA countries increase both in value and volume, that is trade creation, while trade diversion may result in the countries outside the agreement (table 2, second column). What is worth noting is that the rate of exports increase is higher for the MENA countries than for the EU, backing the premise that EU-Med agreements can boost their agricultural exports.

As the results of this scenario reflect to a large extent the expected changes due to completion of Barcelona roadmap, the results in table 3 shed light on the mix of agricultural trade. The main agricultural production of the MENA countries, that is fruits, vegetables and olive oil, benefit substantially from liberalizing agricultural trade.

MENA Special Provisions (Scenario -3)

The Uruguay Round Agreements contain special provisions for developing countries that can be granted to them, among other things, long and gradual phase-in periods for their commitments and fewer obligations in some sectors (The World Bank, 2003). As a consequence, and in order to glean the full benefits from trade openness or trade reforms, MENA countries can take advantage of the market access provisions and adopt adjustments designed to improve their supply response. The trade impact of the reductions in tariff levels on the exports of any one individual MENA country depends on the treatment granted to its products by the importing countries [START_REF] Michalek | Comparative Analysis of Technical Barriers to Trade (TBT) for Central and Eastern European Countries' and Mediterranean Partner Countries' Exports to the EU, FEMISE[END_REF].

In order to evaluate the impact of export potentials of MENA countries, under this scenario, it is allowed MENA countries' exports to face duty-free trade status. As expected, agricultural exports of MENA countries' are grown. Results (table 2) 12 expected, the welfare gains are minimal, though that a slight increase in the world trade was registered (Table 1). Furthermore, a glance at the results in table 3 reveals that for the most of agricultural products a significant export upsurge was recorded. Having in mind that table 3 reflects percentage changes, in can be deducted that increase in vegetables and fruits exports will induce the expansion of the whole sector. This increase will be very conducive for the whole economy of the MENA countries, as agricultural activities, in particular fruits and vegetables, employ the vast majority of the population.

Concluding remarks

Over the last years, MENA countries have attracted the focus of the EU as they present a close and growing market and they maintain huge population reserves, not to mention energy reserves too, that can serve to the EU's consumption and production engine.

Thus, EU very early (1995) engaged in negotiations with Mediterranean countries and signed the EU-Med partnership agreements, aiming at strengthening the bilateral trade and developing the whole Mediterranean area. in line with this agreement, enhancing the EU-MENA countries trade is the larger strategic and prior perspective; but the particular importance of agriculture in retaining the labor surpluses and boosting the development of close related sectors, food and processing industry must not be undermined. The agriculture trade between MENA countries and the EU is of a great importance, as more The trade liberalization process either through EU-Med agreement or WTO could cause substantial changes upon MENA countries trade, growth and welfare. With particular focus on agricultural trade, a Global CGE model was employed to trace out probable impacts of alternative trade liberalization options on MENA countries. Since trade relations through time can be influenced by several non-measurable and unforeseen factors, the assessment of trade liberalization impacts, beyond being a laborious task, is a complex process. Nevertheless, the application of this modeling process can be proven conducive to the discussion of future policy reforms on the base of sound empirical results.

The modeling proceeds first by developing a "base" scenario to represent a stylized view of agricultural production and trade in the world under the current trade policies.

Then, the modeling carries on by building up a global trade reform scenario, where all tariffs and export subsidies on agricultural commodities are eliminated, a full liberalization scenario. Same process is repeated by simulating a scenario that represents somehow the EU-Med agreements. According to this scenario, all trade barriers upon agricultural commodities between the EU and MENA countries are removed. Finally, the last scenario roughly represents WTO agreements, where preferential provisions are granted to developing countries, exports from MENA countries face no trade barriers while barriers to imports remain intact. As expected, results reveal that trade volume, welfare and production of MENA countries are substantially influenced by the liberalization scenarios. The direction of ensuing changes, and changes magnitude to a certain extent, varies significantly among the three formulated scenarios. Thus, results can signal valuable indications to EU and MENA countries policy makers on the pursued trade and integration policies. In particular, scenario -3, MENA countries were granted special provisions, proved to be the most beneficial for them, followed by scenario -2, reflecting the EU-Med agreements. Regarding the gains in terms of exports, both scenario -2 and scenario -3 could induce a substantial export growth.

Finally, export/import flows of agricultural commodities are cast by the model in each and every assumed scenario. Again, scenario -2 and scenario -3 yield the most favorable outcome for vegetables and fruits, and olive oil products, the most vital agricultural commodities for the MENA countries. Exports increase in the aforementioned commodities would foster the agricultural economy and could bring about spillover effects on the local food and processing industry. Sequentially, the drafted roadmap by the EU-Med agreements seems to be beneficial for both the EU and MENA countries and can serve as a safe pace towards further trade liberalization.

  , summit of Barcelona 3 , and a new roadmap was drafted in November 2005. According to this roadmap, a detailed evolution plan of trade liberalization will be negotiated, started in 2006, between the EU and Mediterranean Partner Countries.

  agricultural trade. Several trade reform options are studied using a global trade model. A global analysis of this type provides insights into various agricultural trade options and indicates what the potential effects on MENA countries would be, both positive and negative.

  to cast the trade reform impacts. An overview of the GTAP model is analytically described by

  imports are eliminated, while other agricultural policies remain unchanged. Second, EU and MENA countries trade agreement (scenario -2), a bilateral regional trade agreement, removing all trade barriers on agricultural products is established. In this scenario, the benefit/loss with global trade reform is contrasted.Finally, MENA countries are treated specially and differentially under the provisions in the Agreement on Agriculture for developing countries (scenario-3). In this scenario, trade policies of MENA countries remain in place while MENA countries do not make tariff concessions, given the region's economic development status. Implementation of special and differential treatment would help industrial development and export promotion efforts of MENA countries. Four indicators were used to assess the effects of agricultural liberalization on each country/region. These are: (a) changes in world agricultural prices, (b) changes in the volume of world agricultural trade, (c) changes in exports and imports, and (d) changes

  time are measured by summing the discounted value of this measure over time.

  Since trade liberalization enhances trade, growth in agriculture trade is expected worldwide. Indeed, model results indicate that world agricultural trade increases substantially after liberalization. Removal of all agricultural trade protection worldwide increase in the value of world agricultural trade by almost forty per cent (Table

  expected to gain from the trade-creation effects of a regional agreement with the EU. Factors of production would be reallocated within the MENA countries economy towards the more competitive sectors, as producers take advantage of the new export opportunities and as imports rise to challenge the less competitive sectors. The less competitive sectors of the MENA countries would decline, but gains in the competitive sectors would offset those losses.

  This change is expected to have spillover effects on the whole economy of MENA particular on the food processing industry, associated with the aforementioned commodities.

  reveal that MENA countries' trade flows increase the most. MENA countries' exports (value) would increase by 35.10 per cent and imports by 25.19 per cent. In other words, special provisions treatment induces MENA countries' trade more than in case of scenario 1 or 2. It is clear from the results of this scenario, that MENA countries benefit the most as they are granted a preferential treatment in the world trade. As it could have been

  their exports are destined to the EU and a significant part of their imports originates from EU too.

Table 1 .

 1 MENA trade liberalization options: impacts on world production, GDP and

	Notes 1 Turkey, Morocco, Algeria, Tunisia, Egypt, Lebanon, Libya, Syria, Israel, Gaza Strip
	and West Bank, and Jordan.

Table 2 .

 2 MENA trade liberalization options: impacts on trade flows by

		Scenario -1	Scenario -2	Scenario -3
		Global trade reform	Regional trade	MENA: special
			agreement	provisions
		% change from the base year (1998)
	Exports in value
	Australia/New Zealand	57.69	-0.21	2.26
	Japan	35.65	-0.09	0.70
	United States	44.78	-0.18	2.02
	EU	44.00	5.10	4.32
	MENA	29.57	16.26	35.10
	Rest of the Americas	32.08	-0.23	1.65
	Rest of the World	34.50	-0.28	1.48
	Exports in volume
	Australia/New Zealand	35.14	-0.15	1.51
	Japan	28.18	-0.06	0.47
	United States	27.32	-0.12	1.29
	EU	20.70	3.25	2.77
	MENA	20.06	10.25	22.17
	Rest of the Americas	18.97	-0.15	1.00
	Rest of the World	21.34	-0.19	0.94
	Imports in value
	Australia/New Zealand	18.69	-0.03	0.53
	Japan	58.85	0.00	0.86
	United States	16.66	0.00	0.40
	EU	26.68	2.06	2.01
	MENA	67.16	10.68	25.19
	Rest of the Americas	31.84	-0.02	0.35
	Rest of the World	45.95	0.01	1.16
	Imports in volume
	Australia/New Zealand	11.72	-0.02	0.36
	Japan	34.41	0.00	0.49
	United States	10.23	0.00	0.24
	EU	17.03	1.29	1.28
	MENA	36.20	6.69	16.05
	Rest of the Americas	19.70	-0.01	0.25
	Rest of the World	27.01	0.01	0.75

Table 3 .

 3 MENA trade liberalization options: commodity trade flows Oilseeds includes: soybeans, olives and other oil producing crops.

	Scenario -1 Global trade reform % change from the base year (1998) Scenario -2 Regional trade agreement Scenario -3 MENA: special provisions 92.96 39.07 137.97 35.71 26.67 46.38 27.89 17.64 32.28 10.17 0.27 18.87 14.89 1.95 16.81 42.81 34.01 54.65 15.88 8.65 17.48 70.40 50.85 79.79 107.55 21.16 128.84 14.24 13.24 18.66 41.40 23.59 49.04 16.18 5.12 22.18 1.05 0.24 0.37 1.23 0.15 0.36 53.29 3.18 22.38 33.22 3.53 16.74 148.43 9.89 46.25 30.80 0.37 14.43 23.04 1.57 11.00 76.18 15.55 29.80 21.40 6.14 10.95 145.49 18.04 47.17 117.34 19.61 37.29 28.30 2.09 13.03 63.54 17.44 26.57 210.95 25.58 56.55 -5.21 -0.45 -1.03 -4.72 -0.39 -1.04 *Page Exports value Food grains Feed grains Vegetables&Fruits Oilseeds* Other crops Cattle Other animals Processed meat Other processed meat Processed oil Other processed food Beverages&Tobacco Fishering Non agricultural Imports value Food grains Feed grains Vegetables&Fruits Oilseeds* Other crops Cattle Other animals Processed meat Other processed meat Processed oil Other processed food Beverages&Tobacco Fishering Non agricultural

22 of 22 Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK

	Submitted Manuscript
	F o
	r
	P
	e e r
	R e v i e w

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript

The Average Agriculture's share for MENA countries in the GDP and Exports exceeds the 10 per cent.

The so-called Barcelona Agreement is implemented through Association Agreements with each one of the Mediterranean Partner Countries, and the aim of this Agreement is the formation of a Free Trade Area after 2010.

Model is documented in[START_REF] Diao | Developing Country Interests in Agricultural Reforms Under the World Trade Organization[END_REF] and short description is provided in the Appendix.

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript

Appendix:

The model developed in this study is based on the neoclassical growth theory, and is a static CGE model with a multi-regional and multi-sectoral specification.

Consumption/savings: In each region the representative household owns land, labor and all financial wealth (defined below) to seek to maximize utility. For reasons of simplicity, we assume no independent government investment. Government spends all its tax revenues on consumption or on transfers to households and, hence, fiscal deficits are ignored. The household's utility is:

TC n , which is the aggregate consumption generated from final goods, is as follows:

where C ni is final good i in region n, and P i a ni = 1. The household in each region maximizes (1) subject to a budget constraint:

where Ptc n is consumer price index; wld n is the land rental rate, wlb n is the wage rate; TI n is the lump sum transfer of government revenues; ldt n and lbt n are household land and labor income tax, respectively.

Households allocate their total income flows, including financial and non-financial, between consumption and savings. The current budget constraint for the household is:

where SAV n is n-th region's household savings; wk n is the current capital rental price and kt n is the capital income tax rate; and Ptc nt TC nt are total consumption expenditures.

The traditional Armington functions are all specified. For consumers or investors, goods imported from abroad or produced domestically are not identical. This imperfect substitution relation is reflected with an Armingtonian constant elasticity substitution (CES) function. Furthermore, it can be assumed that the goods consumed by consumers and used for investment are different and there are different substitution elasticities for goods produced at home and imported from abroad. To simplify the analysis, we assume that the composite goods used for consumption or for investment are same goods.

Composite goods are also used as intermediate inputs in each production sector.


Marius-F Danca
email: danca@rist.ro

Qingyun Wang

Synthesizing attractors of Hindmarsh-Rose neuronal systems

Keywords: Mathematics Subject Classification (2000) 34A34, 34D45, 37B25 Hindmarsh-Rose model, chaotic attractors, local attractors, global attractors

In this paper a periodic parameter switching scheme is applied to the Hindmarsh-Rose neuronal system to synthesize certain attractors. Results show numerically, via computer graphic simulations, that the obtained synthesized attractor belongs to the class of all admissible attractors for the Hindmarsh-Rose neuronal system and matches the averaged attractor obtained with the control parameter replaced with the averaged switched parameter values. This feature allows us to imagine that living beings are able to maintain vital behavior while the control parameter switches so that their dynamical behavior is suitable for the given environment.

Introduction

Based on the theory of dynamical systems, Hindmarsh and Rose proposed the phenomenological neuron model, which may be seen either as a generalization of the Fitzhugh equations or as a simplification of the physiologically realistic model proposed by Hodgkin and Huxley [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF][START_REF] Ma | Dynamics and double Hopf bifurcations of the Rose-Hindmarsh system with time delay[END_REF]. The Hindmarsh-Rose (HR) model of neuronal activity is aimed to study the spiking-bursting behavior of the membrane potential observed in experiments of a cell in the brain of the pond snail [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF]. The dynamics of a single HR neuron has been studied well, and it is illustrated that it can exhibit complex dynamical behavior including periodic and chaotic spiking (bursting) motion when certain control parameters of nervous cell models are varied [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF][START_REF] Ma | Dynamics and double Hopf bifurcations of the Rose-Hindmarsh system with time delay[END_REF][START_REF] Holden | From simple to simple bursting oscillatory behavior via chaos in the Rose-Hindmarsh model for neuronal activity[END_REF].

It is well accepted that the HR neuron model is an alternative candidate for studying the dynamics of neuronal systems since it has simpler mathematical forms than Chay's and Hodgkin and Huxley's neuron models [START_REF] Chay | Chaos in a three-variable model of an excitable cell[END_REF][START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. Therefore, this model has been used to study different aspects of neuronal dynamics such as transitions between different firing regimes [START_REF] Gonzalez-Miranda | Observation of a continuous interior crisis in the Hindmarsh-Rose neuron model[END_REF], relations between block structured dynamics and neuronal coding [START_REF] Gonzalez-Miranda | Block structured dynamics and neuronal coding[END_REF], the effect of noise on neuronal signal transduction [START_REF] Reinker | Stochastic and suprathreshold resonances in a stochastic Hindmarsh-Rose model of thalamic neurons[END_REF] and on synchronization [START_REF] He | Stone, Noise-induced synchronization in realistic models[END_REF], the collective dynamics of neuronal networks [START_REF] Zheng | Spatiotemporal patterns and chaotic burst synchronization in a small-world neuronal network[END_REF] and the evolution of spiral waves in coupled Hindmarsh-Rose neurons [START_REF] Ma | Breakup of spiral waves in coupled Hindmarsh-Rose neurons[END_REF].

However, an important phenomenon in neuron activity is the transition between different firing patterns. From the viewpoint of dynamical systems, it is crucial to investigate the transition mechanism of firing patterns for understanding realistic neuronal activities. The local stability and the numerical asymptotic analysis of the Hindmarsh-Rose model are used to investigate the firing evolution of a single Hindmarsh-Rose neuron [START_REF] Corson | Asymptotic Dynamics of the slow-fast Hindmarsh-Rose neuronal system[END_REF]. A simple one-dimensional map method has been applied to the HR neurons to convert irregular spiking and chaotic bursting to regular beating and periodic bursting [START_REF] Sabbagh | Control of chaotic solutions of the Hindmarsh-Rose equations[END_REF].

Recently, we developed a parameter switching method to synthesize attractors of a class of dynamical systems, called hereafter attractors synthesis (AS) algorithm. This method in fact starts from a set of given parameter values, and allows us, via periodic or random parameter-switching, to generate any of the set of all possible attractors of a class of continuous time dynamical systems of integer order [START_REF] Danca | Random parameter-switching synthesis of a class of hyperbolic attractors[END_REF], [START_REF] Danca | A switching scheme for synthesizing attractors of dissipative chaotic systems[END_REF], or of fractional order [START_REF] Danca | Fractional-order attractors synthesis via parameter switchings[END_REF], depending linearly on the control parameter. It has been applied successfully to several dynamical systems such as Lorenz, Chen, Rössler, Lotka-Volterra, minimal networks, fractional Lü systems and so on. Extending this subject, we will synthesize attractors of the HR neuronal model, which can exhibit burst dynamics, via the AS method contributing to a better understanding of neuronal firing transition.

The paper is organized as follows: the HR model will be described in Sec. 2; the AS method will be introduced in Sec. 3, and attractors will be synthesized in Sec. 4. Finally, the paper will end in Sec. 5 with some comments and conclusions.

Description of Hindmarsh-Rose model

The dynamics of an isolated HR neuron is governed by the following set of differential equations [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF]:

• x 1 = bx 2 1 -ax 3 1 + x 2 -x 3 + I, • x 2 = c -dx 2 1 -x 2 , • x 3 = p[s(x 1 -x 1 ) -x 3 ], (1)
where x 1 is the membrane potential, x 2 is associated with the fast current, Na + , or K + and x 3 with the slow current, for example, Ca 2+ . The parameters are [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF]:

a = 1, b = 3, c = 1, d = 5, s = 4, x 1 = -1/2(1 + 5)
-1.6. I and p are the control parameters, and I is a slow parameter while p is a fast parameter. I mimics the membrane input current for biological neurons; p controls the speed of variation of the slow variable.

To find the fixed points we have to solve numerically the following equations: x 2 = -5x 2 1 + 1, x 3 = 4x 1 + 6.4, and x 3 1 + 2x 2 1 + 4x 1 + 2 = 0 which gives the single real solution x 1 = -0.639, x 2 = -1.041, x 3 = 3.8442. Therefore the system (1) has a single equilibrium point: X * = (-0.639, -1.041, 3.844) T Parameter p is the ratio of time scales between fast dynamics and slow dynamics. Therefore, it controls the difference between the slow and the fast dynamics of HR neuron model corresponding to the difference between fast fluxes of ions across the membrane and slow ones. Therefore, it is really interesting to investigate the dynamics of the HR neuron as the parameter p is changed. However, with I as control parameter, the underlying dynamical system does not belong to the class of systems where AS can be applied [see Section 3].

The bifurcation diagram of the first component x 1 versus p is shown in Fig. 1a,b,c for I = 3.4, while that for I = 3.5 in Fig. 1d. We are interested in the case I = 3.4 because this case represents interesting dynamics, where chaos appears in a narrow range of p with a width of about 3.5 × 10 -3 . Fig. 1a gives the bifurcation diagram with respect to the control parameter p in the range [0.0001, 0.15]. It is shown that as the parameter p is changed, the HR neuron firstly bifurcates from period-doubling to chaos, and then, it is stopped via the inverse period-doubling. For a clear vision, the enlargement of Fig. 1a is shown in Fig. 1b, which further confirms the above observation. Further enlargement can clearly guide our assessment as illustrated in Fig. 1c. More details on the dynamics of the HR model can be found in [START_REF] Corson | Asymptotic Dynamics of the slow-fast Hindmarsh-Rose neuronal system[END_REF] 3 Attractors synthesis algorithm

The AS algorithm can be applied to the following general class of continuous-time autonomous and dissipative dynamical systems, modeled by the Initial Value Problem (IVP) [START_REF] Danca | A switching scheme for synthesizing attractors of dissipative chaotic systems[END_REF]

S : ẋ = f p (x), x(0) = x 0 , ( 2
)
where f p is an R n -valued function with a bifurcation parameter p ∈ R, n ≥ 3, and has the expression

f p (x) =g(x) + pAx, (3)
g : R n -→ R n is a continuous-time nonlinear function, A is a real constant n × n matrix, x 0 ∈ R n , and t ∈ [0, T ).
This class of dynamical systems contains the bestknown systems such as Lorenz, Chen, Rössler, Lotka-Volterra, minimal networks, fractional Lü systems and so on (see [START_REF] Danca | Random parameter-switching synthesis of a class of hyperbolic attractors[END_REF][15] [START_REF] Danca | Fractional-order attractors synthesis via parameter switchings[END_REF]).

In [START_REF] Danca | A switching scheme for synthesizing attractors of dissipative chaotic systems[END_REF] it has been shown, via numerical analysis and computer graphic simulations, that the AS algorithm allows the synthesis of any attractor of S by parameter switching following some designed rule.

The AS algorithm can explain what happens with a system when, intentionally or not, the parameter value switches quickly through a set of values. Thus, when p is switched following some designed deterministic [START_REF] Danca | A switching scheme for synthesizing attractors of dissipative chaotic systems[END_REF] or random [START_REF] Danca | Random parameter-switching synthesis of a class of hyperbolic attractors[END_REF] rule, while the system evolves in time, an attractor belonging to the set of attractors is generated (synthesized).

The AS algorithm can be useful in the cases where some desired control parameter value cannot be directly accessed and we want to obtain the corresponding attractor. Let us next suppose the following assumptions hold. Assumption A.1. Throughout, the existence and uniqueness of solutions of the IVP (2)-(3) are assumed.

As known, the computer graphic simulations of the numerical integration results of ( 2)-( 3) can give excellent approximations to the orbits within the invariant sets [START_REF] Stuart | Dynamical Systems and Numerical Analysis[END_REF]. Thus, the orbits which start near a hyperbolic attractor will stay near and they will be shadowed by orbits within the attractor. This happens because at-tractors arise as the limiting behavior of orbits. Therefore, the shadowing property of hyperbolic sets [START_REF] Coombes | Rigorous computational shadowing of orbits of ordinary differential equations[END_REF] enables us to recover long time approximation properties of numerical orbits such as HR's case.

The AS algorithm consists in using a time varying, or more precisely, a periodically switching parameter, according to some design rule. It will be shown, empirically by various experiments, that a desired targeted attractor can be duly obtained by the proposed switching scheme.

Remark 1

The algorithm is robust to some extent: the switching timing and switching parameter values both allow flexibility. Therefore, they do not need to be very precisely determined.

Hereafter the following notations will be used Notation 1 Let us denote by A the set of all global attractors of S, including attractive stable fixed points, limit cycles and chaotic (possibly strange) attractors; P ⊂R the ordered set of the corresponding admissible values of p and P N = {p 1 , p 2 , . . . , p N } ⊂ P a finite ordered subset of P containing N different values p, which determines the set of attractors A N = {A p1 , A p2 , . . . , A pN } ⊂ A. Assumptions A.2. P is considered, as in most cases, to be composed of a single real interval and all the values of P N = {p 1 , p 2 , . . . , p N } for which the system behaves stable and/or chaotic is assumed to be accessible. A.3. S is dissipative i.e.

• f < 0, where [START_REF] Ott | Chaos in dynamical systems[END_REF]).

• f ≡ n i=1 ∂f (x 1 , x 2 , . . . , x n ) /∂x i (see e.g.
Due to the assumed dissipativity, A is non-empty 1 and it follows naturally that for the considered class of systems, a bijection may be defined between the sets P and A. Thus, giving any p ∈ P, a unique global attractor is specified, and vice versa.

Remark 2 Because in this paper computer simulations are used as the major analytical tool, the ω-limit set (actually, its approximation [START_REF] Foias | On the numerical algebraic approximation of global attractors[END_REF]) is considered after neglecting a sufficiently long period of transients. Therefore, by attractors (background on the notion of attractor can be found in [START_REF] Milnor | On the concept of attractor[END_REF]) it is appropriate to understand in this paper the ω-limit set obtained by a numerical method for ODEs with fixed step size h after the transients were neglected.

Let P N = {p 1 , p 2 , . . . , p N }. The AS algorithm relies on the following deterministic time rule applied repeatedly on I

[(m 1 h)p 1 , (m 2 h)p 2 , . . . , (m N h) p N ], (4)
where the weights m i are some positive integers. The algorithm acts as follow: in the first time subinterval of length m 1 h, p will have the value p 1 (i.e. the IVP(2)-(3) will be integrated m 1 steps for p = p 1 ), for the next m 2 integration steps, p = p 2 and so on until the N -th time subinterval of length m N h where p = p N and then the algorithm repeats. In order to simplify the notation in (4), for a fixed step size h, the scheme (4) will be denoted next 1 Attractor sets can exist only for dissipative systems because shrinking of the volume in phase space for conservative systems is ruled out by Liouville theorem

[m 1 p 1 , m 2 p 2 , . . . , m N p N ]. (5)
For example, the scheme [p 3 , 3p 1 , 2p 2 ] represents the infinite sequence of p : p 3 , 3p 1 , 2p 2 , p 3 , 3p 1 , 2p 2 , . . . which means that while the considered numerical method integrates ( 2)-( 3), p switches in each m i h time subinterval. In other words, the numerical method will integrate (2)-(3) one step with p = p 3 , then three times with p = p 1 , then two steps with p = p 2 and so on.

The AS algorithm for the scheme ( 5) is presented in Fig. 2 Notation 2 Let us denote the synthesized attractor obtained with AS algorithm, via [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], by A * and by A p * the averaged attractor corresponding to

p * = N k=1 p k m k N k=1 m k . ( 6
)
Next, the following notion is introduced Definition 1 Two attractors are considered to be identical if their underlying ω-limit sets coincide in the phase space, the identity being considered from a geometric point of view in the phase space, aided by computer graphic analysis.

This identity in the case of chaotic attractors will be considered only asymptotically since they are fully depicted only after an infinite time.

Proposition 1

The synthesized attractor A * and the averaged attractor A p * are identical Proof Let us consider some subset P N . If we denote

α k = m k / N i=1 m i , it is easy to see that p * is a con- vex combination p * = N k=1 α k p k , since N k=1 α k = 1.
Therefore p * belongs within the interval (p 1 , . . . , p N ), whatever the values p i are chosen. Also, taking into account the bijection between P and A, we are entitled to consider that the same convex structure is preserved into A. Therefore, for whatever switched values of p in some subset P N , A * will belong within the set A N considered to be ordered by the mentioned bijection i.e. A p * ∈ (A p1 , A pN ). Next, aided by numerical analysis and via computer graphics, it can be showed that A * is identical (in the sense defined above) to

A p * .Therefore A * ∈ (A p1 , A pN ) . t = 0 repeat for k = 1 to N do p = p k for i = 1 to m k do integrate (2) -(3) t = t + h end end until t ≥ T

Next, we can formulate the main property of the AS algorithm

Proposition 2 For whatever considered set P N , the AS algorithm generates an attractor A * which belongs to (A p1 , A pN ) .

Remark 3

The AS algorithm is useful especially when we want to obtain some attractor A p even the underlying value p cannot be accessible.

In this case A can be synthesized by choosing a corresponding set (p 1 , . . . , p N ) p but p / ∈ P N = {p 1 , . . . , p N }, and a corresponding scheme [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF].

The initial conditions play an important role since for a specific value p ∈ P there is a single global attractor but which could be composed by several local attractors (see e.g. [START_REF] Milnor | On the concept of attractor[END_REF], [START_REF] Hirsch | Stable manifolds and hyperbolic sets[END_REF], [START_REF] Hirsch | Differential Equations, Dynamical Systems and An Introduction to Chaos[END_REF], [START_REF] Kapitanski | Shape and Morse theory of attractors[END_REF]). For example, for the Lorenz system for the control parameter r = 2.5 there are three local attractors: the origin and two symmetrical fixed points X 1 , 2 (±2, ∓ 2, 1.5) while for r = 28 there is a single local attractor which is global too, the known Lorenz strange attractor. To avoid these possible difficulties, all the computer simulations for A * and A p * for a particular case, start from the same initial conditions.

In order to see how the AS algorithm works, let us consider that we want to synthesize the attractor A p .

Then we must choose N and the set P N such that p ∈ (p 1 , . . . , p N ) ( p can be equal or not to one of the elements p i , for i = 2, . . . N -1). Next, choosing empirically the scheme [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], such that the right-hand side of (6) gives p, the initial value problem is integrated and finally A p is obtained. To underline the identity between A * and A p * histograms and Poincaré sections besides the phase plots were drawn after the transients were neglected.

Synthesis of HR attractors

Real neurons often display high nonlinearity, which has been shown in many experiments and confirmed by numerical simulation of many neuron models, including the HR and Chay neuron models [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF][START_REF] Chay | Chaos in a three-variable model of an excitable cell[END_REF]. Chaos is a universal phenomenon in certain neurons, such as those in the human brain where, often, the information is decoded and transduced by means of chaos. It is very important to study the chaos of neuron from different aspects. Next, we will focus on the neuronal firing pattern when p is switched. Firstly, it is easy to prove the following Proposition 3 The HR system [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF] with p considered as control parameter, belongs to the class of systems modeled by the IVP ( 2)- [START_REF] Holden | From simple to simple bursting oscillatory behavior via chaos in the Rose-Hindmarsh model for neuronal activity[END_REF].

Proof Choosing the following substitution

y 1 = x 1 -x 1
and replacing finally y 1 with x 1 , (1) becomes .

x 1 = a 1 x 3 1 + b 1 x 2 1 + c 1 x 1 + d 1 x 2 + e 1 x 3 + f 1 + I, (7)
.

x 2 = a 2 x 2 1 + b 2 x 1 + c 2 x 2 + d 2 , . x 3 = p(sx 1 -x 3 ),
where

a 1 = -a, b 1 = b -3ax 1 , c 1 = x 1 (2b -3ax 1 ) , d 1 = 1, e 1 = -1, f 1 = -ax 3 1 + bx 2 1 , a 2 = -d, b 2 = -2dx 1 , c 2 = -1, d 2 = dx 2 1 .
Thus, the right hand side of ( 7) can be written following [START_REF] Holden | From simple to simple bursting oscillatory behavior via chaos in the Rose-Hindmarsh model for neuronal activity[END_REF] with Because, from a numerical point of view, the mathematical model (1) presents a more accessible form than [START_REF] Gonzalez-Miranda | Block structured dynamics and neuronal coding[END_REF], next we will work with (1).

g(x) = ⎛ ⎝ a 1 x 3 1 + b 1 x 2 1 + c 1 x 1 + d 1 x 2 + e 1 x 3 + f 1 + I a 2 x 2 1 + b 2 x 1 + c 2 x 2 + d 2 0 ⎞ ⎠
The HR system verifies the assumptions A1, A2 and also A3 for large parameters range taking into account that

• f = -3ax 2 1 + 2bx 1 -1p. Therefore, one can apply the AS algorithm to the HR system [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF]. For this purpose we considered the Standard Runge-Kutta numerical method for ODEs with the fixed step size h = 0.005 and, as stated above, I = 3.4.

For the sake of the brevity, in what follows, only the most relevant cases are considered. Let us firstly consider the case N = 2 and P N = {0.004, 0.01} (see Fig. 1b) with m 1 = m 2 = 1. Then, using the scheme [1p 1 , 1p 2 ] while integrating the initial value problem (1), finally, even the attractors A 0.004 and A 0.01 are stable limit cycles, the obtained synthesized attractor A * is chaotic (Fig. 3a). As can be seen from Fig. 3a, the synthesized attractor A * (blue) and the averaged attractor A p * (red) for p * = (0.004 + 0.01)/2 = 0.007 given by (6) coincide. It can be seen that p * ∈ (p 1 , p 2 ) (Fig. 1b) and also A * is situated between the corresponding attractors A p1 and A p2 .

Next, if we use N = 10 values for p defined as follows: p i = 0.0002 + i × 0.0001 and m i = 1, for i = 1, 2, . . . , 10, the synthesized attractor A * is a stable limit cycle which coincides with the averaged attractor A p * with p * = 0.00075 (Fig. 4).

An interesting case appears for N = 2 and P N = {0.0082, 0.008765} and m 1 = m 2 = 1, where p 1 , and p 2 are chosen in a very narrow periodic window (its width is about 1E -4 , see Fig. 1c, Fig. 5b and Fig. 5c). The obtained attractor, A * , is identical to A p * for p * = 0.0084825 (Fig. 5 a). We found numerically that the fixed point X * is unstable for p = p * since the three eigenvalues are: -6.266, 0.179 and 0.02 and, as showed in Fig. 6, the Lyapunov spectrum has two negative exponents while the maximum one is positive (approximately zero). All these lead us to consider that A p * (and consequently A * ) is a stable limit cycle. The ap- parent difference between the attractors in this case is due to the very small difference between p 1 and p 2 .

Conclusions

In summary, we have investigated the synthesis of attractors of the HR neuronal system by means of the AS method. It is shown that when we choose the slow parameter p as the control parameter, the HR neuronal system belongs to a class of systems, where we can apply the proposed switching method. Consequently, we concluded that every attractor can be synthesized by the proposed periodic parameter switching scheme in the HR neuron model. It is accepted that neurons can code information by means of firing patterns, which depends on the variation of key parameters of neuronal systems. Hence, the results presented in this paper may be instructive in understanding the implications of realistic neuronal dynamics.

Fig. 1

 1 Fig. 1 (a)Bifurcation diagram: x 1 vs p with I = 3.4. (b)Detail of Fig. 1a. (c) Consecutive detail of Fig. 1b. (d) Bifurcation diagram: x 1 vs p with I = 3.5.

Fig. 2

 2 Fig. 2 Pseudocode of AS algorithm.

Fig. 3

 3 Fig. 3 AS algorithm with scheme [1p 1 , 1p 2 ] with p 1 = 0.004 and p 2 = 0.01; a) The synthesized and averaged attractors A * and A p * respectively, superimposed for p * = 0.007 ; b) Attractor A 0.004 ; c) Attractor A 0.01 .

Fig. 4

 4 Fig. 4 AS algorithm for N = 10, p i = 0.0002 + i × 0.0001 and m i = 1, for i = 1, 2, . . . , 10; A * and A p * , with p * = 0.00075, are superimposed.

Fig. 6

 6 Fig. 6 Temporal evolution of Lyapunov exponents for p = 0.0084825.

Acknowledgements This work was supported by the National Science Foundation of China (Fund Nos. 10972001, 10702023 and 10832006).


Gabriel Peyré
email: gabriel.peyre@ceremade.dauphine.fr

Jalal Fadili
email: jalal.fadili@greyc.ensicaen.frjulien.rabin@cmla.ens-cachan.fr

Julien Rabin

Wasserstein Active Contours

In this paper, we propose a novel and rigorous framework for region-based active contours that combines the Wasserstein distance between statistical distributions in arbitrary dimension and shape derivative tools. To the best of our knowledge, this is the first variational image segmentation algorithm that involves region-dependent multi-dimensional descriptors based on the optimal transport theory. The distributions are represented owing to non-parametric kernel density estimators (e.g. Parzen), and the exact evolution speed corresponding to the Wasserstein-based segmentation energy is provided. To speed-up the computation and be able to handle high-dimensional features and large-scale data, we introduce a sliced Wasserstein approximation of the original Wasserstein distance. The framework is flexible enough to allow either minimization of the Wasserstein distance to known fixed distributions, or maximization of the distance between the distributions of the regions to be segmented (region competition). Numerical results reported to show the advantages of the proposed optimal transport distance with respect to alternative metrics (such as the Kullback-Leibler divergence). These traditional metrics cannot deal properly with distributions having localized supports, and do not take into account the distance between the modes of the histograms. Additionally, our framework handles distributions in arbitrary dimension, which is crucial to segment color images.

Introduction

Overview of the Literature

Contour-based vs. region-based segmentation methods. Active contours for image segmentation methods can be broadly classified as being either edge-based or regionbased. Starting from the seminal work on the snakes model [START_REF] Kass | Snakes: Active contour models[END_REF], edge-based actives contours are driven towards image edges through the minimization of a boundary integral of functions of features depending on edges. They can be viewed as the computation of a locally [START_REF] Caselles | A geometric model for active contours in image processing[END_REF][START_REF] Caselles | Geodesic active contours[END_REF] or globally [START_REF] Cohen | Global Minimum for Active Contour models: A Minimal Path Approach[END_REF] minimal path for a Riemannian metric usually derived from the image gradient magnitude.

Region-based active contours (RBAC) were proposed later. In these approaches, region-based terms can be advantageously combined with boundary-based ones. The evolution equation is generally deduced from a general criterion that includes both region integrals and boundary integrals; see e.g. [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Ronfard | Region-based strategies for active contour models[END_REF][START_REF] Chan | Active contours without edges[END_REF][START_REF] Paragios | Geodesic active regions: A new framework to deal with frame partition problems in computer vision[END_REF][START_REF] Zhu | Region competition: unifying snakes, region growing, and bayes/mdl for multiband image segmentation[END_REF]. The main issue when dealing with RBAC models is the computation of the velocity vector in the evolution equation from the energy functional, especially when the descriptors are region-dependent, as will be the case in this work. This is mostly due to the fact that the set of image regions does not have the structure of a vector space, preventing us from using in a straightforward way gradient descent methods. To circumvent this difficulty, we here take benefit from the shape derivation principles, see [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF][START_REF] Aubert | Image segmentation using active contours: Calculus of variations or shape gradients[END_REF][START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an eulerian accurate minimization method for image and video segmentation[END_REF].

Statistical segmentation A number of authors have proposed RBAC energy functionals involving statistical region-based terms. These are typically functions of the distribution of some image attributes within the region. The distribution can be either parametric or non-parametric, see for instance [START_REF] Martin | Influence of the noise model on level set active contour segmentation[END_REF][START_REF] Zhu | Region competition: unifying snakes, region growing, and bayes/mdl for multiband image segmentation[END_REF][START_REF] Herbulot | Segmentation of vectorial image features using shape gradients and information measures[END_REF][START_REF] Unal | Information-theoretic active polygons for unsupervised texture segmentation[END_REF][START_REF] Heiler | Natural image statistics for natural image segmentation[END_REF]. In the non-parametric approach, the energy functional usually involves a distance (or divergence), e.g. the Kullback-Leibler divergence, between non-parametric kernel estimates (e.g. Parzen kernel) of the underlying densities. To be able to handle localized distributions, this requires a proper choice of the smoothing kernel bandwidth, see for instance [START_REF] Kim | A nonparametric statistical method for image segmentation using information theory and curve evolution[END_REF].

Optimal transport for imaging. To avoid the drawbacks faced with traditional statistical distances, the authors in [START_REF] Ni | Local histogram based segmentation using the wasserstein distance[END_REF] propose to use the Wasserstein distance. In their formulation, the distribution estimated from patches, but does not take into account the explicit dependence of the Wasserstein distance to the region.

The Wasserstein distance originates from the theory of optimal transport [START_REF] Villani | Topics in Optimal Transportation[END_REF]. It defines a natural metric between probability distributions, that takes into account the relative distance between the different modes. In the case of unimodal distributions, this distance reduces to the distance between the modes, thus enabling a precise discrimination of the features.

This Wasserstein distance have found a wide range of applications such as the comparison of histogram features for image retrieval [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF][START_REF] Rabin | A statistical approach to the matching of local features[END_REF], histogram equalization [START_REF] Delon | Midway image equalization[END_REF] and color transfer [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF].

Relation to Previous Work

To the best of our knowledge, the work of [START_REF] Ni | Local histogram based segmentation using the wasserstein distance[END_REF] is the first, and so far the only one, to clearly address the statistical segmentation problem using a Wasserstein distance. Their work clearly emphasizes the usefulness of optimal transport methods to deal with statistically localized features. Our work, however, departs significantly from theirs in many important ways. First, unlike their work which focused on scalar features with the L 1 Wasserstein distance, we consider a general setting in arbitrary dimension. Secondly, the distance in our energy functional is explicitly regiondependent. Using shape derivative tools, we also provide the exact derivative of the Wasserstein distance with respect to the domain boundary and then deduce the active contour velocity field. In contrast, the work [START_REF] Ni | Local histogram based segmentation using the wasserstein distance[END_REF] use a patch-based local histogram estimation to avoid taking into account the dependence of the statistics to the region. While this solution appears appealing by its simplicity, and provides good results on synthetic and natural images, it faces an important dilemma when choosing the patch size. On the one hand, using large patches is important to get consistent estimates in the setting where the distributions have large variances, which is typically the case for textures or noisy images. But this comes at the price of loss of accuracy near the boundaries as the patches overlap between the inside and the outside regions.

Contributions

Our contributions are threefold. (i) We propose a novel rigorous framework for regionbased active contours that combines the Wasserstein distance in arbitrary dimension and shape derivative tools. (ii) We state theoretical results regarding variational minimization of the Wasserstein distance with respect to a domain. Such results might be of independent interest. (iii) We propose an approximate transport distance to speed up the computations for large-scale images carrying multi-dimensional features.

Notations

In the sequel, we consider a feature map I : u ∈ Σ → R d , where u ∈ Σ indexes the pixel location, and d is the dimension of the feature (for instance d = 3 for a color image). We consider a histogram binning grid Ω ⊂ R d .

In the following, we consider Σ as a continuous domain (equipped with the Lebesgue measure) and Ω as a finite discrete grid (equipped with the counting measure). We thus define the Hilbert spaces L 2 (Ω) and L 2 (Σ) endowed with the inner products

A, B Ω = x∈Ω A(x)B(x) and f, g Σ = Σ f (u)g(u)du.
The set of statistical distributions on Ω is

D(Ω) = A \ A(x) 0 and x∈Ω A(x) = A, 1 Ω = 1 ⊂ L 2 (Ω).
where 1(x) = 1, ∀x ∈ Ω.

Non-parametric Statistical Segmentation

Before detailing our Wasserstein distance segmentation approach, we first introduce in this section the tools needed for the estimation of the region distribution at points taking values in the finite discrete grid Ω, together with the required tools from the theory of shape gradients. Although this is standard material, we detail the exact derivatives since this is crucial for our numerical scheme.

Kernel Density Estimator

Given a fixed feature map I : Σ → R d , and a non-negative weight function w ∈ L 2 (Σ), the kernel density estimator of the distribution underlying I is given by

∀ x ∈ Ω, P (w) = Ψ (w) C, w Σ ∈ D(Ω) , (1)
where the mapping Ψ : L 2 (Σ) → L 2 (Ω) and its adjoint are defined as

Ψ (w) : x → Σ w(u)ψ s (x -I(u))du, Ψ * (f ) : u → x∈Ω f (x)ψ s (x -I(u)),
and

C = Ψ * (1) ∈ L 2 (Σ) i.e. C(u) = x∈Ω ψ s (x -I(u)).
Formally, P is a mapping P : L 2 (Σ) → D(Ω). ψ s is a non-negative symmetric smooth localized window called the kernel, and s is its bandwidth. A common kernel function is a Gaussian kernel, and the corresponding estimator is termed (with an abuse of terminology), the Parzen estimator. There are a plenty of other choices such as the Epanechnikov kernel. The choice of s is even more crucial and results from a traditional bias-variance tradeoff, and should be adapted to the discretization grid Ω and smoothness of the underlying density.

As we will show in the next section, in order to optimize energies with respect to a domain, by the chain rule, we will end up needing to to derive Ψ and its adjoint w.r.t. w. This is precisely the purpose of the following proposition.

Proposition 1 We have DP (w) : δ ∈ L 2 (Σ) → 1 C, w Σ Ψ (δ) -P (w) C, δ Σ ∈ L 2 (Ω) , (2)
DP (w) * : µ ∈ L 2 (Ω) → 1 C, w Σ Ψ * (µ) -C P (w), µ Ω ∈ L 2 (Σ). (3)

Statistical Distance-based Segmentation

Let's consider the problem of variational segmentation of the image domain in two regions Σ = Γ ∪ Γ c , where Γ is a regular bounded open set. Γ and its complement Γ c share the same boundary ∂Γ (denoted also C for short), with normals pointing in opposite directions. The goal is to find a (local) minimizer of an energy including both region (Wasserstein fidelity) and boundary (regularity) functionals. The key principle is to construct a PDE from the energy criterion that changes the shape of the current boundary curve according to some velocity field which can be thought of as a direction of descent of the energy criterion.

Shape derivatives of statistical distances. Let's define the region functional

E(Γ, B) = W (P (χ Γ ), B) (4)
for any fixed

B ∈ D(Ω), where χ Γ (u) is the characteristic function of Γ , i.e. χ Γ (u) = 1 if u ∈ Γ ,
∂ C(m, τ ) ∂τ = v Γ (C(m, τ )) and C(•, 0) = C 0 . (5)
Here v Γ is the so-called shape gradient. It ensures that C(m, τ ) converges when τ increases to a stationary point (hopefully a local minimium) of E(Γ, B). The following proposition gives the expression for the shape gradient.

Proposition 2

The shape gradient v Γ ensuring that (5) converges to a stationary point of E(Γ, B) is

v Γ (u) = G Γ,B (u) N u where G Γ,B = [DP (χ Γ ) * ] (∇ 1 W (P (χ Γ ), B)) , (6)
where N u is the unit inward normal to ∂Γ at u, and ∇ 1 W is the (sub)gradient of W with respect to its first variable, and DP * is given in (3).

Proof. See for instance[1, Theorem 6.1][11, Theorem 2].

Level set implementation. For the numerical implementation of the minimization of the energy with respect to the domain Γ , we use the level set method applied to active contours. The key idea is to introduce an auxiliary ϕ : Σ → R, which is often chosen to be the signed distance to ∂Γ . Thus Γ is represented as

Γ = {u ∈ Σ \ ϕ(u) > 0} and ∂Γ = {u ∈ Σ \ ϕ(u) = 0} .
The energy ( 4) is rewritten as W (P (H(ϕ)), B) where H = χ [0,+∞) is the Heaviside function. Introducing an artificial time τ , the evolution equation ( 5) associated to the energy (4) then becomes

∂ϕ(u, τ ) ∂τ = |∇ϕ(u, τ )|G Γ,B (u) where Γ = {u \ ϕ(u, τ ) > 0} . (7)
Note that the velocity function G Γ,B is computed only on the curve C(•, τ ), but we can extend its expression to the whole image domain Σ. However, the signed distance function ϕ is not a solution of the PDE [START_REF] Delon | Midway image equalization[END_REF], and in practice, it must be periodically reinitialized so that it remains a distance function. This is important to ensure numerical stability of the method.

Statistical Segmentation by Region Competition

In the same vein as [START_REF] Jehan-Besson | Shape gradients for histogram segmentation using active contours[END_REF][START_REF] Herbulot | Segmentation of vectorial image features using shape gradients and information measures[END_REF], we restrict our attention in this paper to a non-parametric variational segmentation method that seeks the maximization of the distance between the respective distributions in Γ and Γ c , i.e. region competition. Of course, our approach can be applied to other energy functionals just as well, e.g. those with terms that favor region homogeneity. The energy functional to be minimized reads

min Γ E(Γ ) = -W (P (χ Γ ), P (χ Γ c )) + λr(C). ( 8
)
where r(C) is a boundary regularity term, e.g. the curve length. Written using the level set formalism, this corresponds to the solution of min ϕ -W (P (H(ϕ)), P (H(-ϕ))) + λR(ϕ).

where R(ϕ) is a suitable regularization associated to r(C). For instance, if r(C) is the length, then R(ϕ) is the TV regularization R(ϕ) = |∇ϕ(u)|du.

The equivalent level set evolution PDE [START_REF] Delon | Midway image equalization[END_REF] that drives an initial contour to a stationary point (hopefully a local minimizer) of ( 8) is

∂ϕ(u, τ ) ∂τ = |∇ϕ(u, τ )| (G Γ c ,B Γ (u) -G Γ,B Γ c (u) + λκ) (10)
where Γ = {u \ ϕ(u, τ ) > 0} is the domain at time τ (we have dropped the dependency on τ for the sake of clarity), where the histogram inside and outside Γ are

B Γ = P (H(ϕ(•, τ )) and B Γ c = P (H(-ϕ(•, τ )),
where G Γ,B Γ is the velocity defined in [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF], and κ = div ∇ϕ |∇ϕ| is the mean curvature of the boundary. In practice, this PDE is discretized with a sufficiently small time step (for instance computed using a line search).

Wasserstein Distance

Previous work, including [START_REF] Jehan-Besson | Shape gradients for histogram segmentation using active contours[END_REF][START_REF] Herbulot | Segmentation of vectorial image features using shape gradients and information measures[END_REF][START_REF] Unal | Information-theoretic active polygons for unsupervised texture segmentation[END_REF], have specialized the above statistical framework to the case of point-wise statistical metrics. Among the most popular is the Kullback-Leibler distance (i.e. symmetrized divergence)

W (A, B) = W (A, B) + W (B, A) where W (A, B) = x∈Ω A(x) log A(x) B(x) . (11)
While this distance enjoys many appealing statistical/information theoretic properties, it suffers from some drawbacks. The most prominent one is that it assumes that both distributions share the same support. In particular W (A, B) = +∞ if A(x) = 0 and B(x) = 0 for some x ∈ Ω, and vice versa. Furthermore, such a metric does not care about the relative positions of highly localized modes in the distribution. Put differently, this means that W (A, B) = W (σ(A), σ(B)) where σ(A)(x) = A(σ(x)) for any permutation σ of the grid points. Both issues make such a metric not very robust when comparing localized distributions. A practical but artificial way to somehow alleviate these difficulties is to use a large smoothing bandwidth s in (1), but this is likely to yield oversmooth distribution estimates.

Our proposal is to use the Wasserstein distance which is not prone to the above drawbacks, since it takes into account the relative distances between the grid points.

Discrete Optimal Transport

We consider two distributions A, B ∈ D(Ω). An optimal transport cost on D(Ω) is defined as

W (A, B) = min P ∈P(A,B) C, P Ω×Ω = (x,y)∈Ω 2 C(x, y)P (x, y) (12)
where C is a fixed cost matrix and P(A, B) is the polytope of stochastic matrices with marginals A and B:

P(A, B) =    P \ P (x, y) 0, y∈Ω P (x, y) = A(x), x∈Ω P (x, y) = B(y)    .
When the cost is C(x, y) = ||x -y|| p , W (A, B) 1/p is often called the L p Wasserstein distance, and is a distance on D(Ω). For the sake of simplicity, we will refer to W as being the Wasserstein distance, and will not use the exponent 1/p.

Wasserstein Subdifferential

The Wasserstein distance can be re-written using the dual problem to ( 12)

W (A, B) = max (u,v)∈Q u, A Ω + v, B Ω = x∈Ω u(x)A(x) + v(x)B(x) (13)
where

Q = (u, v) \ ∀ (x, y) ∈ Ω 2 , u(x) + v(y) C(x, y) ,
see [START_REF] Villani | Topics in Optimal Transportation[END_REF]. Note that Q is independent of A and B. Rephrased in the language of convex analysis, the Wasserstein distance W (A, B) is nothing but the support function of the closed convex set Q. Thus, using properties of Legendre-Fenchel conjugacy, closedness, convexity and subdiffentiability properties of W (A, B) can be established.

Proposition 3

The function A ∈ P(Ω) → W (A, B) ∈ R is closed and convex for any fixed B ∈ P. Its subdifferential at some A ∈ P(Ω) is such that

u ∈ ∂ 1 W (A, B) ⇐⇒ ∃v s.t. (u, v) ∈ Argmax (q,p)∈Q q, A Ω + p, B Ω .
Proof. Closeness and convexity follow from [14, Theorem X. In plain words, a subgradient of W at A, for a given B, belongs to the set of (global) maximizers of [START_REF] Kim | A nonparametric statistical method for image segmentation using information theory and curve evolution[END_REF]. A maximizer can be computed by solving the linear program (13) using dedicated methods in roughly O(|Ω| 3 ) operations [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF].

It is worth noting that these subgradients are defined up to an additive constant. Indeed if (u, v) is a solution to (13), then so is (u + α, vα). Nevertheless, this is of no importance for the computation of the velocity vector [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF] which is invariant to constant perturbations of u, since for any w ∈ L 2 (Ω), [DP (w) * ](1) = 0.

When the subdifferential is a singleton ∂ 1 W (A, B) = {u} (meaning that (13) has a strict global maximizer), then A → W (A, B) is Fréchet differentiable at A and ∇ 1 W (A, B) = u. In the other cases, one can take any subgradient of W (A, B) in lieu of ∇ 1 W (A, B) in [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF].

Sliced Wasserstein Distance

The computation of W (A, B) and ∇ 1 W (A, B) for the Wasserstein distance ( 12) is however demanding for large-scale and high-dimensional histograms. Numerically, the computation time is acceptable (less than a second on a standard laptop computer) for only |Ω| ≪ 10 3 . For the segmentation application we are targeting, typically 3-D color histograms with at least 10 4 grid points, solving this linear program is not an option.

To speed up the computation, we follow [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF] and consider an alternative distance that mimics the properties of the Wasserstein distance, but is faster to compute. In the following, we restrict our attention to the L p Wasserstein distance, meaning that C(x, y) = ||x -y|| p .

Sliced Distance Approximation

The sliced Wasserstein distance reads

SW (A, B) = θ∈Θ W (A θ , B θ ) ( 14
)
where Θ is a finite subset of the unit sphere in R d , and A θ ∈ D(Ω θ ) is the projected distribution in the direction θ, defined on 1-D grid points Ω θ = {x θ = x, θ } x∈Ω , that has the same values as A, i.e. ∀ x ∈ Ω, A θ (x θ ) = A(x). The sliced Wasserstein distance is thus a sum of 1-D Wasserstein distances between the projected distributions. Note that the grid Ω θ that supports these distributions is non-uniform, so that some care is needed to perform the computation, as detailed in the following section.

The method described in [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF] considered probability distributions A(x) defined on a varying grid x ∈ Ω (that depends on A) with constant weights A(x) = 1/|Ω|. But this entails non-convexity of the function A → W (A, B). In contrast, we consider here distributions A ∈ D(Ω) on a fixed grid Ω but with varying values A(x), which thus preserves convexity of the sliced Wasserstein distance, whose subdifferential ∂ 1 SW (A, B) can be easily computed as we will show shortly.

Proposition 4

The function A → SW (A, B) is closed and convex and its subdifferential at A is such that

∀ x ∈ Ω, ∂ 1 SW (A, B)(x) = θ∈Θ ∂ 1 W (A θ , B θ )(x θ )θ.
Proof. Closedness and convexity are preserved under the sum. The formula is a consequence of subdifferential properties and calculus rules [START_REF] Lemaréchal | Convex Analysis and Minimization Algorithms I and II[END_REF] and linearity of the projection A → A θ .

Wasserstein Distance on a Non-Uniform 1-D Grid

The ability to compute [START_REF] Lemaréchal | Convex Analysis and Minimization Algorithms I and II[END_REF] and to (sub)differentiate it, is conditioned to that of computing quickly the Wasserstein distance between two distributions A, B ∈ D(Ω) defined on a 1-D grid Ω ⊂ R. We assume that the grid points Ω = {x i } N i=1 are sorted in increasing order, x i x i+1 . This is achieved in O(|Ω| log(|Ω|)) operations, which will turn out to be main computational cost of the discrete 1-D Wasserstein distance and its (sub)gradients.

L p Wasserstein Distance. The L p Wasserstein distance on the real line for any p 1 reads

W (A, B) = 1 0 |R -1 A (t) -R -1 B (t)| p dt , R A (s) = s -∞ A(x)dx, is the cumulative distribution function (CDF) and R -1 A (t) = inf {s \ R A (s) t} its pseudo-inverse. The latter is well defined as the CDF is non- decreasing. When A ∈ D(Ω) is discrete, the CDF is equal to R A (s) = xi s A(x i ).
Wasserstein distance sub-differential. The following proposition gives the sub-gradients of the Wasserstein distance for discrete distributions.

Proposition 5 Let A, B ∈ D(Ω). For p 1, the sub-gradients of A → W (A, B) are written as

∇ 1 W (A, B) : x i → j i |x j -xj | p -|x j+1 -xj | p (15)
where

xj = x k if R B (x k-1 ) < R A (x j ) < R B (x k ), xj ∈ [x k , x k+1 ] if R A (x j ) = R B (x k ).
Proof. For simplicity, we treat the case where the set of values of R A and R B are disjoints, and A → W (A, B) is thus differentiable. The general case is obtained by replacing R -1 B by a set-valued mapping. We can write

W (A, B) = W (R A , B) where W (R, B) = N -1 i=0 R(xi+1) R(xi) |x i+1 -R -1 B (s)| p ds.
Deriving with respect to R A (x j ) gives

∇ 1 W (R A , B) : x j → |x j -R -1 B (R A (x j )) | p -|x j+1 -R -1 B (R A (x j )) | p .
Using the chain rule, introducing R * the adjoint of A → R A , we obtain

∇ 1 W (A, B) = R * ∇ 1 W (R A , B)
where R * (U ) :

x i → j i U (x j ),
which gives [START_REF] Martin | Influence of the noise model on level set active contour segmentation[END_REF].

5 Numerical Results

Comparison of Wasserstein and Kullback-Leibler Distances

We first illustrate the difficulty of point-wise statistical metrics such as the Kullback-Leibler symmetrized divergence (11) on a simple instructive example with localized 1-D (d = 1) distributions. Figure 1, left, shows an example of a gray-scale image with three concentric regions delimited by two circles of increasing radii r 0 < r 1 .

The distribution of the intensity values within each region Gaussian with a mean value in {0.05, 0.8, 0.95} of small variances, so that the resulting image is a mixture of three localized and barely overlapping Gaussians. As can be seen from the histogram in Figure 1 top-right, the correct segmentation should group together the two outer regions which have close means. Figure 2, shows the energy landscape for circular regions of radius r. The leftmost figure shows this energy for the Wasserstein distance, which provides the proper segmentation. Indeed, circle of radius r 0 is the only local, and hence global, maximum of the L 2 (p = 2) Wasserstein distance between the inside and outside regions. The situation is radically different with the Kullback-Leibler divergence (rightmost figure), where a spurious local minimum for r = r 1 persists, unless an extremely large kernel bandwidth is used. The weakness of point-wise statistical metrics is thus apparent when localized feature histograms come into play, and the smoothing impacts significantly the spatial localization. Note that a similar conclusion is arrived at in [START_REF] Ni | Local histogram based segmentation using the wasserstein distance[END_REF] using a different segmentation method and a L 1 Wasserstein metric.

Synthetic Image Segmentation

We consider a synthetic segmentation example with a realization of two 2-D Gaussian mixtures (hence d = 2)

I(x) ∼ 1 2 N (µ 1 , ρ 2 ) + 1 2 N (ν 1 , ρ 2 ) if x ∈ Γ 0 , 1 2 N (µ 2 , ρ 2 ) + 1 2 N (ν 2 , ρ 2 ) if x / ∈ Γ 0 .
where the standard deviation of the Gaussians is ρ = 0.05, and the means µ i , ν i ∈ R 2 are located on the four corners of a square, see Figure 3, last row. Note that the distributions are 2-D, hence the image colors are located in the red/green 2-D plane.

We solve the region competition problem (8) using the sliced L 2 (hence p = 2) Wasserstein distance [START_REF] Lemaréchal | Convex Analysis and Minimization Algorithms I and II[END_REF]

Natural Images Segmentation

We finally report results on complex color natural images each of size 200×200 pixels.

The results are displayed in Figure 4. For these examples we use |Θ| = 12 random directions in 3-D for the computation of the sliced Wasserstein distance [START_REF] Lemaréchal | Convex Analysis and Minimization Algorithms I and II[END_REF]. In this experiments, we use the L 2 Wasserstein distance, i.e. p = 2.

In each case, the initial contour is a circle of radius 0.4 (assuming an image defined over Σ = [0, 1] 2 ). The figure also depicts 2-D histograms to show how the global color distribution is split by the segmentation algorithm. Only a 2-D slice of the full 3-D histogram is displayed, along the two dominant colors provided by a principal component analysis. The Wasserstein active contour is able to split properly the color space consistently with what an observer would do visually.

Conclusion

We have proposed a mathematically grounded way to handle the statistical segmentation problem in arbitrary dimension. Our framework combines wisely the Wasserstein statistical distance and shape derivative tools. It can handle localized distributions, while this property seems to be out of reach for traditional metrics unless a severe smoothing is applied. Such a smoothing would clearly harm the segmentation precision. This approach is quite general and paves the way to many applications. One may for instance think of more advanced features beside simple colors, typically coefficients in the domain of a multi-scale transform or joint statistics between groups of pixels.

  and 0 otherwise. If we introduce an artificial time τ for the evolution, and consider m ∈ [0, 1] → C(m, τ ) to be a parametric representation of the boundary ∂Γ at time τ , a gradient flow of this boundary reads

1 . 1 . 2 ]

 112 . The subdifferential characterization is a consequence of Fenchel identity [14, Theorem X.1.4.1].

Fig. 1 .

 1 Fig.1. Left: gray-scale image example with three concentric regions with radii r0 < r1. Each region has its intensity values Gaussian-distributed with means {0.05, 0.8, 0.95} and the same variance. Right: estimated densities using the P (χ Γ ) with the optimal partition Γ = Γr 0 and using the Gaussian kernel with bandwidth s = 10 -2 (top) and s = 0.2 (bottom).

Fig. 2 .

 2 Fig. 2. Energy W (P (χ Γr ), P (χ Γ c r )) for a centered circle Γr of radius r as a function of r. Each curve corresponds to a different smoothing bandwidth s in the Parzen kernel estimator.

Figure ( 3 ))Fig. 3 .

 33 Fig. 3. Top row: evolution of the segmentation, with ℓ indexing the iterations, for |Θ| = 8 directions. Middle row: same evolution with |Θ| = 2 axis-aligned orientations. Bottom row: 2-D histograms for the whole domain Σ (which is P (χ Σ )) and inside and outside the central disk Γ0.

PFig. 4 .

 4 Fig. 4. Left: example of natural image segmentation. Below each image is displayed the 2-D histogram of the image (in the space of the two dominant color eigenvectors as provided by a PCA) for the whole domain Σ (which is P (χ Σ )) and over the inside and outside the segmented region.