Pierre Mandel
email: pierre.mandel@veolia.com

Marie Maurel

Cyrille Lemoine

Pascal Roche

Dominique Wolbert

How Bromate and Ozone Concentrations can be Modelled at Full Scale based on Lab-Scale Experiments -A Case Study

Keywords: Ozone, Bromate, Natural Organic Matter, Modelling, Hydraulics, Chemical Kinetics, Drinking Water, Hydroxyl Radicals, Ozonation Plant

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Besides disinfection, management of ozonation processes in potable water treatment has to deal with new constraints: on one hand, regulations on disinfection by-products became very restrictive (10 µg.L -1 for bromate ions since 2008 in the EU); on the other hand, the occurrence and fate of micropollutants in aquatic environments represent matters of increasing concern [START_REF] Huber | Oxidation of Pharmaceuticals during Ozonation and Advanced Oxidation Processes[END_REF][START_REF] Mompelat | Occurrence and Fate of Pharmaceutical Products and By-products, from Resource to Drinking Water[END_REF]). Modelling appears to be one possible solution to face the challenge of managing the ozonation process [START_REF] Rietveld | Good modelling practice in drinking water treatment, applied to Weesperkarspel plant of Waternet[END_REF]. Modelling of ozonation units, however, has to adapt to site-specific conditions. Considering natural water, Natural Organic Matter (NOM) interactions with ozone may vary seasonally, when studying the same water source throughout the year; or locally, when studying water samples from different sources. In addition, flow conditions impact for a large part on disinfection efficiencies and bromate formation levels (Roustan et al., 1996a;Do-Quang et al., 1999). Being adaptive, models for industrial ozonation units therefore need to be calibrated. While calibration can rely on data collected on-site -for chemical kinetics [START_REF] Audenaert | Full-scale modelling of an ozone reactor for drinking water treatment[END_REF] or hydraulics [START_REF] Gujer | A stochastic model of an ozonation reactor[END_REF] -it is of advantage to use instead lab-scale experiments or numerical studies (perfectly known process conditions, experimental ease). In the latter case, modellers nevertheless have to deal with upscaling. In both cases, the model is finally validated at full-scale with data that were not used for calibration purposes. Many modelling studies of industrial ozonation units can be found in the literature, but a model validated on a large experimental dataset predicting ozone and bromate concentrations at full-scale remains, to our knowledge, unpublished. Most studies have dealt with disinfection efficiency, especially since the popularisation of Computational Fluid Dynamics (CFD), [START_REF] Murrer | The development of an ozone contact tank simulation model[END_REF][START_REF] Do-Quang | Study of full scale reactor performance for drinking water advanced oxidation treatment[END_REF][START_REF] Cockx | Use of Computational Fluid Dynamics for Simulating Hydrodynamics and Mass Transfer in Industrial Ozonation Towers[END_REF]; fewer studies have dealt with bromate formation. Besides, modelling studies of full-scale ozonation processes often rely on empirical or semi-empirical chemical modelling. As a consequence, guidelines for the use of quasi-mechanistic predictive models integrating chemical pathways for ozone self-decomposition, influence of inorganic carbon, formation of bromate etc. are rarely published. Aiming at practical application, we asked whether a model calibrated at labscale (chemical kinetics) and with numerical data (hydraulics) can be validated with full-scale data (ozone and bromate concentrations). Since hydraulic and chemical phenomena can be completely decoupled, modellers typically develop two models: one for hydraulics and another for chemical kinetics. Moreover, attempts to model ozonation with a unique model -e.g. with statistical models [START_REF] Sohn | Disinfectant decay and disinfection byproducts formation model development: chlorination and ozonation by-products[END_REF] or artificial neural network models [START_REF] Heck | Modeling the effectiveness of ozone as a water disinfectant using an artificial neural network[END_REF] -were not validated at full-scale. On one hand, different modelling frameworks have been developed for hydraulics. From the less to the most detailed, they are: Axial Dispersion Models (ADMs), systematic networks (patterns of ideal reactors), stochastic models (systematic networks where individual particles are followed through the tank) and CFD. On the other hand, different models for the chemistry of aqueous ozone have been developed in controlled laboratory conditions. Ranking these model from the least detailed to the most detailed, they are: correlations, simplified mechanisms (with less than 10 chemical reactions) and quasi-mechanistic. Different combinations of hydraulic and chemical models have been studied to model industrial ozonation units; we present in Figure 1 a selection of papers proposing full-scale modelling. Some studies that were done using CFD on disinfection efficiency were not included as they are only marginally related to the present study, see [START_REF] Dong | A multiscale mass transfer model for gas-solid riser flows: Part II-Sub-grid simulation of ozone decomposition[END_REF][START_REF] Bolaños | Applicability of computational fluid dynamics to simulate ozonation processes[END_REF][START_REF] Wols | Evaluation of different disinfection calculation methods using CFD[END_REF]. The modelling framework designed for this study combines a hydraulic systematic network with a quasi-mechanistic model for chemical kinetics. Aiming at providing operators with an operational modelling procedure, the following constraints were taken into account: (i) to use no or very few full-scale experimental data to calibrate the models; (ii) to explore the models robustness to seasonal variations and changes in operational conditions.

Model for hydraulics

CFD

Stochastic

Systematic

ADM

THIS STUDY

Correlations Simplified Quasi-mechanistic Model for chemical Mechanisms kinetics [START_REF] Zhang | An integrated design approach for improving drinking water ozone disinfection treatment based on computational fluid dynamics[END_REF][START_REF] Bartrand | High resolution experimental studies and numerical analysis of fine bubble ozone disinfection contactors[END_REF]) [START_REF] Gujer | A stochastic model of an ozonation reactor[END_REF]) [START_REF] Audenaert | Full-scale modelling of an ozone reactor for drinking water treatment[END_REF]) [START_REF] Kim | Modeling Cryptosporidium parvum oocyst inactivation and bromate in a flow-through ozone contactor treating natural water[END_REF] Degree of detail Degree of detail [START_REF] Tang | Modeling Cryptosporidium parvum Oocyst inactivation and bromate formation in a full-scale ozone contactor[END_REF]) [START_REF] Siddiqui | Modeling Dissolved Ozone and Bromate Ion Formation in Ozone Contactors[END_REF][START_REF] Do-Quang | Mathematical Modeling of Theoretical Cryptosporidium Inactivation in Full-Scale Ozonation Reactors[END_REF][START_REF] Smeets | Inactivation of Escherichia coli by ozone under bench-scale plug flow and fullscale hydraulic conditions[END_REF] Legend Unpublished yet

Pilot-scale experiments

Material and Methods

Chemical Models

The models used in this study originate from the literature, whether it be for ozone selfdecomposition [START_REF] Mizuno | Development of Ozone Self-Decomposition Model for Engineering Design[END_REF], or inorganic carbon influence [START_REF] Westerhoff | Application of Ozone Decomposition Models[END_REF]. These models are mechanistic pathways with Arrhenius-law temperature dependence. The model for the influence of NOM on ozone and hydroxyl radicals was proposed in [START_REF] Mandel | A modelling procedure for on-site ozonation steps in potable water treatment[END_REF]; three fractions of NOM-reactive entities are distinguished: ozone consumers, hydroxyl radical initiators and hydroxyl radical promoters. Each fraction is then distributed over two species, with acidic or basic character. The model for the influence of NOM on ozone and hydroxyl radicals is parameterised, so that it requires the calibration of 12 parameters before being implemented. More details on the model can be found in [START_REF] Mandel | A modelling procedure for on-site ozonation steps in potable water treatment[END_REF]. For the formation of bromate, a mechanistic pathway has been derived from previously published results. It is graphically depicted in Figure 2 and detailed in Table 1 and2. In total, 65 chemical reactions were considered. Numerous reactions may be involved in bromate formation and the most commonly encountered in published studies are represented in Figure 2. However, some reactions are not likely to occur due to their relative low kinetics or the weak probability for their reactants to be formed. As a result, they were neglected and are shown in grey in Figure 2. These assumptions were confirmed by numerical tests comparing the modelling results with and without the grey reactions of Figure 2 [START_REF] Mandel | Modelling Ozonation Processes for Disinfection By-Product Control in Potable Water Treatment: From Laboratory to Industrial Units[END_REF]. NOM is known to alter bromate formation by lowering the amounts formed, see e.g. [START_REF] Westerhoff | Numerical Kinetic Models for Bromide Oxidation to Bromine and Bromate[END_REF]. However, the presence of NOM can be the principle stumbling block in producing accurate bromate formation models [START_REF] Jarvis | Bromate formation during drinking water treatment: an assessment of modelling approaches[END_REF][START_REF] Westerhoff | Numerical Kinetic Models for Bromide Oxidation to Bromine and Bromate[END_REF][START_REF] Sohn | Disinfectant decay and disinfection byproducts formation model development: chlorination and ozonation by-products[END_REF]. Consequently, there is no unified model for bromate formation in natural waters yet. In order to take into account the effect of NOM on bromate formation, the isothermal kinetics (i.e. the frequency factor) of the last molecular oxidation step leading to bromate was optimised (Figure 2). Most authors consider instead (i) possible reduction of hypobromous acid or hypobromite ions to bromide ions by NOM reactive species [START_REF] Westerhoff | Numerical Kinetic Models for Bromide Oxidation to Bromine and Bromate[END_REF], and/or (ii) substitution reactions leading to organo-bromine, e.g. [START_REF] Siddiqui | Factors Affecting DBP Formation during Ozone-Bromide Reactions[END_REF][START_REF] Huang | Effect of NOM characteristics on brominated organics formation by ozonation[END_REF]. Implementing these reactions may be chemically more accurate, initial concentrations of the reactive NOM species remain nevertheless to be determined. Seeking a model easy to calibrate, we selected one particular reaction, which drives the formation of bromate (results of a sensitivity analysis): the molecular oxidation

