Nicolas Tchertchian

Dominique Millet

Pierre Alain Yvars

The influence of the level of definition of functional specifications on the environmental performances of a complex system. EcoCSP approach

Keywords: Life cycle assessment, Constraint satisfaction problem, Ecodesign, Complex system, ship design, sustainable development

The tendency towards a homogenous mode of development modeled on that of Western countries means that sustainable development has become increasingly urgent. It is necessary to thoroughly redefine products and their expected performances in such a way that the consequences are compatible with sustainable development.

In the domain of product design, this means that it is no longer sufficient to use assessment tools « after the fact » to check the impact of products whose functional unit was definedprior to production; it is now necessary to rethink the definition of the functional unit itself. This article aims to present an approach basedon a combination of life cycle analysis methods (LCA) and problem solving by constraint satisfaction (CSP). This original approach makes it possible to vary the design of the different dimensions of the functional units of a complex system and thus to make it easier to identify the best architecture along with the best functional definition of the system.

In this study, the EcoCSP approach is applied to define the functional performances of an ecological passenger ferry.

The complexity of couplings between sub-systems and the sheer number of those sub-systems mean that the designer has to use « intelligent » tools. These simulate a great number of scenarios and help him/her to fine-tune the system and make the right technological choices with regard to the right functional specifications.

Introduction

Industrial growth has affected the conditions of life on earth as witnessed by the increase in ecological problems that threaten future generations. In response to this situation, the concept of sustainable development that began in the 80's, has gradually become mainstream. Sustainable development calls on all actors involved in the evolution of developed societies to balance the economic social and environmental dimensions of their activity.This means setting a new paradigm where reasonable consumption/production has become essential. In order to reach this objective, products and systems must be designed to be sustainable. Sustainable design is different from Eco design and Design for Environment because it goes beyond the environmental optimization of goods and services (Van Weenen1995): it attempts to incorporate the considerations demanded by the three pillarsof sustainable development, namely social, economic and environmental factors.

Various authors [START_REF] Daly | Towards a steady-state economy[END_REF][START_REF] Simonis | Preventative environmental policy: Prerequisites, trends and prospects[END_REF]Williams et al. 1987;[START_REF] Herman | In technology and environment[END_REF]Ayres and Kneese 1990;[START_REF] Freeman | The economics of hope: Essays on technical change, economic growth, and the environment[END_REF]) have mentioned the radical nature of the technological transformation that needs to be effected in order to improve the environmental performance of a product or system: they have recommended reducing the proportion of material in the economy using expressions such as X Factor, eco-efficiency, industrial ecology, functional economy, dematerialization, product service-system etc.

Today, traditional Eco design approaches either carry out curative environmental assessments (LCA) (Hauschield 2005), or lead designers towards improved solutions by providing guidelines (Wimmer 2001). Both of these approaches, used in the design of complex systems, most often result in global under-optimizations that are unsuited to the design of complex systems. It thus appears necessary to implement new Eco design practices that are better suited to designing such systems.

A complex system can be defined [START_REF] Krob | Eléments d'architecture des systèmes complexes, Gestion de la complexité et de l'information dans les grands systèmes critiques[END_REF], [START_REF] Cilliers | Complexity and postmodernism: Understanding complex systems[END_REF] as a system comprising numerous sub-sets that are interdependent; each of these sub-sets has several possible alternative solutions. In general, complex systems have different ways of functioning with performances that change according to conditions of use. Finally, the long-term life cycle of a complex system is not easy to predict during the design phase; this is particularly true for life duration, maintenance, component upgrading and end of life.

Any increasein complexity results in the multiplication of technical solutions and thus of possible alternatives. In such cases, design becomes a long process of negotiation within the design team. This negotiation is generally based on an initial definition of the system's specificationsspecifications that arerarely questioned during the design process. In this article, we focus on the necessity for the actors concerned to generate a functional negotiation, that is to say, to select the « right » functions, then the characteristics of the system in order to optimizeits environmental performance.

In the following section of this article we deal with the problem of defining functional units in assessing the life cycle of complex systems. In the third part of the article we give a theoretical introduction to the EcoCSP method; this method combines CSP and LCA in order toidentify the optimal architecture of a system by negotiating the functional unit so that an environmental optimum may emerge. In section 4, the EcoCSP approach is applied in the context of designing a new passenger ferry with hybrid technology. Finally in section 5 we discuss the results and estimate the environmental improvement compared to a system with fixed functionalities. Our conclusion in section 6 summarizes the contributions of this article to a global approach, part of which includes EcoCSP.

2. The innovative Eco designapproach:reassessing functionalities 2.1. Improving environmental performances by reassessing product functions

As underlined by [START_REF] Lagerstedt | Functional And Evironmental Factors In Early Phases Of Product Development-Eco Functional Matrix[END_REF], environmental performance generally depends on product functionalities.

However, from another point of view, the commercial success of a product depends on the functions it offers to users. [START_REF] Lagerstedt | Functional And Evironmental Factors In Early Phases Of Product Development-Eco Functional Matrix[END_REF] mentions the balance that must be found between the « environmental cost » and the « functional gain ». Few methodological supports exist in the domain of tools/methods for environmental improvement or hirarchization, and amongst those that do exist, even fewerenable early intervention in the design process. Current methods of Eco design such as life cycle analysis and other assessment methods derived from this, such as environmental guidelines and checklists, merely identify the causes of environmental problems in order to redesign the product while keeping its functionalities unchanged; this is in contradiction to strategies of radical environmental improvement (X Factor) that necessitate a complete reassessment of product functionalities. Achieving a higher degree of sustainable development requires finding a balance between acceptable impacts and necessary functions. [START_REF] Luttropp | Strategy and material flow in ecodesign[END_REF] presents different ways of reaching this balance: he favors reducing environmental impacts while increasing the level of the product's functional performance -a win-win situation that eliminates all unnecessary functions. On the other hand, he is critical of the « green fix » strategy (using new materials while keeping all the functions) that result in short term, temporary optimizations; he also judges inefficient the « linear down » strategy (improving environmental impact by downgrading or eliminating functions).

2.2. The problem of defining the functional unit in the LCA method.

Life cycle analysis (LCA) is a method of environmental assessment of a product or service over the whole of its life cycle, that is, from the phase of extracting the raw materials and manufacturing the product until the end of life (discharge, recycling, reuse etc.), including distribution, use and maintenance. The methodological framework of LCA is governed by ISO 14040; this distinguishes 4 phasesdefining the objectives and the perimeter of the study, taking the inventory of the life cycle, assessing the impactsof the life cycle and interpreting the results). The phase of defining the objectives and perimeter of the study requires the definition of a functional unit. The functional unit is the « quantified performance of a system of products to be used as the unit of reference in a life cycle analysis » (ISO 14044, 2006).The definition of this functional unit is crucial.

Indeed, in cases where the LCA study aims to analyze the potential impacts of different options, it is imperative that all the options assessed fulfill the same function in order to be comparable [START_REF] Jolliet | Analyse du cycle de vie (ACV), Comprendre et réaliser un écobilan[END_REF]. Now, by constraining the designer to reason by iso-functionality, the LCA methods and its derivations naturally hinder thinking about products that might have a better balance between environmental cost and functional gain (Luttropp and Lagerstedt2006). In general, the available tools, amongst them LCA, are based on a single criterion: the main function expressed in the form of a functional unit [START_REF] Lagerstedt | Functional And Evironmental Factors In Early Phases Of Product Development-Eco Functional Matrix[END_REF]. This means that very different products or concepts can be compared.

Consequently, when comparative LCA's are undertaken for these types of products (i.e. that have several functions), it is important to consider the other sub-functions. If these functions are not identified, broken down, specified and/or prioritized in the right way with regard to the objectives and perimeter of the study, it could result in a functional unit that does not reflect reality. As underlined by [START_REF] Reap | A survey of unresolved problems in life cycle assessment[END_REF], these are important questions, for they can downgrade the precision of the reference flows associated with the chosen functional unit and thus decrease confidence in LCA results.

2.3. State of the art review of environmental assessments (LCA) on systems of transport and vehicles.

In the literature, numerous products and systems have been environmentally assessed using LCA. In order to highlight the problems related to the definition of a functional unit, about thirty LCA studies published by 3 scientific publishers (Springer, Taylor & Francis, Elsevier) between 2003-2013 as well as a few LCA studies presented in the context of doctoral theses were evaluated in the domain of transport systems.

Each of these LCA studies is characterized by the system (mainly vehicles), the functional unit attributed by the author and the parameters that were modified during the sensitivity test. The results are given in Table 1. This non-exhaustive state of the art on new technologies in the transport sector (especially cars) shows that the functional unit is very often assimilated to a principle function: that of transporting a person or an object from A to B over a distance of x thousand kilometers. As stressedin the previous paragraph, and highlighted by [START_REF] Reap | A survey of unresolved problems in life cycle assessment[END_REF], vehicles are complex systems with sub-functions that must be taken into account when two systems are being compared.

Table1:State of the artreview of LCA studies in the domain of transport.

In thirteen studies the variability of functional unit is not taken account. Among the functional modifications suggested, the modification of the driving cycle (annual mileage, vehicle lifetime, use), reduction of the vehicle's mass, the number of battery charges, etc. All these functional modifications are integrated into the process of sensitivity analysis required by the ISO norm. Now, in actual fact, the designer only looks at the extra gains or impacts generated by the modification of one or several parameters ; he/she takes no account of the consequences on the whole of the system. These modifications are not completed by a redefinition of the design parameters. Indeed, if we take the example of the reduction of a vehicle's mass, this could potentially call for a different distribution of the masses of the vehicle, and thus a modification of the system's aerodynamic performances that might generate a resizing of the propulsion.

Presentation of the EcoCSP approach

The EcoCSP approach is a further development of the CSP/ LCA approach proposed by [START_REF] Tchertchian | Benefits and limits of a Constraint Satisfaction Problem/Life Cycle Assessment approach for the ecodesign of complex systems: a case applied to a hybrid passenger ferry[END_REF]). This approach is based on a combination of 2 methods « Constraint Satisfaction Problem»/Life Cycle Assessment.

Definition of a CSP

A CSP (Constraint Satisfaction Problem) is defined by (Montanary, 1974): • Extensive constraints:a constraint in extension describes an explicit and exhaustive list of possible, or on the contrary, impossible combinations -of values (m-tuples) between the m variables at play within the constraint.

