Lhassane Idoumghar
email: lhassane.idoumghar@uha.fr

Mahmoud Melkemi

René Schott

Maha Idrissi Aouad

Hybrid PSO-SA Type Algorithms for Multimodal Function Optimization and Reducing Energy Consumption in Embedded Systems

The paper presents a novel hybrid evolutionary algorithm that combines Particle Swarm Optimization (PSO) and Simulated Annealing (SA) algorithms. When a local optimal solution is reached with PSO, all particles gather around it, and escaping from this local optima becomes difficult. To avoid premature convergence of PSO, we present a new hybrid evolutionary algorithm, called HPSO-SA, based on the idea that PSO ensures fast convergence, while SA brings the search out of local optima because of its strong local-search ability. The proposed HPSO-SA algorithm is validated on ten standard benchmark multimodal functions for which we obtained significant improvements. The results are compared with these obtained by existing hybrid PSO-SA algorithms. In this paper, we provide also two versions of HPSO-SA (sequential and distributed) for minimizing the energy consumption in embedded systems memories. The two versions, of HPSO-SA, reduce the energy consumption in memories from 76% up to 98% as compared to Tabu Search (TS). Moreover, the distributed version of HPSO-SA provides execution time saving of about 73% up to 84% on a cluster of 4 PCs.

Introduction

Several optimization algorithms have been developed over last few decades for solving real-world optimization problems. Among them, we have many heuristics like Simulated Annealing (SA) [START_REF] Locatelli | Simulated annealing algorithms for continuous global optimization[END_REF] and optimization algorithms that make use of social or evolutionary behaviors like Particle Swarm Optimization (PSO) [START_REF] Pant | Particle swarm based metaheuristics for function optimization and engineering applications[END_REF][START_REF] Kennedy | Swarm intelligence[END_REF]. SA and PSO are quite popular heuristics for solving complex optimization problems, but they have some strengths and limitations.

Particle Swarm Optimization (PSO) is based on the social behavior of individuals living together in groups. Each individual tries to improve itself by observing other group members and imitating the better ones. This way, the group members are performing an optimization procedure which is described in [START_REF] Kennedy | Swarm intelligence[END_REF]. The performance of the algorithm depends on how the particles (i.e., potential solutions to an optimization problem) move in the search space, given that the velocity is updated iteratively. Large research body is therefore devoted to the analysis and proposal of different motion rules (see [START_REF]Convergence behavior of the fully informed particle swarm optimization algorithm[END_REF][START_REF] Engelbrecht | Fundamentals of Computational Swarm Intelligence[END_REF][START_REF] Poli | Particle swarm optimization. An overview[END_REF] for recent accounts of PSO research). To avoid premature convergence of PSO, we combine it with SA: PSO contributes to the hybrid approach in a way to ensure that the search converges faster, while SA makes the search jump out of local optima due to its strong local-search ability. In this paper, we present a hybrid optimization algorithm, called HPSO-SA, which exploits intuitively the positive features of PSO and SA. We also validate HPSO-SA using ten benchmark functions given in [START_REF] Suganthan | Problem definitions and evaluation criteria for the cec 2005 special session on real-parameter optimization[END_REF] and compare the results with classical PSO, ATREPSO, QIPSO, and GMPSO algorithms described in [START_REF] Pant | Particle swarm based metaheuristics for function optimization and engineering applications[END_REF], TL-PSO [START_REF] Nakano | Consideration of particle swarm optimization combined with tabu search[END_REF], PSO-SA [START_REF] Che N | Ane wh y b ridalgorithmof particle swarm optimization[END_REF], SAPSO and SUPER-SAPSO presented in [START_REF] Bahrepour | Super-sapso: a new sa-based pso algorithm[END_REF]. We provide also two versions of HPSO-SA (sequential and distributed) for minimizing the energy consumption in embedded systems memories. The two versions, of HPSO-SA, reduce the energy consumption in memories from 76% up to 98% as compared to Tabu Search (TS). Moreover, the distributed version of HPSO-SA provides execution time saving of about 73% up to 84% on a cluster of 4 PCs.

The rest of the paper is organized as follows. Section 2 introduces, briefly, PSO and SA algorithms. Section 3 is devoted to detailed description of HPSO-SA. In Section 4, benchmark functions are applied on HPSO-SA. In Section 5, HPSO-SA is used to solve the energy consumption problem in memory. In addition, simulation results are provided a n dc o m p a r e dw i t ht h o s eo f [START_REF] Aouad | A tabu search heuristic for scratch-pad memory management[END_REF]. Conclusions and further research aspects are given in Section 6.

Background

2.1. Simulated Annealing Algorithm. SA [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] is a probabilistic variant of the local search method, which can, in contrast to PSO, escape local optima. SA is based on an analogy taken from thermodynamics which is as follows: In order to grow a crystal, we start by heating material until it reaches its molten state. Then, we reduce the temperature of this crystal melt gradually, until the crystal structure is formed. A standard SA procedure begins by generating an initial solution at random. At initial stages, a small random change is made in the current solution s c . Then the objective function value of the new solution s n is calculated and compared with that of the current solution. A move is made to the new solution if it has better value or if the probability function implemented in SA has a higher value than a randomly generated number.

Otherwise a new solution is generated and evaluated. The probability of accepting a new solution is given as follows:

p = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1, if f (s n ) -f (s c ) < 0, exp -f (s n ) -f (s c ) T ,o t h e r w i s e . ( 1
)
The calculation of this probability relies on a parameter T, which is referred to as temperature, since it plays a similar role as the temperature in the physical annealing process. To avoid getting trapped at a local minimum point, the rate of reduction should be slow. In our problem we use the following method to reduce the temperature: γ = 0.99 and

T i+1 = γT i,( 2 )
where i = 0, 1, .... T h u s ,a tt h es t a r to fS Am o s tw o r s e n i n gm o v e sm a yb e accepted, but in the end only improving ones are likely to be allowed, which can help the procedure jump out of a local minimum. The algorithm may be terminated after a certain volume fraction of the structure has been reached or after a prespecified runtime. [START_REF] Kennedy | Particle swarm optimization[END_REF], inspired by social behavior patterns of organisms that live and interact within large groups. In particular, it incorporates swarming behaviors observed in flocks of birds, schools of fish, or swarms of bees, and even human social behavior. PSO algorithm is based on an idea that particles move through the search space with velocities that are dynamically adjusted according to their historical behaviors. Therefore, the particles have the tendency to move towards the better and better search area over the course of search process. PSO algorithm starts with a group of random (or not) particles (solutions) and then searches for optima by updating each generation. Each particle is treated as a volume-less particle (a point) in the n-dimensional search space. The ith particle is represented as X i = (x i1 , x i2 , ..., x in ). At each generation, the particles are updated by using following two best values.

Particle Swarm Optimization. PSO is a population based stochastic optimization technique developed by

V i (t) V i+1 (t) g b e sti ( t) -Xi ( t) pbest i (t) -X i (t) gbest i (t) pbest i (t) X i+1 (t) X i (t)
(i) The first value is the best solution (fitness) a particle has achieved so far (the fitness value is also stored). This value is called pbest.

(ii) The second value is the best value tracked by the particle swarm optimizer so far (by any particle) in the population. This best value is a global best and is called gbest.W h e nap a r t i c l et a k e sp a r to fa population as its topological neighbors, the best value is a local best and is called lbest.

At each iteration, these two best values are combined to adjust the velocity along each dimension, which is then used to compute a new iteration step for the particle. A portion of adjustment to the velocity is influenced by the individual's previous best position (pbest), considered as the cognition component, and this portion is influenced by the best in the neighborhood (lbest or gbest), the social component (see Figure 1). With the addition of the inertia factor, ω,by [START_REF] Shi | A modified particle swarm optimizer[END_REF] (for balancing the global and the local search), the equations for velocity adjustment are

v i+1 = ωv i + c 1 * random(0, 1) * pbest i -x i +c 2 * random(0, 1) * gbest i -x i , ( 3
)
x i+1 = x i + v i+1 , (4)
where random(0, 1) is a random number independently generated within the range of [0, 1] and c 1 and c 2 are two learning factors which control the influence of the social and cognitive components (usually, c 1 = c 2 = 2, see [START_REF] Shi | Empirical study of particle swarm optimization[END_REF]). In [START_REF] Kennedy | Swarm intelligence[END_REF], if the sum on the right side exceeds a constant value, then the velocity on that dimension is assigned to be V i min or V i max. Thus, particle velocities are clamped to the range of [V i min; V i max], which serves as a constraint to control the global exploration ability of PSO algorithm. This also reduces the likelihood of particles leaving the search space. Note that the values of x i are not restricted to the range [V i min; V i max]; it only limits the maximum distance that a particle will move during one iteration.

HPSO-SA Hybrid Algorithm