2 3 3 2 O BrO O BrO + → + - -
, which is usually considered as unlimitedly fast. A drastic drop in the kinetics of this reaction may yield abnormally high concentrations of BrO 2 -in the system. Such accumulation however should not destabilise other kinetics since BrO 2 -only reacts to form bromate ions (directly or indirectly, via • O 2 Br). A similar approach was proposed by [START_REF] Westerhoff | Numerical Kinetic Models for Bromide Oxidation to Bromine and Bromate[END_REF], with molecular oxidation

2 3 O BrO O Br + → + - -
. The values for kinetic constants are presented in Table 1 andTable 2, for radical and molecular reactions respectively. .s -1 [START_REF] Buxton | The Radiolysis of Aqueous Solutions of Oxybromine Compounds; the Spectra and Reactions of BrO and BrO 2[END_REF])

O H OBr OH BrOH 2 + → + • • 2.10 9 M -1
.s -1 [START_REF] Buxton | The Radiolysis of Aqueous Solutions of Oxybromine Compounds; the Spectra and Reactions of BrO and BrO 2[END_REF])

- • • - + → + HO Br O OH BrO 2 2 2.10 9 M -1
.s -1 [START_REF] Buxton | The Radiolysis of Aqueous Solutions of Oxybromine Compounds; the Spectra and Reactions of BrO and BrO 2[END_REF])

+ - - • + + → + H BrO BrO O H OBr 2 2 2 2 5.10 9 M -1
.s

-1 [START_REF] Buxton | The Radiolysis of Aqueous Solutions of Oxybromine Compounds; the Spectra and Reactions of BrO and BrO 2[END_REF])

+ - • • + → + H BrO OH Br O 3 2 2.10 9 M -1
.s -1 [START_REF] Buxton | Critical Review of Reactions of Hydrated Electrons[END_REF])

- • • - + → + HO Br OH Br 1.1.10 9 M -1 .s -1 (von Gunten & Hoigné, 1996) - • • - → + HOBr OH Br 10 10 M -1 .s -1 (Zehavi & Rabani, 1972) - • - • + → HO Br HOBr 4.2.10 6 s -1 (Zehavi & Rabani, 1972) 2 3 O OBr O Br + → + • • 1.5.10 8 M -1 .s -1 (von Gunten & Hoigné, 1994) - • - • + → + Br OBr BrO Br 4.1.10 9 M -1 .s -1 (Kläning & Wolff, 1985) OBr CO BrO CO • - - - • + → + 2 3 3 4.3.10 7 M -1 .s -1 (Kläning & Wolff, 1985) Br O CO BrO CO 2 2 3 2 3 • - - - • + → + 1.1.10 8 M -1
.s -1 [START_REF] Buxton | The Radiolysis of Aqueous Solutions of Oxybromine Compounds; the Spectra and Reactions of BrO and BrO 2[END_REF])

O BrO O Br + → + - - 1.6.10 2 M -1
.s -1 [START_REF] Haag | Ozonation of bromide-containing waters: kinetics of formation of hypobromous acid and bromate[END_REF])

2 3 2O Br O BrO + → + - - 3.3.10 2 M -1
.s -1 [START_REF] Haag | Ozonation of bromide-containing waters: kinetics of formation of hypobromous acid and bromate[END_REF])

2 2 3 O BrO O BrO + → + - - 10 2 M -1
.s -1 [START_REF] Haag | Ozonation of bromide-containing waters: kinetics of formation of hypobromous acid and bromate[END_REF])

> 10 5 M -1 .s -1 (von Gunten & Hoigné, 1994) 2 3 3 2 O BrO O BrO + → + - - 4.6.10 4 M -1 .s -1
This study, fitted value

BrOH H BrO → + + - 5.10 10 M -1 .s -1 (Westerhoff, 1995) + -+ → H BrO BrOH 7.9.10 1 s -1
Calculated with pK a = 8.8 and protonation reaction kinetics Laboratory Experiments Two separate chemical models were calibrated based on experimental results: (i) a model for the influence of NOM on ozone and hydroxyl radicals and (ii) a model for the influence of NOM on bromate formation. The experiments were done separately, on two different water samples -both originating from the water works studied at full-scale (Annet-sur-Marne). The 500 mL required for bromate concentration measurement was not compatible with the experimental set-up used for characterising the influence of NOM on ozone and hydroxyl radicals (reactor volume=100 mL). Calibration of the model for the influence of NOM on ozone and hydroxyl radicals was based on the results of fourteen experiments. The experiments are presented in Table 3. Five experimental conditions were explored: changes in temperature, pH, ozone dose, nature of NOM -with and without radical-scavenger -and concentration of NOM. Ozone and parachlorobenzoic acid (pCBA) concentrations were measured for seven contact times with ozone (30''; 1 '; 2'30''; 5'; 15'; 35' and 60'). pCBA was chosen as hydroxyl radical "tracer", due to its kinetics. No pCBA was added when the radical-scavenging species tert-butanol had already been added, because pCBA concentration had been found constant in that case (results of previous experiments not shown here). A specially conceived single-phase set-up was used. The set-up is detailed -including chemical analyses -with a calibration procedure for the model for the influence of NOM on ozone and hydroxyl radicals in [START_REF] Mandel | A modelling procedure for on-site ozonation steps in potable water treatment[END_REF]. The water sample studied was withdrawn at the inlet of the ozonation unit of Annetsur-Marne water works (sand-filtered water). The experiments were done in February 2009. The characteristics of the water sample were as follows: pH=7.4; alkalinity (A T )=2.82 meq.L -1 ; TOC=1.3 mgC.L -1 ; SUVA(254 nm)=1.43 L.m -1 .mgC -1 .

(10mM tert-butanol) X X X X X X pH drop (7.4→6.4) X X X X X X X X Large O 3 dose (1.4 mg.L -1 →2.3 mg.L -1 ) X X X X Temperature drop (20°C →5°C) X X X X X NOM dilution (with ultra-pure water) X X X X
Calibration of the model for the formation of bromate ions was based on the experimental results of fifteen experiments: bromate formation was followed under changing initial bromide concentrations: five levels, from 48 µg.L -1 (natural concentration) to 128 µg.L -1 ; and ozone doses: three levels, 0.9 mg.L -1 ; 1.6 mg.L -1 ; 2.3 mg.L -1 . A full factorial approach was chosen for the design of experiments. Contrary to the experiments presented in Table 3, the contact time with ozone was kept constant at 10 min. Classical single-phase batch experiments were done in 500 mL glass bottles. The solution of aqueous ozone injected in the bottles was prepared as in [START_REF] Mandel | A modelling procedure for on-site ozonation steps in potable water treatment[END_REF]. Samples were analysed using ion chromatography and conductimetric detection; quantification limit was 2 µgBrO 3 -.L -1 . The water sample studied was withdrawn at the inlet of the ozonation unit (sand-filtered water). The experiments were done in the first days of July 2009. The characteristics of the water sample were: pH=8.05; A T =4.24 meq.L -1 ; TOC=1.6 mgC.L -1 ; SUVA(254 nm)=1.03 L.m -1 .mgC -1 ; [Br -]=48 µg.L -1 .

Full-scale Experiments at the Ozonation Unit of Annet-sur-Marne Water Works Annet-sur-Marne ozonation unit consists of two tanks linked by a pipe of diameter 800 mm as depicted in Figure 3. The transfer chambers are shown in yellow and the contact chambers in blue. The two linked tanks have approximately the same volume: 235 m 3 and 252 m 3 , respectively for the first and the second tank. Six sampling points were installed along the tank baffles to follow the profiles for ozone and bromate concentrations. The location of each point was chosen in order that the sampling took place in the main flow, where water velocity was high. Sampling points are also represented on Figure 3; for each sampling point, the time indication corresponds to the hydraulic retention time under nominal flow conditions (1600 m 3 .h -1 ).

Figure 3

View of the ozonation unit investigated in this study [START_REF] Mandel | Modelling Ozonation Processes for Disinfection By-Product Control in Potable Water Treatment: From Laboratory to Industrial Units[END_REF] Aqueous ozone concentration was measured with on-line analysers directly connected to the sampling points. Two types of analysers were used: Solu Comp II (Rosemount Analytical, Chanhassen, USA) and Q 45H/64 (ATI, Collegeville, USA). Gaseous ozone concentration in the exhaust gas at Transfer sampling point was measured with a 964 BT analyser (BMT, Berlin, Germany) connected to a blind flange located on the roof of the unit. Bromate concentration was analysed with ion chromatography and conductimetric detection; quantification limit was 1 µgBrO 3 -.L -1 . For each experiment, the water quality was characterised by the following parameters: pH, A T , TOC, SUVA(254 nm), temperature and bromide concentration. Bromide concentrations were measured using ion chromatography and conductimetric detection. pH, A T , UVA(254 nm) and temperature were directly measured at the ozonation unit during the experiments; TOC was analysed by combustion catalytic oxidation and non-dispersive infrared method. Finally, the following operational process conditions were recorded by the Supervisory Control And Data Acquisition (SCADA) system of the plant: water flow rate (Q l ), ozonated air flow rate (Q g ), aqueous ozone concentrations and ozone concentration in the ozonated air ([O 3 ] g ). Operational process conditions have been divided into two types: the main conditions, water flow rate and setpoint for concentration of residual ozone (called "setpoint" herein), could be manually changed via the SCADA system; secondary conditions, ozonated air flow rate and ozone concentration in the ozonated air, were set by the Programmable Logic Controller. An experimental planning was designed, based on the main operational process conditions: 2 levels were defined for the water flow rate and 3 levels for the setpoint. Accordingly, 17 experiments were scheduled from July to October 2009. Moreover, the experiments were done during two experimental periods (the first period lasting from July 21 st to August 14 th 2009; the second period lasting from September 2 nd to October 10 th 2009). Whereas water characteristics did not vary significantly during the study, many operational process conditions were tested. These conditions and the water characteristics are summarised in Table 4.

