Salwa Belaqziz
email: s.belaqziz@uca.ma

S Mangiarotti

M Le Page

S Khabba

S Er-Raki

T Agouti

L Drapeau

M H Kharrou

M El Adnani

L Jarlan

Irrigation scheduling of a classical gravity network based on the Covariance Matrix Adaptation -Evolutionary Strategy algorithm

Keywords: Irrigation scheduling, Evolutionary algorithm, optimization, gravity irrigation network, irrigation round

Irrigation scheduling is an important task that significantly influences water conservation and crop production. For most gravity irrigation networks located in semi-arid areas, irrigation scheduling is based on available water regardless of crops water needs. The objective of this study is to propose a new approach for optimizing irrigation scheduling taking into account crops water demand, based on the Covariance Matrix Adaptation -Evolution Strategy (CMA-ES) evolutionary strategy algorithm. The objective function, that must be minimized, of this optimization problem is defined as the sum of two terms. The first term is the Irrigation Priority Index (IPI) which characterizes the degree of imbalance between water stress and irrigation timing of the plot. The second one takes into account the various constraints that relate to canals capacity, tasks timing, geographical distances and canal flow rate variations. The approach was applied to an agricultural sector located at 40 km from the city of Marrakech (Morocco). Optimal schedule for the third irrigation, of the 2011-2012 agricultural season, is provided and the comparisons between schedules before and after optimization are made. The obtained results demonstrate that such approach allows reducing the proportion of late irrigated plots (from 22% to 8%) and increasing the proportion of plots irrigated at an appropriate time (from 28% to 40%). We conclude that this approach can be considered as an efficient tool for planning irrigation schedules by considering crops water needs.

Introduction

Sustainable development has become a world concern major element aiming to use and manage rationally human, natural and economic resources in order to meet population's basic needs. Among natural resources, water is the most important and its sustainable management is a necessity, particularly in arid and semi-arid areas where crops development can be characterized by a limited predictability [START_REF] Mangiarotti | Polynomial search and Global modelling: two algorithms for modeling chaos[END_REF]Jarlan et al., 2014). Semi-arid areas occupy currently more than 30% of the globe [START_REF] Schlesinger | Biological feedbacks in global Desertification[END_REF] and have irregular water resources whilst water demand is increasing. Among the different sectors, irrigated agriculture consumes the largest volumes of water resources. In Morocco, it was pointed that uncontrolled irrigation management leads to an alarming situation (PAPNEEI, 2009). Therefore, it is necessary to develop effective solutions to save water and to adopt efficient strategies for the distribution of this resource in the most reasonable way. To ensure agriculture water use efficiency, it is necessary to improve its profitability by providing an irrigation schedule which optimizes the water supply according to real needs of each crop, taking into account environmental conditions, irrigation network constraints and agronomic conditions [START_REF] Kharrou | Water use efficiency and yield of winter wheat under different irrigation regimes in a semi-arid region[END_REF].

The scheduling and management of irrigation is essential. Several delivery methods are used in irrigated agriculture throughout the world, some of these approaches allocate water to different crops at farm level [START_REF] Santhia | A new planning model for canal scheduling of rotational irrigation[END_REF][START_REF] Nixon | A genetic algorithm for optimizing offfarm irrigation scheduling[END_REF][START_REF] Playan | Modernization and optimization of irrigation systems to increase water productivity[END_REF]. Other studies cover larger scales, involving the development of indices and models [START_REF] Al-Kaisi | Evaluation of irrigation scheduling program and spring wheat yield response in southwestern Colorado[END_REF][START_REF] Shang | Application of a simulation based optimization model for winter wheat irrigation scheduling in North China[END_REF][START_REF] Gontia | Development of crop water stress index of wheat crop for scheduling irrigation using infrared thermometry[END_REF][START_REF] O'shaughnessy | A crop water stress index and time threshold for automatic irrigation scheduling of grain sorghum[END_REF][START_REF] Belaqziz | A new irrigation priority index based on remote sensing data for assessing the networks irrigation scheduling[END_REF]. Specific packages and algorithms related to irrigation management have also been addressed in some delivery scheduling models [START_REF] Hess | A microcomputer scheduling program for supplementary irrigation[END_REF][START_REF] Mohan | Expert system applications in irrigation management: an overview[END_REF]. Most of the irrigation models previously developed are based on simulation techniques and also use optimization or scheduling algorithms. Such approaches were applied to different types of irrigation systems, but do not apply to the present problem, that is irrigation rounds in a gravity irrigation system.

The irrigation scheduling is thus the focus of this study. It consists in defining the sequence of irrigated plots and the associated water amount to be applied in order to satisfy all crops water requirements. We focus on the Haouz plain near Marrakech (Morocco) taken as a typical example of irrigated agriculture in semi-arid areas, and especially in the southern Mediterranean countries. In this region, gravity irrigation is the most widely used irrigation delivery method because it is less complex to implement and less demanding in terms of material resources at the plot scale compared to other systems [START_REF] Taky | Diagnostic des pratiques d'irrigation gravitaire et possibilités d'amélioration dans le Gharb au Maroc[END_REF]. However, the implementation of this technique faces several problems: (1) the large variety of situations relating to the type of crop or ploughing, sowing dates, soil and climatic conditions that are difficult to take into account in the irrigation planning; (2) The constraints related to the irrigation network which have effect on the water delivery timing; (3) The decided numbers of water rotations before the agricultural season begins are based on the total available water without considering rain forecasts, and the water amounts allocated to the plots area without taking into account the type of crops; (4) the timing of irrigation may be based only on visual inspection of the crops by the farmers, which may cause either a water excesses or a water shortages; and (5) the bad irrigation scheduling: indeed, irrigation may sometimes coincides with a rainfall event.