X = {x1, x2, x3, …,
• Intensive constraints: a constraint in intension is an explicit equivalence (or non equivalence) linking two variables to each other (equality or inequality). It brings linear and/or non-linear operators into play.

• Logical constraints: conditional constraints (IF....THEN), conjunction of constraints (AND), disjunction of constraints (OR), obtain logical combinations of constraints. In the case of designing a complex system, logical constraints establish composition relationships among the system's components and define « components » whose state of functioning varies over the life cycle.

Environmental criteria come into the CSP approach as constraints to be satisfied in order to respond to an objective of environmental optimization. Depending on the objectives required by the designer, the algorithm of resolution indicates system architectures and operation modes (if these exist) that respect these constraints.

Figure 1 shows the resolution of CSP.

CSP resolution

Resolution by constraint satisfaction results in a complete set of solutions that enable the design team to choose one suited to a design problem according to specific performance variables and constraints. A CSP is typically solved by reducing the domains. The objective of propagating constraints is to replace an initial CSP by an equivalent one which has a more restricted research space. The most basic way of reducing the area for solution is to proceed by trial and error or by dichotomy. A more systematic method is to use filtering techniques which rely on arithmetic of intervals [START_REF] Moore | Life-cycle assessment of a Solar Assist Plug-in Hybrid electric Tractor (SAPHT) in comparison with a conventional tractor[END_REF] and propagation of constraints [START_REF] Mackworth | Consistency in Networks of Relations[END_REF].The most commonly used domain filtering techniques are Arc-Consistency (Mackworth, 1977, Debruyne and[START_REF] Debruyne | Domain Filtering Consistencies[END_REF] for discrete CSP's, Hull-Consistency [START_REF] Lhomme | Consistency Techniques for Numeric CSPs[END_REF][START_REF] Benhamou | Interval constraint logic programming[END_REF][START_REF] Benhamou | Applying Interval Arithmetic to Real, Integer and Boolean Constraints[END_REF] and

Box-Consistency [START_REF] Chiriaev | Interval arithmetic specification[END_REF]Walster, 1997, Benhamou et al., 1999) for discrete and continuous CSP's.As [START_REF] Chenouard | Résolution par satisfaction de contraintes appliquée à l'aide à la décision en conception architecturale[END_REF] points out, using CSP in preliminary design has the advantage of great flexibility for expressing knowledge and modifying models; it resolves generic problems. This is a sought after characteristic in design, for it expresses knowledge without defining how it should be dealt with. CSP makes it easier to manipulate and reuse such knowledge.

The EcoCSP approach : a development of the CSP/ LCA approach

The methodology of life cycle analysis uses a normalized functional unit (UF n ),to facilitate comparisons among systems that show unequal performances.

Our state of the art review and the compilation of [START_REF] Reap | A survey of unresolved problems in life cycle assessment[END_REF]of the main problems posed by LCA, show that defining a functional unit is not sufficient for the radical improvement of the environmental performances of

a complex system. × + × B A C A1 B1 C1 A2 B2 C2 C3 A11 A21 C31 C21 A2B1 A B C Design parameter Xa, Ya, Za, Xb, Yb, Xc, Yc, Zc Constraints Constraint A2=0 Constraint C3=1 A B C A1 C1 B2 C2 C3 × A11 A21 C31 C21 A B C A1 B2 C3 A11 A21 C31

Modelisation of system

Generation of alternatives

Constraint Satisfaction relationships of composition + : AND X : OR The works of [START_REF] Tchertchian | Benefits and limits of a Constraint Satisfaction Problem/Life Cycle Assessment approach for the ecodesign of complex systems: a case applied to a hybrid passenger ferry[END_REF] demonstrated the relevance of an approach combining Constraint Satisfaction Problem solving approach (CSP) and Life Cycle Analysis method to define the best architecture and operation modes of a complex system with respect of environmental constraints. This fruitful research pointed to a way forward. Indeed, the CSP approach allows us to modelise functional requirements as constraints; by exploring these as such, it is then possible to simulate the various architectural alternatives of a complex system while at the same time varying the specifications of that system. We have called this approach EcoCSP.

The general framwork of the EcoCSP approach, represented by a flowchart Figure 2, contains a phase of predefinition of the functional specifications of the system, in which are listed the Negotiable Functions (NF), these areflexible functions in the specification of the " Performance level " expected by the system, and the Non Negotiable Functions (NF̅ ), these are functions whose performances are specified from the beginning(ex: functionsrelated withsecurity).

NF i ∈ [P imin , P imax ], whereP imin is the minimum level of performance of the function i and P i max the maximum level of performance.

In parallel with the definition of overall performance of the system, the design team defines the range of variation of the Functional Unit (FU). The functional unit defines the quantification of performance characteritics of the productintended to be used as reference. In EcoCSP approach the FUis sought to minimize the environmental burden of the system archuitecture.The globally optimized performance of the system (FU GO ) may vary in the domain [P s min , P s max ].

Significant Negociable Function (SNF) which variation of the performance level led to a significant fluctuation of the environmental impact of the system is identified through design of experiments via Taguchi orthogonal array L z [START_REF] Taguchi | Quality Engineering In Production System[END_REF]. L z is an array where z is the number of simulations required to represent all the effect of NF on environmental criteria.

A design of experiments (DoE) is an efficient way to identify the contribution of each function to the environmental impact. It reduces the number of trials (CSP simulations). Without DoE the number of simulations required to test all functional mixes is m n with n the number of functions and m the number of performance level for each function.

In this first methodological proposal, the functions are not supposed to interact among themselves and each of them have m = 2 levels. Each function NF i whose performance P i is in the domain [P i min , P i max ] the two levels are the upper P i max and lower P i min bounds; for functions NF i whose performance is in discrete domain{P i1 , …, P ik } the two levels are the minimum and maximum values.

For example, a functional mix composed of n = 7 NF and with m = 2 performance levels would require 2 7 trials to test all variants of functional mix. The appropriate Taguchi orthogonal arrayis L8 (2 7 ), reducing the number of CSP simulations toz = 8 instead of 128.

Each variant of the functional mix is characterized by a line in DoE (cf. Table2).The results are used to draw the graph of main effect on the average of Environmental Impact (see Figure 3), i.e the effect of each NF i on environmental performance.

In Figure 3, for example, the transition from performance level P 1 min to P 1 max of function F1causes a significant change in the average environmental impact EI Avg , the function NF1 is a significant negotiable function. The three most "significant" NF are retained in the next step to specify their optimal performance level. While the previous approach deals with a functional unit normalised (FU n ) and functional performance defined (P i n ) according to system specifications defined at the beginning of design process ; the EcoCSP approach allows to model the three significant Negotiable Functions as constraints that operate on domains of

F 1 max EI Avg F 3 min Environmental Impact (Pts) F n max F n min F 1 min F 2 max F 2 min F 3 max
performance: [P i min , P i max ]. Similarly, the functional unit FU GO is defined as a constraint in domain [P s min , P s max ].

A limited variation in the level of performance of FU GO compared to the baseline performance of the FU n allows a radical improvement in environmental performance system.

The optimal solution is generated by constraint satisfaction (CSP).Design parameters defining technological choicesfor each subsystem (eg motor sizing), the way of components operate according to the use phase sequences (eg componant activate or not) and the performance level of 3 Functions Negotiable are specified by propagating the constraints.

A CSP solver is used to instantiate the design variables and the performance P i of Negociable Functionsthat minimize the performance variables based on environmental criteria. The constraint solver used is ILOG Solver, developed by IBM. ILOG Solver is a C ++ library. In the flowchart, the solver is involved in 3activity box labelled CSP to satisfy the constraints; generating technological choices, operation modes of components and the performance level of 3 Negotiable Functions.

Finally, we propose an aera of improvement to reduce the impact caused by overproduction of a non-suitable components (oversizing, non-mature technology, etc.) towards the Negotiable Functions performed. Reducing the performance level of the functions of the system creates a need for appropriately sized components to meet performance adapted downwards or upwards. The number of components constituing such databases are not exhaustive, they depend on their availability on the market. Extending the database of components to meet the appropriate performance level should achieve additional environmental benefits. This observation leads to imagine (re) design specific component to achieve the optimal performance level.

EcoCSP allows judgments and choices to be made about functions on the basis of those functions that are deemed negotiable; the approach makes it possible to vary the system's performance in order to reduce environmental impacts.

4. Case Study : designing an eco compatible hybrid passenger ferry 4.1. Simplified modelisation for a complex system -a maritime ferry

The system under study is a maritime passenger ferry that crosses the bay of Toulon. The ferry to be redesigned is equipped with an aluminum hull, two diesel motors and an electric generator to power the auxiliaries.The ferry can transport 100 passengers. The Toulon ferry runs three lines 7 days a week over 300 days per year. Each ferry makes 24 bay crossings daily. The ferry has a lifetime of 20 years. The diesel motors are replaced approximately every 12500 hours (or about 500 000 km). The design project aims to define the architecture and the state of the systems for each sequence of the use cycle according to pre-established conditions of use and assessment criteria (performance variables, environmental, technical and economic criteria respectively). In the case of intercity passenger transport, ship performances are strongly related to conditions of use.

For practical reasons, in this article we have deliberately simplified the system. The passenger ferry is thus broken down into 5 main sub-sets (shown in Figure 4): Hull & Superstructure, Power Generation, Propulsion andSteering, Energy and Auxiliaries.The main relationships governing the system are shown in the appendices.

Defining negotiable functional parameters

To simplify the problem, we identified 6 Negotiable Functions: maximum cruise speed, maximum passenger capacity, number of daily trips, thermal insulation, air conditioning system and the number of charging / day. In the classical approach the performance of Negotiable Functions are defined in the specifications, each function is caracterised by a nominal performance level P i nom (Table3). The functional unit chosen to compare the various systems and functional scenarios is the number of passengers transported per day. The performance corresponding to the nominal Functional Unit (FU n ) is 2400 passengers transported a day.

Table3 shows for each Negotiable Functions the acceptable performance level ; eg the performance for NF1 "maximum speed" is : 11.5, 12]. The Functional Mix 0 (FM0) is the nominal functional mix (maximum speed 12 knots, 100 passengers per day, 24 missions daily, no insulation and air conditioning system, 1 charge per day).