Results and Discussion

Calibration of the Chemical Models Using Lab-Scale Experiments

The model for the influence of NOM on ozone and hydroxyl radicals was calibrated based on the agreement of experimental (observed) and simulated (predicted) concentrations for ozone and pCBA. The comparison of experimental and simulated concentrations is given in Figure 4 and Figure 5, for ozone and pCBA respectively. The coefficient of determination R 2 and the number of points comprised within experimental uncertainty, i.e. for which prediction does not significantly differ from observation, are displayed on each figure. The quality of the results obtained with Annet-sur-Marne water is comparable with that obtained with other water samples (Méry-sur-Oise and Vitré, see [START_REF] Mandel | A modelling procedure for on-site ozonation steps in potable water treatment[END_REF]). The good agreement observed for all those water samples confirms the soundness of the chemical model. The experimental uncertainty was determined by calculating the standard deviation of concentration measurements for replicate experiments; it takes therefore into account the experimental uncertainty of the whole experimental setup, not only the analytical uncertainty. A constant error, whatever the concentration value, has been hypothesised since the analytical methods employed (colorimetry for ozone, mass detection for pCBA, ion chromatography for bromates) rely on the homogeneity of the variance of the measurements over the concentration range for their own calibration (hypothesis of homoscedasticity). Finally, at lab-scale, the experimental uncertainties were found to be: ±0.1 mgO 3 .L -1 ; ±4 µgpCBA.L -1 and ±3 µgBrO 3 -.L -1 . The values of the 12 parameters of the model describing the influence of NOM on ozone and hydroxyl radicals are given in Table 5. These values should not be considered per se since numerical tests showed that, even leading to excellent simulation results, they may be far from the optimum [START_REF] Mandel | Modelling Ozonation Processes for Disinfection By-Product Control in Potable Water Treatment: From Laboratory to Industrial Units[END_REF]. As a result, physicochemical interpretation based on these values should be avoided. Similar reservations about parameter values are generally expressed by process modellers, see e.g. [START_REF] Jeppsson | Modelling Aspects of Wastewater Treatment Processes[END_REF]. The model for the influence of NOM on bromate formation was calibrated based on the agreement of experimental and simulated bromate concentrations. The results are presented in Figure 6. The bromate concentrations increased with initial bromide concentrations and ozone doses. Note however that the model tends to underestimate the impact of an increase in ozone dose and to overestimate the impact of an increase in initial bromide concentration (compare the slope of the points with the same ozone dose to that of the points with the same initial bromide concentration). The model for bromate formation is found to be robust, given that, for a wide range of experimental conditions, predictions are close to observations.

Figure 6 Comparison of experimental and simulated bromate concentrations, lab-scale

The final value for the kinetic constant of reaction

2 3 3 2 O BrO O BrO + → + - -
was set at 4.6.10 4 M -1 .s -1 at 293K, assuming energy of activation of 42 kJ.mol -1 . This value is slightly smaller than what is observed in NOM-free water (≥10 5 M -1 .s -1 ). Such negative influence of NOM on bromate formation is in good agreement with already reported observations, see e.g. [START_REF] Siddiqui | Factors Affecting DBP Formation during Ozone-Bromide Reactions[END_REF].

Calibration of the Hydraulic Model Using CFD Using the results of a previous CFD study of the ozonation unit, a systematic network was proposed to describe the hydraulics of the unit. A systematic network is a pattern of ideal reactors: for this study, single-phase Plug Flow Reactors (PFRs) and Continuously Stirred Tank Reactors (CSTRs) were used; two-phase CSTRs were also considered (two-phase PFRs were modelled by co-or counter-current cascades of two-phase CSTRs). Residence Time Distributions (RTDs) and turbulence analysis by means of CFD were used to calibrate the network, according to the procedure proposed by [START_REF] De Traversay | Hydrodynamics of ozonation tanks: Definition of systemic models from CFD[END_REF]. The two RTDs, for the CFD study and the systematic network, are drawn on Figure 7. The systematic network was designed in two steps by fitting successively the RTD of the first and the second tanks. The resulting systematic network is schematised in Figure 8; it comprises 46 reactors or blocks of reactors. Validation at Full-Scale (with Lab-Scale Calibration) A classical two-film approach was implemented to model mass transfer, see e.g. (Langlais et al., 1991). Values for the mass transfer coefficient k L a generally range from 5.10 -4 s -1 to 0.02 s -1 for similar reactors or tanks [START_REF] Laplanche | Simulation of Ozone Transfer in Water. Comparison With a Pilot Unit[END_REF][START_REF] Roustan | Modeling Hydrodynamics and Mass Transfer Parameters in a Continuous Ozone Bubble Column[END_REF][START_REF] Siddiqui | Modeling Dissolved Ozone and Bromate Ion Formation in Ozone Contactors[END_REF]. Given that ozone dissolution is mostly kinetics-driven, a large value for the k L a was selected: 0.04 s -1 . Preliminary simulation runs showed that a lower value for the k L a would have resulted in worse predictions for ozone concentrations. Full-scale validation results of the two models are gathered in Figure 9 and Figure 10, for ozone and bromate respectively. As for lab-scale experiments, the experimental uncertainty was determined with replicate experiments on the unit; experimental uncertainties were ±0.09 mgO 3 .L -1 and ±1.9 µgBrO 3 -.L -1 . The 17 experiments (complete dataset) are presented in Figure 9; only the 8 experiments from the first period are presented in Figure 10. Bromate concentrations had namely become considerably lower in the second period; this was assessed by comparing replicate experiments from one period to another: ozone profiles were found identical at ±0.04 mg.L -1 , whereas bromate had decreased by more than 4 µg.L -1 . Obviously, the results of the full-scale validation presented in Figure 9 and Figure 10 are not as good as for lab-scale calibration (Figure 4 to Figure 6). However the major part of the simulated concentrations is in close agreement with experimental measurements: 84 % of the concentrations for ozone, 60% for bromate. Besides, the simulation results show that the model may slightly under-predict ozone concentrations and over-predict bromate formation.

Considering the experiments individually, and comparing experimental and simulated data, one can notice two phenomena: -a constant prediction quality for ozone concentration. Experimentally, ozone decomposition rates were found to be stable, so that the modelling results adequately matched the experimental measurements over both experimental periods; -a systematic over-prediction of bromate concentrations. Already observed at the beginning of the study, the phenomenon was enhanced in the meantime separating the first from the second period. Ozone decomposition and hydroxyl radical generation appear to be stable throughout the year as suggested by the fact that a model calibrated with experimental results obtained in February 2009 remains valid with data collected in July 2009. Moreover, other seasonal comparisons with water samples from Annet-sur-Marne and other water works showed the relative stability of the ozone demand [START_REF] Mandel | Modelling Ozonation Processes for Disinfection By-Product Control in Potable Water Treatment: From Laboratory to Industrial Units[END_REF]. Inversely, bromate formation appears more sensitive to seasonal changes. Already reported in [START_REF] Westerhoff | Numerical Kinetic Models for Bromide Oxidation to Bromine and Bromate[END_REF][START_REF] Kim | Modeling Cryptosporidium parvum oocyst inactivation and bromate in a flow-through ozone contactor treating natural water[END_REF], NOM strongly influences bromate formation. Consequently, the overprediction in bromate concentrations may be explained by changes in the nature of NOM that would have occurred during the time interval between lab-scale calibration for the bromate formation model (first days of July) and full-scale experiments, which began on July 21 st . Indeed, it appears from the recordings of concentration measurements stored in the SCADA system of the water works that bromate formation has a seasonal character: bromate formation levels generally increase by the beginning of May; they reach their maximum at the end of June and decrease till November. This variation can be repeatedly observed and is not solely due to changes in temperature, since July and August are generally the two warmest months in Annet-sur-Marne. Moreover, physicochemical properties or operational process conditions are generally constant (pH, A T , [Br -], TOC, flow rates, ozone dose…). Calibrated in August, it is therefore likely that the model would have given more accurate results (see also the change in the SUVA values: 1.03 L.m -1 for calibration; from 1.17 L.m -1 to 1.53 L.m -1 during the case study). Considering that (i) seasonal variations altered much more predictions for bromate ions than for ozone; (ii) the predictions for ozone concentration were matching observations for the two experimental periods; (iii) no molecule could be used to assess hydroxyl radical concentration; it was decided not to modify the values of the parameters presented in Table 5 (model for the influence of NOM on ozone and hydroxyl radicals).