This section presents a new hybrid HPSO-SA algorithm which combines the advantages of both PSO (that has a strong global-search ability) and SA (that has a strong localsearch ability). Other applications of hybrid PSO and SA algorithm can be found [START_REF] Che N | Ane wh y b ridalgorithmof particle swarm optimization[END_REF][START_REF] Bahrepour | Super-sapso: a new sa-based pso algorithm[END_REF][START_REF] Xia | A hybrid particle swarm optimization approach for the job-shop scheduling problem[END_REF][START_REF] Chaojun | Particle swarm optimization algorithm based on the idea of simulated annealing[END_REF][START_REF] Fang | Particle swarm optimization with simulated annealing for tsp[END_REF][START_REF] Wang | Hybrid particle swarm optimization with simulated annealing[END_REF].

This hybrid approach makes full use of the exploration capability of both PSO and SA and offsets the weaknesses of each. Consequently, through application of SA to PSO, the proposed algorithm is capable of escaping from a local optimum. However, if SA is applied to PSO at each iteration, the computational cost will increase sharply and at the same time the fast convergence ability of PSO may be weakened. In order to flexibly integrate PSO with SA, SA is applied to PSO every K iterations if no improvement of the global best solution does occur. Therefore, the hybrid HPSO-SA approach is able to keep fast convergence (most of the time) thanks to PSO, and to escape from a local optimum with the aid of SA. In order to allow PSO jump out of a local optimum, SA is applied to the best solution in the swarm foundsofar ,eachK iterations that is predefined to be 270 ∼ 500 (based on our experimentations).

The hybrid HPSO-SA algorithm works as illustrated in Algorithm 1, where one has the following.

(i) Description of a Particle. Each particle (solution) X ∈ S is represented by its n>0 components, that is, X = (x 1 , x 2 , ..., x n ), where i = 1, 2, ..., n and n represents the dimension of the optimization problem to solve.

(ii) Initial Swarm. Initial Swarm corresponds to population of particles that will evolve. Each particle x i is initialized with uniform random value between the lower and upper boundaries of the interval defining the optimization problem.

(iii) Evaluate Function. Evaluate (or fitness) function in HPSO-SA algorithm is typically the objective function that we want to minimize in the problem.

It serves for each solution to be tested for suitability to the environment under consideration.

(iv) SA Algorithm. If no improvement of the global best solution occur during the last K iterations, then it means that the algorithm is trapped in a local optimum point. To escape out from local optimum, we apply SA algorithm to global best solution. The performance of SA depends on the definition of the several control parameters.

(a) Initial Temperature.K i r k p a t r i c k [START_REF] Kirkpatrick | Optimization by simulated annealing: quantitative studies[END_REF] suggested that a suitable initial temperature is one that resultsinanaverageprobabilityχ 0 of a solution that increases f being accepted of about 0.8. The value of T 0 will clearly depend on the scaling of f and, hence, be problem-specific. It can be estimated by conducting an initial search (100 iterations in next simulations) in which all increases in f are accepted and calculating the average objective increase observed δf. T 0 is then given by

T 0 =- δf ln χ 0 . ( 5
)
(b) Accept Function.FunctionAccept (current solution, Neighbor, T) is decided by the acceptance probability given by ( 1), which is the probability of accepting configuration Neighbor. (c) Generate Function. The neighborhood of each solution x is generated by using the following Equation:

x ←-x + dσr, [START_REF] Poli | Particle swarm optimization. An overview[END_REF] where d is the direction of the new neighborhood and takes either 1 or -1, σ is random number with Gaussian (0,1) distribution and r is a constant that correspond to the radius of neighborhood generator.

Experiments Results

Benchmark Functions.

In order to compare the performance of HPSO-SA hybrid algorithm with those described in [2, 8-10], we use benchmark functions [START_REF] Suganthan | Problem definitions and evaluation criteria for the cec 2005 special session on real-parameter optimization[END_REF] described in Ta b l e 1 . These functions provide a good starting point for testing the credibility of an optimization algorithm. For each of these functions, there are many local optima in their solution spaces. The number of local optima increases with increasing complexity of the functions, that is, with increasing dimension. In the following experiments, we used 10-, 20-and 30-dimensional functions except in the case of Himmelblau and Shubert functions that are twodimensional by definition (see Figure 2 for 3D representation).

Simulation Results and Discussions.

To v e r i f y t h e e fficiency and effectiveness of HPSO-SA hybrid algorithm, the experimental results of HPSO-SA approach are compared with those obtained by [START_REF] Pant | Particle swarm based metaheuristics for function optimization and engineering applications[END_REF][START_REF] Nakano | Consideration of particle swarm optimization combined with tabu search[END_REF][START_REF] Che N | Ane wh y b ridalgorithmof particle swarm optimization[END_REF][START_REF] Bahrepour | Super-sapso: a new sa-based pso algorithm[END_REF]. Our HPSO-SA hybrid

(1) iter ←0, cpt ← 0, Initialize swarm size particles (2) stop criterion←maximum number of function evaluations or Optimal solution is not attained (3) while Not stop criterion do (4) for each particle i ← 1 to swarm size do [START_REF] Engelbrecht | Fundamentals of Computational Swarm Intelligence[END_REF] Ú ÐÙ Ø (particle(i)) if the fitness value is better than the best fitness value (cbest) in history then [START_REF] Poli | Particle swarm optimization. An overview[END_REF] Update current value as the new cbest. [START_REF] Suganthan | Problem definitions and evaluation criteria for the cec 2005 special session on real-parameter optimization[END_REF] end [START_REF] Nakano | Consideration of particle swarm optimization combined with tabu search[END_REF] end [START_REF] Che N | Ane wh y b ridalgorithmof particle swarm optimization[END_REF] Choose the particle with the best fitness value in the neighborhood (gbest) [START_REF] Bahrepour | Super-sapso: a new sa-based pso algorithm[END_REF] for each particle i ← 1 to swarm size do [START_REF] Aouad | A tabu search heuristic for scratch-pad memory management[END_REF] Update particle velocity according to Equation (3) [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] Enforce velocity bounds [START_REF] Kennedy | Particle swarm optimization[END_REF] Update particle position according to Equation (4) ( 14)

Enforce particle bounds ( 15) end [START_REF] Xia | A hybrid particle swarm optimization approach for the job-shop scheduling problem[END_REF] if there is no improvement of global best solution then [START_REF] Chaojun | Particle swarm optimization algorithm based on the idea of simulated annealing[END_REF] cpt ← cpt +1 ( 18) end [START_REF] Wang | Hybrid particle swarm optimization with simulated annealing[END_REF] Update global best solution [START_REF] Kirkpatrick | Optimization by simulated annealing: quantitative studies[END_REF] cpt ← 0 ( 21)

if cpt = K then (22) cpt ← 0 (23)
//Apply SA to global best solution [START_REF] Fradj | Energy aware memory architecture configuration[END_REF] iterSA ← 0, Initialize T according to Equation (5) [START_REF] Gendreau | An Introduction to Tabu Search[END_REF] current solution ← global best solution [START_REF] Tanenbaum | Architecture de l'ordinateur[END_REF] current cost ← Ú ÐÙ Ø (current solution) [START_REF] Kellerer | Knapsack Problems[END_REF] while Not SA stop criterion do [START_REF]Benchmarks[END_REF] while inner-loop stop criterion do (29)