P 1 ∈ [P 1 min , P 1 max ] = [
For each mix, the CSP model is solved using ILOG solver generating architectures and operation modes minimizing environmental impact with Eco Indicator 99 (EI 99) scores.

For each Functional mix FM j , Table5shows the main elements of Architecures A j (of the simplified model) selected from the component libraries.

The function « objective » is to minimize the environmental impact over the life cycle (Raw Materials + Manufacturing phase, Use phase and Maintenance phase, the End of Life phase is not included).

The diagram (Figure 5) shows the distribution of impacts over the three life cycle phases that are assessed. The predominant phase is the Use phase with 87% of impacts, then the Maintenance phase with 8% and the Raw Materials +Manufacturing phase with 5%.

The eight Functional Mix scenarios described above are assessed environmentally using the indicator of a single EI99 score in order to make the results clearer. It is understood that a multicriteria assessment is recommended in order that the study be robust. We therefore provide a summary with the results of the multicriteria assessment below in annex.

Table5 : Specifications of main elements of the system. Each architecture generated by functional mix described by design of experiment (Table4) is environmentally assessed (Table7). The graph of main effects (Figure 6) illustrates the influence of Negotiable Functions on the Raw materials and Manufacturing impact. Functions from 1 to 5 have little impact on Raw materials and Manufacturing phase. The number of battery charge (NF6), explains largely environmental gains measured for Functional Mix scenarios 1, 2, 5 and 6. In fact, recharging batteries more often, the amount of energy to be stored becomes less important and the need for batteries is less. Batteries, with the superstructure of the ship, are the main contributors to the impacts in Raw materials and Manufacturing phase.

• Assessment of the Use phase

The use phase represents over 85% of the impacts generated by the passenger ferry. Moreover, as shown by the various assessments of scenarios 1 to 7 (Table7), the environmental impact of use is sensitive to the variation of negotiable functionalities.

The best functional mixes allow more than 10% gains. The performance level of NFs leads to a variation of the average environmental impact generated by the eight functional mix scenarios of the DoE of :

 4.5% for NF 3 ,  3.3% for NF 1 ,  2.7% forNF 2 ,  2.5% for NF 4 ,  1.3% for NF 6 ,  and 0.9% for NF 5 .

The analysis of different scenarios associated with a functional mix explains in part the impacts resulting of the use phase. Reducing the maximum speed or the number of missions per day, reduces the fuel consumption. The increase of the number of batteries chargeleads to reducethe mass of batteries and thus reducing the need of propelling poweras system presents less drag. In the same way reducing passenger capacity per mission also reduces the propulsive power. The combination of air conditioning and insulation of system is more difficult to predict, better insulation improves thermal efficiency but weighed down the system increasing fuel consumption while air conditioning increases the energy demand of the system.

 Assessment of the maintenance phase

The ferry undergoes maintenance throughout its entire life cycle: motors are changed (lifetime of 12 500 running hours) and batteries are charged and uncharged (600 cycles for this study). Initially, 5 motors are used over the ship's life time involving 4 changes of motor. For a single charge/uncharge cycle per day, batteries must be replaced 9 times. The batteries are the components generating the highest impact. The number of missions per day leads to a variation of 5% of the average environmental impact of maintenance phase. The use of the system is due to the number of mission it performs by day, which generates wear of these main components.

 Global assessment of scenarios

The overall impact (Raw materials, Use and maintenance) of the various scenarios associated with a Functional Mix of DoE follows the trend drawn by the use phase, in fact 87% of the environmental impacts caused by the use phase.The three significant functions are determined by the functions generating the greatest variation on the overall impact, they are characterized by observation on the graph of main effects (Figure 6). The performance level of NF1, NF2 and NF3 leads, respectively, to a variation of 3.3%, 2.7% and 4.7% of the average environmental impact generated by the eight functional mix scenarios of the DoE.

In the following, the three Significant Negotiable Functions are modeled as constraints and the three non significant are set at their nominal values. The performance level of NF is identified by CSP.

The performance of three SNF to optimize the environmental performance of the system is characterized by a 11.5 knots as maximum speed, 98 passengers as a maximum capacity and 23 crossings per day (Table8). The main components of the system satisfying the requirements of functional mix globally optimized (FM GO ) are defined Table9. The Functional Unit Globally Optimized(FU GO ) is 2300 passengers per day. This corresponds to a reduction of about 5% of the nominal performance. However, in line with the trend observed on the effects of SNF (Figure 6), environmental impacts compared to the reference system are reduced by approximately 13%. Table9 : Specifications of main elements for system with FMGO.

Reducing the system's performances or eliminating certain functions raises the question of outcomes for the passenger. For example, in this type of intercity transport, the number of passengers is not constant throughout the day. It fluctuates, and there are more people during rush hours. In the above simulations, the environmental gain is achieved to the detriment of « social » considerations; this is in contradiction to the concept of sustainable development. Reducing the amount of space on the ship results in constraints for the user. In parallel with initiatives to define a coherent functional mix (that we suggest with the EcoCSP tool), it is therefore necessary to set up measures to make sure that the system's ecological performances are not achieved at the expense of the users experience. This means, for example, setting up incentives to obtain a more regular flow of passengers throughout the day, such as preferential tariffs for certain time bands etc.

4.5. Influence of the exhaustiveness of technological solutions on environmental performance.

In this section, we propose to modify the component database "Engine". The database used to model the system comes from the manufacturers catalogs. In Figure 7the speed of the nominal functional mix is reduced by 0.5 kts.

In the first case the database "Engine" (DB 0 ) is not modified. In the second case, a new engine of 70kW is added to the database (DB 1 ) (Figure 7).The integration of a 70kW diesel engine in the component database of CSP model has allowed an environmental gain of 8% compared to the use of the initial database (DB 0 ) using diesel engines from manufacturers catalogs. More the alternatives are important more opportunities to generate better solutions for environment are high.

Functional negotiations must lead a questioning about the components to use in the system, which could otherwise limit the benefits of environmental performance.

Discussions and Conclusions

The CSP / LCA approach relies on the CSP solver makes it possible to instantiate the design variables of the system that optimize the performance variables, in our case it is environmental criteria (but it could be also economic criteria). It allows to identify early in design process what are the best combinations of technologies, among many alternatives, for each subsystem whose functions are already specified in the specifications. The Philosophy behind the use of CSP is to allow the designer to model and identify viable concepts reconciling environmental and economic aspects (CSP / LCA), but also the social aspects by acting on the definition of functional performance to get closer to a model of sustainable development.

In this article we have therefore enriched the CSP/LCA approach by constructing the EcoCSP approach. This enables us to anticipate the configuration of a system's architecture by adapting the performances of negotiable functions. A complex system such as a passenger ferry has numerous sub-functions. A slight downgrading of the performances related to these functions can generate substantial environmental gains. The complexity of couplings among sub-systems and their sheer number obliges the user to make use of « intelligent » tools, that by simulating many different scenarios, help the designer to fine-tune and choose the right technologies for sustainable systems. This EcoCSP approach breaks with traditional design conventions, and allows us to define firstly combinations of technologies and secondly the functional mix that will significantly reduce environmental impacts. Life Cycle Analysis (LCA) is thus no longer used as a tool for system assessment and comparison, but as a tool for ecodesign.

In the experimentation phase we showed the significance of the number of alternatives suggested by the CSP model for reducing environmental impacts. The greater the number of alternative techniques, the higher the number of possibilities for generating better environmental solutions. Finally, modifying the functional performances of a system results in new social, economic and environmental constraints. In parallel therefore, it is necessary to reflect on all the consequences of such modifications in order to avoid destabilizing the three mainstays of sustainable development.

  xn}, a set of variables, n being the number of variables of the problem. To keep the generic element, we say that these variables may relate to design, performance or state. Design parameters structure the design and their values distinguish between two design configurations. The instantiation of all the design parameters defines the complete potential design solution. Performance parameters translate the state or the quality of a design alternative and compare it to a reference from the specifications or one related to the state of the art of the company or sector concerned. These characteristics are linked to the translation of a given configuration in physical terms and are generally directly linked to the design parameters. D = {d1, d2, d3, …, dn}, a set of domains. Each domain, associated to a variable, can be discrete or continuous. C = {c1, c2, c3, …, cp}, a set of constraints, p being the number of constraints of the problem. The constraints translate how the structuring functions are carried out by the system during the life situation in question. The constraints take the form of explicit relationships between several variables. These relationships impose restrictions on the domains of possible values for the variables of the problem. More precisely, it can be a logical combination of several elementary constraints, among the following:

Figure 1 :

 1 Figure 1: Process of CSP solution.

Figure 2 :

 2 Figure 2: Flowchart of the EcoCSP approach.

Figure 4 :

 4 Figure 4: Model of hybrid passenger ferry.

Figure 5 :

 5 Figure 5 : Average distribution of impacts on life cycle phases.

Figure 6 :

 6 Figure 6 : Main effect of Negotiable Function on Environmental Impact.

  Table3 : Definition of nominal performence Pi nom and performance level [Pi min, Pi max] of NFi.

								Parameters
		P1: Maximal	P2: Number of	P3: Number			P5:	P6 : Nombre	Ps : Unité
	Level	Speed	passengers	of missions	P4: isolation	Climatisation	de charges	fonctionnelle
	Pi nom		12		100	24 1 (insulation)	1 (AC)	1	2400
	Pimin		11.5		97	23 0 (no insulation) 0 (No AC)	1	2231
	Pimax		12		100	24 1		1	12	2400
	4.3. Defining significant negotiable functional parameters
	They are 2 6 possible functional mix for 6 Negotiable Functions with two performance levels. In order to reduce
	the number of CSP simulation, the appropriate Taguchi orthogonal array is L 8 (2 6 ) is implemented Table4,
	representing the main effects of each NF on the environmental impact.
	Functional Mix	Performance	Performance	Performance	Performance	Performance	Performance
	(FM)		level forNF1	level for NF2	level for NF3	level for NF4	level for NF5	level for NF6
	FM0		12		100	24		0	0	1
	FM1		12		100	24		1	1	12
	FM2		12		97	23		0	0	12
	FM3		12		97	23		1	1	1
	FM4		11.5		100	23		0	1	1
	FM5		11.5		100	23		1	0	12
	FM6		11.5		97	24		0	1	12
	FM7		11.5		97	24		1	0	1