The difficulties of scheduling irrigation may result in crop water stress or, on the contrary, water losses due to excess intake. However, applying the right amount of water at the correct time with a lower cost and lower water losses, taking human and technical constraints into account, remains a complex problem. This is increasingly the case where crops, soils and climatic conditions are contrasted and when the total volume of water is limited due to restrictions or scarcity. However with help of optimization techniques and models, it has been made possible to deal with such situation. Considerable research works have been done to develop mathematical models optimizing irrigation water management for different irrigation systems. [START_REF] Suryavanshi | Optimal operation schedule of irrigation distribution systems[END_REF] for the first time proposed a 0-1 linear programming (LP) model for preparing the optimal operation schedule for irrigation canal outlets later improved by [START_REF] Wang | Improved 0-1 programming model for optimal flow scheduling in irrigation canals[END_REF]; [START_REF] Reddy | Optimal scheduling of irrigation for lateral canals[END_REF] and Anwar and Clark (2001). Other works by [START_REF] Garg | Two level optimization model for lower Indus basin[END_REF] and [START_REF] Tzimopoulos | Irrigation network planning using linear programming[END_REF] were developed to solve a similar problem using linear programming. [START_REF] Naadimuthua | A Heuristic Dynamic Optimization Algorithm for Irrigation Scheduling[END_REF]; [START_REF] Anwar | Irrigation scheduling II: Heuristics approach[END_REF] and [START_REF] Almiñana | WISCHE: A DSS for water irrigation scheduling[END_REF] introduced heuristic solutions which were shown to be more efficient compared with linear programming approaches. [START_REF] Wardlaw | Comparison of Genetic Algorithm and Linear Programming Approaches for Lateral Canal Scheduling[END_REF] applied a genetic algorithm to the problem described by Anwar and Clark (2001) and claimed better solution quality by scheduling supplies as close as possible to the times requested by farmer. Ul [START_REF] Ul Haq | Evaluation of a Genetic Algorithm for the Irrigation Scheduling Problem[END_REF] demonstrated that heuristics such as genetic algorithms is computationally more efficient than the Integer Programming (IP) and appears to have considerable scope as a tool to solve the problem of preparing an irrigation schedule. Several other studies [START_REF] Chen | Management of water resources using improved genetic algorithms[END_REF][START_REF] Srinivasa Raju | Irrigation Planning using Genetic Algorithms[END_REF][START_REF] Mathur | Optimal Operation Scheduling of Irrigation Canals Using Genetic Algorithm[END_REF]2009;Elferchichia et al., 2009;[START_REF] Sharif | Multireservoir systems optimization using genetic algorithms: case study[END_REF][START_REF] Kuo | Decision support for irrigation project planning using a genetic algorithm[END_REF] demonstrated the efficiency and the strength of genetic algorithm approach as an optimization tool to provide good solutions for an irrigation scheduling problem. Paly and Zell (2009) compared five Evolutionary Algorithms (Real Valued Genetic Algorithm, Particle Swarm Optimization, Differential Evolution, and two Evolution Strategy-based Algorithms) on the problem of irrigation optimization and showed that Differential Evolution and Particle Swarm Optimization are able to optimize irrigation schedules achieving results which are extremely close to the theoretical optimum.

A new irrigation scheduling approach is presented in this paper aiming to rationalize the water distribution in a classical gravity network based on an ad hoc index recently introduced to define the irrigation priority [START_REF] Belaqziz | A new irrigation priority index based on remote sensing data for assessing the networks irrigation scheduling[END_REF]. This index requires information about the crop water need which is based and estimated on a spatialized version of the FAO-56 method [START_REF] Allen | Crop Evapotranspiration: Guidelines for computing crop water requirements[END_REF] guided by time series of remote sensing imagery (Simmoneaux et al., 2008;Er-Raki et al., 2010a). The rationalization is formulated as an optimization problem taking into account the main technical and human constraints of the irrigation. It is performed based on an evolutionary algorithm. The approach is applied to an irrigated sector located at 40 km from the city of Marrakech (Morocco) in the eastern part of the Tensift plain. The paper is organized as follows: the second section gives an outline of the study area and the irrigation method currently applied. The spatialized method used to estimate the soil water depletion is also presented. In the third section, the CMA-ES algorithm is introduced and the optimization approach is presented with a detailed description of the objective function formulation. Results are presented and discussed in the fourth section. Finally, the conclusions of this work are drawn in the last section.

Presentation of the case study

Study area

The study area is located 40 km east of the city of Marrakech (Morocco) in the eastern part of the Tensift plain. Centred on the geographical coordinates 7°39'4.16"W and 31°41'33.90"N, this area, named "R3", has a 2800 ha surface and is mainly used for cereal crops (45% in 2011-2012 and the remainder of the irrigated area corresponds to annual crops and tree crops). The climate of the region is semi-arid, characterized by high temperature in summer (38 °C on average, in July) and mild temperature in winter (5 °C on average, in February) (Er-Raki et al., 2010b). The annual precipitation average is about 250 mm, whereas the evaporative demand is about 1500 mm/yr [START_REF] Duchemin | Monitoring wheat phenology and irrigation in Central Morocco: on the use of relationships between evapotranspiration, crops coefficients, leaf area index and remotely-sensed vegetation indices[END_REF].

The irrigation network

The study area is irrigated by a classical gravity network. Water is provided by the Hassan 1 dam located in the neighboring catchment of Oum Erabia through the Rocade canal, feeding two primary canals: P1 on the right bank of the sector and P2 on the left bank. In this work, we focus on the irrigation scheduling of the second primary canal R3P2 presented in Figure 1. Irrigation water is conveyed to each plot by a network of secondary and tertiary open air canals. The area is divided into irrigation blocks generally supplied by a tertiary canal. This later feeds a variable number (from 1 to 10) of blocks. Some of them are directly irrigated by the primary or secondary canals. The plot irrigation is done by flooding method.

The detailed data on irrigation distribution is provided by the Regional Office of Agricultural Development of Haouz (ORMVAH). The volume of water made available and the number of irrigation rounds (from 3 to 7) depends on the available water in early autumn and the volume accumulated during the winter season. The applied volume for irrigation water is calculated based on the duration of canal opening. The number of plots watered can vary at each irrigation round. This number relies principally on the farmer's choice and funds, but also on available water in the dam.

The current irrigation method

In collaboration with three farmers associations, the local ORMVAH centre is in charge of preparing the irrigation schedule according to the sowing dates and water availability in the dam; and managing the irrigation rounds. In practice, a global amount of water is assigned to the irrigation sector at the beginning of the agricultural season. The number of irrigation rounds, irrigation rounds starting dates and the volumes assigned to each round are discussed and fixed by managers and users.

At each round, the farmers receive an amount of water according to the owned area. This amount does not take into account the type of crops and its water requirement, even if some plots are uncultivated. As a consequence, an amount of water may be lost for the crop and the irrigation distribution is not always adapted to the crop needs and water stress level which may have indirect significant impact on growth and yields. At the plots scale, the level of complexity is compounded by the spatial heterogeneity of the sowing dates (Er-Raki et al., 2010a), the size of tilled plots, and the quantity of nitrogen used [START_REF] Hadria | Calibration and validation of the shoot growth module of STICS crop model: application to manage water irrigation in the Haouz plain, Marrakech plain[END_REF]. A high spatial and temporal heterogeneity of irrigation water results from this situation which complexity is not usually considered in water distribution. In such complex system, finding the optimal combination of irrigation starting times is difficult to formulate and obtain.