Validation at Full-Scale (with a Readjusted Model for Bromate Formation) In order to test the possibility of a calibration bias for the bromate formation model, a smaller value for the kinetic constant rate of the molecular oxidation

2 3 3 2 O BrO O BrO + → + - -
was determined (see Table 6). The value was optimised based on the bromate concentrations of the first full-scale experiment (chosen so as to reproduce the use of a simulation software tool by the operators at the water works), then injected in the model and finally used to model and to be validated by the 7 remaining experiments of the first period. A similar approach was applied for the second experimental period. In this case however, the model for bromate formation was calibrated based on a single measurement: the outlet bromate concentration of the first experiment of the second period. Validation examples are given in Figure 11 to Figure 14. In these figures, temporal concentration profiles are given for ozone and bromate, for the different sampling points of the tank and for different operational process conditions. The global results for ozone and bromate concentrations are gathered in Figure 15 and Figure 16, respectively. Second period experiment As it appears from the comparison of the concentration profiles obtained with lab-scale calibration or with on-site readjustment, the formation of bromate and the decomposition of ozone may be simultaneously predicted. This observation confirms that the model for bromate formation may be applied to simulate industrial ozonation units. However, it should periodically be readjusted (this was done here by considering only isolated concentration measurements) on the basis of a regular sampling (monthly for example). The frequency of the sampling / readjustment has to be determined. Regarding ozone decomposition, the model appears to be more stable, but this has to be confirmed over a longer period of testing.

([O 3 ] g =10.2 g.Nm -3 ; Q g =80.2 Nm 3 .h -1 ; Q l =990 m 3 .h -1 ; setpoint=0.2 mg.L -1 ) Figure 13 First period experiment ([O 3 ] g =13.9 g.Nm -3 ; Q g =149 Nm 3 .h -1 ; Q l =1851 m 3 .h -1 ; setpoint=0.4 mg.L -1 ) Figure 14 Second period experiment ([O 3 ] g =13.1 g.Nm -3 ; Q g =81.6 Nm 3 .h -1 ; Q l =1100 m 3 .h -1 ; setpoint=0.3 mg.L -1 )
The modelling procedure (systematic network associated to quasi-mechanistic chemical models) developed in this study presents several advantages over other approaches used for modelling full-scale ozonation units: One possible limitation of this study is the poor quality of the modelling of gaseous ozone concentrations. Experimental gaseous ozone concentration measurements were used to assess the transfer efficiency of the unit. Experimental efficiencies were high: for the four first experiments, a value of 83 % was calculated. The simulations overestimated this parameter: the calculations gave an efficiency of 97 %. This discrepancy may partly arise from the geometry of the unit, since gaseous ozone is sampled in a headspace common to several chambers (chambers D and E in Figure 8) whereas modelling is specific to individual chambers [START_REF] Siddiqui | Modeling Dissolved Ozone and Bromate Ion Formation in Ozone Contactors[END_REF]. This discrepancy may also arise from an overconsumption of ozone in the model. Three main explanations can be given to this overconsumption: -The ozone doses used for calibrating the model for the influence of NOM on ozone and hydroxyl radicals are much larger than those encountered on-site; -

The CFD study used to calibrate the systematic network was calculated for a single-phase problem (no gaseous ozone). Gaseous flow rate may impact on water hydraulics and mass transfer. -Differences in the behaviour of NOM towards ozone may have caused the model not to be exactly calibrated for the experimental conditions of the case study. If the predictions for ozone concentrations are generally well predicted (84% of predictions matching observations) over a practical working domain (0-0.5 mg.L -1 ), the experimental values observed for bromate ions could not be used to fully check the validity of the model. A reduced domain (0-6 µg.L -1 ) could be considered. Still, the amount of data remains appreciable with respect to other full-scale modelling studies.

Conclusion

The modelling framework of this study relies on the combination of a quasi-mechanistic chemical reaction pathway with a hydraulic systematic network. The reaction pathway consists of reaction sets from previous studies for: ozone self-decomposition, influence of inorganic carbon and influence of NOM on ozone and hydroxyl radicals. In addition, a simplified reaction pathway comprising 18 reactions was proposed to model bromate formation. In total, 65 chemical reactions were considered. Based on lab-scale experiments, kinetic parameters were adjusted in order to account for the influence of NOM on ozone and hydroxyl radical (12 parameters) and on bromate concentrations (1 parameter). The modelling results adequately reproduced the phenomena observed at lab-scale (concentration profiles) or generated numerically (RTDs). Chemical models were applied to simulate full-scale experiments. The majority of the modelled concentrations matched measurements (±experimental error): 84% and 60% of the concentrations for ozone and bromates, respectively. Moreover, the model can adapt to different operational process conditions: changes in liquid flow rate, gaseous ozone concentration and gaseous ozone flow rate. The validity of the model was confirmed in the case of bromate ions by a local adjustment, which led to a considerable improvement of the simulation results: more than 90% of the modelled bromate concentrations were found in agreement with measurements. The modelling approach presented in this paper was conceived as a procedure for practical application. Calibration is based on the results of classical lab experiments (chemical models) and of CFD studies (hydraulic model), which are now commonly performed. If a continuous use of the models throughout the year is sought, it is likely that the chemical model may have to be regularly readjusted. This should represent, however, only a limited task (one sampling every two months in this study). Prospects for this work involve the extension of the full-scale study over a longer time period in order to assess the robustness of the model. Moreover, the modelling framework developed in this study may be applied to simulate disinfection efficiencies and micropollutant removal. The possibilities of model parameter readjustments shall also be considered in future work (e.g. sensitivity analysis of the model parameters).

Figure 1

 1 Figure 1Modelling frameworks of previous full-scale studies

Figure 2

 2 Figure 2Bromate formation pathway used in this study

Figure 7

 7 Figure 7RTDs obtained with CFD and with the systematic network developed

Figure 9

 9 Figure 9Comparison of experimental and simulated ozone concentrations, validation at fullscale

  Figure 11First period experiment([O 3 ] g =13.5 g.Nm -3 ; Q g =84.4 Nm 3 .h -1 ; Q l =1481 m 3 .h -1 ; setpoint=0.2 mg.L -1 ) Figure 12 Second period experiment ([O 3 ] g =10.2 g.Nm -3 ; Q g =80.2 Nm

Figure 16

 16 Figure 16Experimental and simulated bromate concentrations, full-scale

Table 4

 4 Summary of the experimental conditions of the case study

			Water characteristics
		pH	A T (meq.L	-1 )	TOC (mg.L -1 )	SUVA (L.m -1 .mgC	-1 )	T (°C)	[Br -] (µg.L -1 )
	Minimum	7.2	3.58		1.5	1.17		16.6	30
	Median	7.5	3.67		1.7	1.35		20.6	32
	Maximum	7.5	3.78		1.8	1.53		22.8	41

Table 6 Kinetic

 6

parameters for bromate formation, first and second experimental periods

Acknowledgements

This research is partially supported by the European Commission in frame of the integrated project TECHNEAU. We thank Ms. Delphine Bourdin and Mr. Huy-Hoang Pham for their technical assistance. We thank Messrs Eric Cablan and Eric Pujol, from Veolia Water, for their cooperation. We thank Dr. N.J. Maximilian Macaluso for refining the English in the manuscript and the referees for their helpful comments.

-

Contrary to studies coupling simple chemical kinetics with accurate hydraulic description by CFD [START_REF] Zhang | An integrated design approach for improving drinking water ozone disinfection treatment based on computational fluid dynamics[END_REF][START_REF] Bartrand | High resolution experimental studies and numerical analysis of fine bubble ozone disinfection contactors[END_REF][START_REF] Wols | Evaluation of different disinfection calculation methods using CFD[END_REF], calculations are quite fast (~1 s) and parameter values can easily be changed, for example water flow rates. The long computational time required for simulation explains for a large part that very few studies consider the effects of water quality or operational conditions on CFD models [START_REF] Zhang | An integrated design approach for improving drinking water ozone disinfection treatment based on computational fluid dynamics[END_REF].

-Hydraulic flow conditions can also be modelled by stochastic models and coupled with chemical models, for example with pseudo first-order ozone decomposition and with disinfection kinetics with lag-phase, as in [START_REF] Gujer | A stochastic model of an ozonation reactor[END_REF]. Such models are well suited to reproduce tank behaviours, both globally -RTDs -and locally -by tracking particles. However, the high computational cost here also prohibits the refinement of chemical models [START_REF] Gujer | A stochastic model of an ozonation reactor[END_REF]. Additionally, this approach has not yet been supported by full-scale experimental validation.

-Some authors couple one-dimensional ADMs with chemical correlations [START_REF] Tang | Modeling Cryptosporidium parvum Oocyst inactivation and bromate formation in a full-scale ozone contactor[END_REF] or with quasi-mechanistic models [START_REF] Kim | Modeling Cryptosporidium parvum oocyst inactivation and bromate in a flow-through ozone contactor treating natural water[END_REF]. ADMs describe non-ideal tubular reactors and can be applied to model bubble columns. Extrapolating to other geometries is difficult [START_REF] Cockx | Use of Computational Fluid Dynamics for Simulating Hydrodynamics and Mass Transfer in Industrial Ozonation Towers[END_REF], but results were promising in terms of predictions for inactivation of Cryptosporidium parvum oocysts and gave interesting insights in bromate formation pathway for NOM-containing water. These results need however to be extended, either by considering larger ozonation units [START_REF] Kim | Modeling Cryptosporidium parvum oocyst inactivation and bromate in a flow-through ozone contactor treating natural water[END_REF], or longer experimental periods [START_REF] Tang | Modeling Cryptosporidium parvum Oocyst inactivation and bromate formation in a full-scale ozone contactor[END_REF].

-Systematic networks were generally coupled with simplified chemical approaches: chemical correlations [START_REF] Siddiqui | Modeling Dissolved Ozone and Bromate Ion Formation in Ozone Contactors[END_REF][START_REF] Do-Quang | Mathematical Modeling of Theoretical Cryptosporidium Inactivation in Full-Scale Ozonation Reactors[END_REF][START_REF] Smeets | Inactivation of Escherichia coli by ozone under bench-scale plug flow and fullscale hydraulic conditions[END_REF] or simple mechanisms [START_REF] Audenaert | Full-scale modelling of an ozone reactor for drinking water treatment[END_REF]. Bromate formation was modelled by multiple linear correlation [START_REF] Siddiqui | Modeling Dissolved Ozone and Bromate Ion Formation in Ozone Contactors[END_REF] or by a single chemical reaction [START_REF] Audenaert | Full-scale modelling of an ozone reactor for drinking water treatment[END_REF]. The results presented at full-scale were promising in terms of predictions for disinfection and ozone concentrations, but remain limited for bromate concentrations.

Figure 15

Experimental and simulated ozone concentrations, validation at full-scale


C Chaux
email: chaux@univ-mlv.fr

M El-Gheche
email: el.gheche@univ-mlv.fr

J Farah

J.-C Pesquet
email: pesquet@univ-mlv.fr

B Pesquet-Popescu

A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation

Keywords: Stereo, disparity, color images, illumination variation, proximity operator, total variation, tight frame, convex optimization, parallel proximal algorithms

Proximal splitting algorithms play a central role in finding the numerical solution of convex optimization problems. This paper addresses the problem of stereo matching of multi-component images by jointly estimating the disparity and the illumination variation. The global formulation being non-convex, the problem is addressed by solving a sequence of convex relaxations. Each convex relaxation is non trivial and involves many constraints aiming at imposing some regularity on the solution. Experiments demonstrate that the method is efficient and provides better results compared with other approaches.

Introduction

Depth estimation has become an increasingly important problem in computer vision as it constitutes a fundamental task in the reconstruction of realistic 3D scenes. A large number of algorithms have been proposed to perform this estimation. The problem being ill-posed, prior information (e.g. smoothness of the desired solution) needs to be taken into account in order to obtain reliable results. In recent years, much progress has been made in stereo matching with the development of i) powerful global combinatorial optimization methods such as graph-cuts [START_REF] Kolmogorov | Computing visual correspondence with occlusions using graph cuts[END_REF][START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] and belief propagation [START_REF] Klaus | Segment-based stereo matching using belief propagation and a self-adapting dis-a) Left image b) Right image c) True disparity d) True illumination e) MAE=1[END_REF][START_REF] Yang | Realtime global stereo matching using hierarchical belief propagation[END_REF]; ii) variational approaches which have proven to be also very effective for globally solving the matching problem [START_REF] Deriche | Optical-flow estimation while preserving its discontinuities: A variational approach[END_REF][START_REF] Miled | Disparity map estimation using a total variation bound[END_REF][START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF].