Table4 : Design of Experiment -Taguchi orthogonal array L8 (2 6 ).

  Table6 : Life Cycle Assessment of Architecture Aj.

								Architectures
				Functional Mix	Functional Mix	Functional Mix	Functional Mix	Functional Mix	Functional Mix	Functional Mix	Functional Mix
				0		1	2	3	4		5	6	7
		Power		80 000 W	80 000 W	70 000 W	80 000 W	70 000 W	70 000 W	70 000 W	70 000 W
	Engine
		Weight		800 Kg	800 Kg	750 Kg	800 Kg	750 Kg		750 Kg	750 Kg	750 Kg
	•	Power Assessment of the Raw Materials + Manufacturing phase 24 000 W 24 000 W 22 000 W 24 000 W	24 000 W	24 000 W	22 000 W	24 000 W
	Motor
		Weight		94 Kg	94 Kg	85 Kg	82 Kg	94 Kg		94 Kg	70 Kg	94 Kg
		Power		10000	10000	10000	10000	10000		10000	10000	10000
	Generator
		Weight		200	200	200	200	250		250	250	250
		Energy mission (Winter/ Summer)	15 kWh/ 12,4 13,7 kWh/ 17,7 14,6 kWh/ 11,9 13,7 kWh/ 17,7 15,3 Wh/18,0 14,0 Wh/12,6 14,8 Wh/17,5 14,0 Wh/12,6
	Batteries	Weight		501 Kg	53 Kg	42 Kg	565 Kg	577 Kg		40 Kg	53 Kg	464 Kg
		Energy Density	50 Wh/Kg	50 Wh/Kg	50 Wh/Kg	50 Wh/Kg	50 Wh/Kg	50 Wh/Kg	50 Wh/Kg	50 Wh/Kg
	Hull	Weight		15 000 Kg	15 000 Kg	15 000 Kg	15 000 Kg	15 000 Kg	15 000 Kg	15 000 Kg	15 000 Kg
		Weight		27 539 kg	27 648 Kg	25 504 Kg	27 769 Kg	27 587 Pts	27 520 Kg	25 526 Kg	27 462 Kg
	Ship	Impact of Use phase / Mission Phase Impact of Manufacturing	2,03 Pts 16 325 Pts	1,82 Pts 13 901 Pts	1,79 Pts 13 562 Pts	2,04 Pts 16 670 Pts	1,82 Pts 16 690 Pts	1,87 Pts 13 789 Pts	1,67 Pts 13 848 Pts	1,81 Pts 16 079 Pts
		Impact of Maintenance	27134 Pts	31250 Pts	24403 Pts	30244 Pts	30804 Pts	32177 Pts	30839 Pts	25134 Pts
						8%
							5%
								Life Cycle Phases
								Raw Materials +
								Manufacturing
								Use
								Maintenance
							87%
	Functional Mix (FM)	Raw Material + Manufacturing	Use	Maintenance	Life Cycle	Performance level PSforFU
	FM 0		16325	318077	27134	361536	2400
	FM 1		13901	331917	31250	377068	2231
	FM 2		13562	259115	24403	297080	2231
	FM 3		16670	296337	30244	343251	2400
	FM 4		16690	268935	30804	316428	2400
	FM 5		13789	269761	32177	315727	2231
	FM 6		13848	270843	30839	315530	2231
	FM 7		16079	285645	25134	326858	2400
	Average Impact 15108	287579	27936	330623

Table7 : Impact of architectures for functional mix with Design of Experiment L8 (2 6 ).

Table :

 : Multicriteria assessment according CML Method 1 .

						Diesel Engine DB 1
						Type	Power (kW)	Weight (kg)
		Diesel Engine DB 0		1	60	550
		Type	Power (kW)	Weight (kg)	2	70	750
		1		60	550	2	80	800
		2		80	800	4	90	850
		3		100	900	3	100	900
			Functional Mix normalized	Functional Mix globally optimized/ DB0	Functional Mix globally optimized/ DB1
		Power	80 000 W	80 000 W		70 000 W
	Engine
		Weight	800 Kg	800 Kg		750 Kg
		Power	24 000 W	22 000 W		22 000 W
	Motor
		Weight	94 Kg	85 Kg		85 Kg
		Power	10000	10000		10000
	Generator
		Weight	200		200		200
		Energy mission	15 kWh/ 12,4	15,2 Wh/12,4		15,2 Wh/12,6
	Batteries	Weight	501 Kg	488 Kg		460 Kg
		Energy Density	50 Wh/Kg	50 Wh/Kg		50 Wh/Kg
	Hull	Weight	15 000 Kg	15 000 Kg		15 000 Kg
		Impact of Use phase (Pts)	292320	279460		251160
		Impact of Manufacturing Phase (Pts)	16 325	16 079		16 125
	Ship	Impact of Maintenance (Pts)	27134,0	25103,9		25250,8
		Global	335 779	320 643 (-4.7%)		292 536 (-12,9%)
		Functional Unit	2400		2300		2300

1

Guinée, J.B., Gorrée, M., Heijungs, R., Huppes, G., et al. 2001

. "CML-Guide to Life Cycle Assessment." Institute of Environmental Sciences (CML), Leiden University, NL

The EcoCSP tool allows us to make functional judgments and choices to optimize negotiable functional performances and thereby reduce environmental impacts. Nevertheless, this does not mean we should neglect consideration of the social consequences that these choices have on the system's use.

Abiotic_depletion Eutrophisation Acidification

GWP Eco


Shi-Qing Xin

Bruno Lévy

Zhonggui Chen
email: chenzhonggui@xmu.edu.cn

Yaohui Yu

Changhe Tu

Wenping Wang

Centroidal Power Diagrams with Capacity Constraints: Computation, Applications, and Extension

Keywords: centroidal power diagram, displacement interpolation, convex decomposition, blue noise Concepts:, Computing methodologies → Computer graphics;

This article presents a new method to optimally partition a geometric domain with capacity constraints on the partitioned regions. It is an important problem in many fields, ranging from engineering to economics. It is known that a capacity-constrained partition can be obtained as a power diagram with the squared L2 metric. We present a method with super-linear convergence for computing optimal partition with capacity constraints that outperforms the state-of-the-art in an order of magnitude. We demonstrate the efficiency of our method in the context of three different applications in computer graphics and geometric processing: displacement interpolation of function distribution, blue-noise point sampling, and optimal convex decomposition of 2D domains. Furthermore, the proposed method is extended to capacity-constrained optimal partition with respect to general cost functions beyond the squared Euclidean distance.

Introduction

The central theme of the present paper is to apply and extend the notion of power diagram to computing optimal partition of a geometric domain with capacity constraints. A partition of a geometric domain is a collection of disjoint constituent parts whose union is equal to the whole domain. A partition is said to be optimal if all its partitioned regions, also called cells, have compact shapes with respect to some specified cost kernel that defines the cost function.

Here "compactness" can be roughly measured by the isoperimetric quotient. Capacity constraints on a partition mean that each cell of the partition has a specified size, measured in volume (or area in 2D) or mass. Hence, an optimal partition with capacity constraints is a partition such that its cells have compact shapes as well as sizes equal to some specified quantities, unlike the optimal mass transport [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]] that concentrates on capacity constraints regardless of shape compactness or the centroidal Voronoi tessellation [START_REF] Du | Centroidal Voronoi tessellations: applications and algorithms[END_REF][START_REF] Liu | On centroidal Voronoi tessellation: energy smoothness and fast computation[END_REF]] that considers the opposite.

Finding an optimal partition by some criteria is a central issue in many applications. Optimal partitions, with or without capacity constraints, are used in computer graphics and geometric computing for point sampling [START_REF] De Goes | Blue noise through optimal transport[END_REF], information visualization [START_REF] Balzer | Voronoi treemaps for the visualization of software metrics[END_REF], shape segmentation [START_REF] Lu | Variational 3D shape segmentation for bounding volume computation[END_REF]], shape matching [START_REF] Su | Optimal mass transport for shape matching and comparison[END_REF], area-preservation mapping [START_REF] Zhao | Area-preservation mapping using optimal mass transport[END_REF]] and mesh generation [START_REF] Du | Centroidal Voronoi tessellations: applications and algorithms[END_REF][START_REF] Liu | On centroidal Voronoi tessellation: energy smoothness and fast computation[END_REF]]. The optimal partitioning problem with capacity constraints is also studied in engineering and economics for optimal resource location or provision of services. For example, it appears as the optimal coverage problem in robotics and sensor networks [Cortés 2010;[START_REF] Parke | Centroidal areaconstrained partitioning for robotic networks[END_REF], the pattern formation problem in material science [Bourne et al. 2014], and the land consolidation problem in agriculture [START_REF] Borgwardt | Geometric clustering for the consolidation of farmland and woodland[END_REF].

The contribution of this paper is twofold.

1. An L-BFGS method with empirically superlinear convergence is developed for computing capacity-constrained centroidal power diagram, assuming the squared Euclidean distance as the cost kernel. We show its superiority to previous method and validate it in two applications: displacement interpolation of function distribution and blue-noise point sampling. Furthermore, the method is extended to computing optimal capacity-constrained partition with a general cost function;

2. Based on the centroidal power diagram, an effective method is developed for decomposing a 2D polygonal region into convex and compact cells with capacity constraints.

Related Work Voronoi Diagram & Power Diagram

The Voronoi diagram is a partition of a given domain Ω ⊂ E d into regions based on the straight-line distance to a set of points {xi ∈ Ω} n i=1 (called sites or generators), where xi dominates the subregion

Ωi : {x ∈ Ω | x -xi ≤ x -xj , j = i},
which is called the Voronoi cell of the site xi. A power diagram is defined as an extension of Voronoi diagrams, where each site xi has an associated coefficient (called a "weight") and dominates the subregion

Ωi : {x ∈ Ω | x -xi 2 -wi ≤ x -xj 2 -wj, j = i}.
Ωi is called the power cell of the site xi. The power diagram reduces to the Voronoi diagram when all the weights are equal [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF][START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF].

Centroidal Voronoi Tessellation (CVT)

A centroidal Voronoi tessellation (CVT) of a set of sites X = {xi} n i=1 is a special Voronoi diagram for which each site xi coincides with the centroid (i.e. the center of mass) of its Voronoi cell Ωi [START_REF] Du | Centroidal Voronoi tessellations: applications and algorithms[END_REF]]. Alternatively, a CVT can be also defined as a stationary point of the following objective function:

Q(X) = n i=1 Ω i ||x -xi|| 2 ρ(x)dσ, (1)
where ρ(x) is a C 1 -smooth density function on the whole domain Ω. Domain partition using CVT yields convex polyhedral cells (barring those abutting boundaries) that have compact shapes and similar sizes, assuming that ρ(•) is uniform. Due to these properties, CVT is widely used for computing optimal partitions in various applications [START_REF] Lu | Build-to-last: strength to weight 3D printed objects[END_REF][START_REF] Zhou | Multiscale superpixels and supervoxels based on hierarchical edge-weighted centroidal voronoi tessellation[END_REF], and in some cases even for optimal partition with capacity constraints [START_REF] Balzer | Capacity-constrained Voronoi diagrams in finite spaces[END_REF][START_REF] Chen | Variational blue noise sampling[END_REF][START_REF] Li | Fast capacity constrained Voronoi tessellation[END_REF]]. However, it is known that the class of partitions given by the Voronoi diagram is too restrictive to properly accommodate capacity constraints [START_REF] Zhang | Capacity constrained blue-noise sampling on surfaces[END_REF].

A prevailing method for computing CVT is Lloyd's method, which behaves as a gradient-descent method and thus converges linearly.

A quasi-Newton method for computing CVT is proposed in [START_REF] Liu | On centroidal Voronoi tessellation: energy smoothness and fast computation[END_REF]] that is based on the limited-memory BFGS technique and is shown to demonstrate super-linear convergence in practice.

Optimal Mass Transport (OMT)

The optimal mass transport (OMT) was formalized by the French mathematician Gaspard Monge [1781] who studied the most economical way of moving soil from one area to the other and introduced the theory in terms of minimizing the L1 norm of the distance transported. The problem of minimizing L2 cost whose kernel is the squared distance (note that the square root of the minimum L2 cost is called the 2-Wasserstein metric), however, takes on additional structure and can be expressed as the fully nonlinear degenerate elliptic Monge-Ampère equation. In 1949, Kantorovich proved the existence and uniqueness of the optimal transport plan and showed how Monge's problem is connected to linear programming [START_REF] Kantorovich | Application of mathematical methods to problems of analysis of freight flows[END_REF]. The Monge-Kantorovich optimization technique has been studied in a wide range of research fields in the past two decades, including imaging, adaptive meshing, geophysical fluid dynamics, and cosmology.

It is shown in [START_REF] Pogorelov | Generalized solutions of Monge-Ampère equations of elliptic type. In A tribute to Ilya Bakelman[END_REF]], [START_REF] Gangbo | The geometry of optimal transportation[END_REF] (example 1.6) and [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF]] (or as a direct consequence of [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF]) that for the semi-discrete OMT (from a discrete probability distribution function to a continuous one), the class of partitions parameterized by the power diagram provides the proper representation to accommodate capacity constraints in the following sense. Given a set of fixed sites in a domain, there exists a unique set of weights (up to an additive constant) such that the power cells, which are determined by the given sites and the associated weights, satisfy the specified capacity constraints and minimize the sum of least squares costs of all the cells. Here, the least squares cost of the cell Ωi associated with a site xi is defined to be

Ω i x -xi 2 ρ(x)dσ.
This means that we just need to look for the capacity-constrained partition within the set of power diagrams with different weighting schemes, rather than having to search in the set of arbitrary partitions of the domain, which obviously cannot be parameterized by a finite number of parameters.

Centroidal Power Diagram (CPD)

In this paper, we are interested in computing the optimal partition of a geometric domain with specified capacity constraints. Similar to the conventional OMT problem, the solution to such a capacity-constrained optimal partition problem with the least squares cost function can be represented by a power diagram [START_REF] Balzer | Capacity-constrained point distributions: A variant of Lloyd's method[END_REF]]. More specially, the solution is given by a centroidal power diagram (CPD) which is a special power diagram in which each site coincides with the centroid of its power cell. Interestingly, the weights of the power diagram exactly correspond to the degrees of freedom required to account for capacity constraints. However, one should notice that the relation between the weights and the capacities of the cells is non-trivial (in this context, the term "weight" is misleading, they are completely different from the masses of the cells).

CPD was used for blue-noise sampling in computer graphics. In numerical implementation, [START_REF] Balzer | Capacity-constrained point distributions: A variant of Lloyd's method[END_REF] proposed several iterative schemes to compute capacity-constrained optimal partitions, including the method based on pixel-based clustering and the methods that enforce capacity constraints on the Voronoi diagram as well as the power diagram. These methods are easy to implement. But because they update only one single pixel or site in each iteration, the methods have very slow convergence and low efficiency. Later, a Lloyd-like method that repeatedly moves each site to the mass center of its cell is proposed to compute CPD for blue-noise sampling [START_REF] De Goes | Blue noise through optimal transport[END_REF]], which has only linear convergence.

The above discussions naturally lead to two questions. First, is it possible to develop an efficient method that outperforms these existing methods? Second, remarking that the least squares cost function is in fact the squared L2 metric, we want to know whether the new numerical method is able to deal with the CPD problem with a general cost. We will give positive answers to both questions in this paper.

Capacity-Constrained Optimal Partition

Problem Formulation

We first recall the notions used for a power diagram in E d , ddimensional Euclidean space. Let Ω ⊂ E d be a convex, closed, bounded and connected domain. Let X = {xi} n i=1 ⊂ Ω be a set of points, called sites, with its associated weights {wi} n i=1 ⊂ R. Then the power diagram induced by the weighted points (xi, wi) is the collection of convex polytopes, or cells, P D = {P Di} n i=1 , where P Di = {x| x -xi 2 -wi ≤ x -xj 2 -wj, ∀j, j = i}. The restriction of the power diagram P D to the domain Ω induces a partition of Ω, denoted by P DΩ = {P Di ∩ Ω} n i=1 .

Now consider partitioning a domain Ω in Euclidean space E d into a set of n regions R = {Ri} n i=1 in a capacity-constrained manner. Let X = {xi} n i=1 ⊂ Ω denote n given points, also called sites, with associated capacity constraints ci > 0, i = 1, 2, . . . , n. We define the least squares cost of each region Ri, given by xi, to be R i

x -xi 2 ρ(x)dσ, where ρ(x) is a C 1 -smooth density function on Ω. The capacity constraint here requires that |Ri| ≡ R i ρ(x)dσ = ci, i = 1, 2, ..., n, that is, the capacity (i.e. area or volume) of each region Ri matches the constraint ci. (Of course, i ci is assumed to be equal to |Ω| ≡ Ω ρ(x)dσ.) Then the problem of capacity-constrained optimal partition is to find the partitions Ri and the sites xi such that the total cost

Q(X, R) = n i=1 R i ||x -xi|| 2 ρ(x)dσ (2)
is minimized, subject to the constraints |Ri| = ci, i = 1, 2, ..., n.

In the case of the L2 cost, it follows from [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF]] that the solution (X * , R * ) = arg min Q(X, R) exists and is given by a power diagram; See [START_REF] Pogorelov | Generalized solutions of Monge-Ampère equations of elliptic type. In A tribute to Ilya Bakelman[END_REF]] and [START_REF] Gangbo | The geometry of optimal transportation[END_REF] for a proof, as well as the derivations in [START_REF] Gu | Variational principles for Minkowski type problems, discrete optimal transport, and discrete Monge-Ampère equations[END_REF]. [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF] further showed that the optimal power diagram can be found by extremizing

F (X, W ) = n i=1 P D(W ) i ||x-xi|| 2 ρ(x)dσ- n i=1 wi(|P D(W )i|-ci),
(3) where P D(W )i (or P Di) is a power cell dominated by xi. A key observation made in [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF]] is that, when the sites X in F (X, W ) are fixed, their weights can be computed as a maximizer of a concave function to make the cells of the resulting power diagram meet the capacity constraints |P Di| = ci, i = 1, 2, ..., n. This observation leads to a Lloyd-like method [Cortés 2010;[START_REF] De Goes | Blue noise through optimal transport[END_REF] for minimizing the objective function F (X, W ) in Eqn. (3), i.e., alternatively optimizing the weights to meet the given capacity constraints and moving the sites to their respective mass centers. We will discuss these numerical techniques in detail in the next section, and introduce a new method that better exploits the derivatives of the objective function, resulting in a faster convergence. Let µ denote a measure with non-vanishing density ρ(.), supported by the domain Ω. Suppose that each site xi has capacity constraint ci > 0, such that n i ci = |Ω| ≡ Ω dµ, which is the total mass of the measure µ. Again, let d(x, y) denote a distance cost between any two points x and y in Ω. Let R = {Ri} n i=1 denote a partition of Ω, that is, ∪Ri = Ω and Ri ∩ Rj = ∅, ∀i = j. Define the cost of the partition R with respect to the sites X = {xi} n i=1 as

Partition with General Cost Functions

Q(X, R) = n i=1 R i d(x, xi)dµ. (4)
In what follows, we support that the distance d(x, y) satisfies the twist condition, i.e. ∇xd(x0, y) is one-to-one ∀x0, which ensures the existence and uniqueness of the optimal transport plan (see e.g. [McCann and Guillen 2010] p. 10 and [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]). Then the capacity-constrained optimal partition problem is to find the sites X * and partition R * that minimize the function Q(X, R) subject to the capacity constraints ci = |Ri| ≡ R i dµ, i = 1, 2, . . . , n. The difference between this definition of the general setting and the previous simple setting of the squared L2 cost is that here Ω may be a curved surface and that d(x, y) is a general cost kernel.

Next we shall show that the optimal partition R * with capacity constraints assuming a general cost function Q(X, R)) can be found by only considering the candidates of extended power diagrams where each region (possibly with curved boundaries) is associated with a weight. This is a natural extension from [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF][START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]]. Its proof is given in the appendix.

Theorem 1: For any fixed sites X = {xi} n i=1 , the capacityconstrained partition R = {Ri} n i=1 of the domain Ω that minimizes the cost function in Eqn. ( 4) is the one such that the boundary of any two adjacent regions Ra and R b , which are associated with the sites a and b in X respectively, is defined by the equation d(x, a) -wa = d(x, b) -wa for some constants wa and w b .

Similar to the treatment in [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF]] for the case of the least squares cost function, we build a new objective function to accommodate the optimal partition according to the Lagrangian Multiplier Theorem:

F (X, W ) = n i=1 GD(W ) i d(x, xi)dµ - n i=1 wi(|GD(W )i| -ci).
(5) It can be shown that F (X, W ) is concave w.r.t. W and therefore has a unique maximizer (up to an additional constant), denoted

W * = {w * i } n i=1 .
When the sites X = {xi} n i=1 are fixed, the optimal partition exists uniquely and is given by the generalized Voronoi diagram GD(W * ). When the capacity constraints are met, F (X, W ) is exactly equal to the cost function in Eqn. (4). Furthermore, the gradient of F (X, W ) w.r.t. W can be shown to be [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF][START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]]

∇W F (X, W ) = -(|GD(W )1|-c1, . . . , |GD(W )n|-cn), (6)
which is a normal vector to the tangent plane of F (X, W ). Hence, to seek the maximizer W * of F (X, W ), we have ∇W F (X, W * ) = 0, which means that all the capacity constraints are satisfied. To summarize, we optimize over all the sites X and weights W to find the generalized weight Voronoi diagram that gives an optimal partition with capacity constraints. Because there is a unique set of weights W (up to an additional constant) for every set of sites X to make the resulting weighted Voronoi diagram satisfy the capacity constraints, all the capacity-constrained partitions can be represented by a surface W = W * (X), which will be called the constraint surface (see the green curve in

Algorithm for Computing Centroidal Power Diagrams

In this section we propose an efficient method for computing CPDs, as well as the implementation details.

Overall Algorithm

As discussed above, if (X * , W * ) defines the optimal partition, the cost function F (X, W ) (see Eqn. ( 5)) is minimized by X * while maximized by W * . In other words, finding (X * , W * ) is neither a pure maximization problem nor a minimization problem. de Goes et al. [2012] suggested updating X and W in an alternative style, where the two key operations are: (1) finding the optimal weighting scheme W to meet the given capacity constraints using Newton's method while fixing X [START_REF] Benamou | Numerical solution of the optimal transportation problem using the Monge-Ampère equation[END_REF]; and (2) moving each site xi to the mass center of its corresponding power cell. However, this iterative scheme has only linear convergence due to the ineffectiveness of the latter operation. As de Goes et al. [2012] observed, developing an efficient CPD algorithm with super-linear convergence is non-trivial.

Based on the discussion in Section 3, the optimal weights W * maximizing the cost function in Eqn. ( 5) is uniquely determined (up to an additional constant) for a given site collection X. Therefore the cost function can be reconsidered as a function of X :

Minimize F (X, W * (X)) = n i=1 GD(W * ) i d(x, xi)dµ - n i=1 w * i (|GD(W * )i| -ci), (7)
where W * depends on X so that the specified capacity constraints are met. Based on this formulation, we will use the combination of the reduced-gradient technique and the L-BFGS solver [START_REF] Lasdon | Nonlinear optimization using the generalized reduced gradient method[END_REF][START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF] to move (X, W ) on the constraint surface W = W * (X) towards the minimum of F (X, W * (X)).

Our method works as follows. We start from a feasible point (X0, W0), i.e. a point on the constraint surface W0 = W * (X0).

Without loss of generality, we assume that (X k , W k ), the position in the k-th iteration, is a feasible point as well. Next, we compute the reduced-gradient vector ∇x i F (X, W * (X)) and feed it to the LBFGS solver to update the sites X k to X k+1 . Note that the weights are computed according to the up-to-date sites before we evaluate the objective function F . Thus the (k + 1)-th iteration yields the next position (X k+1 , W k+1 ). Repeat the process until some stopping criteria is met.

We now consider the derivation of ∇x i F (X, W * (X)). Based on the envelope theorem [START_REF] Milgrom | Envelope Theorems for Arbitrary Choice Sets[END_REF], we know

∇x i F (X, W * (X)) = ∇x i F (X, W ). ( 8
)
According to Reynolds' transport theorem [START_REF] Leal | Advanced transport phenomena: fluid mechanics and convective transport processes[END_REF]] that recasts derivatives of integrated quantities, we can further show that when the cost kernel d is differentiable everywhere in the domain, the gradients of F (X, W * (X)) also have a simple form:

∇x i F (X, W * (X)) = GD(W * ) i ∇x i d(x, xi)dµ. (9)
The pseudocode of our method is given in Algorithm 1. It is worth noting that for general cost kernel that is different from the 2-Wasserstein metric whose kernel is the squared distance, instead of using Newton's method, we use the L-BFGS method to update W , because the Hessian matrix is much more complicated in this case than that in the case.

X k+1 ← X k 11 k ← k + 1 12 until ||∇X F (X, W * (X))|| < 13 Output: (X k , W k )
Algorithm 1: Computing CPDs using L-BFGS.

Optimization Options

The computation of W * is crucial to accurate computation of ∇x i F (X, W * (X)). As de Goes et al. [2012] observed, it usually takes 3 to 5 iterations to bring the residual of capacity constraints to be within an accuracy of 10 -12 . Therefore we require ||∇W F (X, W )|| < 10 -12 during the step of updating W . In implementation, rather than always resetting W to 0 at the very beginning of Newton's method, it is better to use the previous values of W to initialize Newton's method. This is helpful for reducing the computation cost especially when ||∇W F (X, W )|| becomes very small. We observed an improvement of about 10% due to this "trick" in our experiments. When empty cells occur, we reduce the step size by half repeatedly until each site dominates a nonempty cell, as in [START_REF] De Goes | Blue noise through optimal transport[END_REF][START_REF] Zhao | Area-preservation mapping using optimal mass transport[END_REF]]. Note that the convergence of this Newton step-control scheme was recently proved [START_REF] Kitagawa | A Newton algorithm for semi-discrete optimal transport[END_REF]].

For the 2-Wasserstein metric, the optimal partition is given by power diagrams. In this case, the Hessian matrix of F (X, W ) w.r.t. W has an elegant form [START_REF] De Goes | Blue noise through optimal transport[END_REF]. But for a general cost kernel different from the squared Euclidean distance, the coefficients of the Hessian matrix (gradient of the distance function integrated over the facets of the power diagram) are much more complicated, and in most cases not known in closed form. For this reason we use the L-BFGS method that only needs the first-order derivatives (that always correspond to the difference between the capacity constraints and the measures of the Laguerre cells).

Numerical Integration and Boundary Treatment

When the density function ρ is constant and the cost kernel is the 2-Wasserstein metric, the integrand in Eqn. ( 7) has an explicit representation by using the Green formula [START_REF] Dobrovolsky | Sur l'histoire de la classification des points singuliers des équations différentielles[END_REF]] to transform the computation of the cost function, as well as its gradients, into contour integration. But in general cases, the cost function cannot be reduced to an explicit representation and thus discretization of the domain, typically a sufficiently refined triangulation, is required. Therefore, at the very beginning of the algorithm, we randomly distribute a set of points inside the domain and optimize their positions [START_REF] Liu | On centroidal Voronoi tessellation: energy smoothness and fast computation[END_REF]]. After that, we project the outer-ring points (whose Voronoi cells intersect the boundary) onto the boundary and re-optimize the other points. Finally we compute the constrained Delaunay triangulation w.r.t. these points as discretization of the domain; See the inset figure .   x

i x i x j x j x k p 1 v 1 p p 1 p 2 v 2 v 3 v 3 v 2 v 1 p 3 p 2
Besides, the optimal partition of the domain for a general cost kernel d(•, •) may not be the power diagram. We have to develop an effective technique to build the partition with the help of the base triangulation. First of all, we take the sites in X as sources and compute the weighted discrete distance field using the fast sweeping technique so that each vertex v of the base triangulation finds its nearest site xi in X, i. Another implementation detail we need to mention is boundary treatment. During the optimization process, some sites may jump out of the domain Ω and we need to pull them back to ∂Ω or the interior side of Ω. Suppose that in the j-th iteration, the i-th site x (j) i

is outside Ω and we wish to pull x (j) i back. We find the intersection point p between ∂Ω and the connecting line segment of x (j) i and its previous position x (j-1) i , and take p as the corrected position of x (j) i . Note that the horizontal axis in (g-l) denotes the timing (s).

Performance Evaluation

example of n = 500 sites (see the middle column), the density function at x = (x (1) , x (2) ) is set to ρ(x) = 0.1 + x (1) . For the example of n = 1, 000 sites (see the right column), the density is set to ρ(x) = e -8(x (1) -0.5) 2 -8(x (2) -0.5) 2 . Fig. 2(j-l) shows that our algorithm exhibits a faster convergence. In detail, we observe that (1) the resultant cells by our algorithm are of more compact shapes than [START_REF] De Goes | Blue noise through optimal transport[END_REF] after the same number of iterations; (2) the cost by our algorithm reaches the stable state earlier than [START_REF] De Goes | Blue noise through optimal transport[END_REF]; and (3) the gradient of the cost function decreases more sharply. Note that the horizontal axis of Fig. 2(g-l) denotes the timing.

Since the operation of building power cells is the most timeconsuming part, Table 2 shows the number of building power cells to achieve the same specified threshold 10 -8 , as well as the total computation time. From the statistics we can see that our reduced gradient method has empirical super-linear convergence and significantly outperforms [START_REF] De Goes | Blue noise through optimal transport[END_REF].

General Cost Functions

Many real-life application occasions can be formulated as a CPD problem. Consider deploying some service centers (e.g. postal offices or medical centers) to serve the residents of a city, where each center has a capacity constraint. It's better to take the Euclidean distance, rather than the squared distance, as the service cost kernel. Another example is when the domain is a curved surface where we often take geodesic distances as the cost kernel. As pointed out in Section 3, our algorithmic framework supports general kernels. Here we use six typical kernels (see Fig. 3) to exhibit the extensibility of our algorithm, including

d(x, xi) = (x (1) -x (1) i ) 2 + (x (2) -x (2) i ) 2 in R 2 ; d(x, xi) = (x (1) -x (1) i ) 2 + (x (2) -x (2) i ) 2 in R 2 ; d(x, xi) = (x (1) -x (1) i ) 4 + (x (2) -x (2) i ) 4 in R 2 ;
Geodesic distance on mesh surfaces; Squared geodesic distance on mesh surfaces.

Cubic geodesic distance on mesh surfaces.