Spatialized estimates of soil water depletion

The water stress coefficient (K s ) is a key factor that affects the growth, development and crop production. It allows characterizing the water needs of crops by calculating the water amount in the root zone. Based on the dual approach of FAO-56 method [START_REF] Allen | Crop Evapotranspiration: Guidelines for computing crop water requirements[END_REF], the concept is to consider a fraction RA W of readily available water, of the total available soil water at the root zone TA W , to be extracted by the crop without suffering water stress. K s = 1 (no stress) when there is no limitation to water uptake and K s < 1 (severe stress) when the available soil water is a limiting factor. K s is expressed as:

1 RA TA r TA     W W D W K s (1)
where D r is the root zone depletion in (mm), i.e. the water shortage relative to field capacity, when D r ≤ RA W , K s = 1. TA W is expressed by the difference between the water content at field capacity and wilting point:

r WP FC TA ) ( 1000 Z θ θ W   (2)
where FC θ is the water content at field capacity in (m 3 m -3 ), WP θ the water content at wilting point in (m 3 .m -3 ) and Z r the rooting depth in (m) which varies according to the plant development. RA W , the fraction of TA W is expressed by: TA RA W ρ W 

(3) where  is the average fraction of TA W that can be depleted from the root zone before moisture stress occurs.

Before each irrigation round, a map of the plot water stress i K S corresponding to plot i, can be estimated. The SAMIR software (SAtellite Monitoring of IRrigation) [START_REF] Simonneaux | Estimation spatialisée de l'Evapotranspiration des cultures irriguées par télédétection. Application à la gestion de l'Irrigation dans la plaine du Haouz (Marrakech, Maroc)[END_REF] is directly used for this purpose. This software uses satellite images to compute estimates of evapotranspiration and soil water budget over large areas.

The optimization problem

Decision variables

The object of this work is to rationalize water distribution taking into account the crops water demand in regard to available volumes, while respecting some constraints resulting from practical implementation of water distribution. Such problem can be viewed as an optimization problem. Optimization algorithms aim at finding the values of decision variables allowing for optimizing (minimizing or maximizing) a function called the objective function. This function has to be defined cautiously in order to properly characterize the optimization problem. The present optimization problem can be formulated as follows: given a gravity irrigation system subject to times constraints, flow rates constraints, irrigation network constraints and duration constraints of each canal regarding the crops needs, how to obtain a feasible (all the tasks can be performed) and optimal (a rational distribution of water) irrigation schedule?

During an irrigation round, it is assumed that the irrigation duration of a plot does not rely on the opening time, the closing time of each zone thus directly results from the opening time. The problem therefore relies on scheduling the dates of canals opening (the closing date being directly deduced from the crops needs). The opening times i δt can thus be taken as the only decision variables in the optimization process. 26 decision variables have to be considered when the constraints of the irrigation network are taken into account. 116 decision variables are necessary when assuming that irrigation constraints can be relaxed, that is, each plot can be irrigated independently from any others.

Objective function

Any constraint violation causes a schedule to become infeasible. Indeed, by precluding the estimation of the objective function where constraints are not satisfied, strict constraints may impede the exploration of decision variables space. In practice, a penalty method can be introduced to control infeasibility. The penalty method consists in turning a strict constraint into a smoother one to suit the application of evolutionary computing techniques. It is thus chosen to formulate constraints as a penalty rather than a prohibitive constraint, and to include this penalty in the objective function. The objective function F will thus include two parts:

2 1 F F F   (4)
The role of the first term F 1 is to rationalize the water distribution by considering the need of water; an irrigation criterion named Irrigation Priority Index (IPI) presented in [START_REF] Belaqziz | A new irrigation priority index based on remote sensing data for assessing the networks irrigation scheduling[END_REF] is used to satisfy this objective. The role of the second term F 2 is to take the practical constraints into account in a smooth manner in order to guaranty the exploration of decision variables space in an effective way. Based on this formulation the aim of the optimization problem can be stated as follows:

  F N k k x    .. 1 min arg (5) where   N k k x .. 1 
corresponds to the optimized schedule of N opening/closing tasks to be performed. Optimal IPI index opt  and optimal values of opt 1 F and opt 2 F also rely on Eq. 4.

The IPI criterion  is an ad hoc index that characterizes the priority irrigation of each plot at each irrigation round for a specific area [START_REF] Belaqziz | A new irrigation priority index based on remote sensing data for assessing the networks irrigation scheduling[END_REF]. It is ad hoc because it can only refer to the chosen area and does not allow any comparison between different areas. It is also a relative criterion since it allows distinguishing the heterogeneities inside the area of study in terms of irrigation priority. It applies to heterogeneous areas which plants may undergo any level of water stress from none to very high. Once the area of study is chosen, this index can be calculated for each plot i and expressed as a linear function of two terms, one accounting for the water stress of the crop, another one accounting for the irrigation delay after the irrigation round is started. The index of plot i is expressed as follows:

(6) where i K S is the water stress at the beginning of the irrigation round, and min In practice, irrigation scheduling relies on canals and irrigation cannot be scheduled plot by plot individually but canal by canal. To consider the optimization of all the plots, the  index should thus be considered on the average of plot groups irrigated by the same canal.

It should be also noted that the dynamical response to water stress is, by nature, dissymmetrical: under stress conditions, a late irrigation may lead to severe damages for the crop impacting its development and production than an early one. In other words, a moderately negative index may have severe impact on the crops whereas a moderately positive index may just require a moderate over-consumption of water. Such a dissymmetry may be advantageously introduced into the objective function. To do so, the first term of the objective function is expressed as follows:

      N k d k γ γ F 1 2 1 (7)
where k is the canal index, N the total number of canals, and d γ the index shift to take into account for the dissymmetrical response to stress.

T δt K K K K γ i i i     min S max S min S S

Constraints

Five constraints are identified which present a decreasing level of priority. (1) The capacity constraint ensures that supplies can never exceed the total capacity of the canal. (2) The interval constraint ensures that all the irrigation tasks can be scheduled during the irrigation round which dates of earlier possible beginning and later possible ending are fixed.

(3) The overlap constraint ensures that all the practical actions can be applied consistently (all actions can be performed but several actions cannot be performed simultaneously), taking into account (4) the geographical distance between the locations where the actions must be applied and the irrigation time span required for all the plots of a same canal. (5) The daily working time of operators. These constraints are formulated as a weighted sum of functions:

   5 1 2 j j j P r F (8)
where r j is the weight associated to each constraint j, chosen arbitrarily in order to give more or less priority to the constraints; and P j the corresponding penalty. The penalty is expressed as:

∑ 1 , N k k j j υ P   ( 9
)
where N is the total number of tasks to be applied and the k th task is such as 0

,  k j  if the task is satisfying the constraints, 1 ,  k j  otherwise.
In practice, all the tasks are considered at the end of each simulation from which the P j are estimated. More specific descriptions of the constraints are given hereafter.

The flow rate constraint

The flow rate constraint ensures that, at any time, the supply did never exceed the discharge capacity Θ , Q and q of the primary, secondary and tertiary canals respectively (in m 3 s -1 ). The constraint can be expressed as:

    t Θ t Q J j j     1 (10)     t Q t q j J I j i ij   ≤ ∑ , 1 ,  (11)
where i and j denote the tertiary and secondary canals' number, respectively, opened at the same time δt .

Irrigation interval constraint

The irrigation interval constraint ensures that all the tasks were scheduled within the specified irrigation round interval, that is:

E t t i start ≤ ≤  (12) E D t i i ≤   (13)
where start t  refers to the starting time of the irrigation round; i t  the opening time scheduled for canal i; E the total duration allowed for the irrigation round, and D i the total irrigation duration scheduled for the plots supplied by canal i.

Overlap, travel time and time working constraints

An operator cannot perform several tasks simultaneously, considering that each canal opening or closing task requires a time span shorter than 15 min. (traveling time not included). The traveling time required for the operator to travel by moped from one canal gate to another to perform the opening/closing actions obviously depends on the distance between canals which can be significant. The travel time between two gates is estimated linearly considering the spatial distance between two points and assuming a moving speed of 30 km h - 1 . Finally, each task (open or close a canal) is to be scheduled within the specified working time that is between 8h00 and 18h00.

Optimization algorithm

The present optimization problem can be seen as the imbrications of three coupled problems: a basic scheduling, a locally linear optimization problem and a traveling salesman problem. This nesting makes the problem difficult to formulate simply. The Covariance Matrix Adaptation -Evolution Strategy (CMA-ES) algorithm [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF] is used to perform the optimization of such problem. Evolution Strategy algorithms belong to the class of the metaheuristic techniques and distinguish from Genetic algorithm in the adaptation of the algorithm parameters that control the search of optimal solution, as a function of the objective function under exploration. At present, the CMA-ES algorithm is one of the most powerful techniques for optimization of single-objective problems, where the multivariate normal distribution has a mean and a covariance matrix continually updated during the optimization process. Benchmarking of several algorithms has established that CMA-ES is one of the most efficient metaheuristics for dealing with difficult numerical optimization problems [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking bbob-2009[END_REF]. Another reason to prefer such algorithm is that a modification of the constraints may not require a reformulation of the problem, which, in practice, may be of considerable usefulness. This choice may thus also have some operational advantages.

In order to avoid local solutions by providing practical a "a priori" solutions which may bias the optimization approach, it was preferred to simulate the initial parameterization of the search algorithm randomly, and to reiterate the algorithm 10 times. Such ensemble of run will also allow checking that the best results are independent from the initial seed of decision variables which are chosen randomly to initialize the algorithm. The solution corresponding to the lower value of F among final solutions is considered as the optimal solution (see Eq. 5). This search process is summarized in Figure 2. Five input parameters are required to run the optimization algorithm: (1) the total irrigation duration corresponding to each canal. This time span is deduced from the volume of water required for each plot to meet the needs of the crops, as provided by the SAMIR software. (2) The map of water stress before the irrigation round, which is required to estimate the  index, is also obtained from the SAMIR software. (3) The irrigation network, which is used to estimate the traveling times between canal gates. And the parameterization of the CMA-ES optimization itself: the number of iterations (4) has been chosen as 200 and the population (5) size chosen such as λ = 600. Finally (6) the CMA-ES also requires one initial solution associated with an initial variance of the distribution, used to generate the initial population set.

Results and discussion

Focus is put on the third irrigation round of the 2011-2012 agricultural season because crops are well developed and the water needs are significant. This round lasted from 9th to 23th of February 2012. The map of the K S stress coefficient before the irrigation round is displayed in Figure 3a (heterogeneity at higher scale is discussed below). At this date, an important contrast can be observed in the area of study with, broadly speaking, a low level of stress observed in the southern part, a high level in the east and west zones which are very patchy, and a moderate level of stress elsewhere. As a reminder, the achieved scheduling at this time was based on the water amount available into the dam, such that at each irrigation round, farmers receive a water volume according to the plot area. A map of the  index that was practically applied can thus be reconstructed [START_REF] Belaqziz | A new irrigation priority index based on remote sensing data for assessing the networks irrigation scheduling[END_REF]. This map is shown in Figure 3b from which three types of behavior can be distinguished: areas of stress (accounting about 19% of the total surface) for which irrigation was applied too late, that is with  indices lower than -0.2 (in black); areas out of stress (54%) for which irrigation was applied too early, with  indices greater than +0.2 (dark gray), and areas for which irrigation was applied appropriately (27%), that is with  indices between -0.2 and +0.2 (light gray). These results are summarized in Table 1.

By applying the optimization procedure presented above with a dissymmetrical shift set to zero ( 0

 d γ
), a first schedule was obtained. The convergence of the algorithm is illustrated in Figure 4 for one of the runs. A small range of control parameters is explored during the first hundred iterations of the optimization (see Figure 4b). The necessity to enlarge this range is detected by the CMA-ES algorithm, the covariance is then adapted which allows reducing rapidly the objective function (Figure . 4a). The optimization process is then pursued which allows refining the parameter search till the convergence is reached. To evaluate the robustness of the results, it was checked that the randomly chosen initial values of the decision variable (see Figure . 5a) would lead to similar optimized values. The optimal values obtained for the 10 runs are often undistinguishable from one to another (see Figure . 5b, plain grey lines). A sensitivity analysis was also performed in order to test the robustness of the results. Numerous factors may introduce uncertainty in the results. These factors include satellite data error, modeling approximation, imprecise parameterization model, etc. These errors are difficult to characterize and to estimate. The error of the water stress K S resulting from the combinations of these factors was assumed to be equivalent to an error level of 5%. Based on this hypothesis, an additive Gaussian noise of 5% was added to the K S before the optimization was processed. These perturbations led to moderately modify the optimal decision variables (see Figure . 5b, dotted grey lines) confirming the robustness of the results.

All the constraints could be satisfied except for six actions for which the last constraint of working time was not compatible with the working day, leading to a minimum value of F 2 different from 0 (Figure 4a; dark gray). This means that exceptional tasks should be managed, or that a preemptive scheduling for all the tasks may be tried. A map of the  index was reconstructed from this optimized schedule (Figure 3c). The resulting indices clearly illustrate the efficiency of the optimization approach by obtaining values of  close to zero everywhere in the area of study. This efficiency is fully salient by considering the distributions of the  indices by canal, before and after optimization (Figure 6a). The distribution is very dispersed and covers a very large range of the IPI before optimization whereas all the indices are concentrated between -0.2 and +0.2, after.

However, as shown in Figure 7a, an important spatial heterogeneity of water stress is clearly observed between the plots irrigated by the same canal. This heterogeneity results from the choice of crops variety (durum or bread wheat, barley, olives or beet), the sowing dates, the size of tilled plots, the quantity of nitrogen used, and also from the pedology (soil type and thickness) which exhibits slight variations in the area of study [START_REF] Hadria | Calibration and validation of the shoot growth module of STICS crop model: application to manage water irrigation in the Haouz plain, Marrakech plain[END_REF].

Estimates of the  indices can thus also be mapped at the plot scale. Values corresponding to the practical case are plotted in Figure 7b. As expected, the map clearly exhibits a high number of plots (about 50%) under low stress for which an early irrigation was applied ( ). These results are summarized in Table 2. At the plots scale, it is found that the optimization is only marginally able to improve the results (see Figure 6b). This limited improvement does not directly arise from the definition of the γ index which aim is to characterize the heterogeneity in terms of priority. Such difficulty resultsin firstfrom the practical constraints of the irrigation network which do not allow for independent irrigation of plots alimented by the same canal. Indeed, as mentioned before, the optimization approach perfectly allowed for centering the distribution around zero, but with large dispersion resulting from the important heterogeneity of the plot indices which cannot be optimized for a group of plots corresponding to the same canal. This limited improvement also resultsin secondfrom the symmetrical and linear formulation of the IPI criterion, which assumes that the balance between the water stress and the irrigation time can be formulated mathematically in a symmetrical and proportional way. This difficulty was solved by introducing a dissymmetrical shift ). In other words an early irrigation is systematically preferred, which is favorable for a yield increase [START_REF] Belaqziz | A new irrigation priority index based on remote sensing data for assessing the networks irrigation scheduling[END_REF].

This only moderate improvement of IPI values is mainly due to the irrigation constraints, especially those related to irrigation network which does not permit to take the spatial heterogeneity of water stress. If the plots could be irrigated independently, the results would be significantly different. This situation was simulated by performing an optimization assuming that the 116 plots could be irrigated separately, even when depending on the same canal. The obtained results after optimization show that the proportion of irrigated plots in time is obviously improved (89% with ). A comparison of the water quantities that would be consumed when using the present scheduling approach with the amounts which were presently allocated for irrigation was also performed. This comparison could be achieved for each tertiary canal (see Fig. 8). Since based on independent evaluations of the crop needs, applied and estimated water amounts can differ significantly. Although not systematic, a generally lower amount of water was prescribed by the SAMIR tool. While improving water scheduling with regards to crops water need, this comparison also shows that an optimized scheduling may achieve a reduction of more than 25% of consumption of water resources ranging from an amount of 580 280 m 3 to 429 440 m 3 . Such reduction would be considerable.

2 . 0 2 . 0      , see

Conclusion

The rationalization of irrigation in a classical gravity network is considered in this study. The problem is viewed as an optimization problem in which control is ensured by irrigation scheduling. An approach is developed to obtain rationalized solutions accounting for the practical monitoring of irrigation. The ad hoc Irrigation Priority Index is used to balance irrigation time with water stress. Crop water stress is derived from the SAMIR software controlled by satellite images. An evolutionary algorithm, the CMA-ES algorithm, is used to perform the optimization.

The approach is applied to the third irrigation round of the 2011-2012 agricultural season in a gravity irrigation network located at the eastern part of the Tensift Al-Haouz plain.

The approach shows a very good efficiency and allows providing an optimized schedule of the irrigation. Nonetheless, because the practical monitoring of the network does not allow for rescheduling freely the plots irrigation supplied by a same canal, some of the heterogeneity cannot be solved by the approach. To take this heterogeneity into account, a dissymmetrical term is introduced in the objective function in order to rationally foster early irrigation rather than water stress. The approach shows a high efficiency by allowing rationalization of irrigation scheduling, amount of water to be provided, and by avoiding high levels of water stress.

The optimization approach could allow significant decreasing proportion of late irrigated plots (from 22% to 8%) and increasing proportion of plots irrigated at appropriate time (from 28% to 40%). Therefore, we conclude that the proposed approach is computationally efficient and robust to optimize irrigation scheduling in a gravity irrigation network.

In this semi-arid area, water restrictions leads to common appearance of hydric stress, making the use of K s relevant for irrigation scheduling. Nevertheless, further investigations should be performed by considering the root zone depletion D r instead of the water stress coefficient K s . However, this would require the definition of another priority criterion.

An optimized irrigation scheduling may also achieve potentially a significant 25% reduction of water consumption. Such reduction may be investigated under real condition in a future work in order to check its practical interest and, especially, its potential consequences in terms of yields that will, eventually, be the only guarantee of the approach operational efficiency.

Acknowledgements

This research was conducted within the International Joint Laboratory TREMA (Télédétection et Ressources en Eau en Méditerranée semi-Aride, http://trema.ucam.ac.ma/) as well as the RS/09/2011 project financed by the CNRST of Morocco. The authors would like to thank the MISTRALS/SICMED program, CESBIO and the IRD for financial support. CNES is also acknowledged for access to the SPOT5 images through the ISIS program. The staff of the Regional Office of Agricultural Development ORMVAH is also kindly acknowledged for their assistance during the course of this study. The implementation of CMA-ES used here was adapted from the original code by Nikolaus Hansen from TAO team in INRIA. We are grateful to Marc Schoenauer for helpful interactions.

Table caption

Table 1: Obtained results at the aggregated scale (grouped plots). The ranges results from the K s perturbations.

Table 2: Obtained results at the plots scale. The ranges results from the K s perturbations.

  minimum and maximum of the S K map, that is the most and less stressed plots, respectively. i δt is the time delay (in days) between the beginning time of the irrigation round and the irrigation time of plot i and T the irrigation day of the last plot, that is T= max i t δ . Values of i δt range between 0 and T and are chosen as the decision variables of the optimization process. Equation (6) ranges the IPI values between -1 and +1. Extreme values reflect a non rational water monitoring, 1    corresponding to the most stressed plot irrigated the last day of the irrigation round; 1    to the less stressed plot irrigated in first. Contrarily, 0   corresponds to a reasonable tradeoff between level of stress and irrigation delay. Therefore, value of  close to 0 is an indicator of optimal water distribution during an irrigation round.

  of plots under high stress are irrigated lately (

  promote moderately earlier irrigation rather than water stress. By shifting the distribution of the γ indices on the right, this modification allows avoiding strong stresses.Applying this shift, the optimized schedule (see Figure7c) leads to a narrower range of plots irrigated lately (about 8% in time regarding their relative levels of water stress (

Figure captions Figure 1 :

 captions1 Figure captionsFigure1: The hierarchical irrigation network of the studied area. The primary canal R3P2 (thick lines) distribute water to the 217 plots through the secondary (dark lines) and the tertiary canals (thin lines).

Figure 2 :

 2 Figure 2: Flowchart of the optimization procedure.

Figure 3 :

 3 Figure 3: Map of the water stress coefficient Ks averaged by canal areas just before the third irrigation round, February 8th 2012 (a); Map of the  index averaged by canal areas without optimization, same date (b); and

Figure 4 :

 4 Figure 4: Convergence of the function F1 (light gray), F2 (dark gray) and resulting objective function F = F1 + F2 (black) during the optimization process (a). Convergence of the control parameters is also provided (b).

Figure 5 :

 5 Figure 5: Random initial (a) and final optimal solution (b) of the decision variables. The very best solutions corresponding to the minimal objective function F (black bold); the best results obtained for each optimized run based on the randomly chosen initial values of the decision variables presented in (a) (plain grey lines); and the best results for each optimized run under perturbed K S conditions (dashed grey lines).

Figure 6 :

 6 Figure 6: Distribution of the  indexes intervals for canals (a) and for plots (b) when no optimization (solid line) and with the centered (

Figure 7 :

 7 Figure 7: Same as Figure 3, except that Ks (a),  index (b) and opt  (c) index are not averaged by

Figure 8 :

 8 Figure 8: Comparison between the applied and estimated water amounts (in m3) at each of the 26 tertiaries.

  Figure 1

Figure 6 Figure 7 Figure 8

 678 Figure 6

Table 2

 2

	) and the proportion of

Table 1

 1

		aggregated plots
			with optimization
		without optimization
			optimal solution	range
	IPI < -0.2	19%	0%	[0%, 0%]
	-0.2 < IPI < + 0.2	27%	100%	[69%, 100%]
	IPI > + 0.2	54%	0%	[0%, 19%]

Table 2

 2

			Non-aggregated plots
				with optimization
		without optimization	with constraints	without constraints
			Optimal solution	Range	Optimal solution	Range
	IPI < -0.2	22%	8%	[8%, 9%]	0%	[0%, 0%]
	-0.2 < IPI < + 0.2	28%	40%	[39%, 40%]	89%	[64%, 89%]
	IPI > + 0.2	50%	52%	[52%, 53%]	11%	[11%, 36%]


David Picard

Ngoc-Son Vu

Inbar Fijalkow

PHOTOGRAPHIC PAPER TEXTURE CLASSIFICATION USING MODEL DEVIATION OF LOCAL VISUAL DESCRIPTORS

Keywords: Image texture, Image classification, Image texture analysis

This paper investigates the classification of photographic paper textures using visual descriptors. Such classification is called fine grain due to the very low inter-class variability. We propose a novel image representation for photographic paper texture categorization, relying on the incorporation of a powerful local descriptor into an efficient higher-order model deviation where texture is represented by computing statistics on the occurrences of specific local visual patterns. We perform an evaluation on two different challenging datasets of photographic paper textures and show such advanced methods indeed outperforms existing descriptors.

INTRODUCTION

Texture classification is an active research topic in the field of computer vision and image processing, playing a significant role in many applications such as medical image analysis, remote sensing, object recognition, document analysis, environment modeling, contentbased image retrieval and many more. As the demand of such applications increases, texture classification has received considerable attention over the last decades and numerous novel methods have been proposed [START_REF] Oana | Compact representation of bidirectional texture functions[END_REF][START_REF] Permuter | A study of gaussian mixture models of color and texture features for image classification and segmentation[END_REF][START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF][START_REF] Chen | Wld: A robust local image descriptor[END_REF][START_REF] Melendez | Multi-level pixel-based texture classification through efficient prototype selection via normalized cut[END_REF].

Recently, the orderless Bag-of-Words (BoW) approach [START_REF] Lazebnik | A sparse texture representation using local affine regions[END_REF], representing texture images statistically as histograms over a discrete texton dictionary, has proven successful in texture classification tasks. The strength of the BoW for texture classification comes from its two complementary components: the powerful local texture descriptors and global statistical histogram characterization. The former attempts to extract a set of robust and discriminative features from local patches to encode local patterns; while the latter relies on the fact that the texture image is a combination of those local patterns viewed as a random sampling process.

In this paper, we tackle the identification of photographic paper textures with papers having similar physical properties (same sheet, same package, same manufacturer, etc). In such context, we must consider features allowing to discriminate between small local differences in textures. To that end, we investigate the recent advances in BoW approaches, which we argue are able to capture such tiny local variations. More formally, we aim at incorporating a recent powerful local descriptor, Patterns of Oriented Edge Magnitudes (POEM) [START_REF] Vu | Enhanced patterns of oriented edge magnitudes for face recognition and image matching[END_REF], into an efficient second order model deviation method, Vectors of Locally Aggregated Tensors (VLAT) [START_REF] Picard | Efficient image signatures and similarities using tensor products of local descriptors[END_REF], and then evaluate the gain of such approaches compared to global descriptors alone. By experiments carried out on two very challenging photographic paper texture datasets [START_REF] Messier | Automated surface texture classification of inkjet and photographic media[END_REF] with very low interclass variability, we show significant improvements over using the descriptors alone.

The rest of the paper is organized as follows: in the next section we present the related work on texture descriptors and BoW approaches. We then detail in Section 3 the proposed combination of descriptors and signatures. In Section 4, we evaluate these approaches and analyze the results before we conclude.

RELATED WORK

Most of earlier work on texture analysis focused on the development of filter banks and on characterizing the statistical distributions of their responses. Davis [START_REF] Davis | Polarograms -a new tool for image texture analysis[END_REF] exploited polarograms and generalized co-occurrence matrices to obtain rotation invariant statistical features. In [START_REF] Duvernoy | Optical digital processing of directional terrain textures invar-iant under translation, rotation, and change of scale[END_REF], Duvernoy proposed Fourier based descriptors to extract texture feature on the frequency domain. In [START_REF] Eichmann | Topologically invariant texture descriptors[END_REF] Eichmann and Kasparis presented texture descriptors based on line structures extracted by the Hough transform. Although many efforts have been carried out along this direction, the supremacy of filter bank-based descriptors for texture analysis has been challenged by several authors [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF][START_REF] Varma | A statistical approach to material classification using image patch exemplars[END_REF] who showed that it is possible to discriminate between textures using the intensities or differences of pixel within small scale neighborhoods. They demonstrated that despite the global structure of the textures, very good discrimination could be achieved by exploiting the distributions of such pixel neighborhoods. Two particularly important works along these lines are the VZ-Joint classifier [START_REF] Varma | A statistical approach to material classification using image patch exemplars[END_REF] and the Local Binary Pattern (LBP) method [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF]. The efficient texture descriptor LBP may be the preferable choice over the VZ-Joint classifier due to its simplicity.

Due to its computational efficiency and good discriminative property, the LBP descriptor [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] has gained considerable attention, and has already been used in many other applications, e.g., face recognition [START_REF] Ahonen | Face description with local binary patterns: Application to face recognition[END_REF][START_REF] Vu | Exploring patterns of gradient orientations and magnitudes for face recognition[END_REF]. Despite its great success, the conventional LBP operator comes with limitations, such as small spatial support regions, loss of local textural information, and sensitivities to rotation and noise. A lot of effort has been devoted to overcome these limitations, and many different LBP like descriptors have been presented such as the completed LBP model of Guo et al. [START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF] including both the magnitudes of local differences and the pixel intensity itself, the Dominant LBP in [START_REF] Liao | Dominant local binary patterns for texture classification[END_REF], the Local Ternary Patterns of Tan and Triggs [START_REF] Tan | Enhanced local texture feature sets for face recognition under difficult lighting conditions[END_REF], the extended LBP model of Liu et al. in [START_REF] Liu | Extended local binary patterns for texture classification[END_REF], the circular symmetric LBP of Heikkila et al. [START_REF] Heikkila | Description of interest regions with local binary patterns[END_REF], as well as the Patterns of Oriented Edge Magnitudes (POEM) features of Vu and Caplier [START_REF] Vu | Enhanced patterns of oriented edge magnitudes for face recognition and image matching[END_REF].

More recently, several authors proposed a texture model based on visual dictionaries [START_REF] Lazebnik | A sparse texture representation using local affine regions[END_REF][START_REF] Varma | A statistical approach to texture classification from single images[END_REF][START_REF] Varma | A statistical approach to material classification using image patch exemplars[END_REF]. In these approaches, texture are considered to be sampled from a distribution of characteristic visual patterns or textons (named visual words in this context). This distribution is estimated by computing the histogram of occurrences of the visual words, hence the name Bag of Words. The similarity between textures is then measured with standard metric on histograms (e.g., χ 2 distance or Canberra distance). Robust and discriminative local texture descriptors and global statistical histogram characterization In fact, the BoW approach is originally borrowed from the object recognition field [START_REF] Sivic | Video google: A text retrieval approach to object matching in videos[END_REF], where local distinctive features are extracted from regions of interest, quantized in a dictionary (or visual codebook) of visual words. The histogram of occurrences of the visual words is then computed and used as a signature. Recently in object recognition, the BoW approach has been superseded by methods encoding more efficiently the statistical information in the set of local descriptors. In those approaches, the visual codebook is considered as universal model of the descriptors space, and the images are modeled as deviations from this model. Such methods include the popular Fisher Vectors [START_REF] Perronnin | Improving the fisher kernel for large-scale image classification[END_REF] that achieve state of the art performances in object recognition challenges, the first order method [START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF] as well as the higher order deviations using a tensor framework [START_REF] Picard | Improving image similarity with vectors of locally aggregated tensors[END_REF].

PROPOSED METHOD

In this paper, we propose to combine the powerful local descriptor, Patterns of Oriented Edge Magnitudes (POEM) [START_REF] Vu | Enhanced patterns of oriented edge magnitudes for face recognition and image matching[END_REF], with an efficient second order model deviation method, Vectors of Locally Aggregated Tensors (VLAT) [START_REF] Picard | Efficient image signatures and similarities using tensor products of local descriptors[END_REF] for texture classification. For the sake of completeness, this section details these two components of the proposed texture classification pipeline.

POEM feature extraction

The key idea of the POEM algorithm is to characterize object appearance by the relations between the distributions of local gradients of neighboring patches. The POEM feature extraction process consists of the three following steps:

(s1) Gradient computation and orientation quantization: let θ(p) and m(p) denote the gradient orientation and magnitude at pixel p. The gradient orientation at each pixel is evenly discretized over 0-π (unsigned) or 0-2π (signed). A pixel feature is represented as a d-dimensional vector with only one non-null element [ m1(p), ..., md (p)] where d is the number of discretized orientations, the i-th element of this vector takes the original magnitude ( mi(p) = m(p)) if the discretized orientation of the current pixel belongs to this bin, and mi(p) = 0 otherwise (see Fig. 3 where only one non-null vector emits from pixel p).

(s2) Magnitude accumulation: to incorporate information from neighboring pixels (Fig. 3), a local histogram of orientations over all pixels within a local image patch (also called "cell") is computed. As a result, a pixel conveys now richer information, i.e., the distribution of edge direction of a local image structure. At pixel p, the feature vector is [ m1(p), ..., md (p)] where mi(p) = p j ∈C mi(pj) and C refers to cell centered on the considered pixel p.

(s3) Self-similarity calculation: the accumulated magnitudes are encoded using the self-similarity LBP-based operator within a more extended image patch (called "block"). The original LBP operator labels the pixels of an image by thresholding the neighborhood surrounding the pixel with the intensity value of central pixel, and considering the sequence of 8 result bits as a number. The similar procedure is applied on the accumulated magnitudes and across dif- ferent directions to build the POEM features. Formally, at pixel p, for the discretized direction i, the POEM feature is calculated as:

P OEM i L,w,n (p) = n j=1 f ( mi(p) -mi(cj))2 j ( 1
)
where cj are surrounding pixels; L, w refer to the size of blocks and cells, respectively; n is number of pixels surrounding the considered pixel p; and f is defined as: f (x) = 1 if x ≥ τ, f (x) = 0 otherwise where the value τ is slightly larger than zero to provide some stability in uniform regions. Then, the POEM feature set at that pixel p is the concatenation of these d unidirectional POEM values: P OEML,w,n(p) = [P OEM 1 , ..., P OEM d ].

VLAT Signature

The POEM feature extraction is performed on a dense grid basis, using several scale of regions (usually in the order of 5% to 15% of the whole image). The descriptors dri extracted from an image i are stored in an unordered set Ai = {dri}. Using a large set of descriptors A = {dri} sampled from the entire texture collection, we compute a clustering of the descriptors space using the k-means algorithm. For each cluster k, we compute the first and second order moments using the partition C k of A corresponding to the descriptors belonging to cluster k:

µ k = 1 |C k | d ri ∈C k dri (2) Σ k = 1 |C k | d ri ∈C k (dri -µ k )(dri -µ k ) ⊤ (3)
µ k is thus the mean of all descriptors belonging to cluster k and Σ k is its covariance matrix. The set of (µ k , Σ k ) for all k is a model of the whole descriptors space. The set of {µ k } alone is akin to the classical visual words dictionary of the BoW approach.

For each texture i, we use the same clusters to compute the partitions C ki of its set of descriptors Ai. For each cluster, the VLAT aggregation corresponds to the deviation between the texture second order moments centered on µ k and the model covariance matrix Σ k :

τ ki = d kri ∈C ki (d kri -µ k )(d kri -µ k ) ⊤ -Σ k (4)
The VLAT signature is the flattening and concatenation of matrices τ ki for all k. We argue this second order deviation is well suited for fine grain identification, since not only the pattern prototypes are identified as it would in standard BoW approaches, but also the slight variations that might differentiate very closely related textures (e.g., paper from the same manufacturer, but not from the same sheet).

EXPERIMENTS

We evaluate our pipeline on two different photographic paper texture datasets (bw and inkjet) [START_REF] Messier | Automated surface texture classification of inkjet and photographic media[END_REF] containing each 120 images of photographic paper sheets. Each dataset contains 9 categories consisting of ten images with similar but varying physical properties (same sheet, same package, same manufacturer, etc), and 30 random images not related to the other categories. Images from both the bw and inkjet datasets are presented in Figure 1. As we can see, it is very challenging for non-experts to distinguish the categories.

In our evaluation, we compared three different descriptors, namely HOG, LBP and POEM, as well as their aggregation using the VLAT signature model. For each texture, we compute the euclidean distance between its feature and the ones of every other texture. In Figure 2, we show the different distance matrices obtained for all the 6 possible combinations. The ideal distance matrix is block diagonal for the first 9 blocks of 10 rows/columns (each block corresponding to a category), and random for the last 30 rows/columns (corresponding to the distractors). As we can see, the bw set is far more difficult than the inkjet, although it is not obvious simply by looking to the images. Using the VLAT aggregation, we are able to capture more of the block structure expected in the ideal distance matrix.

To obtain more quantitative results, we computed the precision/recall curves associated with these distance matrices. Each Figure 4 shows these precision/recall curves for all possible combinations of descriptors on the inkjet dataset. As we can see, single descriptors are systematically outperformed by their aggregated counterpart. The precision curve has a very slow decrease, which means that most categories are almost perfectly sorted. Figure 5 shows the precision/recall curves for the bw dataset. Contrarily to the inkjet dataset, the results are much lower, which is consistent with the observed distance matrices. The aggregation also outperforms related single descriptors in all tested combinations. It would be worth noting that on these photographic paper databases, the texture LBP descriptor surprisingly performs worst than HOG and POEM. This can be explained by the fact that these images are of very high resolution whereas LBP has too small spatial support regions. In other words, LBP cannot capture global enough information of those images.

To sum up the retrieval performances of each combination, Table 1 shows the mean Average Precision (mAp -area under the precision/recall curve). As we can see, the aggregation allows to boost the performances of the descriptors by a fair 3% to 6% of mAp, which is considerable. The VLAT+POEM combination always obtains the best performance on both datasets, which means that the model by a distribution of patterns of varying intensity is the best among those tested.

CONCLUSION

In this paper, we evaluated several descriptors and aggregation model in the case of fine grain texture classification. More specifically, we tackle the identification of photographic paper textures, where paper with the same physical properties have to be retrieved among a collection of many paper textures. Unlike common global approaches, we propose to model the texture as a distribution of specific visual pattern sampled from a dictionary. In that sense, we proposed to adapt efficient local descriptors with a recent model deviation based approach from the object classification community to our texture identification problem. We performed the comparison of these descriptors with their aggregated counterparts on two photographic paper texture datasets, and show that the aggregation method performs best.

  (a) Images from the bw dataset. (b) Images from the inkjet dataset.

Fig. 1 :

 1 Fig. 1: Images extracted from the texture datasets. The first row contains images from the bw dataset, while the second row contains images from the inkjet dataset. Each image is from a different category.

Fig. 3 :

 3 Fig.3: POEM feature extraction for one pixel p. In (s3), the algorithm considers the length of vectors.

Fig. 2 :

 2 Fig. 2: Distance matrix for the different descriptors and signatures on both datasets: first row for the bw dataset, second row for the inkjet dataset. Dark colors indicate small distances, while lighter colors indicate greater distances.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Precision/Recall curves on the inkjet set.

Table 1 :

 1 mAP for all tested combinations. The first row corresponds to the inkjet dataset, while the second row is for the bw dataset image acts as a query, and the precision/recall is measured by the rank of the other 9 images of the same category within the corresponding sorted row of the distance matrix.

		HOG LBP POEM VLAT+HOG VLAT+LBP VLAT+POEM
	inkjet	63.8	45.7	66.3	69.7	67.8	73.6
	bw	26.9	24.1	26.3	31.0	26.4	31.3