In a majority of dense disparity estimation techniques, the scene is assumed to be Lambertian. However, in the presence of illumination variation (that often occurs in practice), this assumption is violated. Several techniques have thus been proposed to either derive illumination invariant features or to model the illumination changes. Likewise, Cox et al. [START_REF] Cox | Dynamic histogram warping of image pairs for constant image brightness[END_REF] proposed a dynamic histogram warping that consists of directly matching histogram values and performing a global optimization via dynamic programming. However, this approach is restricted to the case of a spatially invariant relationship between the intensities of the two images. Otherwise, methods based on non-parametric local transforms followed by normalized cross correla-tion (NCC) and rank-transform [START_REF] Zabih | Non-parametric local transforms for computing visual correspondence[END_REF] have been used. In [START_REF] Davis | BRDF invariant stereo using light transport constancy[END_REF], a new invariant for stereo reconstruction called Light Transport Constancy was introduced to formulate a rank constraint and enable correspondence of non-Lambertian surfaces. The work in [START_REF] Gennert | Brightness-based stereo matching[END_REF] shows that the intensities of two corresponding points are related by a spatially varying multiplicative term. Following this model, a stereo algorithm based on a nonconvex cost function was developed. The minimization of the cost function from the associated Euler-Lagrange equations becomes a difficult task. In the same context, based on this multiplicative model, dense disparity map estimation is formulated in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF] as a constrained optimization problem in which a strictly convex quadratic objective function is minimized under various convex constraints. The resulting optimization problem is solved via a parallel block iterative algorithm involving subgradient projections.

In this paper, we present an alternative strategy relying on a parallel proximal optimization approach [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF]. Similarly to the subgradient based methods, we use the multiplicative model by Gennert which allows us to define a global similarity measure which is minimized over an intersection of convex constraint sets. The original contributions of our paper not only consist of the ability to consider multicomponent images in the presence of illumination variation, but also in the flexibility in minimizing various (possibly nonsmooth) convex similarity measures and combining them with various convex constraints.

The remaining of the paper is organized as follows. The notation, background and considered model are presented in section 2. Then, in section 3, we describe the adopted parallel proximal algorithm that allows us to solve the derived convex minimization problem. Some simulation results are shown in section 4 on the MiddleBury data set and other real images. Finally, some conclusions are drawn in section 5.

Stereo matching model

We consider two multicomponent images of the same scene acquired by a stereoscopic camera. The left (resp. right) view is given by the function

I L : R 2 → R K s → I L (s) = (I (1) L (s), I (2)
L (s), ..., I (K) L (s)) ⊤ (1)
(resp. I R : R 2 → R K s → I R (s) = (I (1)
R (s), I (2)
R (s), ..., I

(K) R (s)) ⊤ ).
For every k ∈ {1, . . . , K}, I

(k) R (s) (resp. I (k) L (s)) rep- resents the k-th component of image I R (resp. I L ) at position s ∈ R 2 .
Particular cases of this framework include color images by taking K = 3. Each component then corresponds to one of the color channels of a specific color system (RGB, YUV, ...). Throughout this paper, it will be assumed that I R and I L are differentiable functions.