Neighbor ← Ò Ö Ø (current solution) (30) Neighbor cost ← Ú ÐÙ Ø (

Comparison with Results

Obtained by Using TL-PSO Algorithm [START_REF] Nakano | Consideration of particle swarm optimization combined with tabu search[END_REF]. In this section we compare HPSO-SA approach with TL-PSO method [START_REF] Nakano | Consideration of particle swarm optimization combined with tabu search[END_REF]t h a ti sb a s e do n combining the excellence of both PSO and Tabu Search. As described in [START_REF] Nakano | Consideration of particle swarm optimization combined with tabu search[END_REF], we apply HPSO-SA algorithm to the four following benchmark problems: Rastrigin, Schwefel, Griewank and Rosenbrock function. Here, the number of particles in the swarm is 30. The number of dimension of searching n = 20 and the number of objective function evaluations is 60000 (i.e., 2000 × 30). The results obtained after numerical simulations are shown in Ta b l e 2 .T h e s e results indicate the Mean, Best, and Worst values obtained under the same condition over 50 trials. By analyzing Ta b l e 2 , we conclude that the results obtained by HPSO-SA algorithm are preferable in comparison with those obtained by TL-PSO algorithm. [START_REF] Pant | Particle swarm based metaheuristics for function optimization and engineering applications[END_REF]. Performance of four Particle Swarm Optimization algorithms, namely classical PSO, Attraction-Repulsion based PSO (AT RE P S O), Quadratic Interpolation based PSO (QIPSO) and Gaussian Mutation based PSO (GMPSO)i s evaluated in [START_REF] Pant | Particle swarm based metaheuristics for function optimization and engineering applications[END_REF]. The algorithms presented in this paper are guided by the diversity of the population to search the global optimal solution of a given optimization problem, where as GMPSO uses the concept of mutation and QIPSO uses the reproduction operator to generate a new member of the swarm.

Comparison with Other PSO Algorithms Described in

20 + e -20e -0.2((1/n) n i=1 x 2 i ) 1/2 -e (1/n) n i=1 cos(2πxi) [-30.0; 30.0] 0 Multimodal Michalewicz - n i=1 sin(x i )sin 2m ((i -x 2 i )/π)[ -π; π]- M u l t i m o d a l Himmelblau (x 2 + x 2 1 -11) 2 +(x 1 + x 2 2 -7) 2 +
In order to make a fair comparison between classical PSO, ATREPSO, QIPSO, GMPSO and HPSO-SA approach, we fixed, as indicated in [START_REF] Pant | Particle swarm based metaheuristics for function optimization and engineering applications[END_REF] ,t h es a m es e e df o rr a n d o m number generation so that the initial swarm population is same for all five algorithms. The number of particles in the swarm is 30. The algorithms use a linearly decreasing inertia weight ω which starts at 0.9 and ends at 0.4, with the user defined parameters c 1 = c 2 = 2.0. For each algorithm, the number of objective function evaluations is 300000. A total of 30 runs for each experimental setting were conducted and the average fitness of the best solutions throughout the run was recorded. The mean solution and the standard deviation (note that the standard deviation indicates the stability of the algorithms), found by the five algorithms, is listed in Ta b l e 3 . The numerical results given in Ta b l e 3 show the following.

(i) All the algorithms outperform the classical Particle Swarm Optimization.

(ii) HPSO-SA algorithm gives much better performances in comparison to PSO, QIPSO, ATREPSO, and GMPSO, out of the Sphere's and Ackley's functions.

(iii) On Sphere's function, QIPSO obtains better results than those obtained by HPSO-SA approach. But when the maximum number of iterations is fixed to 1.5 × 10 6 , HPSO-SA obtains the optimal value.

(iv) The analysis of the results, obtained for Ackley's function, shows that QIPSO obtains better mean result than HPSO-SA algorithm. However, HPSO-SA has a much smaller standard deviation.

Comparison with

Other PSO Algorithms Described in [START_REF] Bahrepour | Super-sapso: a new sa-based pso algorithm[END_REF]. In this section four benchmark functions are used to compare the relative performance of HPSO-SA algorithm with SUPER-SAPSO, SAPSO, and PSO algorithms described in [START_REF] Bahrepour | Super-sapso: a new sa-based pso algorithm[END_REF]. For all comparisons, the number of particles was set to 30. HPSO-SA algorithm uses a linearly decreasing inertia weight ω which starts at 0.9 and ends at 0.4, with the user defined parameters c 1 = c 2 = 2.0, 20 runs are conducted for each experimental setting and, for each algorithm, the average value is given in Ta b l e 4 .

In all above experiments, HPSO-SA algorithm obtains better results in comparison to those obtained by both the standard PSO and SAPSO algorithm [START_REF] Bahrepour | Super-sapso: a new sa-based pso algorithm[END_REF]. A comparison of HPSO-SA algorithm and SUPER-SAPSO [START_REF] Bahrepour | Super-sapso: a new sa-based pso algorithm[END_REF], shows that the last one converges faster than HPSO-SA.

SUPER-SAPSO uses an expression for the particle movements (x t+1 = (x t + v t+1 )T where T ≪ 1) which is welladapted to the case where the global optimum is 0. This is the reason why SUPER-SAPSO needs a very small number of iterations in this case. [START_REF] Che N | Ane wh y b ridalgorithmof particle swarm optimization[END_REF]. In this section performances of HPSO-SA are compared with these of PSOSA [START_REF] Che N | Ane wh y b ridalgorithmof particle swarm optimization[END_REF], Genetic Algorithm and hybrid algorithm [START_REF] Morten | Hybrid Particle Swarm Optimization with Breeding and Subpopulations,S p r i n g e r[END_REF].

Comparison with PSO-SA Algorithm Described in

Ta b l e 5 lists different results obtained for three different dimensions of each function. The optimum value of Sphere, Rastrigrin and Griewank was set to be 1e -10 and the goal value of Rosenbrock function was set to be 1e -06 (as indicated in [START_REF] Che N | Ane wh y b ridalgorithmof particle swarm optimization[END_REF]).

To make a fair comparison, the maximum number of function evaluations allowed was set to 20000, 30000 and 40000 for HPSO-SA and PSOSA algorithms when the number of particle was set to 20. HPSO-SA algorithm uses a linearly decreasing inertia weight ω which starts at 0.9 and ends at 0.4, with the user defined parameters c 1 = c 2 = 2.0.

The numerical results given in Ta b l e 5 show that:

(i) Over four benchmark functions, HPSO-SA and PSOSA do better than standard GA and hybrid algorithm [START_REF] Morten | Hybrid Particle Swarm Optimization with Breeding and Subpopulations,S p r i n g e r[END_REF].

(ii) For Sphere, Rastrigin and Griewank functions, HPSO-SA and PSOSA algorithms obtain optimal solutions within specified constrains (number of objective function evaluations).

(iii) For Rosenbrock function, PSOSA obtains better results than HPSO-SA for dimension 20, but for dimensions 10 and 30, HPSO-SA does better and has smaller standard deviation.

Reducing Memory Energy Consumption in Embedded Systems

5.1. Description of the Memory Problem. According to trends in [START_REF]ITRS, System Drivers[END_REF], memory will become the major energy consumer in an embedded system. Indeed, embedded systems must integrate multiple complex functionalities which needs bigger battery and memory. Hence, reducing memory energy consumption of these systems has never been as topical.

In this paper, we will focus on software techniques for the memory management. In order to reduce memory energy consumption, most authors rely on Scratch-Pad Memories (SPMs) rather than caches [START_REF] Aouad | A survey of scratch-pad memory management techniques low-power and -energy[END_REF]. Although cache memory helps a lot with program speed, it is not the appropriate for most of the embedded systems. In fact, cache increases the system size and its energy cost (cache area plus managing logic). Like cache, SPM consists of small, fast SRAM. The main difference is that SPM is directly and explicitly managed at the software level, either by the developer or by the compiler which makes it more predictable. SPM requires up to 40% less energy and 34% less area than cache [START_REF] Fradj | Energy aware memory architecture configuration[END_REF]. In this paper, we will therefore use an SPM in our memory architecture. Due to the reduced SPM size, we allocate space for interesting data only whereas, the remaining is placed in main memory (DRAM). In order to determine interesting data, we use data profiling to gather memory access frequency information. The Tabu Search (TS) approach consists of allocating space for data in SPM based on TS principles [START_REF] Gendreau | An Introduction to Tabu Search[END_REF]. More details about how TS is implemented can be found in [START_REF] Aouad | A tabu search heuristic for scratch-pad memory management[END_REF].

In order to compute energy cost of the system, we propose an energy consumption estimation model, for our memory architecture composed by an SPM, an instruction cache and a DRAM. Equation [START_REF] Suganthan | Problem definitions and evaluation criteria for the cec 2005 special session on real-parameter optimization[END_REF] gives the energy model where the three terms refer to the total energy consumed, respectively, in SPM, in instruction cache and in DRAM.

E = E tspm + E tic + E tdram . ( 7
)
In this model, we distinguish between the two cache write policies: Write-Through (WT)a n dW r i t e -B a c k( WB). In a WT cache, every write to the cache causes a synchronous write to DRAM. Alternatively, in a WB cache, writes are not immediately mirrored to DRAM. Instead, the cache tracks which of its locations have been written over and then, it marks these locations as dirty. The data in these locations is w r i t t e nb a c kt oD R A Mw h e nt h a td a t ai se v i c t e df r o mt h e cache [START_REF] Tanenbaum | Architecture de l'ordinateur[END_REF]. In this paper, the aim is to minimize the energy for the detailed estimation model presented as follows:

E = N spmr * E spmr (8) + N spmw * E spmw (9)
h ik * E icr + 1 -h ik * [E dramr + E icw + (1 -WP i ) * DB ik * (E icr + E dramw ) (10)
+ Nicw k=1 WP i * E dramw + h ik * E icw + (1 -WP i ) * 1 -h ik * E icw + DB ik * (E icr + E dramw ) (11)
+ N dramr * E dramr [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] + N dramw * E dramw. [START_REF] Kennedy | Particle swarm optimization[END_REF] Equations ( 8)a n d( 9) represent, respectively, the total energy consumed during reading and writing from/to SPM. Equations (10)a n d( 11) represent, respectively, the total energy consumed during reading and writing from/to instruction cache. When, Equations ( 12)and( 13)represent, respectively, the total energy consumed during reading and writing from/to DRAM. The various terms used in this energy model are explained in Ta b l e 6 .

As SPM has got a lot of advantages, it is clearly preferable to put as much data as possible in it. In other words, we must maximize terms N spmr and N spmw in the model. Hence, the problem becomes to maximize the number of accesses to the SPM. It is therefore a combinatorial optimization problem like knapsack problem [START_REF] Kellerer | Knapsack Problems[END_REF]. We want to fill SPM that can hold a maximum capacity of C with some combination of data from a list of N possible data each with size i and access number i so that the access number of the data allocated into SPM is maximized. This problem has a single linear constraint, a linear objective function which sums Reading access number to SPM.

N spmw

Writing access number to SPM.

E icr

Energy consumed during a reading from instruction cache.

E icw

Energy consumed during a writing into instruction cache.

N icr

Reading access number to instruction cache.

N icw

Writing access number to instruction cache.

E dramr

Energy consumed during a reading from DRAM.

E dramw

Energy consumed during a writing into DRAM.

N dramr

Reading access number to DRAM. N dramw Writing access number to DRAM.

WP i

The considered cache write policy: WT or WB. In case of WT, WP i = 1 else in case of WB then WP i = 0.

DB i k

Dirty Bit used in case of to indicate during the access k if the instruction cache line has been modified before (DB i = 1) or not (DB i = 0).

h i k

Type of the access k to the instruction cache. In case of cache hit, h i k = 1. In case of cache miss,

h i k = 0.
the sizes of the data allocated into SPM, and the added restriction that each data will be in the SPM or not. If N is the total number of data, then a solution is just a finite sequence s of N terms such that s[n]i se i t h e r0o rt h es i z eo ft h en th data. s[n] = 0 if and only if the n th data is not selected in the solution. This solution must satisfy the constraint of not exceeding the maximum SPM capacity (i.e., N i=1 s[i] ≤ C).

Discrete Sequential

Hybrid HPSO-SA Algorithm. This section should be considered as an attempt to use hybrid evolutionary algorithms for reducing energy consumption in embedded systems. Here, the focus is on the use of HPSO-SA algorithm designed in previous sections. Since the problem under consideration is dicrete and has specific features, HPSO-SA needs changes.

solution (Particle).

A solution can be represented by an array having size equal to the number of the data. Each element from this array denotes whether a data is included in the SPM ("1") or not ("0"). The HPSO-SA algorithm starts with an initial swarm which is randomly initialized.

Evaluate Function. It is the objective function that we want to minimize in the problem. It serves for each solution to be tested for suitability to the environment under consideration

Evaluate(solution) = Total Number Access all data

-Number Access(solution). [START_REF] Shi | A modified particle swarm optimizer[END_REF] Position Update Equation. Each dimension j of the particle i is updated by using [START_REF] Shi | Empirical study of particle swarm optimization[END_REF]:

x ij = ⎧ ⎪ ⎨ ⎪ ⎩ 1, rand < sigm v ij , 0, otherwise , (15)
where sigm(v ij ) is the sigmoid function, used to scale the velocities between 0 and 1, defined as:

sigm v ij = 1 1+exp -v ij ( 16
)
Generate Function. SA uses a notion of neighborhood relation. Let S be the set of all feasible solutions to the problem and f : S → Ê the objective function to be minimized. A neighborhood relation is a binary relation N ⊆ S × S with some desired properties. The interpretation of N (s, s ′ )isthat solution s is a neighbor of solution s ′ in the search space of all solutions S. A neighbor heuristic proceeds in steps. Starting search at some initial solution s 0 and then each step moves from the current solution to some neighbor according to rules specific to the heuristic. At each iteration, SA algorithm generates a random neighbor of the current solution (line 10). The neighborhood relation is defined as follows:

(1) with probability equal to 0.03, the value of each element of current solution is flipped from 1 to 0 or from 0 to 1; (2) validate solution:

while Not feasible (current solution must satisfy the constraint of not exceeding the maximum SPM capacity.) (current solution)do Remove the data j having a low number of access from current solution:(current solution[ j] ← 0).

Accept Function. the key idea in the SA approach is the function Accept which specifies the probability of accepting the move from current solution to a Neighbor solution, which also depends on so called temperature (T). The function Accept should satisfy the following conditions:

( Neighbor is worse than current solution the value of p decreases with time and tends to 0.

The function Accept(C cost , N cost , T) is decided by the probability of accepting configuration Neighbor. This probability is given by the following formula:

p = ⎧ ⎪ ⎨ ⎪ ⎩ 1, N cost <C cost , 1 2 rand * 1+e (Ccost-Ncost)/T ,o t h e r w i s e , ( 17
)
where T is the temperature and rand is a random number independently generated within the range of [0, 1].

Discrete Cooperative Distributed

Hybrid HPSO-SA Algorithm. For Distributed hybrid HPSO-SA (HPSO-SA Dist) algorithm, we use independent subswarms of particles with their own fitness functions which evolve in isolation, except for an exchange of some particles (migration). A set of m = 3 0p a r t i c l e si sa s s i g n e dt oe a c ho ft h eP processors, for a total population size of m × P. The set assigned to each processor is its subswarm. The processors are connected by an interconnection network with a ring topology. Initial subswarms consist of a randomly constructed assignment created at each processor. Each processor, disjointly and in parallel, executes the HPSO-SA Seq algorithm on its subswarm for a certain number of generations. Afterwards, each subswarm exchanges (HPSO -SA Dist runs in Asynchronous mode: at the 100th iteration, each processor sends its best solution and continues the improvement of its subswarm and verifies if it does not receive a solution from its neighbor.) its best particle (migrant) with its neighbors. We exchange the particles themselves (i.e., the migrant is removed from one subswarm and added to another). Hence, the size of the subswarm remains the same after migration (the worst particle is removed). The process continues with the separate improvement of each current solution for a maximum number of iterations. At the end of the process the best solution that exists constitutes the final assignment.

Experimental Results.

In order to compute the energy cost of studied memory architecture composed by an SPM, an instruction cache and a DRAM, we proposed an energy consumption estimation model which is explained in [START_REF] Aouad | A tabu search heuristic for scratch-pad memory management[END_REF]. Hybrid HPSO-SA algorithms and TS have been implemented on a cluster of PCs running under Windows XP Professional version 2002. The cluster is composed by 4 Pentium (D) machines running at 3 GHz. Each processor has 1 Gbyte of memory. Ta b l e 7 gives a description of the benchmarks used and they also can be downloaded from [START_REF]Benchmarks[END_REF].

In experiments, 30 different executions for each heuristic are performed and the best and average results obtained on these 30 executions are recorded. In this case, the best and the average solutions give similar results. Figure 3 shows that both HPSO-SA Seq and HPSO-SA Dist achieve better performances than TS on energy savings. In fact, hybrid HPSO-SA heuristics consume from 76.23% (StatemateCE) to 98.92% (ShaCE) less energy than TS.

As HPSO-SA Seq and HPSO-SA Dist give similar results, we decide to experiment their behavior when considering their execution time. We recorded the average execution times needed by HPSO-SA Seq and HPSO-SA Dist (running on a cluster of 4 PCs) to achieve the 30 executions. requires 73.16% (AdpcmCE) to 84.65% (CntCE) less execution time than HPSO-SA Seq. The Distributed HPSO-SA versionisalwaysfasterthantheSequentialHPSO-SAversion.

Conclusion and Perspectives

In this paper, we have designed a hybrid algorithm (HPSO-SA) that combines the exploration ability of PSO with the exploitation ability of SA, and is capable of preventing premature convergence. Compared with QIPSO, ATREPSO and GMPSO [START_REF] Pant | Particle swarm based metaheuristics for function optimization and engineering applications[END_REF], Tl-PSO [START_REF] Nakano | Consideration of particle swarm optimization combined with tabu search[END_REF], PSO-SA [START_REF] Che N | Ane wh y b ridalgorithmof particle swarm optimization[END_REF] and SUPER-SAPSO [START_REF] Bahrepour | Super-sapso: a new sa-based pso algorithm[END_REF] on well-known benchmark functions and for the problem of reducing energy consumption in embedded systems memories, it has been shown that HPSO-SA performs well in terms of accuracy, convergence rate, stability and robustness. In future, we will also compare the performances of HPSO-SA with the above mentioned algorithms on the embedded systems memory saving problem.

In addition, we will compare HPSO-SA algorithm with other hybrid algorithms (PSO-GA, PSO-MDP, PSO-TS) whose design is in progress by the authors. Comparison will also be done on additional benchmark functions and more complex problems including functions with dimensionality larger than 30.

Figure 1 :

 1 Figure 1: Movement of each particle.

  (d) SA Stop Criterion. The stopping criterion of SA algorithm defines when the system has reached 3000 function evaluations or maximum number of functions evaluations or Optimal solution are not attained. (e) Decrementing the Temperature. The most commonly used temperature reducing function is geometric (see (2)). In next simulations γ = 0.99. (f) Inner-Loop. The length of each temperature level determines the number q = 150 of solutions generated at each temperature, T.

Figure 2 :

 2 Figure 2: 3D mathematical benchmark functions.

Figure 3 :

 3 Figure 3: Energy consumed by benchmarks studied in this work.

Figure 4 :

 4 Figure 4: Execution time used by HPSO-SA algorithms on benchmarks studied in this work.

Table 1 :

 1 Standard benchmark functions adopted in this work.

	Function	Problem	Range	f (x * )	Classification
	Sphere	n i=1 x 2 i	[-100; 100]	0	Unimodal
	Rastrigin Griewank	n i=1 (x 2 i -10 cos(2πx i ) + 10) (1/4000) n i=1 x 2 i -n i=1 cos(x i / √ i)+1	[-5.12; 5.12] [-600; 600]	0 0	Multimodal Multimodal
	Rosenbrock	n-1 i=1 (100(x i+1 -x 2 i ) 2 +(x i -1) 2 )[ -2.048; 2.048]	0	Unimodal
	Quartic	( n i=1 ix 4 i )+rand[0,1]	[-1.28; 1.28]	0	Noisy
	Schwefel	420.9687n -n i=1 (x i sin( |x i |)[ -500.0; 500.0]	0	Multimodal
	Ackley

Table 2 :

 2 Comparison of hybrid HPSO-SA algorithm with TL-PSO approach[START_REF] Nakano | Consideration of particle swarm optimization combined with tabu search[END_REF].

	x 1	[-5.0; 5.0]	-3.78396	Multimodal

Table 3 :

 3 Comparison of mean/standard deviation of solutions obtained by using hybrid HPSO-SA algorithm and approaches described in[START_REF] Pant | Particle swarm based metaheuristics for function optimization and engineering applications[END_REF].

	Function	PSO	QIPSO	ATREPSO	GMPSO	HPSO-SA
	Rastrigin	22.339158 15.932042	11.946888 9.161526	19.425979 14.349046	20.079185 13.700202	0 0
	Sphere Griewank Rosenbrock	1.167749e -45 5.222331e -46 0.031646 0.025322 22.191725 1.615544e +04	0.000000 0.000000 0.011580 0.012850 8.939011 3.106359	4, 000289e -17 0.000246 0.025158 0.028140 19.490820 3.964335e +04	7, 263579e -17 6.188854e -17 0.024462 0.039304 14.159547 4.335439e +04	5.3656e -32 2.98492e -31 3, 32255e -20 2.68415e -20 0, 227048188 0.243978057
	Noisy	8.681602 9.001534	0.451109 0.328623	8.046617 8.862385	7.160675 7.665802	0, 002019998 0.000650347
	Schwefel Ackley Michalewizc	-6178.559896 4.893329e + 02 3.483903e -18 8.359535e -19 -18.159400 1.051050	-6355.586640 477.532584 2.461811e -24 0.014425 -18.469600 0.092966	-6183.677600 469.611104 0.018493 0.014747 -18.982900 0.272579	-6047.670898 482.926738 1.474933e -18 1.153709e -08 -18.399800 0.403722	-8379.66 2.20425e -19 7.43546e -16 1.09382e -15 -19, 62555806 0.00926944
	Himmelblau Shubert	-3.331488 1.243290 -186.730941 1.424154e -05	-3.783961 0.190394 -186.730942 0.000000	-3.751458 0.174460 -186.730941 1.424154e -05	-3.460233 0.457820 -186.730942 1.525879e -05	-3, 78396 3.16001e -15 -186, 730942 8.66746e -14

Table 4 :

 4 Performance results of HPSO-SA, SUPER-SAPSO, SAPSO, and PSO algorithms on benchmark functions.

	Function	Algorithm	Number of iterations	Average Error
		PSO	2989	0.097814
	Rastrigin	SAPSO SUPER-SAPSO	2168 5	0.07877 0.0
		HPSO-SA	2022.43	0.0
		PSO	805	0.094367
	Sphere	SAPSO SUPER-SAPSO	503 4	0.085026 0.0
		HPSO-SA PSO	501.15 2004	7, 92e -32 0.082172
	Griewank	SAPSO SUPER-SAPSO	1517 3	0.075321 0.0
		HPSO-SA PSO	1496.35 4909	3, 61e -20 0.099742
	Ackley	SAPSO SUPER-SAPSO	3041 5	0.099461 0.0
		HPSO-SA	3020.65	8, 92e -16

Table 5 :

 5 Performance comparison between PSOSA, Sid.GA, Hybrid and HPSO-SA for benchmark functions[START_REF] Che N | Ane wh y b ridalgorithmof particle swarm optimization[END_REF].

	Function	Dimension	PSOSA	Std.GA	H ybrid	HPSO-SA
	Sphere Rastrigin Griewank Rosenbrock	10 20 30 10 20 20 10 20 30 10 20 30	0 ± 02 .43e -04 ± 1.14e -05 0 ± 00 .00145 ± 6.22e -05 0 ± 00 .00442 ± 1.78e -04 0 ± 03 .1667 ± 0.2237 0 ± 01 6 .8732 ± 0.6007 0 ± 04 9 .3212 ± 1.1204 0 ± 0 283.251 ± 1.812 0 ± 0 611.266 ± 3.572 0 ± 0 889.537 ± 3.939 0.17856 ± 1.25988 109.810 ± 6.212 0.00043 ± 0.00111 146.912 ± 10.951 0.57431 ± 4.05976 199.730 ± 16.285	2.42e -04 ± 2.17e -05 0.00212 ± 2.75e -04 0.01203 ± 6.33e -04 3.0599 ± 0.1535 11.6590 ± 0.3602 27.8119 ± 0.8059 0.09078 ± 0.03306 0.00459 ± 0.01209 0.09911 ± 0.00106 43.521 ± 16.047 169.112 ± 21.535 187.033 ± 22.960	0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0.146 ± 0.224193 0.246897 ± 0.231982 0.439149 ± 0.304347
		Nicr
	+
		k=1

Table 6 :

 6 List of terms.

	Te r m	M e a n i n g
	E spmr	Energy consumed during a reading from SPM.
	E spmw	Energy consumed during a writing into SPM.
	N spmr

Table 7 :

 7 List of Benchmarks. algorithm that produces a 160-bit message digest for a given input.BitcountCEMiBench Tests the bit manipulation abilities of a processor by counting the number of bits in an array of integers.

	Benchmark	Suite	Description
	ShaCE The secure hash FirCE MiBench SNU-RT Finite impulse response filter (signal processing over a 700 items long sample).
	JfdctintCE AdpcmCE	SNU-RT Mälardalen	Discrete-cosine transformation on 8 × 8pixelblock. Adaptive pulse code modulation algorithm.
	CntCE	Mälardalen	Counts nonnegative numbers in a matrix.
	CompressCE	Mälardalen	Data compression using lzw.
	DjpegCE	Mediabenchs	JPEG decoding.
	GzipCE	Spec 2000	Compression.
	NsichneuCE	Wcet Benchs	Simulate an extended Petri net. Automatically generated code with more than 250 if-statements.
	StatemateCE	Wcet Benchs	Automatically generated code.

Acknowledgments

The authors are grateful to anonymous referees for their pertinent comments and suggestions. Dawood Khan helped the authors with the intricaties of the english language. The work of M. Idrissi Aouad is supported by the French national research agency (ANR) in the Future Architectures program.


Gabriele Facciolo
email: gabriele.facciolo@upf.edu

Rida Sadek
email: rida.sadek@upf.edu

Aurelie Bugeau
email: aurelie.bugeau@labri.fr

Vicent Caselles
email: vicent.caselles@upf.edu

Temporally consistent gradient domain video editing

Keywords: Video editing, Poisson editing, temporal consistency, total variation, convective derivative, numerical methods

In the context of video editing, enforcing spatio-temporal consistency is an important issue. With that purpose, the current variational models for gradient domain video editing include space and time regularization terms. The spatial terms are based on the usual space derivatives, the temporal ones are based on the convective derivative, and both are balanced by a parameter β. However, the usual discretizations of the convective derivative limit the value of β to a certain range, thus limiting these models from achieving their full potential. In this paper, we propose a new numerical scheme to compute the convective derivative, the deblurring convective derivative, which allows us to lift this constraint. Moreover, the proposed scheme introduces less errors than other discretization schemes without adding computational complexity. We use this scheme in the implementation of two gradient domain models for temporally consistent video editing, based on Poisson and total variation type formulations, respectively. We apply these models to three video editing tasks: inpainting correction, object insertion and object removal.

Introduction

In the context of static image editing, the insertion/removal of content in an image is generally performed using gradient domain methods. These methods allow editing an image without introducing artifacts at the boundaries of the edited regions. Consequently, gradient domain methods are widely used in image processing for: seamless cloning and compositing [START_REF] Pérez | Poisson image editing[END_REF][START_REF] Georgiev | Image reconstruction invariant to relighting[END_REF], shadow removal [START_REF] Finlayson | On the removal of shadows from images[END_REF], HDR compression [START_REF] Fattal | Gradient domain high dynamic range compression[END_REF], image inpainting [START_REF] Arias | A variational framework for exemplar-based image inpainting[END_REF][START_REF] Komodakis | Image completion using efficient belief propagation via priority scheduling and dynamic pruning[END_REF][START_REF] Kwatra | Texture optimization for examplebased synthesis[END_REF], and matting [START_REF] Sun | Poisson matting[END_REF] among others. Essentially, gradient domain image editing is based on the manipulation of the gradients of an image instead of its graylevels. The modified gradients are then integrated to recover the resulting image. This procedure prevents the appearance of seams at the boundaries of the edited region. For a more detailed introduction to gradient-domain methods, the reader is referred to [START_REF] Agrawal | Gradient domain manipulation techniques in vision and graphics[END_REF]. In particular, Poisson image editing [START_REF] Pérez | Poisson image editing[END_REF] is one of such techniques; it formulates the problem variationally as min u O⊂Ω ∇u -g 2 dx;

with u| ∂O = u 0 ,

where Ω ⊂ R 2 is the image domain, O ⊂ Ω is the region to be edited, g : O → R 2 is the guidance vector field (e.g. gradient of the image to be composed), u : O → R is the solution which best approximates the field g, and u 0 : Ω → R is the original image which provides the boundary conditions needed for reconstructing the solution u. The solution of this problem is computed solving the Poisson equation with Dirichlet boundary conditions u| ∂O = u 0 . Since a video is nothing but a stack of images captured at evenly spaced times, it is natural to apply image editing techniques for video editing tasks. However, video editing poses some challenges that were not present in still image editing such as preserving the temporal consistency of the video. To cope with these challenges, video editing tasks are generally broken down into three steps: a) tracking the region where the editing is being performed [START_REF] Bai | A geodesic framework for fast interactive image and video segmentation and matting[END_REF][START_REF] Grundmann | Efficient hierarchical graph-based video segmentation[END_REF], a task that may also require detecting occlusions of the region; b) computing the new content of the region; and c) blending the new content with the original video so that the modification is not noticeable. The first two steps are usually handled by tracking and video inpainting techniques.

The focus of this work lies in computing the missing content in the editing region and enforcing the temporal consistency during the blending step. We consider three application scenarios: the correction of artifacts due to illumination changes in inpainted video sequences; the insertion of objects in a video; and the removal of objects. We restrict our study to the case where the inserted or removed object is affixed to a surface (for instance a sign on a wall). However, we do not restrict the motions of the camera nor of the surface in the scene. Our method also copes with occlusions and disocclusions of the inserted/removed object.

Temporal inconsistencies are sometimes more conspicuous than spatial ones. They tend to manifest in the edited video by introducing an annoying flickering. A simple solution to avoid the flickering is to consider the video as a threedimensional volume and manipulate the spatio-temporal gradient to perform the editing operations. In [START_REF] Wang | Videoshop: A new framework for spatiotemporal video editing in gradient domain[END_REF], it is shown that setting the spatio-temporal gradient as the guidance field eliminates the flickering artifact. However, this approach is only valid when the edited region is not moving. That is, if the target region is moving, then the temporal consistency must be enforced by following its movement in the scene. This may be the case of inserting new content in a video, like the advertisement in Figure 1. In [START_REF] Bugeau | Coherent Background Video Inpainting through Kalman Smoothing along Trajectories[END_REF] the temporal consistency is imposed by a Kalman smoothing process applied along the trajectories of the pixels, which are computed by integrating the optical flow. Other variational approaches, such as [START_REF] Bhat | Gradientshop: A gradient-domain optimization framework for image and video filtering[END_REF][START_REF] Bhat | Using photographs to enhance videos of a static scene[END_REF], enforce temporal consistency by adding a spatial regularization term to a temporal one, based on a derivative along the direction of the movement (the convective derivative). Both terms are then balanced by a parameter β > 0. However, due to numerical problems that arise in the discretization of the convective derivative, β is restricted to a certain range of values. This restriction manifests itself in practice by requiring the guidance field to be present in all frames.

In this paper, we propose a novel numerical scheme for computing the convective derivative that allows β to take any desired value and, by that, allowing the model to achieve its full potential. In particular, in some applications, we may avoid computing the guidance field in all frames. Furthermore, the proposed scheme introduces fewer errors than other discretization schemes. We use this scheme in the definition of two gradient domain video editing formulations that use L 2 -and L 1 -norms, which lead to Poisson [START_REF] Pérez | Poisson image editing[END_REF] and total variation type models, respectively. While the L 2 -norm forces smooth solutions, the L 1 -norm allows sharp transitions which may be desirable to better preserve textures.

Let us describe the organization of the paper. Section 2 presents our formulation for gradient domain video editing enforcing temporal consistency. In Section 3 we discuss the discretization of the convective derivative. In Section 4 we present the proposed new discretization scheme. Section 5 shows applications of the presented models to inpainting correction, object insertion, and object removal. The experiments show the ability of our model to handle occlusions automatically. Finally, in Section 6 we include some concluding remarks and future work.

Video editing model (continuous setting)

We propose a functional to perform gradient domain editing in a video while enforcing the temporal consistency. For that we minimize an energy defined on the spatio-temporal volume

Π := Ω × T E p (u) = O⊂Π |∂ v u(x, t)| p + β ∇ x u(x, t) -g(x, t) p dx dt, (1)
where Ω ⊂ R 2 is the rectangular image domain, T = [0, T ] is the temporal domain, O ⊂ Ω × T is the spatio-temporal domain where the editing is performed (Figure 2 illustrates these domains), u : Ω × T → R is a scalar function representing the video, and v : Ω × T → R 2 is a known velocity field obtained from the original video u 0 . We consider p ∈ {1, 2} and β ≥ 0. Equation ( 1) is solved with Dirichlet boundary conditions u| ∂O\∂Π = u 0 at ∂O \ ∂Π and homogeneous Neuman boundary conditions at O ∩ ∂Π.

The second term in (1) is similar to other gradient domain image editing models [START_REF] Pérez | Poisson image editing[END_REF][START_REF] Georgiev | Image reconstruction invariant to relighting[END_REF] and it is responsible for image editing at each frame. The spatial gradient with respect to the two spatial dimensions is denoted by ∇ x and the field g : O → R 2 is the guidance field which determines the new content of the frames.

The first term in (1) imposes the temporal consistency on the video. We denote by ∂ v u(x, t) the convective derivative of u along v (a derivative of the video u along the optical flow field v) which we describe next. Please, note that throughout the text, we will refer to temporal boundary as the set of points (x, t) ∈ ∂O such that the scalar product [v 1 (x, t), v 2 (x, t), 1], n(x, t) = 0, where n(x, t) is the normal to ∂O, and v i (•) denotes the i-th component of v.

Spatial boundary condition

Temporal boundary condition

Fig. 2: Illustration of the spatio-temporal domain.

Temporal consistency and the convective derivative. A particle moving in the real world describes a trajectory s : T → Ω when seen in a video. The velocity field v(s(t), t) = d dt s(t) characterizes the motion of all the particles in the video. Then, given any particle and the associated trajectory s, the temporal consistency means that along s the graylevel of u(s(t), t) is constant, which can be stated as d dt u(s(t), t) = 0.

Applying chain's rule to this equation we obtain

∂ v u(x, t) := ∇ x u(x, t) • v(x, t) + ∂ t u(x, t) = 0, ∀(x, t) ∈ Ω × T. (2)
Here, ∂ t denotes the derivative with respect to time ( ∂ ∂t ). The left hand side of (2) corresponds to the definition of the convective derivative of u, which is nothing else than the directional derivative in the direction of the flow v. Therefore the temporal consistency can be stated in terms of the convective derivative:

∂ v u(x, t) = 0, ∀(x, t) ∈ Ω × T.
Analysis of the energy. In order to highlight the roles of the terms of (1), let us re-write it as

T Ot |∂ v u(x, t)| p dx dt + β T Ot ∇ x u(x, t) -g(x, t) p dx dt, (3)
where O t denotes the slice of O at time t. Observe that the rightmost term is a stack of independent gradient domain image editing problems, each spatial term attaches to its guidance field g(•, t) and to the corresponding spatial boundary condition. The leftmost term enforces the temporal consistency by introducing a dependence between consecutive instants. In particular, it penalizes the changes of intensity along the trajectories specified by v, while attaching to the temporal boundary conditions. The parameter β ≥ 0 determines the mixing between the temporal consistency and the spatial term. When β → ∞, the temporal consistency is irrelevant and the solution is equivalent to solving many independent image editing problems. On the other hand, when β = 0 the functional only enforces the temporal consistency. In this case, the information at the temporal boundary is transported along the trajectories specified by v; however, we are not able to deal with the illumination changes that are incorporated by the space boundary conditions. Moreover, in the case β = 0, if we consider only the boundary condition at time t = 0, then solving (1) amounts to solve the advection equation [START_REF] Georgiev | Image reconstruction invariant to relighting[END_REF]. A common choice of β for enforcing the temporal consistency is β ∼ 1 10 [15]. However, in some cases it is interesting to consider even lower values of β. This could be the case in the absence of the guidance field, here the information needs to be transported from the temporal boundaries. As we will see in Section 3.1 standard discretizations of (1) do not always behave in a satisfactory way when β → 0.

The choice of p ∈ {1, 2} leads to two models with different characteristics of the solution. Let us discuss these two cases.

Case p = 2: Here, the energy (1) is quadratic and its solution is computed by solving the linear system

(∂ * v ∂ v • +βdiv x ∇ x •) u = βdiv x g, (4)
where div x denotes the spatial divergence and 4) is of Poisson-type, and we solve it by using the conjugate gradient method.