In Fig. 3, the first row shows that CPD has freedom degrees to support user-specified capacity constraints while CVT cannot control cell sizes. The second row gives three CPDs with different cost kernels assuming the same uniform density. The third row shows the case of a Gaussian density function. The last row exhibits three CPDs on curved surfaces. When handling mesh surfaces, each site xi is constrained to move along the tangent plane and projected onto the surface before evaluating the cost function. Furthermore, to find the optimal CPD partition, we have to frequently query the geodesic distances [START_REF] Surazhsky | Fast exact and approximate geodesics on meshes[END_REF][START_REF] Xin | Improving Chen and Han's algorithm on the discrete geodesic problem[END_REF] between two points, which is computationally expensive. Due to this limitation, each of the input meshes in Fig. 3(i-k) contains about 1K triangle faces in our experiments. With the termination condition ∇W F (X) being less than 10 -4 , we need about 15 iterations and 5 minutes to deal with each example. (Note that the partition of the surface is approximate and the evaluation of the cost function is approximate as well. That's why we set the tolerance to be 10 -4 , rather than 10 -8 .) However, this can be greatly accelerated using several recent methods for geodesic computation [START_REF] Crane | Geodesics in heat: A new approach to computing distance based on heat flow[END_REF]Ying et al. 2013b;Ying et al. 2013a]. With a faster geodesic computation querier, we believe that it has a great potential for many geometry processing tasks.

Applications

In this section, we shall use it for three problems: displacement interpolation of function distributions, blue-noise point sampling, and optimal convex decomposition of 2D domains.

Displacement Interpolation

The traditional Wasserstein barycenter problem, i.e. displacement interpolation, provides a generic method for interpolating between distributions or functions based on advection instead of blending [START_REF] Bonneel | Displacement interpolation using Lagrangian mass transport[END_REF]. It has been found very useful in many applications such as shape interpolation [START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF]; Oberman and Ruan 2015] and texture synthesis [START_REF] Bonneel | Sliced and radon Wasserstein barycenters of measures[END_REF]. The displacement interpolation is to find a target probability density ρ defined by Gaussian Density where {ρ k } m k=1 are a collection of given probability densities and {λ k } m k=1 are user-specified weights for computing the average probability density with regard to the 2-Wasserstein distance W(•, •). Recently many research works [START_REF] Cuturi | Fast computation of Wasserstein barycenters[END_REF][START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF][START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF]] suggested an energy term with an additional entropy regularization, bringing convexity to the problem and make it easier to solve. However, such an entropy term would result in an over-smoothed result in exchange for better efficiency. Besides, although lots of works proposed to compute the traditional Wasserstein barycenter distributions, few works aim at the semi-discrete setting. Here we propose a framework that supports multiple input distributions, each of which is represented by a collection of sites associated with capacities or a continuous probability density function. Our algorithm produces adaptively distributed samples to encode the final distribution.

min ρ m k=1 λ k W 2 2 (ρ k , ρ), (10)
(f) 2-Wasserstein (g) 1-Wasserstein (h) (∆x (1) ) 4 + (∆x (2) ) 4 (i) Geodesic distance (j) Squared (k) Cubic
Consider the following semi-discrete minimization problem, denoting R k = {R i k } n i=1 : Ω → {xi} n i=1 as an arbitrary assignment between the support of ρ k and X

         min X m k=1 λ k W 2 2 (ρ k , X, R k ), s.t. R i k ρ k dx = ci, for all i, k. (11)
The constraint term in (11) forces the mass transported to xi under ρ k to be equal to a given ci. Based on the conclusion of [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF]], each assignment R k minimizing the 2-Wasserstein cost under the capacity constraint should be a power diagram, and thus Eqn. ( 11) can be further reduced to the following minimization problem [de Goes et al. 2012]:

Minimize X F b (X) = m k=1 λ k n i=1 GD(W ) i k ρ k ||x -xi|| 2 dx - n i=1 wi(|GD(W ) i k | -ci) , (12)
with

∇x i F b (xi) = m k=1 2ciλ k (xi -b i k ), (13)
where b i k is the centroid of the i-th power cell under the k-th distribution. In implementation, the optimization procedure can follow a similar procedure to computing CPDs, presented in Algorithm 2. Fig. 4 shows the result from two identical circular spots. The ideal distribution when the weights of these two inputs are both 0.5 should be of the same circular shape at the middle position. We compare our result with that by [START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF]] and visualize the color-coded errors in Fig. 4(c&g). The resolutions of both input images are 200 × 200 in this example. We can see that when n = 1.5K, we achieve an error of RMSE=0.078 in about 30 seconds, which is at the same efficiency/accuracy level with [START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF]; But when n increases to 10K, the error of RMSE decreases to 0.0441 at the cost of more computation time -about 30 minutes for this example (the inefficiency is due to the density setting where lots of density values are very close to 0). Note that for comparison purpose, Fig. 4(e-f) are produced by an adaptive kernel estimation as the postprocessing [START_REF] Abramson | On bandwidth variation in kernel estimates-a square root law[END_REF]]. Fig. 5 shows the results of 2-Wasserstein barycenters with multiple inputs.

Blue Noise Sampling

Stochastic point distributions with blue noise properties are widely used for various applications in computer graphics [START_REF] Yan | A survey of blue-noise sampling and its applications[END_REF]]. de [START_REF] De Goes | Blue noise through optimal transport[END_REF] presented a fast algorithm to generate highquality blue noise point distributions by finding a stationary point of the objective function defined in Eqn. (5). However, extremizing the objective function does not prevent the regular patterns from emerging. In fact, the hexagonal grid patterns are the solutions to the extremization problem. de [START_REF] De Goes | Blue noise through optimal transport[END_REF] propose to jitter the points and run the optimization again when local regularities are detected. It is expected that a "shallower" extremum of the objective function which contains no regular patterns is obtained by a local searching algorithm. However, the optimization method of [START_REF] De Goes | Blue noise through optimal transport[END_REF]] used a very loose stopping criteria, which usually terminates in 10 iterations, well before convergence. A more faithful local extremum can be efficiently achieved by our reduced gradient method that optimizes the same function. Experimental results show that "extreme" optimization of the objective function always results in structured point distributions even starting from a random point distribution (see Fig. 6(a)).

It is worth noting that all point samples are enforced to have the same capacity c in the existing methods based on capacityconstrained Voronoi tessellation [START_REF] Balzer | Capacity-constrained point distributions: A variant of Lloyd's method[END_REF][START_REF] Li | Fast capacity constrained Voronoi tessellation[END_REF][START_REF] Chen | Variational blue noise sampling[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF]]. To prevent regular patterns in the point distributions, [START_REF] Zhang | Capacity constrained blue-noise sampling on surfaces[END_REF] suggested using non-uniform capacity constraints and dynamically adjust them during optimization. Inspired by this, we propose to jitter the capacity constraints for the points before optimization. A random noise with a maximum magnitude τ is used to jitter the capacity constraint for each point. The new capacity constraint c * i for each point is now specified by the following equation:

ci = (1 + τ • rand(-1, 1))c, c * i = ci ci nc,
where n is the sample number, and the second formula is to normalize the constraints so that they all sum to the total capacity of the domain.

Fig. 6(a) and 6(c) show point sets generated by extremizing the objective function with uniform capacity constraints and the jittered capacity constraints, respectively. A strict convergence threshold ( ∇X F ≤ 10 -7 ) is used for both cases. We can see that the jittered capacity constraints avoid regularity in the resultant point distribution, producing high-quality blue noise. For the sake of fast computation, we set τ as a small number 0.2 and use a loose stopping condition with ∇X F ≤ 0.05 √ nm 3 in practice. Under this setting, our optimization method is able to converge after a dozen of iterations, and generate high-quality blue noise samplings, as shown in Fig. 6(c).

As we jitter the capacity constraints for sampling points, an immediate concern is that the point distribution may not adapt to the given density function precisely. Fig. 7 shows a sampling result generated by our method for an intensity ramp. We count the point density in each quarter of the ramp, and compare it to the reference density of the ramp. We can see that jittering the capacity constraints break local structured patterns while returns good adaption to the given density function.

Centroidal Convex Decomposition

Another important problem in geometric computing, especially in mesh generation, is the convex decomposition of a geometric domain [START_REF] Greene | The decomposition of polygons into convex parts[END_REF][START_REF] Ghosh | Fast approximate convex decomposition using relative concavity[END_REF], that is, decomposing a general non-convex domain into the union of a number of convex cells. The problem is relevant to numerous applications in pattern recognition, Minkowski sum computation, motion planning, computer graphics, and origami folding [START_REF] Ghosh | Fast approximate convex decomposition using relative concavity[END_REF]. In this section we use the optimal capacity-constrained partition to provide a novel and effective solution to the convex decomposition of a 2D planar region.

Existing algorithms for convex decomposition usually minimize the number of convex components; See Fig. 9(a). In this paper, we propose a different type of convex decomposition, called the optimal convex decomposition that meets the following criteria: (1) each decomposed cell is convex; (2) the decomposed cells have specified areas; and (3) the decomposed cells are of compact shapes. Fig. 9(b) shows such a convex decomposition. Note that the term "optimal" is used to emphasize shape regularity rather than the number of decomposed cells.

Suppose that the boundary of the input 2D domain to be decom- posed is by a polygonal curve. It is observed that simply running a CPD algorithm would cause a concave vertex on the boundary to be located in the interior of a single CPD cell, resulting in a non-convex cell (see Fig. 10(a)). Our key idea is to define all concave boundary vertices as virtual sites that are fixed and associated with a sufficiently small capacity constraint. Thus, the new cost function becomes

F (X, W ) = n+k i=1 R W i ρ(x) x -xi 2 dx (14) subject to ci = R W i ρ(x)dx( mi), 1 ≤ i ≤ n + k, ( 15
)
where n is the number of variable sites and k is the number of virtual sites on the boundary. Here {mi, 1 ≤ i ≤ n} is a set of capacity constraints for the moveable sites, and the number of such sites, n, is usually specified by the user. Let maverage be the average value of {mi, 1 ≤ i ≤ n}. The capacity constraint for the virtual sites, say, mvirtual, is given by

mvirtual = × maverage, (16)
where is set to be 10 -3 in our experiments.