Corresponding points in the images of the stereo pair tend to have similar values. The pixel at position (x, y) ∈ R 2 in the left image I L corresponds to a pixel at position (x ′ , y ′ ) ∈ R 2 in the right image I R , the disparity between these pixels being equal to (xx ′ , yy ′ ). When stereo images are rectified [START_REF] Fusiello | A compact algorithm for rectification of stereo pairs[END_REF], the vertical component of the disparity vector vanishes. The disparity thus reduces to u(x, y) = xx ′ . Finding for each pixel (x, y) in the left image the corresponding pixel (xu(x, y), y) in the right image constitutes the goal of stereo matching. The problem is then equivalent to finding a disparity field u : R 2 → [0, +∞[ which minimizes a similarity measure. However, under varying illumination conditions, the corresponding points in a stereo pair do not have the same component values. The illumination variation model we will employ is grounded on the work by Gennert [START_REF] Gennert | Brightness-based stereo matching[END_REF] and takes the following form:

(∀(x, y) ∈ R 2 ) I R (x -u(x, y), y) ≃ v(x, y)I L (x, y) (2)
where v : R 2 → [0, +∞[ represents the illumination field. The problem now is not only to estimate the disparity u but also to estimate the illumination field v. The next section aims at better formulating the considered problem.

Problem formulation

In what follows, we propose a variational approach to jointly estimate the disparity u and the illumination variation v. Based on Model (2), we can formulate the stereo matching problem as the minimization of a similarity measure J:

J(u, v) = K k=1 s=(x,y)∈A\O φ (k) v(s)I (k) L (s) -I (k) R (x -u(s), y) (3)
where, for every k ∈ {1, . . . , K}, φ (k) is assumed to belong to Γ 0 (R), the class of proper lower semi-continuous convex functions from R to ] -∞, +∞]. In addition, A ⊂ Z 2 is the considered finite image domain, u = (u(s)) s∈A , and v = (v(s)) s∈A . The occlusion areas denoted by O correspond to pixels only visible from one view of the stereo pair, which should not be taken into account in the computation of the similarity measure. Consequently, they have been discarded in the expression of the cost function.

In this paper, we will be mainly concerned with convex optimization approaches. Unfortunately, J is nonconvex with respect to the variable u. To tackle this difficulty, we perform the first-order Taylor expansion of the disparity compensated right image around an initial value of the disparity ū. The latter can be derived from a rough estimation, e.g. by a block-based correlation method [START_REF] Tsai | The evaluation of normalized cross correlations for defect detection[END_REF]. In practice, we can iteratively update our initial value in order to make our final solution weakly dependent on it. When the magnitude of the difference of the fields u and ū = (ū(s)) s∈A is small enough, we obtain the following expression: for every k ∈ {1, . . . , K} and s = (x, y) ∈ A,

I (k) R (x -u(s), y) ≃ I (k) R (x -ū(s), y) -(u(s) -ū(s))∇ (1) I (k) R (x -ū(s), y) (4)
where ∇ (1)

I (k)
R denotes the horizontal gradient of the k-th component of the right image. As a consequence of (3) and (4), we can approximate the cost function by:

J(u, v) = K k=1 s∈A\O φ (k) T (k) 1 (s)u(s) + T (k) 2 (s)v(s) -r (k) (s) (5)
where, for every k ∈ {1, . . . , K} and s = (x, y) ∈ A,

     T (k) 1 (s) = ∇ (1) I (k) R (x -ū(s), y) T (k) 2 (s) = I (k) L (s) r (k) (s) = I (k) R (x -ū(s), y) + ū(s)T (k) 1 (s). (6)
As previously mentioned, our objective is to jointly estimate u and v. Thus, by defining w = (u, v),

(∀s ∈ A) w(s) = u(s) v(s) (7)
and, for every k ∈ {1, . . . , K} and s

∈ R 2 , T (k) (s) = [T (k) 1 (s), T (k)
2 (s)], the above expression can be reexpressed more concisely as

J(w) = K k=1 J (k) (w) (8)
where, for every k ∈ {1, . . . , K},

J (k) (w) = s∈A\O φ (k) T (k) (s)w(s) -r (k) (s) . (9)
Despite the convexity of the function J, optimizing this criterion is an ill-posed problem, since we have two variables to estimate at each point s and the components of T (k) (s) may locally vanish for some k and s. We thus need to incorporate additional prior information on the desired disparity and illumination variation fields.

Introducing prior information

Our objective here is to introduce prior information on our target solution. This can be done either by adding to J some regularization terms or by incorporating some convex constraints to the problem. In this paper, we will follow the second approach. Let (S i ) 1≤i≤m denote the m nonempty closed convex sets modelling the constraints one wants to impose. The resulting optimization problem can then be formulated as:

Find w ∈ S = m i=1 S i such that J(w) = inf J(S). (10)
Constraint sets defined on the Hilbert space H = R |A| × R |A| can be described as lower level sets:

(∀i ∈ {1, ..., m}) S i = {w ∈ H | f i (w) ≤ δ i } , (11)
where, for every i ∈ {1, . . . , m}, f i : H → R is a convex function and δ i ∈ R. One of the potential advantages of this approach with respect to the regularization formulation is that it appears often easier from a physical viewpoint to set the constraint bounds than determining regularization parameters.

We will now review some of the constraints that can be applied to our problem.

Constraints that can be applied on the disparity

Range values As the minimum and maximum values u min and u max of the disparity field are often known, we propose to introduce a constraint related to these bounds. The corresponding constraint can be expressed as

S 1,1 = {(u(s)) s∈A ∈ R |A| | (∀s ∈ A) u min ≤ u(s) ≤ u max }. (12)
It can be noticed that the disparity is always nonnegative, so that u min ≥ 0.

First-order smoothness constraint A piecewise constant behaviour of the disparity field is often expected. So, it is meaningful to introduce a constraint enforcing smoothness in homogeneous areas while preserving edges. One can think of using the Total Variation semi-norm which has played a key role in image recovery problems [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF]. Let ∇ (1) and ∇ (2) denote discrete horizontal and vertical gradients (for example, obtained by cyclic convolutions [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF][START_REF] Pustelnik | Parallel ProXimal algorithm for image restoration using hybrid regularization[END_REF]). Then, a discrete version of the total variation (TV) is the following one:

(∀u ∈ R |A| ) TV(u) = s∈A | ∇ (1) u(s)| 2 + | ∇ (2) u(s)| 2 . (13)
The associated constraint set is

S 1,2 = u ∈ R |A| TV(u) ≤ τ 2 (14)
with τ 2 > 0.

Frame analysis constraint Alternatively, we can adopt a frame analysis approach to construct a smoothness constraint. In particular, it well-known [START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Chambolle | Nonlinear wavelet image processing: variational problems, compression and noise removal through wavelet shrinkage[END_REF] that wavelet frames, through their relationships with Besov spaces, constitute appropriate tools to characterize useful classes of regular signals. Frame representations [START_REF] Han | Frames, bases, and group representations[END_REF] and more precisely tight frame representations have become very popular during the last decade. Such transforms can be described by an analysis frame operator F : R |A| → R Q with Q ≥ |A|, which associates to u ∈ R |A| its frame coefficients (F u) q 1≤q≤Q . The considered frame analysis constraint is expressed as

S ′ 1,2 = u ∈ R |A| Q q=1 η q |(F u) q | ≤ τ ′ 2 (15)
where

(η q ) 1≤q≤Q ∈ [0, +∞[ Q and τ ′ 2 > 0.
F is said to be a tight frame when F ⊤ F = νI, where ν > 0. A simple example of a tight frame is the union of ν orthonormal wavelet bases [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

Second-order constraint Some recent works [START_REF] Lefkimmiatis | Hessian-based norm regularization for image restoration with biomedical applications[END_REF] have shown that considering the second-order derivatives of the target image can be helpful in data recovery problems. Let a discrete version of the Hessian of u ∈ R |A| at s ∈ A be

∇ 2 u(s) = ∇ (1,1) u(s) ∇ (1,2) u(s) ∇ (1,2) u(s) ∇ (2,2) u(s) .
Consequently, a second-order variant of the discrete total variation can be defined as

(∀u ∈ R |A| ) TV 2 (u) = s∈A | ∇ (1,1) u(s)| 2 + | ∇ 2,2) u(s)| 2 + 2| ∇ (1,2) u(s)| 2 (16)
and the associated constraint set is

S 1,3 = u ∈ R |A| TV 2 (u) ≤ τ 3 ( 17
)
with τ 3 > 0.

Constraints that can be applied on the illumination field

Range values Similarly to the disparity field, estimations of the minimum and maximum values v min ≥ 0 and v max of the illumination field are often available.

The corresponding constraint can be expressed as

S 2,1 = { v(s) s∈A ∈ R |A| | (∀s ∈ A) v min ≤ v(s) ≤ v max }. (18)
First-order smoothness constraint Since the illumination field is usually smoothly varying, we can add a quadratic constraint on the discrete gradient of the illumination field. This amounts to considering the following constraint set:

S 2,2 = v ∈ R |A| ∇v 2 ℓ2 ≤ κ 2 (19)
where κ 2 > 0 and

(∀v ∈ R |A| ) ∇v ℓ 2 = s∈A | ∇ (1) v(s)| 2 +| ∇ (2) v(s)| 2 1/2 . ( 20
)
Second-order constraint We can also introduce a constraint on the discrete Hessian of the illumination field, thus yielding

S 2,3 = v ∈ R |A| ∇ 2 v 2 ℓ2 ≤ κ 3 ( 21
)
where κ 3 > 0 and, for every v ∈ R |A| ,

∇ 2 v ℓ 2 = s∈A | ∇ (1,1) v(s)| 2 + | ∇ (2,2) v(s)| 2 + 2| ∇ (1,2) v(s)| 2 1/2 . (22)

Resulting constraints applied on vector w

Based on the constraints presented above, four different constraint sets are considered in this work:

1. The first one concerns the range values taken by the disparity and the illumination fields. It is defined as

S 1 = S 1,1 × S 2,1 . (23)
2. The second set is based on the first-order regularity of the illumination and disparity fields:

S 2 = S 1,2 × S 2,2 (24)
3. The third set involves the frame analysis formulation instead of the total variation for the disparity, and reads:

S ′ 2 = S ′ 1,2 × S 2,2 . (25)
4. The last constraint set models the second-order smoothness constraints:

S 3 = S 1,3 × S 2,3 . (26)
3 Proximal approaches for convex optimization

Proximity operator

In a seminal paper [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF], Moreau proposed an extension of the notion of projection by considering the following optimization problem:

minimize z∈H ϕ(z) + 1 2 z -y 2 (27)
where H is a real Hilbert space, ϕ ∈ Γ 0 (H) and y ∈ H. This problem admits a unique solution which is called the proximity operator prox ϕ y of ϕ at y. Consequently, by taking ϕ = ι C , where the indicator function ι C of a nonempty closed convex set C ⊂ H is defined as

(∀z ∈ H) ι C (z) = 0 if z ∈ C +∞ otherwise, ( 28
)
prox ιC reduces to the projection operator P C onto C [START_REF] Combettes | The foundations of set theoretic estimation[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Youla | Image restoration by the method of convex projections: Part I -theory[END_REF]. Proximity operators have very attractive properties [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] that make them particularly well-suited to design iterative minimization algorithms.

Optimization background

The constraints proposed in section 2.2.3 are separable since each constraint set S i with i ∈ {1, . . . , m} can be expressed as S 1,i × S 2,i . 1 In addition, for every i ∈ {1, . . . , m}, one can express S 1,i (resp. S 2,i ) as

L -1 1,i (C 1,i ) (resp. L -1 2,i (C 2,i )) where C 1,i (resp. C 2,i
) is a nonempty closed convex subset of R N1,i (resp. R N2,i ) and L 1,i (resp. L 2,i ) is a matrix in R N1,i×|A| (resp. R N2,i×|A| ). For example, we have chosen L 1,1 and L 2,1 equal to the |A| × |A| identity matrix. By setting L 1,2 and L 2,2 equal to the concatenation of the horizontal and vertical gradient operators (N 1,2 = N 2,2 = 2|A|), C 1,2 and C 2,2 reduce to ℓ 2,1 balls [START_REF] Van Den | Probing the Pareto frontier for basis pursuit solutions[END_REF][START_REF] Eldar | Robust recovery of signals from a structured union of subspaces[END_REF]. A similar choice can be made for constraint S 3 by considering linear operators related to the second-order derivatives (N 1,3 = N 2,3 = 3|A|). In the considered cases, the projection onto each convex set C 1,i or C 2,i takes a closed form or it can be computed in a finite number of operations [START_REF] Van Den | Probing the Pareto frontier for basis pursuit solutions[END_REF]. The considered optimization problem then takes the following generic form:

minimize L 1,i u ∈ C 1,i , L 2,i v ∈ C 2,i , i ∈ {1, . . . , m} J(w). (29)
Numerical solutions to this problem can be provided by parallel proximal splitting methods. These methods basically consist of iterating computations of proximity operators.

PPXA+ algorithm

In this work, we employ the Parallel Proximal Algorithm (PPXA+) [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF]. This algorithm is closely related to augmented Lagrangian approaches [START_REF] Afonso | An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems[END_REF][START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF]. The existing links have been investigated in [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF]Section 6] where it is shown that the parallel versions of the alternating direction method of multipliers in [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF][START_REF] Attouch | Augmented Lagrangian and proximal alternating direction methods of multipliers in Hilbert spaces -applications to games, PDE's and control[END_REF] constitute a particular case of PPXA+. Alternative algorithms could be employed in our context [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Esser | A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science[END_REF][START_REF] Combettes | Proximity for sums of composite functions[END_REF][START_REF] Briceño-Arias | A monotone + skew splitting model for composite monotone inclusions in duality[END_REF][START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, lipschitzian, and parallel-sum type monotone operators[END_REF]. However, PPXA+ was observed to have a good convergence profile with respect to other competing algorithms.

As shown by Algorithm 1, we proceed by splitting function J in the K terms corresponding to each image component and by exploiting the separability of the closed convex constraint sets (C 1,i × C 2,i ) 1≤i≤m . The algorithm is initialized by setting the positive weights (ω 1 , ω 2 , . . . , ω m+K ) ∈]0, +∞[ m+K as well as the variables (z 1,i,0 ) 1≤i≤m+K and (z 2,i,0 ) 1≤i≤m+K associated 1 Constraint S ′ 2 will be substituted for S 2 in some of our experiments.

with the convex constraints (C 1,1 × C 2,1 , . . . , C 1,m × C 2,m ) and the split data fidelity criterion J → w → (J (1) (w), . . . , J (K) (w)).

• N otation N 1,m+1 = . . . = N 1,m+K = |A| N 2,m+1 = . . . = N 2,m+K = |A| For k = 1, 2, . . . , K do L 1,m+k = I L 2,m+k = I end For • Initialization (ω 1 , ω 2 , . . . , ω m+K ) ∈]0, +∞[ m+K (z 1,i,0 ) 1≤i≤m+K ∈ R N 1,1 × R N 1,2 × • • • × R N 1,m+K (z 2,i,0 ) 1≤i≤m+K ∈ R N 2,1 × R N 2,2 × • • • × R N 2,m+K Q 1 = m+K i=1 ω i (L 1,i ) ⊤ L 1,i -1 Q 2 = m+K i=1 ω i (L 2,i ) ⊤ L 2,i -1 u 0 = Q 1 m+K i=1 ω i (L 1,i ) ⊤ z 1,i,0 v 0 = Q 2 m+K i=1 ω i (L 2,i ) ⊤ z 2,i,0 • M ain loop For n = 0, 1, . . . do • Direct projections For i = 1, 2, . . . , m do p 1,i,n = P C 1,i (z 1,i,n ) and p 2,i,n = P C 2,i (z 2,i,n ) end For • Prox computations For k = 1, 2, . . . , K do (p 1,m+k,n , p 2,m+k,n ) = prox J (k) ω m+k (z 1,m+k,n , z 2,m+k,n ) end For • Averaging c 1,n = Q 1 m+K i=1 ω i (L 1,i ) ⊤ p 1,i,n c 2,n = Q 2 m+K i=1 ω i (L 2,i ) ⊤ p 2,i,n • Updates For i = 1, 2, . . . , m do z 1,i,n+1 = z 1,i,n + λn(L 1,i (2c 1,n -un) -p 1,i,n ) z 2,i,n+1 = z 2,i,n + λn(L 2,i (2c 2,n -vn) -p 2,i,n ) end For For k = 1, 2, . . . , K do z 1,m+k,n+1 = z 1,m+k,n + λn(2c 1,n -un - p 1,m+k,n ) z 2,m+k,n+1 = z 2,m+k,n + λn(2c 2,n -vn - p 2,m+k,n ) end For u n+1 = un + λn(c 1,n -un) v n+1 = vn + λn(c 2,n -vn) end For

Algorithm 1: PPXA+

The main loop first consists of computing the proximity operators of the data fidelity terms, as well as the direct projections onto the different convex sets mod-eling prior information. By using basic properties of proximity operators [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], it can be shown that, for every k ∈ {1, . . . , K}, ω k > 0, and

(z 1 , z 2 ) ∈ R |A| × R |A| , prox J (k) ω k (z 1 , z 2 ) = (z 1 (s) + µ (k) (s)T (k) 1 (s), z 2 (s) + µ (k) (s)T (k) 2 (s)) s∈A (30)
where z 1 = (z 1 (s)) s∈A , z 2 = (z 2 (s)) s∈A , and

µ (k) (s) =              1 γ (k) (s) prox γ (k) (s)φ (k) ω k (τ (k) (s)) -τ (k) (s) if γ (k) (s) = 0 and s ∈ O 0 otherwise, (31)
τ (k) (s) = T (k) 1 (s)z 1 (s) + T (k) 2 (s)z 2 (s) -r (k) (s) (32)
γ (k) (s) = |T (k) 1 (s)| 2 + |T (k) 2 (s)| 2 . ( 33
)
Subsequently, an averaging of the variables resulting from the computation of the projections and proximity operators is performed before the update equations. It must be emphasized that the matrix inversions required in the algorithm can be efficiently performed by making use of discrete Fourier diagonalization techniques [START_REF] Pustelnik | Relaxing tight frame condition in parallel proximal methods for signal restoration[END_REF].

To gain insights in the convergence properties of this algorithm, assume that, for every k ∈ {1, . . . , K}, φ (k) is a finite function. In our problem, the data fidelity term is convex with respect to the variable (u, v) and each individual constraint is also convex with respect to either u or v, as mentioned in sections 2.1 and 2.2. Then, provided that the sequence of relaxation parameters [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF]Prop. 5.2] allows us to guarantee that the sequence (u n , v n ) n∈N generated by the algorithm converges to a solution to problem [START_REF] Davis | BRDF invariant stereo using light transport constancy[END_REF]. The value of this parameter is set to 1.5 for all iterations. As the relaxation parameter constitutes an internal parameter of the algorithm, its choice does not modify the optimality of the solution. However, its value plays a role in the convergence speed as illustrated in Table 1. According to this table, the value of λ n has been chosen so as to maximize the convergence speed. Note finally that PPXA+ can be similarly applied in order to minimize J over S 1 ∩ S 2 or S 1 ∩ S 3 .

(λ n ) n∈N involved in Algorithm 1 is such that (∀n ∈ N) λ ≤ λ n+1 ≤ λ n < 2, where λ ∈]0, 2[,

Results

Error measures

First of all, we would like to define the quality measures we adopt. The measures below allow us to provide quantitative evaluations of our results:

1. Mean Absolute Error (MAE) MAE(z, z) = 1 |A \ O| s∈A\O |z(s) -ẑ(s)|. ( 34
)
2. Average percentage of bad pixels

Err(z, z) = 1 |A \ O| s∈A\O 1 {|z(s)-ẑ(s)|>T } ( 35
)
where T is some given positive threshold value chosen here equal to 2 and, for every s ∈ A,

1 {|z(s)-ẑ(s)|>T } = 1 if |z(s) -ẑ(s)| > T 0 otherwise. ( 36
)

Initialization

As previously mentioned, initial estimates ū and v have been obtained based on a block matching technique. Usually, two criteria are considered in the minimization performed in the block matching procedure: the sum of squared differences and the sum of absolute differences. In this work, we consider another popular measure, namely the Normalized Cross Correlation (NCC) similarity measure [START_REF] Tsai | The evaluation of normalized cross correlations for defect detection[END_REF]. Note that more sophisticated measures could be envisaged such as the one proposed in [START_REF] Yoon | Adaptive support-weight approach for correspondence search[END_REF]. NCC is defined as follows:

(∀(x, y) ∈ A)(∀u ∈ N) ψ x,y (u) = K k=1 ψ (k) x,y (u) (37)
where

ψ (k) x,y (u) = (i,j)∈B I (k) L (x + i, y + j)I (k) R (x -u + i, y + j) (i,j)∈B (I (k) L (x + i, y + j)) 2 (i,j)∈B (I (k) R (x -u + i, y + j)) 2 (38)
and B corresponds to the matching block (here of size 5 × 5) centered at pixel (x, y) in the reference image.

In order to get the corresponding pixel (xū(x, y), y) in the right image, the matching process determines ū(x, y) which maximizes the similarity measure ψ x,y :

ū(x, y) ∈ argmax u ∈ N∩[umin,umax] ψ x,y (u). (39)
Once ū(x, y) is found, v(x, y) is evaluated by the following weighted least squares estimation:

v(x, y) = K k=1 θ (k) (i,j)∈B I (k) L (x + i, y + j)I (k) R (x -ū(x, y) + i, y + j) K k=1 θ (k) (i,j)∈B I (k) L (x + i, y + j) 2 , (40)
where (θ (k) ) 1≤k≤K are nonnegative real weights. In our tests, the illumination variation is the same for all the components and, for YUV images, the best results have been obtained by choosing a weight equal to 1 for the Y component, and 0 for the U and V components. However, with the previous technique, some remaining artefacts can be observed in the disparity maps, e.g. contouring effects and oversmoothing in some areas.

In order to reduce these artifacts and limit the influence of occlusion areas, a bidirectional matching process is performed. More precisely, we firstly apply the NCC method to find the left disparity map denoted by ūL (x, y), by proceeding as described above. Similarly, we apply the NCC method to find the right disparity map denoted by ūR (x, y) by taking now the right view as the reference image. Finally, we derive the initial disparity ū(x, y) as follows:

ū(x, y) = ūR (x -ūL (x, y), y). (41)
The upper and lower bounds in the constrained formulation are estimated using prior knowledge when available or from the initial estimates ū and v, otherwise.

Algorithm implementation

In this section, we summarize the algorithm steps and parameters that are required in our joint disparity and illumination variation estimation method.

Inputs:

-Occluded areas O estimated from the initial disparity map based on the ordering [START_REF] Yuille | A generalized ordering constraint for stereo correspondence[END_REF] and border constraints.

-(T (k) 1 , T (k) 2 , r (k) ) for k ∈ {1, 2, 3}.
-Minimum and maximum estimates of the fields to be estimated. -Upper bounds for the constraints.

Parameter choice:

-For S 1 , ω 1 =100, ω 2 =100. -For S 2 , S 3 or S 4 : ω 3 =200, ω 4 =200. -For J (1) ,J (2) and J (3) : ω 5 =ω 6 =ω 7 =10. λ n ≡ 1.5.

The weights were chosen so that all the terms in the averaging step have a contribution of the same order.

Stopping criterion:

-The algorithm is stopped when the following condition is satisfied for 10 successive iterations: u n+1u n < ǫ u n where ǫ = 10 -5 .

Outputs:

-The final values of the disparity vector u n and the illumination variation vector v n .

Gray level images without variation illumination

Although the purpose of this work is to focus on situations where illumination variation arises, in this first experiment we aim at illustrating the performance of the proposed approach on grayscale images (K = 1) without any variation of the illumination field (for every s ∈ A, v(s) = 1). This allows us to compare our method with some recent state-of-the-art methods referenced in the Middlebury website 2 . In Table 2, we display the results obtained on 3 different stereo pairs (Teddy, Venus and Cones) for which we compare the results obtained by:

the correlation-based approach used to initialize our approach (NCC), -a graph-cut based approach (Gc+Occ) [START_REF] Kolmogorov | Computing visual correspondence with occlusions using graph cuts[END_REF], -a total variation based algorithm (ConvexTV) [START_REF] Pock | A convex formulation of continuous multilabel problems[END_REF], -a method based on the knowledge of Ground Control Points (GlobalGCP) [START_REF] Wang | Global stereo matching leveraged by sparse ground control points[END_REF], -a histogram-based method (HistoAggr) [START_REF] Min | A revisit to cost aggregation in stereo matching: how far can we reduce its computational redundancy?[END_REF], -a curvelet-based approach (CurveletSupWgt) [START_REF] Mukherjee | Stereo matching algorithm based on curvelet decomposition and modified support weights[END_REF], -the approach based on set theoretic estimation (subg. proj.) in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF], -our proposed method based on PPXA+ minimizing an ℓ 1 cost subject to constraints S 1,1 , S 1,2 and S ′ 1,2 . Note that the frame employed to define constraint S ′ The proposed approach appears to be competitive with respect to the other ones, although most of them are not able to deal with illumination variation. Moreover, as already mentioned, the proposed method can be implemented on a GPU architecture, thus reducing drastically the computation time. Indeed, the proximity operators and the projections onto the different convex sets can be computed in parallel. For example, with a GPU implementation, when considering the range values S 1,1 and the total variation constraints S 1,2 , the following executing times per iteration were measured: for Teddy: 21.7ms, for Venus: 5.72ms, and for Cones: 21.98ms.

Let us now turn our attention to more challenging scenarii including illumination variation both for grayscale and color images.

Gray level images with illumination variation

In this section, we evaluate our method on grayscale images (K = 1) in the presence of real and simulated illumination changes. The results are obtained on three stereo pairs depicted in Fig. 1, named "Books", "Dolls" and "Parking meter". In a first experiment, we estimate the disparity map and the illumination variation (real case) on the Books pair downloaded from MiddleBury website. We take here φ (1) = | • | for the cost function and consider the two convex constraints S 1 and S 2 for which the upper bounds were evaluated on the ground truth maps. In Fig. 2, we compare our results with the subgradient projection approach proposed in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF], where a strictly convex quadratic function is used. We also compare our approach with two naive procedures for the disparity field estimation: the correlation based approach ignoring the illumination variation and a histogram equalization technique followed by a variational approach for disparity estimation, in the spirit of the approach we developed. It turns out that the ℓ 1 cost is a good choice, and that the use of second-order constraints leads to a marginal improvement.

In a third experiment (see Fig. 5), we show the results provided by the proposed method (φ (1) = | • |) for a real image stereo pair from JISCT database3 , taken under simulated illumination variation. To introduce a significant illumination variation, the left image is kept untouched and we modified the right image by multiplying it with the Gaussian profile g depicted in Fig. 3 and defined as: where (x 0 , y 0 ) is the image center, ξ = 512 is the illumination deviation, α 1 and α 2 are given constants equal to 1.8 and -0.6 respectively. In this case, the constraint bounds in S 1 and S 2 are not known. Hence, they were firstly computed on the initial disparity and then divided by a factor 2, thus assuming that the expected result should be smoother than the initial estimate.

Conclusion

In this paper, we have investigated the application of a parallel proximal algorithm to the dense disparity estimation problem for multicomponent (e.g. color) stereo matching under illumination variation conditions. The proposed approach is flexible as it allows us to consider various convex objective functions and constraints. It is also able to exploit the potentials offered by multicore/GPU parallel architectures. However, one of the current limitations of the proposed method is that it is limited to convex cost functions, thus requiring a linearization of the original disparity model and making the estimation potentially sensitive to high values of matching errors. In our future work, we therefore plan to consider extensions of this approach to nonconvex cost functions.

1 , 2

 12 simply consists of an overcomplete Haar wavelet representation carried out over 1 resolution level. Other wavelet families have been tested without observing significant improvements.

Fig. 1

 1 Fig. 1 Stereo images: "Books" (left), "Dolls" (center) and "Parking meter" (right).

Fig. 2

 2 Fig. 2 Results for "Books" stereo pair: a)-b) ground truths, c)d) subgradient projection method (c) disparity and d) illumination fields) and e)-f) proposed approach with (u min , umax) = (20, 75),(v min , vmax) = (0.1, 1.1), (τ 2 , κ 2 ) = (74000, 400) (e) disparity and f) illumination fields), g) estimation of the disparity by NCC, h) estimation of the disparity after histogram equalization.