∂ * v is given by ∂ * v f = -∂ t f + div x (vf ). Equation (
The solution of (4) smoothly adapts to the boundary conditions of O. Moreover, any error due to inconsistencies between the boundary conditions and the potential field g is smoothly spread across the whole domain O.

Case p = 1: Here, the energy (1) takes the form of a total variation minimization problem. To solve it, we perform an implicit gradient descent:

u j+1 = arg min u E 1 (u) + 1 2λ u -u j 2 , ( 5
)
where λ is a positive number. Each iteration of the gradient descent entails the resolution of a convex problem similar to the total variation model for denoising [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. The solution of ( 5) is computed by solving its dual problem [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF].Defining the dual variables ψ : O → R and ξ : O → R 2 , we perform the fixed point iteration with time step τ ≤ 1/8:

ψ k+1 = ψ k + τ ∂ v [∂ * v ψ k + βdiv x ξ k + u j /λ] 1 + τ |∂ v [∂ * v ψ k + βdiv x ξ k + u j /λ]| , ξ k+1 = ξ k + τ ∇ x [∂ * v ψ k + βdiv x ξ k + u j /λ] + g/λ 1 + τ ∇ x [∂ * v ψ k + βdiv x ξ k + u j /λ] + g/λ ,
and at convergence the solution is recovered as u = u j + λ(∂ * v ψ + βdiv x ξ). This method allows discontinuities of u in O and at its boundary. Its solution attaches to the boundary conditions reducing the effects of the illumination changes, but as opposed to the case p = 2 the transitions may not be smooth. While setting p = 2 favors smooth transitions, the model with p = 1 produces sharp transitions which may be desirable in some circumstances. It also allows a better preservation of textures. Figure 7e shows the distribution of the error in both cases and it also highlights the fact that texture is better preserved (we refer to Section 5 for details).

Remark 1. The formulation of (1) can also easily accommodate spatial weights for controlling the blending as in [START_REF] Tao | Error-tolerant image compositing[END_REF], and the temporal consistency term can be modified as in [START_REF] Bhat | Using photographs to enhance videos of a static scene[END_REF] to keep the reflectance properties of the objects.

Discretization of the model

This section deals with the discretization of (1). The convective derivative is defined using the velocity field v, so we start by commenting on the use of the optical flow as an approximation for v. We argue that the optical flow is suited for imposing the temporal consistency. Then, we study some numerical schemes for computing the convective derivative and their limitations. In next section, we introduce the deblurring convective derivative, a new numerical scheme for computing it.

In what follows we work in a discrete setting. The spatial domain is now a square lattice in Z 2 , Ω = {0, 1, . . . , N } 2 , and the temporal domain is T := {0, 1, . . . , T }; therefore a video is represented as a stack of T + 1 digital images (frames). The spatial gradient is computed using forward differences and is denoted by ∇ +

x .

About the optical flow: As a discrete approximation of the velocity field, we use the optical flow computed on the original video u 0 : Π → R. We define the forward optical flow v + : Ω ×T → R 2 between two frames u(•, t) and u(•, t+1) as a vector field such that u(x, t) and u(x + v + (x, t), t + 1) correspond to the same point in the scene. Similarly, the backward optical flow v -, relates the frame at t with the one at time t -1. We discretize the temporal consistency constraint (2) using the forward optical flow as

∂ + v u(x, t) := û(x + v + (x, t), t + 1) -u(x, t) = 0, (6)
Interpolated pixel forward optical flow Fig. 3: Results with and without pre-processing of the optical flow. From left to right: a diagram describing the pre-processing, the removed flows, the result obtained using the pre-processed flow and without using it.

where û(x+v + (x, t), t+1) is the bilinear interpolation of u(•, t+1) at x+v + (x, t). In a similar way, we can also discretize the temporal consistency constraint using the backward optical flow. The optical flow used in our experiments is obtained with the algorithm described in [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF], but any optical flow algorithm with subpixel precision and regularization (or better, with edge preserving regularization) could be used.

As a last remark, the optical flow is not defined at the occluded/disoccluded areas of a frame. In these cases there may be no correspondence for a pixel in the next frame. While some optical flow algorithms produce occlusion maps, many others do not, so we pre-process the flow to identify the occlusion/disocclusion areas and remove them from the energy. A simple and direct way to achieve that is to compute |û 0 (x + v + (x, t), t + 1) -u 0 (x, t)| > tol, where u 0 is the original video, and tol is a tolerance for the change in the gray levels. Flow vectors that do not satisfy the tolerance criteria are removed from the energy. The discrete energy becomes:

E p (u) = t∈T x∈Ω Occ(x, t) ∂v u(x, t) p + β t∈T x∈Ω ∇ + x u(x, t) -g(x, t) p , (7)
where ∂v denotes a discretization of the convective derivative, for instance ∂ + v u(x, t), and Occ : Ω × T → {0, 1} indicates if a vector v + (x, t) satisfies the tolerance criterion or not. The inclusion of Occ(•, •) implies that occluded/disoccluded pixels are only influenced by the spatial regularization term. Figure 3 illustrates with an experiment the benefits of pre-processing the optical flow.

Discretization of the convective derivative

The discretization of the convective derivative given in (6) corresponds to an implicit upwind scheme with an adaptive stencil [START_REF] Zhou | Explicit Finite Difference Schemes for the Advection Equation[END_REF]. That is, for computing the bilinear interpolation û(x + v(x, t), t + 1) the convective derivative considers stencil points surrounding (x + v I (x, t), t + 1) (where v I (x, t) is the integer part of the vector v + (x, t)), this results in a potentially different stencil for each x. When used for simulating the advection equation, this adaptive scheme allows to achieve stability with time steps beyond the one prescribed by the Courant-Friedrichs-Lewy condition [START_REF] Zhou | Explicit Finite Difference Schemes for the Advection Equation[END_REF]. Columns from left to right show the: original frame, results obtained with v --scheme (explicit), v + -scheme (implicit), and DCD scheme for discretizating the convective derivative (see Section 4).

Using the forward optical flow v + we define the v + -scheme for computing the convective derivative as

∂ + v u(x, t) := û(x + v + (x, t), t + 1) -u(x, t) if t < T, 0 if t = T or (x + v + (x, t), t + 1) / ∈ O,
where û : Z 2 × T → R is a bilinear interpolation of u(•, t + 1) at x + v + (x, t), modified in order to account for the Neumann boundary conditions at ∂(Ω × T).

Similarly, with the backward optical flow v -we define the discrete v --scheme for computing the convective derivative as

∂ - v u(x, t) := u(x, t) -û(x +v -(x, t), t -1) if t > 0, 0 if t = T or (x + v -(x, t), t -1) / ∈ O.
We implement these operators as sparse matrices, which allow us to easily compute their adjoints.

Evaluation: Let us consider the following experiment. We solve [START_REF] Kwatra | Texture optimization for examplebased synthesis[END_REF] with β = 0 (spatial term disabled) and with a Dirichlet boundary condition only at t = 0. Note that the guidance field g is undefined in this experiment. In the solution, we expect the content of the first frame to be transported according to the known optical flow v from time t = 0 to all subsequent frames. This experiment is related to simulating the advection equation ∂ v u = 0 with initial condition at t = 0. When applied in this setting, the v + -and v --schemes for computing the convective derivative correspond to applying explicit or implicit upwind schemes for solving the advection equation, respectively. Using the v --scheme can be seen as applying an explicit scheme, where a filter is applied to u(•, t -1) to obtain u(•, t), while the v + -scheme behaves as an implicit one where an inverse filter relates u(•, t) with u(•, t + 1). Figure 4 shows the results of this transport experiment. It is not surprising to see that for both, the explicit and implicit discretizations, the solutions are completely distorted only a few frames away from t = 0. The reason being that these schemes introduce a numerical diffusion (or oscillations) that is accumulated over time. These artifacts are due to approximation errors in the discretization of the differential operators. In particular, in the case of the above upwind schemes the approximation is only first-order accurate. Artifacts as those seen in Figure 4 are in general never observed, since models like ( 7) are always applied with β > 0 and with a known guidance field g for all the frames [START_REF] Bhat | Gradientshop: A gradient-domain optimization framework for image and video filtering[END_REF]. However, in some circumstances, like in the case of partial occlusions, the guidance field may not be easily computed. Moreover, as we see in Figure 5, even with β > 0 and with a known g, these discretizations induce noticeable errors in the result. We have seen that due to discretization issues we are unable to cope with a small value of β and, therefore, unable to use [START_REF] Kwatra | Texture optimization for examplebased synthesis[END_REF] to its full potential.

The deblurring convective derivative (DCD)

It is interesting to note that, in the experiment shown in Figure 4, the effects of the explicit and implicit schemes for discretizing the convective derivative are somehow opposite. The explicit scheme introduces blurring in the solution, while the implicit scheme sharpens the solution but also introduces oscillations. This motivates the search for a new scheme which allows the convective derivative to preserve the transported information for longer periods of time. The idea of the deblurring convective derivative (DCD for short) is to attain this objective as the balance of two opposing processes: implicit and explicit schemes. That is, applying both the v + -and v --schemes to moderate each other's effects.

Let us analyze in detail the effect of using the explicit upwind scheme (v -scheme) while solving (7) by considering an illustrative example. We consider a similar problem as before, where only the initial frame (t = 0) is known, and we have a constant optical flow v -= [0. problem reduces to 1D). Taking β = 0 and p = 2 in ( 7) and analyzing the first two frames of the sequence we get:

min u x∈Ω ∂ - v u(x, 1) 2 .
It is easy to see that the values of u(•, 1) which minimize this energy are explicitly determined by applying the filter [0.7, 0.3] (coefficients specific for this flow) to the rows of frame u(•, 0). Similarly, frame u(•, 2) is obtained by filtering u(•, 1), and so on. Denoting the filtering operator as M v -, then we can describe this relation as: u(•, t + 1) = M v -u(•, t). In this case, the minimum of the complete energy (for all the frames) is 0, and the solution is increasingly blurry with t since u(•, t) = (M v -) t u(•, 0), which confirms what we observed in Figure 4.

A similar analysis for the implicit upwind scheme (v + -scheme) reveals that the solution of the problem restricted to the first frame, min

u x∈Ω ∂ + v u(x, 0) 2 , satisfies the relation M v + u(•, 1) = u(•, 0), where M v + applies the filter [0.3, 0.7]
to the rows of u(•, 1). Therefore, u(•, 1) is given by the pseudo-inverse of M v + applied to u(•, 0). The repeated application of the pseudo-inverse acts as an inverse smoothing. This sharpening enhances the high frequencies in the solution but also introduces numerical artifacts, as seen in Figure 4. Let us mention that, by solving the linear system using the conjugate gradient method, this problem is mitigated by its regularization effect.

The DCD takes advantage of the fact that these two schemes have opposite effects on the data (one blurs while the other sharpens) and alternates between them. Shortly, if between t = 0 and t = 1 we apply the v + -scheme (implicit), then from frame t = 1 to t = 2 we apply v --scheme (explicit) and so on. The diagrams in Figure 6 shows how the temporal derivatives are taken.

Note that there is no computational overhead in the use of DCD with respect to implicit or explicit schemes. It can also be applied sequentially on pairs of frames [START_REF] Bhat | Gradientshop: A gradient-domain optimization framework for image and video filtering[END_REF] without the need of solving the complete system associated to [START_REF] Kwatra | Texture optimization for examplebased synthesis[END_REF].

In the DCD, a step with the v + -scheme permits to recover the frequencies smoothed by the previous step (v --scheme). However, it may also introduce other high frequencies, which on the long term will build up as high frequency artifacts. These artifacts can then be removed using a low-pass filter as a postprocessing step. In general the DCD preserves the transported information for much longer periods of time without noticeable decay. The performance of the proposed method is shown in Figures 4 and5, where it is compared with the implicit and explicit upwind schemes. Observe that the texture is preserved for more than 70 frames (which is considered as a very long time).

Applications

What we have presented so far can be used for a variety of applications. Here we present some experiments that illustrate the usefulness of these techniques in three applications:

1. Correcting an inpainted domain by making it temporally consistent. 2. Inserting an object on a surface while handling situations where the object is partially occluded/disoccluded. 3. Removing an object from a surface while handling situations where the object is partially occluded/disoccluded.

In what follows we discuss the experimental results for each application. The full set of results can be viewed at http:\\gpi.upf.edu\static\emmcvpr11.

Inpainting correction: In this scenario, we are given an inpainted video sequence. Thus, we can compute the guidance field in all frames. We would like to correct this inpainted sequence so that it becomes spatio-temporally consistent.

In order to achieve that, we use the inpainted video as our input and solve [START_REF] Kwatra | Texture optimization for examplebased synthesis[END_REF] for p = {1, 2} and β = 1. Figure 7a shows the original sequence. Figure 7b shows the resulting inpainted sequence. The sequence has been inpainted with an algorithm that generates a temporally consistent inpainting, yet not consistent with respect to the illumination changes at each frame. Figures 7c and7d show the result obtained by solving [START_REF] Kwatra | Texture optimization for examplebased synthesis[END_REF] with p = 1 and p = 2, respectively. Notice that the inconsistency which is present in the inpainted sequence has been corrected and the result is temporally and spatially consistent. We also observe that in this case, the results of p = 1 and p = 2 are very similar.

Object insertion: As we have pointed out in Section 3, we pre-process the optical flow for detecting occlusions/disocclusions and the tolerance criterion used in the experiments is tol = 50. In this scenario, we wish to insert a twodimensional object into a video sequence and affix it to a surface. Thanks to the DCD, we are able to coherently transport information present in one frame into subsequent (or previous) frames in the video sequence. Basically, we start inserting the object into a chosen frame that indicates the first appearance of the object in the video. Let us call this the first frame. Then, by solving Equation ( 7) with a small value of β (usually 0.01 for p = 1, and 0.0001 for p = 2) and Fig. 9: Example of object removal. We remove the object only in the first frame (by inpainting) and then transport the information to subsequent frames by solving Equation [START_REF] Kwatra | Texture optimization for examplebased synthesis[END_REF] with p = 1. The result for p = 2 is similar. The complete sequence is available at the webpage mentioned at the beginning of Section 5.

setting the first frame as a Dirichlet temporal boundary condition, we are able to transport the first frame to the others. Since we are affixing an object on a surface, the inserted object inherits the optical flow from this surface. The occlusions are handled automatically by the temporal consistency using the preprocessed optical flow and it suffices to insert the object into the first frame. If we want to handle disocclusions, then we need also to insert the object into a later frame in the sequence and set it as another Dirichlet boundary condition. Let us call this the last frame. We note that in this setting we only need the information in the first and the last frames. For the intermediate frames we just have a hole that we fill-in with the new information.

Figure 8 shows an experiment where we insert a poster on a door replacing an existing one. The area where we insert the poster is being occluded and then disoccluded by a moving man. The total number of frames of the sequence is 31. The results shown in Figure 8 are obtained by solving Equation ( 7) with p = {1, 2}. As we can see, our method handles both occlusions and disocclusions. Notice that sometimes, on the boundary of the occluding object, we can see some small inconsistencies due to the inaccuracy of the optical flow. We will address this issue in a future work. Object removal: This object insertion case can be easily adapted to object removal. We can remove an object from the first and last frames and then fill-in the hole in these frames either by inpainting or by copying a zero gradient field using Poisson editing, for instance. Then we can just transport the information from the temporal boundaries to all the domain in-between. Figure 9 shows an example.

Conclusions and future work

In this work we have presented the deblurring convective derivative, a novel scheme for computing the convective derivative that allows us to maintain information for a long period of time without decay. We have also integrated this derivative in a functional for video editing that imposes both spatial and temporal consistency. We applied this formulation to three video editing applications: object insertion, object removal and inpainting correction.

Our solution relies extensively on the quality of the optical flow and we have seen that inconsistencies in the optical flow affect considerably the quality of the results. Our future work will mainly focus on how to detect and handle inconsistent flows in order to improve the results.

Fig. 1 :

 1 Fig. 1: Temporal consistency in gradient domain video editing. From left to right, a video frame with motion vectors superimposed, the result of video editing imposing the temporal consistency disregarding the motion of the scene, and the result of video editing considering the motion of the scene.

Fig. 4 :

 4 Fig. 4: Temporal propagation from frame t = 0, by solving (7) using β = 0 and p = 2.

Fig. 5 :

 5 Fig. 5: Temporal propagation with β = 0.01, horizontal motion 0.175px/frame, g known and p = 2. After 40 frames the results of v --scheme (explicit scheme) or v +scheme (implicit scheme) are more distorted than with DCD. The plots of the evolution of the error w.r.t. the ground truth (right) confirm it.

Fig. 6 :

 6 Fig.6: On the left, we show a diagram depicting the v + -scheme for discretizing the convective derivative. On the right, we show a diagram corresponding to the DCD: it alternates between v -and v + -schemes.

  (a) Original sequence. (b) Inpainted sequence. (c) Corrected inpainting with p = 1 (d) Corrected inpainting with p = 2 (e) See caption.

Fig. 7 :

 7 Fig. 7: Correcting an inpainted domain to impose spatio-temporal consistency using Equation (7) with p = {1, 2} and β = 1. Figure (e), left (resp. right), shows the modulus of the difference of the magnitude of the gradients corresponding to the first column of images 7b and 7c (resp. 7b and 7d). These images have been jointly scaled to take values in [0, 255] (max value = 16). Note that in the case p = 2 the differences are spread across the domain.

( a ) 2 Fig. 8 :

 a28 Fig. 8: Object insertion experiment with handling of occlusions and disocclusions. Figures (a) show some frames of the original sequence, with the editing region superimposed. Figures (b) and (c) show the results obtained after solving[START_REF] Kwatra | Texture optimization for examplebased synthesis[END_REF] with p = 1 and p = 2, respectively. In this experiment the new object is only provided in the first and last frames of the sequence, and the editing region in-between is unknown. A closer inspection of the results corresponding to p = 1 and p = 2 reveals that the model with p = 1 preserves better the texture.

Acknoledgements: We acknowledge support by MICINN project, reference MTM2009-08171, and by GRC reference 2009 SGR 773. VC also acknowledges partial support by IP project "2020 3D Media: Spatial Sound and Vision", financed by EC, and by "ICREA Acadèmia" prize for excellence in research funded both by the Generalitat de Catalunya. We acknowledge P. Arias and E. Meinhardt for their assistance and useful comments.