Furthermore, we assign a large multiplicative penalty to the cost terms defined by the virtual sites such that those small-sized power cells can be made compact and stay around the virtual sites (see Empirically we take the penalty coefficient to be

λ = mnormal maverage . (18)
After the above cost function is minimized (see Fig. 10(c)), we remove those virtual sites to get a clean convex decomposition in which each concave point is located on the common boundary of two adjacent convex power cells (see Fig. 10(d)). Fig. 11 demonstrates three examples of convex decomposition produced by our method based on CPDs with various density settings. We note that the density value should be larger in narrow regions and around concave points, generally implying a denser site distribution, to achieve a successful decomposition.

Conclusions and Future Works

In this paper, a super-linear convergent method is developed for computing the capacity-constrained centroidal power diagram. Experimental results show that our algorithm outperforms the existing method. Due to its effectiveness, we apply our fast CPD algorithm to three problems in computer graphics and geometric processing: displacement interpolation, blue-noise point sampling, and optimal For the Guitar shape (b) and the Gear-Teeth shape (c), the density is higher in narrow regions or around the user-specified concave points.

convex decomposition of 2D domains. Furthermore, we show that the proposed algorithm can be extended to the capacity-constrained optimal partition with respect to general cost functions.

Our algorithm, in its current form, has several limitations. First, the capacity constraints are required to be met throughout the optimization process, which may be a cause of inefficiency. It is therefore valuable to develop an optimization technique that can accommodate an infeasible start. Second, the optimal mass transport problem in 3D space has been studied [START_REF] Su | Measure controllable volumetric mesh parameterization[END_REF]Lévy 2014]. We are therefore interested in extending our method to 3D or spaces of higher dimensions. Finally, we have only studied optimal convex decomposition of a 2D domain. The optimal convex decomposition of 3D domain is more relevant problem in practice and is much more challenging that the 2D version, so would need further study.

  Given a general distance cost kernel d(x, y), where x, y ∈ Ω ⊂ E d , a generalized power diagram (also called a Laguerre diagram) of Ω with the weighted sites {(xi, wi)} n i=1 is defined to be the collection of cells GD = {GDi} n i=1 , where GDi = {x|d(x, xi)wi ≤ d(x, xj) -wj, ∀j, j = i}. When d(x, y) is the squared distance, the generalized power diagram reduces to the traditional power diagram whose cells are convex polytopes. However, with a general cost, the cells GDi are not necessarily convex and their boundaries are curves or curved surfaces rather than straight lines or planes.

Figure 1 :

 1 Figure 1: An optimal partition (X * , W * (X * )) is a saddle point of the cost function F (X, W ) in Eqn. (5).

  Fig. 1(b)). For each point (Xa, W * (Xa)) on the constraint surface, the restriction F (Xa, W ) of F (X, W ) (Eqn. (5)) to the subspace X = Xa attains its unique maximum at the point (Xa, W * (Xa)) among all the weighting schemes W (see the red curve in Fig.1(b)). Meanwhile, the restriction F (X, W * (X)) of F (X, W ) to the constraint surface attains its local minimum at a point (X * , W * (X * )) that gives an optimal capacity-constrained partition, which is marked with the blue ball on the green curve in Fig.1(b). Hence, an optimal partition (X * , W * (X * )) is a saddle point of the cost function F (X, W ) in Eqn. (5).

128

  Input: domain D, density ρ, cost kernel d, number of points n, capacity constraints {ci} and a threshold as the termination condition. 3 Initialization: set k = 0 and X0 to be n randomly generated sites. 4 repeat //L-BFGS 5 repeat //To meet the capacity constraints 6 Update W k by Newton's method or L-BFGS method 7 until ||∇W F (X, W )|| < 10 -12 Compute the gradients ∇X F (X, W * (X)) in Eqn. (9) 9 Estimate the step-size for updating X by a line-search 10

  e. d(v, xi) -wi ≤ d(v, xj) -wj, ∀j = i. After that, we compute the decomposition of each triangle depending on the following cases: If the site xi gives the minimum weighted distances to v1, v2, v3 at the same time, we simply assign the triangle f to xi. If the site xi gives the minimum weighted distance to v1 and the site xj gives the minimum weighted distances to v2, v3, we come to find two dividing points p1 ∈ v1v2 and p2 ∈ v1v3 such that p1p2 can be taken as a segment delimiting xi's cell and xj's cell. The dividing point p1 can be found from the equation d(xi, v1) -wi + (x -v1) • ∇xd(x, xi) |v 1 = d(xj, v2) -wj + (x -v2) • ∇xd(x, xj) |v 2 . The other dividing point p2 can be computed likewise. See the inset figure.If the minimum weighted distances to v1, v2, v3 are respectively given by xi, xj, x k , then we find three dividing points p1 ∈ v1v2, p2 ∈ v2v3, p3 ∈ v3v1 as mentioned above. We denote the center point by p = (p1 + p2 + p3)/3. Finally, we assign the subcell v1p1pp3 to xi, v2p2pp1 to xj, and v3p3pp2 to x k , as the inset figure shows.

Figure 2 :

 2 Figure 2: Performance comparison with [de Goes et al., 2012].The first row: snapshots of our algorithm. The second row: snapshots of[START_REF] De Goes | Blue noise through optimal transport[END_REF]. The third row: the plots of the cost function. The last row: the plots of the gradient of the cost function. Note that the horizontal axis in (g-l) denotes the timing (s).

  1-Wasserstein (e) (∆x (1) ) 4 + (∆x (2) ) 4

Figure 3 :

 3 Figure 3: Our algorithm can be extended to general cost kernels. The first row shows that CPD can generate compact cells with user-specified capacity constraints while CVT cannot. For the other examples, the capacity constraints divide the domain equally.

12

  Input: densities {ρ k } m k=1 , weights {λ k } m k=1 , capacities {ci} n i=1 , and a threshold . 3 Output: {xi} n i=1 . 4 repeat //L-BFGS 5 for k = 1, 2, 3, • • • , m 6 //for the k-th distribution 7 Optimize W k until each power cell amounts to the given capacity 8 Compute the centroids of the power diagram under k-th distribution 9 end 10 Compute the cost function F b (X) in Eqn. (12) 11 Compute the gradient w.r.t. X; See Eqn. (13) 12 Update X by step search 13 until ||∇X F b (X)|| < Algorithm 2: Capacity Constrained Wasserstein Barycenter with LBFGS.

Figure 4 :

 4 Figure 4: Displacement interpolation of two identical circular spots. (a)&(d): the two inputs of circular spots; (b)&(c): interpolated distribution by[START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF]] and the corresponding color-coded point-wise error (RMSE=0.0794), generated in 30 seconds (using the Matlab code provided by the authors of[START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF]). (e-g) our results for n=1.5K and n=10K. When n = 1.5K, we achieve RMSE=0.078 in about 30 seconds; When n increases to 10K, the error of RMSE decreases to 0.0441 at the cost of increased computation time -about 30 minutes for this example.

Figure 6 :

 6 Figure 6: Sampling results with different τ 's over a periodic domain with constant density. The blue noise properties can be analyzed through comparing the averaged periodograms, radial mean power (red curve), and anisotropy (blue curve).

Figure 5 :

 5 Figure 5: The interpolation by our method of the four inputs (500×500 in size) at the four corners of the array (highlighted in red). Each output frame is generated by 10K sites in about 15 minutes.

Fig. 8

 8 Fig.8shows two image stippling results generated by our optimization framework. The sampling points {xi} n i=1 are initialized randomly according to the piecewise-constant density function defined by the gray scale of the input image. Integrations involved in our formulation are exactly computed by constructing the pixel-cell intersections, as done in[START_REF] De Goes | Blue noise through optimal transport[END_REF]].

Figure 7 :

 7 Figure 7: Sampling of a quadratic density function with 1000 points, the percentages indicate point density in each quarter: (a) input quadratic density function; (b) sampling result by [de Goes et al. 2012]; (c) our result.

Figure 8 :Figure 9 :

 89 Figure 8: Image stippling. (a)&(c) Input images; (b)&(d) stippling results with 10K points generated by our optimization framework in 129 seconds and 146 seconds, respectively.

Figure 10 :

 10 Figure10: Computation of the optimal convex decomposition: (a) Generally the concave points may be covered totally by a single CPD cell, giving rise to a non-convex cell. (b) Even if we enforce a small capacity for the virtual site at the concave point, the power cell stay non-convex after optimization. (c) After adding a large penalty coefficient λ for the virtual sites, the minimization makes all the normal power cells almost convex. (d) By removing the virtual sites from the results in (c), we obtain a convex decomposition of the domain.

Figure 11 :

 11 Figure 11: Convex decomposition of complex 2D shapes. The blue dots denote the user-specific concave points on the given boundary. For the Tortoise shape, we set a uniform density distribution. For the Guitar shape (b) and the Gear-Teeth shape (c), the density is higher in narrow regions or around the user-specified concave points.

Table 1 :

 1 Statistics of the number of operations of building power cells for the three examples shown in Fig.2. Note that we take ∇F ≤ 10 -8 as the termination condition.

	Examples	# of building power cells [de Goes et al. 2012] Ours	Timing (s) [de Goes et al. 2012] Ours
	100 sites	4091	279	20.8	1.6
	500 sites	2944	464	246.7	39.5
	1000 sites	4307	471	556.0	58.4
	We implemented and experimented with our CPD algorithm on a
	computer with a 64-bit version of Win7 system, a 3.07 GHz In-

tel(R) Core(TM) i7 CPU and 6 GB memory. The coding language is C++. Now we shall compare our CPD algorithm with the state of the art

[START_REF] De Goes | Blue noise through optimal transport[END_REF]

. For the examples shown in Fig.

2

, we take ∇F ≤ 10 -8 as the termination condition. For the example of n = 100 sites (see the left column), we assume that the domain [0, 1] × [0, 1] is associated with a uniform density. For the

Acknowledgements

We are grateful to the anonymous reviewers for their helpful comments. This work is supported by NSF of China (61332015, 61300168, 61472332, 61100105, 61572021, 61272019), GRF of RGC ( 17263316), National 973 Program (2015CB352501) and SZfund (SIRI0404201431). It is also partially supported by ANR grants BECASIM and MAGA.

Appendix

Theorem 1: For any fixed sites X = {xi}, the capacityconstrained partition R = {Ri} n i=1 of the domain Ω that minimizes the cost function in Eqn. ( 4

On the other hand, it follows from the optimality of R, in particular, the optimality of Ra and R b that swapping Bx and By between Ra and R b , while maintaining the capacities of Ra and R b , does not reduce the total cost. Therefore for some constants wa and w b . This completes the proof.