Fig. 3

 3 Fig. 3 Synthetic Gaussian profile.

(

  ∀(x, y) ∈ A) g(x, y) = α 1 exp -(xx 0 ) 2 + (yy 0 ) 2 2ξ 2 + α 2 (42)

4. 6

 6 Color images with illumination variationFor color images, we present tests performed on Dolls pair also downloaded from Middlebury website. Constraint sets S 1 and S 2 were used, the associated bounds being computed on the ground truth fields. The results are provided in Fig.6for YUV color images and the grayscale level images (Y component only) by using, for

Fig. 4

 4 Fig. 4 Results for Dolls stereo pair: proposed approach considering for c)-d)-ℓ 1 -cost and S 1 ∩ S ′ 2 with (u min , umax) = (20, 75), (v min , vmax) = (0.1, 1.1), (τ ′ 2 , κ 2 ) = (92000, 230), e)f)-ℓ 4 cost and S 1 ∩ S ′ 2 , g)-h) ℓ 1 cost and S 1 ∩ S 2 ∩ S 3 with (τ 2 , κ 2 ) = (34000, 230), (τ 3 , κ 3 ) = (180000, 570), i)-j)-ℓ 3/2 cost and S 1 ∩ S 2 ∩ S 3 .

Fig. 5

 5 Fig. 5 Results for "Parking meter" stereo pair. a)-b) Subgradient projection method and c)-d) our approach with (u min , umax) = (1, 8), (v min , vmax) = (0.1, 1.1), (τ 2 , κ 2 ) = (20000, 0.4).

Table 1

 1 Influence of λn on execution time (in seconds) for PPXA+ algorithm. The stopping condition described in section 4.3 is used.

	λn	0.1	0.3	0.5	0.7	0.9	1	1.1	1.3	1.5	1.7	1.9
	Teddy	90.79	30.27	18.16	12.97	10.09	9.09	8.26	6.98	6.29 6.77	7.35
	Cones	98.07	32.75	19.69	14.08	10.96	9.89	9.01	10.00	9.05 10.55	11.07

Let z ∈ R |A| denote an estimate of z ∈ R |A| .

Table 2

 2 Comparative results based on MAE for Teddy, Venus and Cones stereo pairs. The red superscript numbers represent method ranking.

http://vision.middlebury.edu/stereo/

http://vasc.ri.cmu.edu/idb/html/jisct/

Acknowledgements We would like to thank Prof. Wided Miled for providing us with her codes. We would also like to thank Dr. Raffaele Gaetano and Giovanni Chierchia for their implementation of our approach on GPU architectures.